Unexpected behavior when dtchangeable=false

[termi-official]

Describe the bug 🐞

Long-running ODE solves can accumulate round-off errors in the current time (i.e. integrator.t). Having dtchangeable=false behaves in an unexpected ways when hitting tstops which are very slighly off.

This code here errors

using OrdinaryDiffEqCore, OrdinaryDiffEqLowOrderRK
f(x,p,t) = -x
dt = 0.01
tspan = (0.0,1005.0)
prob = ODEProblem(f,[0.0],tspan)
integrator = init(prob, Euler(), dt=dt)
integrator.dtchangeable = false
add_tstop!(integrator, 1000.0)
solve!(integrator)

with

ERROR: The current setup does not allow for changing dt.

While this code here ignores that dtchangeable=false without error

using OrdinaryDiffEqCore, OrdinaryDiffEqLowOrderRK
f(x,p,t) = -x
dt = 0.01
tspan = (0.0,1000.0)
prob = ODEProblem(f,[0.0],tspan)
integrator = init(prob, Euler(), dt=dt)
integrator.dtchangeable = false
solve!(integrator)

(set in this line

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl

Line 19 in 125184e

integrator.dt = integrator.t - integrator.tprev

)

Expected behavior

I would expect that either both error, or that dt is not touched. Changing dt here looks like a dangerous endavour, as some time integration algorithms rely on the fact that dt will not change.

Minimal Reproducible Example 👇

See above. Can reproduce this in master (and older versions).

Additional context

Not sure if this is a bug or if this is just undocumented behavior. I understand why it happens and why the current logic is designed the way it is, but I think we should settle on a consistent behavior here.

[oscardssmith]

@ChrisRackauckas do you have thoughts on the intended behavior here?

[ChrisRackauckas]

The issue here is actually that Euler should be dtchangable 😅. If not dtchangable, definitely makes sense to error in most of these situations. We could allow for fudging dt close to eps or something in such cases, but it's a guardrail that exists for algorithms that cannot change dt easily. For example, if you use a multistep method, just changing dt is wrong: you need to recalculate all of your delay coefficients. Adaptive multistep methods have that built in, general fixed time multistep methods do not. So this is setup to force them to error if you ever try to change dt, as doing so would just give you incorrect answers.

So, the issue is... why did Euler become !dtchangable? That might be a bug from the package split. I'll probably fix this in a second.

[oscardssmith]

In that case, this is a bad reduction. This was originally found from a SplitIntegrator that truly isn't dtchangable. The other place this becomes tricky is that if you shrink dt (e.g. with a callbacack) in a non-adaptive, but dtchangable, method there is no mechanism for dt to go back to the previous value which can lead to dt getting ridiculously small.

[termi-official]

The issue here is actually that Euler should be dtchangable 😅.

It is. Please note that I just forced dtchangeable=false here in an attempt to have a minimal example. The error should also be reproducible with algorithms who have dtchangeable=false by default.

[ChrisRackauckas]

I see. Yeah we should allow for the epsilon changes with non-dtchangable methods.

• https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqAdamsBashforthMoulton/src/alg_utils.jl
• https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqBDF/src/alg_utils.jl
• https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqIMEXMultistep/src/alg_utils.jl

We need to mark those as not dtchangable. 😅 Just noticed that.

But okay how to do this. There is a routine for fixing small nudges around tstops:

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

Lines 289 to 302 in 5509482

	function fixed_t_for_floatingpoint_error!(integrator, ttmp)
	if has_tstop(integrator)
	tstop = integrator.tdir * first_tstop(integrator)
	if abs(ttmp - tstop) <
	100eps(float(max(integrator.t, tstop) / oneunit(integrator.t))) *
	oneunit(integrator.t)
	tstop
	else
	ttmp
	end
	else
	ttmp
	end
	end

What we need to do is allow it to stop to the tstop if the condition is met

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

Lines 47 to 51 in 5509482

	elseif integrator.dtchangeable && !integrator.force_stepfail
	# always try to step! with dtcache, but lower if a tstop
	# however, if force_stepfail then don't set to dtcache, and no tstop worry
	integrator.dt = integrator.tdir *
	min(abs(integrator.dtcache), abs(tdir_tstop - tdir_t)) # step! to the end

and also not error here if the condition is met

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

Lines 395 to 396 in 5509482

	elseif integrator.dt != integrator.dtpropose && !integrator.dtchangeable
	error("The current setup does not allow for changing dt.")

So not difficult to fix.

[oscardssmith]

what do we want to do in the case where we have something like dt = 0.1, tstops=[1.33] or something like that where increments of dt don't take us anywhere near the correct tstop? It feels like the least bad thing might be to take the step, interpolate it back to the tstop, and then change dtprev but not dt or something like that.

[termi-official]

Not sure if that is the full story. There are two more things to consider here IIUC. First, there is a potential dtmin, which needs to be taken into account here. Second, we need to recover the old dt after doing the eps step.

[oscardssmith]

dtmin mostly doesn't exist anymore (unless you specifically ask for it). Also, IIUC Chris's suggestion was to take a slightly bigger step the step before (which we already do for adaptive integrators). Currently for adaptive integrators, we will take a step that is within a tolerance of the requested dt if it lets us hit a tstop exactly, and the argument is that we possibly should do that for non-adaptive ones as well.

[termi-official]

[...] we will take a step that is within a tolerance of the requested dt if it lets us hit a tstop exactly [...]

So we do not actually touch dt, but only adjust the new t?

what do we want to do in the case where we have something like dt = 0.1, tstops=[1.33] or something like that where increments of dt don't take us anywhere near the correct tstop? It feels like the least bad thing might be to take the step, interpolate it back to the tstop, and then change dtprev but not dt or something like that.

I also thought that "overstepping", interpolating back and potentially resetting parts of the cache (e.g. for multipstep methods) should be the actual thing that happens. But that is definitely overkill for the scenario that we miss the tstop due to float point error accumulation.

[ChrisRackauckas]

It feels like the least bad thing might be to take the step, interpolate it back to the tstop, and then change dtprev but not dt or something like that.

I also thought that "overstepping", interpolating back and potentially resetting parts of the cache (e.g. for multipstep methods) should be the actual thing that happens. But that is definitely overkill for the scenario that we miss the tstop due to float point error accumulation.

We do this:

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

Lines 464 to 466 in 5509482

	DiffEqBase.change_t_via_interpolation!(integrator,
	integrator.tdir *
	pop_tstop!(integrator), Val{true})

Though we should probably automate the reinit like:

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/integrators/integrator_interface.jl

Lines 468 to 473 in 5509482

	if alg_extrapolates(integrator.alg) || !isdtchangeable(integrator.alg)
	reinit!(integrator, integrator.u;
	t0 = t,
	reset_dt = false,
	reinit_callbacks = false,
	reinit_cache = false)

non-dtchangable methods are probably the least used and tested so we might want to do a few things here.

First, there is a potential dtmin, which needs to be taken into account here.
Second, we need to recover the old dt after doing the eps step.

The system is already setup to do those checks and tweaks.

Currently for adaptive integrators, we will take a step that is within a tolerance of the requested dt if it lets us hit a tstop exactly, and the argument is that we possibly should do that for non-adaptive ones as well.

Yes exactly. The system already has a thing where if a tstop is nearby, it can just tweak to step to that. isdtchangable just disables all dt tweaks. The point is that if it's just a eps(t) sized thing to correct for floating point error in dt, it should just do the tweak since that doesn't give a major change to the truncation error. It will tweak, and it'll return to the old dt through the dtcache. And it attempts to grow instead of shrink exactly because doing an extra step of size eps() doesn't make sense.