How to use DiscreteCallbacks with sundials.jl IDA solver? How to implement reinit!?

Specific Domains Numerics
, , ,
1

I’m solving large DAE-Systems derived from PDEs for electrochemical systems. I’m using the IDA solver because the ODE solvers (DAE formulated as ODE in mass matrix form) are not working (all of them are unstable) for my DAE System. IDA works well, but unfortunately, I could not figure out yet how to use DiscreteCallbacks in combination with the IDA solver. In another issue @ChrisRackauckas gave the advice to use reinit! :

You can use reinit! as described above, but it hasn’t been integrated into the library for an automatic fix. DABDF2 and DFBDF still need a bit of work though, so I’d just point people to IDA for now. When those get optimized (which should be soon, we have the DAE benchmarks up and running as of this week to start doing it), then we’ll sweep back around to fix this. Until then the solvers aren’t too useful though since they are slow, which puts this at a lower priority.

Originally posted by @ChrisRackauckas in #589 (comment)

I am not sure how to implement it in my code:

param=[3E4;9000]
tstop=[30.0]

function condition(u,t,integrator)    
  t in tstop
end

function affect!(integrator)
  integrator.p[1] = 1E4
end

save_positions = (true,true)
cb = DiscreteCallback(condition, affect!, save_positions=save_positions)
prob_DAE = DAEProblem(electrolyzer_DAE,du0,u0_ss,tspan,param,differential_vars=diff_vars_bool)
sol_DAE = solve(prob_DAE,IDA(use_linesearch_ic=false),callback = cb, tstops=tstop)

like this it’s not working. implementation of reinit! is missing. How should reinit! be implemented?

Many thanks in advance!

2

First of all, this is a repost: Using DiscreteCallback with IDA solver for DAE systems. How to implement reinit! ? · Issue #845 · SciML/DifferentialEquations.jl · GitHub . Please only post your questions in one place.

Secondly, this is very misguided. I did not give advice to use reinit!. I said to another dev that reinit! could be used as a workaround until that issue is solved. That issue is solved, the reinit!s are done internally. You should not use reinit! like that, it will break things, don’t do it. It’s done within the package correctly which is why the issue has a linked PR and was closed.

So third, what is your actual issue? You probably answered that yourself in the first line:

Almost certainly this is pointing to an issue of having a higher index DAE which is only working in IDA because the initialization got a bit lucky (which can happen with Index 2). So first question, what is the index of the DAE? Second question, how are the other methods failing? Is it at initialization? The answer is likely to stabilize the solve.

3

Many thanks for the fast answer. I don’t know the DAE Index, I have to check it. Do you know a good way to check the index of the DAE?

The methods in mass matrix form are all failing with this error:

┌ Warning: Instability detected. Aborting
└ @ SciMLBase C:\Users\maitb\.julia\packages\SciMLBase\DXiE6\src\integrator_interface.jl:351

I think they are failing at initialization.

My main problem is that DiscretCallback is not working with IDA. I always get the error:

[IDA ERROR]  IDACalcIC
  The linesearch algorithm failed: step too small or too many backtracks.

For normal simulation from initial conditions till steady-state, the IDA solver is working well.

I also found that the DImplicitEuler() is working for simulation till steady-state and also works for DiscreteCallback. But when the applied Discontinuity in DiscreteCallback gets too large then it takes longer and longer to get a solution. For large Discontinuities, I have to stop the simulation because it takes too long.

Many thanks for your help!

4

Hi Tobias, do you have an MWE that I can run and reproduce this?

5

BTW, while we currently don’t have a utility function for checking the exact index number, you could check if your system is index 1 easily by checking the matrix pencil. Check out this blog post where it mentions matrix pencils A Survey of Differential-Algebraic Equations (1) .

6

Thanks for your help! What do you mean with MWE?

7

Minimal Working Example. A small program that allows us to reproduce your problem.

8

The model is quite large and contains different files. Maybe I could send them to you in a private message?

9

Well, minimal working example means that you invest some time and effort to reduce your code to the minimal size that allows to reproduce the problem…