SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
523 stars 199 forks source link

`OrdinaryDiffEq` arguably uses the wrong DAE initializealg after callbacks. #2143

Open oscardssmith opened 5 months ago

oscardssmith commented 5 months ago

Describe the bug 🐞 After hitting a callback, daes are re-initialized with their initializealg (https://github.com/SciML/OrdinaryDiffEq.jl/blob/31d91572c490b6eb685d6793d8bfd80389160dae/src/integrators/integrator_interface.jl#L39) rather than BrownFullBasicInit. This means that if the user is using a custom initializealg (or ShampineColocationInit), that initialization will be called after callbacks which is likely not desired. Note, however that fixing this correctly is somewhat difficult (hence the lack of a PR), because if the user is using a custom init for different reasons (e.g. turning off autodiff for initialization), we don't want to use the default BrownFullBasicInit as doing so would try to use autodiff and potentially crash. The code below reproduces this (and the problem with the obvious fix)

using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra, DiffEqBase, ADTypes, NonlinearSolve, Plots
f(u::Vector{Float64}, p, t) = u
fun = ODEFunction(f, mass_matrix=Diagonal([1., 0.]))
callback = PresetTimeCallback(0.5, Returns(nothing))
struct CustomInitialize <: DiffEqBase.DAEInitializationAlgorithm
end
function OrdinaryDiffEq._initialize_dae!(integrator, prob::ODEProblem, alg::CustomInitialize, iip)
    integrator.u .= 1
    OrdinaryDiffEq._initialize_dae!(integrator, prob, BrownFullBasicInit(;nlsolve=RobustMultiNewton(autodiff=AutoFiniteDiff())), iip)
end
prob = ODEProblem(fun, zeros(2), (0., 1.); callback, initializealg=CustomInitialize())
sol = solve(prob, FBDF(autodiff=false))
plot(sol)

image

staticfloat commented 5 months ago

So it sounds to me like one path forward might be to have two initializealg arguments, something like initializealg and reinitializealg, where reinitializealg would just default to BrownFullBasicInit? And reinitializealg would be what gets called after discontinuities?

topolarity commented 5 months ago

What information/flag does the standard initializealg use that prevents it from modifying differential variables when re-initializing after a callback?

oscardssmith commented 5 months ago

the standard initializealg is BrownFullBasicInit.

topolarity commented 5 months ago

Ah, and BrownFullBasicInit does not modify any differential variables.

ChrisRackauckas commented 4 months ago

So it sounds to me like one path forward might be to have two initializealg arguments, something like initializealg and reinitializealg, where reinitializealg would just default to BrownFullBasicInit? And reinitializealg would be what gets called after discontinuities?

In general you may want to tag callbacks with a reinitialization routine, because what variables to keep constant could depend on which ones you just modified.