CVODE_BDF: Event points are not stored

[MartinOtter]

Take the following bouncing ball model

using DifferentialEquations, Sundials

function f(du,u,p,t)
    du[1] = u[2]
    du[2] = -9.81
end

function condition(u,t,integrator)
    z = -u[1]
    return z
end

function affect(integrator)
    integrator.u[2] = -0.8*integrator.u[2]
    #auto_dt_reset!(integrator)
    #set_proposed_dt!(integrator, integrator.dt)
end

cb = ContinuousCallback(condition,affect)

u0 = [1.0,0.0]
tspan = (0.0,2.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,CVODE_BDF(),saveat=0.1,callback=cb)

println("... show time")
for t in sol.t
    @show t
end

If CVODE_BDF is used, then sol.t is:

t = 0.0
t = 0.1
t = 0.2
t = 0.3
t = 0.4
t = 0.5
t = 0.6
t = 0.7
t = 0.8
t = 0.9
t = 1.0
t = 1.1
t = 1.2
t = 1.3
t = 1.4
t = 1.5
t = 1.6
t = 1.7
t = 1.8
t = 1.9
t = 2.0

so no event points are stored in "sol". If integrator Tsit5() is used, the correct output is produced:

t = 0.0
t = 0.1
t = 0.2
t = 0.3
t = 0.4
t = 0.4515236409857255
t = 0.4515236409857255
t = 0.5
t = 0.6
t = 0.7
t = 0.8
t = 0.9
t = 1.0
t = 1.1
t = 1.1739614665628781
t = 1.1739614665628781
t = 1.2
t = 1.3
t = 1.4
t = 1.5
t = 1.6
t = 1.7
t = 1.7519117270245959
t = 1.7519117270245959
t = 1.8
t = 1.9
t = 2.0

[MartinOtter]

The following workaround does not work (savevalues! has no effect):

function affect(integrator)
    savevalues!(integrator)
    integrator.u[2] = -0.8*integrator.u[2]
    savevalues!(integrator)  
end

The following workaround works (is a "hack" - not clear whether it will always work):

function affect(integrator)
    push!(integrator.sol.t, integrator.t)
    push!(integrator.sol.u, integrator.u)    
    integrator.u[2] = -0.8*integrator.u[2]
    push!(integrator.sol.t, integrator.t)
    push!(integrator.sol.u, integrator.u)      
end

and produces the correct output:

t = 0.0
t = 0.1
t = 0.2
t = 0.3
t = 0.4
t = 0.4514751773312246
t = 0.4514751773312246
t = 0.5
t = 0.6
t = 0.7
t = 0.8
t = 0.9
t = 1.0
t = 1.1
t = 1.1738290747976405
t = 1.1738290747976405
t = 1.2
t = 1.3
t = 1.4
t = 1.5
t = 1.6
t = 1.7
t = 1.7517022432564546
t = 1.7517022432564546
t = 1.8
t = 1.9
t = 2.0

[MartinOtter]

The proposed work-around does not work, because the value of "sol.u" is wrong at the event (not clear why). Changing the workaround by first copying integrator.u into a vector and then copying this vector into sol.u, gives a correct solution:

function affect(integrator)
    push!(integrator.sol.t, deepcopy(integrator.t))
    push!(integrator.sol.u, deepcopy(integrator.u))    
    integrator.u[2] = -0.8*integrator.u[2]
    push!(integrator.sol.t, deepcopy(integrator.t))
    push!(integrator.sol.u, deepcopy(integrator.u))      
end

[ChrisRackauckas]

It sounds like the issue is just that savevalues! isn't working with Sundials. Just noticed a small typo in there:

#310

Try that branch and see if that fixes it. If so, we can make this a test case.