OrdinaryDiffEq.jl icon indicating copy to clipboard operation
OrdinaryDiffEq.jl copied to clipboard

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

Open oscardssmith opened this issue 1 year ago • 5 comments

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

oscardssmith avatar Feb 20 '24 18:02 oscardssmith

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?

staticfloat avatar Feb 20 '24 18:02 staticfloat

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

topolarity avatar Feb 20 '24 19:02 topolarity

the standard initializealg is BrownFullBasicInit.

oscardssmith avatar Feb 20 '24 19:02 oscardssmith

Ah, and BrownFullBasicInit does not modify any differential variables.

topolarity avatar Feb 20 '24 20:02 topolarity

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.

ChrisRackauckas avatar Feb 21 '24 14:02 ChrisRackauckas