FunctionCallingCallback not correctly synchronized with events

[MartinOtter]

trafficstars

The FunctionCallingCallback seems to be not correctly synchronized with events as shown with the following test model:

module Test_bouncingBall_withFunctionCallBack

using DifferentialEquations

function f(du,u,p,t) #u[1]: height, u[2]:velocity
    du[1] = u[2]
    du[2] = -9.81
    return nothing
end


# Own result storage (emulated by print statements)
function outputs!(x, t, integrator)::Nothing
    println("... Store result at time = $t: h = ", integrator.u[1], ", v = ", integrator.u[2])
    return nothing
end


# State event
function condition(out,u,t,integrator)
    out[1] = u[1] 
    return nothing
end

function affect!(integrator, event_index)
    println("... State event at ", integrator.t)
    
    # Store result before event
    outputs!(integrator.u, integrator.t, integrator)
    
    # Process event
    integrator.u[2] = -0.8*integrator.u[2]
    auto_dt_reset!(integrator)
    set_proposed_dt!(integrator, integrator.dt) 
    
    # Store result after event
    outputs!(integrator.u, integrator.t, integrator)   
    return nothing    
end

u0     = [1.0,0.0]
tspan  = (0.0,0.6)
tspan2 = 0.0:0.1:0.6 
prob = ODEProblem(f,u0,tspan)

cb1 = FunctionCallingCallback(outputs!, funcat=tspan2)
cb2 = VectorContinuousCallback(condition,affect!,1)
callbackSet = CallbackSet(cb1,cb2)

sol = solve(prob, 
            Tsit5(), 
            save_everystep=false, 
            save_start=false, 
            save_end=false,
            adaptive=true,
            dt=0.1,
            callback=callbackSet)
end

This results in the following output:

... Store result at time = 0.0: h = 1.0, v = 0.0
... Store result at time = 0.1: h = 0.9509500000000003, v = -0.9809999999999998
... State event at 0.4515236409857323
... Store result at time = 0.4515236409857323: h = 1.6653345369377348e-15, v = -4.429446918070022
... Store result at time = 0.4515236409857323: h = 1.6653345369377348e-15, v = 3.543557534456018
... Store result at time = 0.2: h = 1.6653345369377348e-15, v = 3.543557534456018
... Store result at time = 0.3: h = 1.6653345369377348e-15, v = 3.543557534456018
... Store result at time = 0.4: h = 1.6653345369377348e-15, v = 3.543557534456018
... Store result at time = 0.5: h = 0.418002671515621, v = 2.087004452526053
... Store result at time = 0.6: h = 0.418002671515621, v = 2.087004452526053

This means the affect-function is called at time = 0.45, but the FunctionCallingCallback for time = 0.2, 0.3, 0.4 is called after the affect-function is called at time = 0.45.

[MartinOtter]

This above result was generated with DifferentialEquations.jl version v6.15.0.

[MartinOtter]

Unfortunately, the issue is worse than expected:

It does not matter whether results are stored at an event point (I thought this is a work-around): The FunctionCallingCallback is still called for time instants before the event, but using the state information returned from affect!(..). This means that the result storage of saveat points before the event is completely wrong.

[MartinOtter]

Any progress here?

I have the following (bad) work-around: At an event instant, I call outputs!(..) myself for all time instants stored in integrator.sol.t/sol.x from the previous event until the current time instant and store the result in my result data structure and only afterwards process the event.

I have now an application where my current work-around is not working: I am using IDA() to avoid solving large equation systems in the model. Before calling outputs!(..), I need to compute the derivatives by interpolation in the integrator (integrator(du, t, Val{1})). However, interpolation seems to be only possible within the current step:

# Error message:
[IDAS ERROR]  IDAGetDky
  Illegal value for t. t = 0.2 is not between tcur - hu = 0.67939 and tcur = 0.77939.  # event is at 0.75

┌ Warning: IDAGetDky failed with error code =
│   flag = -26
└ @ Sundials ...\packages\Sundials\k9hc3\src\simple.jl:20

and I need it for previous steps since outputs!(..) is called for previous steps.

It seems my only option is currently to call IDA(..) directly (with all the logic for event handling and communication points), and not using DifferentialEquations.