ContinuousCallback with u[1]^3 doesn't find correct event time

[frankschae]

When replacing the condition of the ContinuousCallback used in the bouncing/free falling ball example https://diffeq.sciml.ai/stable/features/callback_functions/#Example-1:-Bouncing-Ball by u[1]^3=0, one expects to get the same event time t for both callbacks. However, one finds:

• t=3.194 for u[1]=0
• t=3.146 for u[1]^3=0.

MWE:

using OrdinaryDiffEq
condition(u,t,integrator) = u[1] # Event when event_f(u,t) == 0
condition2(u,t,integrator) = u[1]^3 # Event when event_f(u,t) == 0
function affect!(integrator)
  # for simplicity the callback doesn't do anything, except saving the state
  integrator.u[2] += zero(integrator.u[2]) #-integrator.p[2]*integrator.u[2] 
end
cb = ContinuousCallback(condition,affect!,save_positions=(true,true))
cb2 = ContinuousCallback(condition2,affect!,save_positions=(true,true))
function fiip(du,u,p,t)
  du[1] = u[2]
  du[2] = -p[1]
end
abstol=1e-4
reltol=1e-4
u0 = [50.0,0.0]
tspan = (0.0,5.0)
p = [9.8, 0.8]
prob = ODEProblem(fiip,u0,tspan,p)
sol1 = solve(prob,Tsit5(),callback=cb, abstol=abstol,reltol=reltol,
            save_everystep=true,save_start=true)
sol2 = solve(prob,Tsit5(),callback=cb2, abstol=abstol,reltol=reltol,
        save_everystep=true,save_start=true)
sol1.t != sol2.t

when dt_max=0.1 is set, the second callback results in an error:

ERROR: Double callback crossing floating pointer reducer errored. Report this issue.
Stacktrace: Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] find_callback_time(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{ContinuousCallback{typeof(condition2), typeof(affect!), typeof(affect!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Nothing, Int64}}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, callback::ContinuousCallback{typeof(condition2), typeof(affect!), typeof(affect!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Nothing, Int64}, counter::Int64)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/lULzQ/src/callbacks.jl:697
  [3] find_first_continuous_callback
    @ ~/.julia/packages/DiffEqBase/lULzQ/src/callbacks.jl:459 [inlined]
  [4] handle_callbacks!
    @ ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/integrators/integrator_utils.jl:249 [inlined]
  [5] _loopfooter!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{ContinuousCallback{typeof(condition2), typeof(affect!), typeof(affect!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Nothing, Int64}}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/integrators/integrator_utils.jl:204
  [6] loopfooter!
    @ ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/integrators/integrator_utils.jl:168 [inlined]
  [7] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{ContinuousCallback{typeof(condition2), typeof(affect!), typeof(affect!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Nothing, Int64}}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:478
  [8] #__solve#404
    @ ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:5 [inlined]
  [9] #solve_call#56
    @ ~/.julia/packages/DiffEqBase/lULzQ/src/solve.jl:61 [inlined]
 [10] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(fiip), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::Tsit5; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:callback, :abstol, :reltol, :save_everystep, :save_start, :dtmax), Tuple{ContinuousCallback{typeof(condition2), typeof(affect!), typeof(affect!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Nothing, Int64}, Float64, Float64, Bool, Bool, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/lULzQ/src/solve.jl:82
 [11] #solve#57
    @ ~/.julia/packages/DiffEqBase/lULzQ/src/solve.jl:70 [inlined]
 [12] top-level scope
    @ none:1

The final zero_func(Θ) in https://github.com/SciML/DiffEqBase.jl/blob/master/src/callbacks.jl#L700 looks relatively stable:

@show Θ
Θ1 = Θ+100eps(Θ)
Θ2 = Θ-100eps(Θ)
@show ODE_DEFAULT_NORM(zero_func(Θ),Θ)
@show ODE_DEFAULT_NORM(zero_func(Θ1),Θ1)
@show ODE_DEFAULT_NORM(zero_func(Θ2),Θ2)

Θ = 3.142840820939239
ODE_DEFAULT_NORM(zero_func(Θ), Θ) = 4.09986205491231
ODE_DEFAULT_NORM(zero_func(Θ1), Θ1) = 4.099862054901799
ODE_DEFAULT_NORM(zero_func(Θ2), Θ2) = 4.099862054922821
integrator.t = 3.142840820939239