# 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` 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 = (sqrt(2)*230*sin(2*pi*50.0*t) - u * p) *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 == 1
u
else
return 1
end
end

function affect_off!(integrator)

integrator.p = 1e10  # total resistance (load + diode)
integrator.p = 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 == 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) + u * 100)) - 1
else
ret =  -1
end

end

function affect_on!(integrator)
# println("applied affect on at t=\$(integrator.t)")
integrator.p = 100  # total resistance (load + diode)
integrator.p = 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  # <- outputs come from here
(u, )
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 = 0.0
integrator.t = 0.0010355149302289026
u = 0.02466443779613447
integrator.t = 0.0010355149302289026
u = 0.09489608425988974
integrator.t = 0.0016705724725223623
u = 0.20457147552261043
integrator.t = 0.002348954342345277
u = 0.3470141746597498
integrator.t = 0.0030822014110460615
u = 0.5151192104287117
integrator.t = 0.0030822014110460615

...truncated...

integrator.t = 0.052876422347891805
u = 0.49957763366843355
integrator.t = 0.053380394153775476
u = -1.0304261248426776e40   # <- this does not seem correct
integrator.t = 0.05444431272273977
u = -3.528390922450684e-8
integrator.t = 0.05444431272273977
u = -2.0560975872248562e-8
integrator.t = 0.06
u = -3.213201469371499e-8
integrator.t = 0.06
u = -3.254155700933845e-8
integrator.t = 0.06
u = -3.2109302361073413e-8
integrator.t = 0.06
u = -3.0880836404674164e-8
integrator.t = 0.06
u = -2.890823031148805e-8
integrator.t = 0.06
u = -2.625004076845552e-8
integrator.t = 0.06
u = -2.2971309978110092e-8
integrator.t = 0.06
u = -1.9143565658578353e-8
integrator.t = 0.06
u = -1.4844821043579996e-8
integrator.t = 0.06
u = -1.0159574882427828e-8
integrator.t = 0.06
u = -5.1788114400275675e-9
integrator.t = 0.06
u = -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 = 0.0
integrator.t = 0.0010355149302289026
u = 0.02466443779613447
integrator.t = 0.0010355149302289026
u = 0.09489608425988974
integrator.t = 0.0016705724725223623
u = 0.20457147552261043
integrator.t = 0.002348954342345277
u = 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. 