Partials of states have different values for non-stiff ODE solvers with a callback at the event location

odedual.jl

using OrdinaryDiffEq, ForwardDiff, LinearAlgebra

y0val  =  0.6
vy0val = -0.1
gyval  = 0.005
#                      val     d/dy0  d/dvy0 d/dgy
y0 = ForwardDiff.Dual( y0val,  1.0,   0.0,   0.0  )
vy0= ForwardDiff.Dual(vy0val,  0.0,   1.0,   0.0  )
gy = ForwardDiff.Dual( gyval,  0.0,   0.0,   1.0  )
u0dual = [y0;       vy0]
u0     = [y0val; vy0val]
Δy, Δvy, Δgy = 1.0e-4 .* (y0val, vy0val, gyval)
tspan = (0.0, 2.0)

function st!(du, u, p, t)
    @inbounds begin
        du[1] = u[2]
        du[2] = (1-u[1]) - p
    end
    return nothing
end

cnd(u, t, integ) = u[1] - 1
CC = ContinuousCallback(cnd, terminate!, save_positions=(true,false),
                        interp_points=100, affect_neg! = nothing, abstol=0, reltol=0)

prbdual    = ODEProblem(st!, u0dual, tspan, gy    )
prbcenter  = ODEProblem(st!, u0    , tspan, gyval )
prb_y_back     = ODEProblem(st!, u0 .- [Δy;   0], tspan, gyval )
prb_y_forward  = ODEProblem(st!, u0 .+ [Δy;   0], tspan, gyval )
prb_vy_back    = ODEProblem(st!, u0 .- [0;  Δvy], tspan, gyval )
prb_vy_forward = ODEProblem(st!, u0 .+ [0;  Δvy], tspan, gyval )
prb_gy_back    = ODEProblem(st!, u0            , tspan, gyval-Δgy)
prb_gy_forward = ODEProblem(st!, u0            , tspan, gyval+Δgy)

sldual1 = solve(prbdual, Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0)
sldual2 = solve(prbdual, Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0, callback=CC)

slcenter      = solve(prbcenter,      Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0)
sl_y_back     = solve(prb_y_back,     Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0)
sl_y_forward  = solve(prb_y_forward,  Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0)
sl_vy_back    = solve(prb_vy_back,    Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0)
sl_vy_forward = solve(prb_vy_forward, Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0)
sl_gy_back    = solve(prb_gy_back,    Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0)
sl_gy_forward = solve(prb_gy_forward, Vern9(), abstol=1e-12, reltol=1e-12, tstops=0.0:0.01:2.0)

#Event time
te = sldual2.t[end].value

#Numerical differentiation du(te)/dy0 without callback
duₜₑ_dy₀_numerical  = (sl_y_forward(te)  .- sl_y_back(te) ) ./ (2Δy)
#Dual differentiation du(te)/dy0 without callback
getindex.(getproperty.(sldual1(te), :partials), 1)
#Dual differentiation du(te)/dy0 with callback
getindex.(getproperty.(sldual2[end], :partials), 1)
"""
julia> duₜₑ_dy₀_numerical  = (sl_y_forward(te)  .- sl_y_back(te) ) ./ (2Δy)
2-element Array{Float64,1}:
 -0.257299187312691
 -0.966331789921393

julia> getindex.(getproperty.(sldual1(te), :partials), 1)
2-element Array{Float64,1}:
 -0.2572991873207459
 -0.9663317899169441

julia> getindex.(getproperty.(sldual2[end], :partials), 1)
2-element Array{Float64,1}:
 -8.765352157211923e-17
 -0.969489370005787
"""

#Numerical differentiation du(te)/dvy0 without callback
duₜₑ_dvy₀_numerical = (sl_vy_forward(te) .- sl_vy_back(te)) ./ (2Δvy)
#Dual differentiation du(te)/dvy0 without callback
getindex.(getproperty.(sldual1(te), :partials), 2)
#Dual differentiation du(te)/dvy0 with callback
getindex.(getproperty.(sldual2[end], :partials), 2)
"""
julia> duₜₑ_dvy₀_numerical = (sl_vy_forward(te) .- sl_vy_back(te)) ./ (2Δvy)
2-element Array{Float64,1}:
  0.9663317899255562
 -0.2572991872951125

julia> getindex.(getproperty.(sldual1(te), :partials), 2)
2-element Array{Float64,1}:
  0.9663317899169441
 -0.2572991873207459

julia> getindex.(getproperty.(sldual2[end], :partials), 2)
2-element Array{Float64,1}:
  5.670180841891999e-16
 -0.24544034683690863
"""

#Numerical differentiation du(te)/dgy without callback
duₜₑ_dgy_numerical  = (sl_gy_forward(te) .- sl_gy_back(te)) ./ (2Δgy)
#Dual differentiation du(te)/dgy without callback
getindex.(getproperty.(sldual1(te), :partials), 3)
#Dual differentiation du(te)/dgy with callback
getindex.(getproperty.(sldual2[end], :partials), 3)
"""
julia> duₜₑ_dgy_numerical  = (sl_gy_forward(te) .- sl_gy_back(te)) ./ (2Δgy)
2-element Array{Float64,1}:
 -1.2572991878734283
 -0.9663317898089828

julia> getindex.(getproperty.(sldual1(te), :partials), 3)
2-element Array{Float64,1}:
 -1.2572991873207455
 -0.9663317899169441

julia> getindex.(getproperty.(sldual2[end], :partials), 3)
2-element Array{Float64,1}:
 -4.651293499446739e-16
 -0.9817613873476325
"""