diff --git a/src/preset_time.jl b/src/preset_time.jl index 3b7e3f95..df9728af 100644 --- a/src/preset_time.jl +++ b/src/preset_time.jl @@ -6,9 +6,9 @@ struct PresetTimeFunction{T, T2, T3} end function (f::PresetTimeFunction)(u, t, integrator) return if hasproperty(integrator, :dt) - insorted(t, f.tstops) && (integrator.t - integrator.dt) != integrator.t + is_preset_time(t, f.tstops) && (integrator.t - integrator.dt) != integrator.t else - insorted(t, f.tstops) + is_preset_time(t, f.tstops) end end @@ -26,11 +26,28 @@ function (f::PresetTimeFunction)(c, u, t, integrator) for tstop in _tstops add_tstop!(integrator, tstop) end - return if insorted(t, tstops) + return if is_preset_time(t, tstops) f.user_affect!(integrator) end end +function is_preset_time(t, tstops) + insorted(t, tstops) && return true + i = searchsortedfirst(tstops, t) + return (i <= length(tstops) && is_close_preset_time(t, tstops[i])) || + (i > firstindex(tstops) && is_close_preset_time(t, tstops[i - 1])) +end + +function is_close_preset_time(t, tstop) + _t, _tstop = promote(t, tstop) + if _t isa AbstractFloat && _tstop isa AbstractFloat && isfinite(_t) && isfinite(_tstop) + scale = max(abs(_t), abs(_tstop)) + return abs(_t - _tstop) <= 100 * eps(scale) + end + + return t == tstop +end + """ ```julia PresetTimeCallback(tstops, user_affect!; diff --git a/test/preset_time.jl b/test/preset_time.jl index c88f7bee..1ec8e717 100644 --- a/test/preset_time.jl +++ b/test/preset_time.jl @@ -34,6 +34,18 @@ sol = integrator.sol @test 0.6 ∈ sol.t @test p != startp +struct PresetTimeRoundoffIntegrator{T} + t::T + dt::T +end + +roundoff_integrator = PresetTimeRoundoffIntegrator(8.999999999999996, 0.1) +for tstops in ([2, 9], range(2, 9; length = 2), [2.0, 9.0]) + roundoff_cb = PresetTimeCallback(tstops, integrator -> nothing) + @test roundoff_cb.condition(nothing, 8.999999999999996, roundoff_integrator) + @test !roundoff_cb.condition(nothing, 9.0 - 1.0e-10, roundoff_integrator) +end + notcalled = true prob = ODEProblem(some_dynamics, u0, tspan, p) cb = PresetTimeCallback([1.2], integrator -> notcalled = false)