DAE (Mass Matrix) Solver Compatibility with Callbacks

[taylormcd]

The following mass matrix DAE solvers fail to solve when used with callbacks that modify the current state: Rodas4, Rodas42, Rodas4P, Rodas4P2, Rodas5.

MWE:

# original problem
using DifferentialEquations
function f(du, u, p, t)
    du[1] = -p[1]*u[1] + p[2]*u[2]*u[3]
    du[2] = p[1]*u[1] - p[2]*u[2]*u[3] - p[3]*u[2]*u[2]
    du[3] = u[1] + u[2] + u[3] - 1.
end
M = [1 0 0; 0 1 0; 0 0 0]
p = [0.04, 10^4, 3e7]
u0 = [1.,0.,0.]
tspan = (0., 1e6)
prob = ODEProblem(ODEFunction(f, mass_matrix = M), u0, tspan, p)
sol = solve(prob, Rodas5()) # solves fine

# problem with callback that resets the state to the initial conditions
tstops = [0.1]
affect!(integrator) = integrator.u .= u0
cb = PresetTimeCallback(tstops, affect!)
cprob = ODEProblem(ODEFunction(f, mass_matrix = M), u0, tspan, p; callback = cb)
csol = solve(cprob, Rodas5()) # aborts due to instability right after callback (at t = 0.1)

Note that while this example uses Rodas5 as the solver, the same result is obtained when using any of the aforementioned solvers.

I tested the same problem using the following solvers and they do not have the same issue: Ros3P, Rodas3, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A, Ros4LStab

[ChrisRackauckas]

This isn't just mass matrix: both need to reinit after an event.