# Wrong interpolated value in `SavingCallback` output

This is a follow-up question on DiffEq skips over event in the first `solve` but works after `reinit!` - #16 by ChrisRackauckas. The code simulates a circuit in which an AC power source is supplying an RL load through a diode. I used a `SavingCallback` to save `u[1]` at fixed intervals, but the result contains data point(s) that are extremely large.

The code below is a bit long, but the suspect issue is related to the `SavingCallback`.

``````using DifferentialEquations
using Plots
using BenchmarkTools

du[1] = (sqrt(2)*230*sin(2*pi*50.0*t) - u[1] * p[1]) *2
end

dt = [];

# test if diode needs to be turned off
function condition_off(u, t, integrator)
# println("test affect off at \$(integrator.t), u=\$(integrator.u)")
# println("In condition_off")
# @show integrator.t integrator.dt

# push!(dt, integrator.dt)
if integrator.p[2] == 1
u[1]
else
return 1
end
end

function affect_off!(integrator)

integrator.p[1] = 1e10  # total resistance (load + diode)
integrator.p[2] = 0
# println("In affect_off!")

end

# test if diode needs to be turned on
function condition_on(u, t, integrator)
# println("In condition_on")

if integrator.p[2] == 0
# voltage across the diode minus the minimum turn-on voltage (1)
ret = (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator)[1] + u[1] * 100)) - 1
else
ret =  -1
end

end

function affect_on!(integrator)
# println("applied affect on at t=\$(integrator.t)")
integrator.p[1] = 100  # total resistance (load + diode)
integrator.p[2] = 1

end

cb_off = ContinuousCallback(condition_off, nothing, affect_off!,
repeat_nudge=1//10)  # only when condition goes negative
cb_on = ContinuousCallback(condition_on, affect_on!, nothing,
repeat_nudge=1//10)  # only when condition goes positive

saved_values = SavedValues(Float64, Tuple{Float64})

function save_u(u, t, integrator)
@show integrator.t u[1]  # <- outputs come from here
(u[1], )
end

cb_save = SavingCallback(save_u, saved_values,
saveat=0:1/2000:0.06,
)

cb_set = CallbackSet(cb_off, cb_on, cb_save);

u0 = [0.0]
p = [100., 1.] # resistance, diode flag
tspan = (0.0, 0.06)

prob = ODEProblem(I_load!, u0, tspan, p)

solver = Rodas5P()

integrator = init(prob, solver,
callback=cb_set,
abstol=1e-5,
reltol = 1e-6,
dtmax=1/60,
saveat=0:1/2000:0.06;
);

sol = solve!(integrator)
plot(sol, label="I Rodas5P")
``````

Output:

``````integrator.t = 9.999999999999999e-5
u[1] = 0.0
integrator.t = 0.0010355149302289026
u[1] = 0.02466443779613447
integrator.t = 0.0010355149302289026
u[1] = 0.09489608425988974
integrator.t = 0.0016705724725223623
u[1] = 0.20457147552261043
integrator.t = 0.002348954342345277
u[1] = 0.3470141746597498
integrator.t = 0.0030822014110460615
u[1] = 0.5151192104287117
integrator.t = 0.0030822014110460615

...truncated...

integrator.t = 0.052876422347891805
u[1] = 0.49957763366843355
integrator.t = 0.053380394153775476
u[1] = -1.0304261248426776e40   # <- this does not seem correct
integrator.t = 0.05444431272273977
u[1] = -3.528390922450684e-8
integrator.t = 0.05444431272273977
u[1] = -2.0560975872248562e-8
integrator.t = 0.06
u[1] = -3.213201469371499e-8
integrator.t = 0.06
u[1] = -3.254155700933845e-8
integrator.t = 0.06
u[1] = -3.2109302361073413e-8
integrator.t = 0.06
u[1] = -3.0880836404674164e-8
integrator.t = 0.06
u[1] = -2.890823031148805e-8
integrator.t = 0.06
u[1] = -2.625004076845552e-8
integrator.t = 0.06
u[1] = -2.2971309978110092e-8
integrator.t = 0.06
u[1] = -1.9143565658578353e-8
integrator.t = 0.06
u[1] = -1.4844821043579996e-8
integrator.t = 0.06
u[1] = -1.0159574882427828e-8
integrator.t = 0.06
u[1] = -5.1788114400275675e-9
integrator.t = 0.06
u[1] = -4.968781781758675e-16
``````

The issue is that it silently fails. If the `abstol` and `reltol` are set to other values, the extremely large values will appear at other times (but all close to discontinuities).

While debugging this issue, I found another issue. `SavingCallback` saves timestamps at the given `saveat`, and the saved values match that in `sol.u`. The issue is that the reported `integrator.t` from within the saving callback differs from the `saveat` times.

``````julia> saved_values.t
121-element Vector{Float64}:
0.0
0.0005
0.001
0.0015
0.002
0.0025
0.003

... truncated ...
``````
``````julia> saved_values.saveval
121-element Vector{Tuple{Float64}}:
(0.0,)
(0.02466443779613447,)
(0.09489608425988974,)
(0.20457147552261043,)
(0.3470141746597498,)
(0.5151192104287117,)
(0.7014919141095118,)

... truncated ...
``````

Output from the saving callback:

``````integrator.t = 9.999999999999999e-5
u[1] = 0.0
integrator.t = 0.0010355149302289026
u[1] = 0.02466443779613447
integrator.t = 0.0010355149302289026
u[1] = 0.09489608425988974
integrator.t = 0.0016705724725223623
u[1] = 0.20457147552261043
integrator.t = 0.002348954342345277
u[1] = 0.3470141746597498
``````

What did I miss here? Is it because `integrator.t` is not modified while interpolating? Thanks.

(edit: compared the saved values with `sol` and found it to be correct, except for the extremely large value).

Could this be a bug in the interpolation routine? I lowered the tolerances and lowered the `repeat_nudge`, but the issue persists.

I tracked the issue down to the use of `first(get_tmp_cache(integrator))` in `SavingAffect`.

The first value saved in the cache seen by the very last `affect!` function is close to zero. But the retrieved first value by the saving callback from the cache is mutated.

The print-out of `get_tmp_cache(integrator)` from the `ContinuousCallBack` looks like

``````curt=0.013, curu=[0.3303721934837517] before giving to integrator
curt=0.013, curu=[0.21664209586008248] inplace update currentu
In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
In Cond off: ([-0.0006075553301955638], [0.0009123193020741382])
In Cond off: ([-1.2564339208389116e-6], [0.0009123193020741382])
In Cond off: ([-2.593250008765789e-9], [0.0009123193020741382])
In Cond off: ([-5.34254518576491e-12], [0.0009123193020741382])
In Cond off: ([-1.865278874804446e-15], [0.0009123193020741382])
In Cond off: ([8.97371217660693e-15], [0.0009123193020741382])
In Cond off: ([3.0533229421397907e-15], [0.0009123193020741382])
In Cond off: ([1.0908325831142319e-15], [0.0009123193020741382])
In Cond off: ([-8.701559074838522e-16], [0.0009123193020741382])
In Cond off: ([9.721148422111273e-17], [0.0009123193020741382])
In Cond off: ([9.721148422111273e-17], [0.0009123193020741382])
In Affect off: ([-1.3433720342587264e-9], [2.674489451237605e-5])  # !!<- first here does not match the one below
curt=0.0132, curu=[-2.4091928323155922e10] before giving to integrator
curt=0.0132, curu=[2.2505962311451964e9] inplace update currentu
In Cond off: ([7.978805348180991e-8], [-0.7438025523740635])
``````

At other normal steps, the value in the cache seen by a `ContinuousCallBack` is what the next `SavingCallBack` retrieves, such as

``````  In Cond off: ([-0.0002850830064895942], [3.039545393851428e-8])
In Cond off: ([-0.00028514796297349465], [3.039545393851428e-8])
In Cond off: ([-0.0002852128703656998], [3.039545393851428e-8])
In Cond off: ([-0.0002852128703656998], [3.039545393851428e-8]) # !!<- first here matches below
curt=0.0134, curu=[-0.0002852128703656998] before giving to integrator
curt=0.0134, curu=[-0.0002850108807624704] inplace update currentu
``````

Iâ€™m stuck here because Iâ€™m not sure where the cache gets mutated. Any help is appreciated!

This is the same issue as SavedValues not stored at events Â· Issue #557 Â· SciML/DifferentialEquations.jl Â· GitHub and FunctionCallingCallback not correctly synchronized with events Â· Issue #126 Â· SciML/DiffEqCallbacks.jl Â· GitHub . I actually plan to fix it next week because I thought of an elegant solution. The issue is that discrete callbacks are applied prior to the continuous callback, so you get values from the â€śpre pullbackâ€ť phase in the SavingCallback/FunctionCallingCallback. What actually needs to happen is that they need to happen after the continuous callback has changed the time interval, or have a choice in where to be applied.

Chris, thanks a lot for your reply. I started debugging into OrdinaryDiffEq.jl but will now wait for a fix.

By the way, do you know if the current Julia Debugger in VSCode can be used to step into DiffEq? I got an error in `FunctionWrappers.jl` regarding `llvmcall` not being supported in an interpreter.

Just use `SciMLBase.FullSpecialize` in the construction to avoid the funciton wrappers.

1 Like

Thanks! Now I recall reading that in the FAQ.

Iâ€™m not sure if thatâ€™s the cause of my case. In `handle_callbacks!`, continuous callbacks are applied prior to discrete callbacks. This may be a separate issue.

The `u` value in the cache suddenly changes to a very large value after calling `reeval_internals_due_to_modification!` of the continuous callback. More specifically, `_ode_addsteps!` updated the cached to a large value.

I stepped into `_ode_addsteps!` invoked by `reeval_internals_due_to_modification!`. In the former, `cache.tmp` keeps diverging (as the norm of k1 through k8 keeps increasing). It might be a numerical issue, but I could be wrong.

How can I force a Jacobian update in `_ode_addsteps!`? I saw a comment in that function saying

``````        ### Jacobian does not need to be re-evaluated after an event
### Since it's unchanged
jacobian2W!(W, mass_matrix, dtgamma, J, true)
``````

The event in my code changed a parameter and thus the Jacobian. Iâ€™m poking aroundâ€¦

No that happens after to it. Whatâ€™s going on is that you start with an interval `(t, t+dt)`. Then you find an event at `t + theta * dt` for `theta in [0,1]`. You then shrink the interval to `(t, t + theta * dt)` and re-evaluate the stepper to get a valid interpolation over that new interval. In this interval youâ€™d then want to run the callbacks. But they should know that the real `dt` is `theta * dt`.

I figured out how your described event timing and time modification work. I suspect that due to a silent inexact error, `t + theta * dt` ends up on the other side of the zero crossing. In this case, shouldnâ€™t the Jacobian be reevaluated? My problem has large parameters `1e10` so the exact error can be possible.

If I reduce that huge parameter from `1e10` to `1e3`, the result is correct. The solution is different due to the parameter change.

If the parameter is changed to `1e4`, the issue appears:

If the parameter is increased to `1e5`, the issue appears at different times as shown below. All these times are in the proximity of events.

Iâ€™ll stop here and wait for your fix. Thank you for your patience.