Rosenbrock solvers fail or interpolate poorly due to (algebraically) inconsistent `u` after callback

[topolarity]

Describe the bug 🐞

Using a DiscreteCallback (specifically a PresetTimeCallback) with a DAE appears to be poorly behaved with a couple solvers.

Expected behavior

I expected each solver to solve from the discontinuity forward "as if" the simulation had been restarted.

Minimal Reproducible Example

julia> using OrdinaryDiffEq, DiffEqCallbacks, LinearAlgebra
julia> f(u, p, t) = [u[2], u[2] - p * sin(t)]
julia> mass_matrix = Diagonal([1.0,0.0])
julia> prob = ODEProblem(ODEFunction(f; mass_matrix), [0.0, 0.0], (0., 2π), 1,
                         callback=PresetTimeCallback(0:.3:2π, integ->(integ.p=-integ.p;)))

julia> sol = solve(prob, Rodas5P())
julia> plot(sol)

julia> sol = solve(prob, Rosenbrock23())
julia> plot(sol; xlim=(0, 2π))

Rodas5P has a bad interpolant for a couple timesteps after each discontinuity:

Rosenbrock23 fails to solve past t = 0.3 (and ROS34PRw fails similarly):

[topolarity]

For reference, FBDF appears to solve this correctly:

[ChrisRackauckas]

Only algebraic variables?

[topolarity]

As far as I've noticed, yeah - Only algebraic discontinuities are an issue

[topolarity]

IDA also appears to struggle to get past the first stop:

julia> using OrdinaryDiffEq, DiffEqCallbacks, Sundials
julia> f(du, u, p, t) = [du[1] - u[2], u[2] - p * sin(t)]
julia> prob = DAEProblem(DAEFunction(f), [0.0, 0.0], [0.0, 0.0], (0., 2π), 1,
                         callback=PresetTimeCallback(0:.3:2π, integ->(integ.p=-integ.p;)),
                         differential_vars=[true,false])
julia> sol = solve(prob, IDA())
julia> plot(sol; xlim=(0, 2π))

[ChrisRackauckas]

Is this all on latest?

[topolarity]

Yes it is - 6.70.0

[gstein3m]

Rosenbrock methods require consistent initial values for DAEs, i.e. algebraic equations should be fulfilled. If you induce discontinuities in the algebraic equations, the current solution is no longer consistent. Therefore, the integrator will usually abort or the interpolation is no longer correct. The following modification should solve the problem:

function affect!(integrator) 
    integrator.p = -integrator.p
    integrator.u[2] = -integrator.u[2]
end

prob = ODEProblem(ODEFunction(f; mass_matrix), [0.0, 0.0], (0., 2π), 1,
                               callback=PresetTimeCallback(0:.3:2π, affect!))

[topolarity]

Yeah, the problem is certainly that u is inconsistent after the discontinuity, although I'm not so sure that it should the responsibility of the callback to find a consistent starting condition (especially since it can require a non-linear solve in general).

To take the analogy with initialization (where we have the same problem):

julia> prob = ODEProblem(ODEFunction(f; mass_matrix), [0.0, 100.0], (0., 2π), 1) # provide inconsistent u0[2]
julia> sol = solve(prob, Rosenbrock23())
retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 12834-element Vector{Float64}:
 0.0
 1.0e-6
 ⋮
 6.283185307179586
u: 12834-element Vector{Vector{Float64}}:
 [0.0, 0.0]
 [4.999999999999791e-13, 9.999999999999288e-7]
 ⋮
 [-3.633302206017213e-11, -2.442385622089885e-16]

The inconsistent conditions are resolved automatically by the solve routine.

[ChrisRackauckas]

We should do a reinitialization after u_modified! on callbacks. We do it on DAEs but I don't think we do it on the mass matrix ODEs right now.

[oscardssmith]

I think this might also explain some of the weird behavior I've seen with FBDF timestepping after discontinuities.

[topolarity]

Do we know why IDA fails to solve? Is it also missing a reinit?

[ChrisRackauckas]

IDA has a reinitialization. It must not be reinitializable with your choice of initialize. Brown basic?