Ordering of isoutofdomain and continuous callbacks in DifferentialEquations.jl

Specific Domains Numerics
1

I am encountering some undesired behavior with a combination of DifferentialEquations.jl’s marvelous isoutofdomain and event handling features. I seem to be having integrations fail because of the ordering of a Continuous Callback and isoutofdomain when evaluating whether a step should be accepted/rejected:

  1. Step passes isoutofdomain condition
  2. Continuous Callback occurs before end of step
  3. Interpolated solution at continuous callback root violates isoutofdomain condition but is allowed to occur
  4. Integration fails from dt <= dtmin because domain violated

(These are my deduction from diagnostics in my code because I can’t find where exactly the source code is :no_mouth:…)

My goal is to get the solver to reject the proposed timestep at point 3 because it violates the domain.

Here is a kind of stupid MWE of this behavior:

using DifferentialEquations

# set up differential equation problem 
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
# make negative numbers out of domain
isoutofdomain(u,p,t) = u < 0. 
# set up continuous callback when u = 1 with the effect of setting u = -1 
# (i.e., out of the domain)
findu1(u,t,integrator) = u - 1.
makeunegative!(integrator) = integrator.u = -1.
cb = ContinuousCallback(findu1,makeunegative!)
# solve problem 
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8,isoutofdomain=isoutofdomain,callback=cb)
# see that solver accepts step from ContinuousCallback 
# allowing out of domain modification of u 
println("last accepted u = $(sol.u[end]) where isoutofdomain = $(isoutofdomain(sol.u[end],nothing,sol.t[end]))")
# prints: last accepted u = -1.0 where isoutofdomain = true

In this scenario, my goal from above translates into the integration failing at (something like) sol.u[end]=0.99999 rather than sol.u[end]=-1.

Is there a way I can use the affect! function input into the ContinuousCallback to have the integrator reject the step the ContinuousCallback wants to make when that step violates isoutofdomain? (In my actual problem, the integration is able to satisfy both the solution domain conditions and the Continuous Callback effect if a smaller timestep is taken unlike in this MWE.)

2

You can try some of these options here:
https://diffeq.sciml.ai/stable/features/callback_functions/

In particular:

  • rootfind = SciMLBase.NoRootFind or interp_points = 2? That could help the integrator…
    (But I’m not sure if that fixes the problem)

Just for understanding, so the solution starts already outside of the domain? So, could it be that the problem is also that the root finding would actually never find a root inside the interval of the timestep?

3

Thanks, just to check that I understand the documentation correctly: rootfind = SciMLBase.NoRootFind basically turns a Continuous Callback into a Discrete Callback that only applies affect! at the end of a timestep without modifying the timestep?

Just for understanding, so the solution starts already outside of the domain? So, could it be that the problem is also that the root finding would actually never find a root inside the interval of the timestep?

There isn’t a problem with the root finding. The interpolated solution at the root is out of the domain, so when the callback is applied, the step taken is to an out of domain point from which the solver can’t move forward and crashes.

I tried putting into the event affect! a conditional like this

if isoutofdomain(integrator.u,integrator.p,integrator.t)
	SciMLBase.set_proposed_dt!(integrator,integrator.t-integrator.tprev)
	integrator.u = integrator.uprev
	integrator.t = integrator.tprev
end

to try to undo the step, but then I still messed up the solution. But I don’t really understand the different save options, which is maybe why that approach didn’t work?

4

Callbacks are only ran on accepted time steps, and if a solution does not pass isoutofdomain then the step is rejected. The ordering is here:

The key to your question seems to be:

It is not guaranteed that the whole interpolation is within the domain, just the end point. We could in theory change the API here so that isoutofdomain has the integrator so that you can check the interpolant. And with the interpolant definition, we could specialize max and min to give back the maximum and minimum value of the interpolant, which would make this cheap and easy. But we do not have all of that machinery implemented. It’s worth opening an issue though.