`PresetTimeCallback` at t=0. triggers `dt <= dtmin` if f(u) = 0.

[fabern]

I'm trying to use PresetTimeCallback and solve an ODE with it. However, there seems to be some kind of instability with a Callback at t=0.0 and when my ODE returns 0.0, i.e. no change.

See below a MWE:

using DifferentialEquations

u0 = [0.02]
p  = (1.0) # unused
tspan = (-0.0001, 150.)

function f(du,u,p,t)
    du[1]=0.0
end
function jump_solution!(integrator)
    integrator.u[1] += 0.0
end

# Without callback it solves without error
sol_0 = solve(ODEProblem(f,u0,tspan,p), progress = true);
# With empty callback it solves without error
jump_cb_1 = PresetTimeCallback([],jump_solution!);
sol_1 = solve(ODEProblem(f,u0,tspan,p,callback = jump_cb_1), progress = true);
# With callback at 0.1 it solves without error
jump_cb_2 = PresetTimeCallback([0.1],jump_solution!);
sol_2 = solve(ODEProblem(f,u0,tspan,p,callback = jump_cb_2), progress = true);

# With callback at 0.0 it throws an error
jump_cb_3 = PresetTimeCallback([0.0],jump_solution!);
sol_3 = solve(ODEProblem(f,u0,tspan,p,callback = jump_cb_3), progress = true);

The third trial with PresetTimeCallback([0.0],jump_solution!) returns:

┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ DiffEqBase ... /.julia/packages/DiffEqBase/TXlS1/src/integrator_interface.jl:343

Does this make sense to somebody? Might this be linked to https://github.com/SciML/DifferentialEquations.jl/issues/601 ?

[fabern]

Note that if the ODE doesn't return 0., the callback at t=0. seems to work.

# With callback at 0.0, but f not returning 0.0, it solves without error
function f_01(du,u,p,t)
    du[1]=0.1
end
sol_3b = solve(ODEProblem(f_01,u0,tspan,p,callback = jump_cb_3), progress = true);

[ChrisRackauckas]

Yeah this is an odd issue. The issue is that the default dt happens to be 1e-20 away from the tstop that is generated for the callback, so it steps dt, then it tries to step 1e-20, realizes that's below dtmin and ends the integration. It should just grow that first dt by 1e-20, which it does in the other steps. The workaround is just to do like dt=1e-3 to give it a different starting dt to have it around, but we should make sure to avoid this issue.

[ChrisRackauckas]

It's another instance of https://github.com/SciML/DifferentialEquations.jl/issues/570

[fabern]

Thanks, @ChrisRackauckas . Glad to hear that you are aware of the cause of this behaviour.

Indeed, dt=1e-3 solves the issue, and should be perfectly fine in my use case.

sol_3c = solve(ODEProblem(f,u0,tspan,p,callback = jump_cb_3), 
    progress = true, dt=1e-3)

Note that any dts down to 1e-15 solve the issue, but at 1e-16 it reappears.

[ChrisRackauckas]

but at 1e-16 it reappears.

1e-16 < eps(Float64), so it's below dtmin.