SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.84k stars 225 forks source link

Using DiscreteCallback with IDA solver for DAE systems. How to implement reinit! ? #845

Open TobiasFranzBRB opened 2 years ago

TobiasFranzBRB commented 2 years ago

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 https://github.com/SciML/DifferentialEquations.jl/issues/589#issuecomment-995248038

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 its not working. implementation of reinit! is missing. How should reinit! be implemented?

Many thanks in advance!

ChrisRackauckas commented 2 years ago

I answered this in https://discourse.julialang.org/t/how-to-use-discretecallbacks-with-sundials-jl-ida-solver-how-to-implement-reinit/75674 . Please don't post in two places.