Wrong interpolated value in `SavingCallback` output

[cuihantao]

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

function I_load!(du, u, p, t)
    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).

[cuihantao]

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).

[cuihantao]

@cuihantao

u[1] = -1.0304261248426776e40   # <- this does not seem correct

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

[cuihantao]

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!

[ChrisRackauckas]

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.

[cuihantao]

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.

[ChrisRackauckas]

@cuihantao

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.

[cuihantao]

Thanks! Now I recall reading that in the FAQ.

[cuihantao]

@ChrisRackauckas

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.

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.

[cuihantao]

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…

[ChrisRackauckas]

@cuihantao

The event in my code changed a parameter and thus the Jacobian.

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.

[cuihantao]

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.

[cuihantao]

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.

[cuihantao]

I’ll stop here and wait for your fix. Thank you for your patience.