Callbacks on matrices with Sundials solver

[Antomek]

I have been running into issues with using callbacks on a solution which looks like a matrix, when I passed an algorithm hint to the solver (e.g. alg_hints = [:stiff]). The solution is found when this hint is not passed. I think I traced the issue down to using a solver from Sundials with a matrix shaped solution. I believe that Sundials maps the matrix onto a vector, which messes up the event finding.

A minimal example that reproduces the error:

using DifferentialEquations, Sundials, Plots

I = 6.

function f(du, u, p, t)
    du[1, 1] = -u[1,1] + I
    du[1, 2] = -u[1, 2]
    du[2, 1] = -0.9 * u[2,1] + I
    du[2, 2] = -u[2, 2]
end

u0 = [0. 0.; 0. 0.]

Th = 5.
Reset = 0.
Reset2 = 1.

function condition(u,t,integrator) # Event when event_f(u,t) == 0
    maximum(u[:, 1]) - Th
end

function affect!(integrator) # What to do upon event
    u = integrator.u
    maxidx = findmax(u[:,1])[2]
    u[maxidx,1] = Reset
    u[maxidx,2] = u[maxidx,2] + Reset2
    nothing
end

callback = ContinuousCallback(condition, affect!, nothing)

tspan = (0.0, 10.0)
prob = ODEProblem(f, u0, tspan)

#@time sol = solve(prob; alg_hints=[:stiff], callback = callback)
@time sol = solve(prob, CVODE_BDF(), callback = callback);
#@time sol = solve(prob, callback = callback);
display(plot(sol))

I will note in passing that Sundials seems to have an issue with VectorContinuousCallback as well, in the same vein. I will explore this further and open an issue when I have found a minimal example that reproduces this.

[ChrisRackauckas]

Thanks for pointing this out. Yeah, we may need to be a bit more careful with the reshapes to deal with this.