Using DiffEq isoutofdomain to enforce positivity

General Usage
1

I’m trying to enforce positivity of my solution using the isoutofdomain solver option in the DifferentialEquations package, yet it seems that an error within my derivative function (which requires positivity at the time step) is still being thrown. I thought isoutofdomain should prevent this from occurring.

Here’s my setup:

rhs(m, par, t) = get_aerosol_growth_3mom(m, par, t, v_up)
prob = ODEProblem(rhs, moments_S_init, tspan, ODE_parameters)
sol = solve(prob, reltol = tol, abstol = tol, isoutofdomain = (m,p,t) -> any(x->x<0, m))

The derivative function “get_aerosol_growth_3mom(m, par, t, v_up)” throws an error and cannot proceed if any element of m is negative. When I run the script, this error is still being thrown (rather than being averted via isoutofdomain). Inserting some print statements for the arguments to the derivative function:

time: 0.0
prognostic moments: [100.0, 10825.08853082, 1.4756905287634078e6, 0.02]

time: 3.5670907415192044e-6
prognostic moments: [100.0, 11008.214660723124, 1.4807861187372427e6, 0.02000000186343498]

time: 0.0
prognostic moments: [100.0, 10825.08853082, 1.4756905287634078e6, 0.02]

time: 1.7835453707596022e-6
prognostic moments: [100.0, 10916.651595771562, 1.4782383237503252e6, 0.020000000931717492]

time: 3.884165474098689e-5
prognostic moments: [100.0, -14624.568385517678, 1.560013303242658e6, 0.02000002029073535]

ERROR: LoadError: all moments need to be nonnegative.

So the second element is negative, despite the isoutofdomain option.

Can anyone see where I went wrong here? Thanks for the help.

2

It will get a negative and then reject that time step. You need to make your derivative function robust to this (for example, just give Infs), and then it’ll pull back when it realizes it goes negative.

3

Thanks Chris, I’ve thought about this and am not sure I understand. Could you clarify what you mean when you say “make your derivative function robust to this (for example, just give Infs)”. Should my derivative function interact somehow with the callback or with the time step? What do you mean by “just give Infs”?

4

On another related note, I’m having a similar issue when I try to use callbacks instead. As a toy example, I’m trying:

`# implement callbacks to halt the integration
  function condition(m,t,integrator) 
    println("Condition checked")
    t>=1e-9
  end

  function affect!(integrator) 
    terminate!(integrator)
  end

This exact callback works for the simple toy “bouncing ball” example from the DiffEq documentation, but it is never even called when I use it in my code. What’s going on?

5

Hi @emily_dj, could you provide a complete working (or non, as the case may be) example with the callbacks? We need to see how you are hooking them into the solver to answer your question.

6

Yes, of course.

  function condition(m,t,integrator) 
    println("Condition checked")
    t>=1e-9
  end

  function affect!(integrator) 
    terminate!(integrator)
  end

  cb=DiscreteCallback(condition, affect!)

  # set up ODE
  rhs(m, par, t) = get_aerosol_growth_3mom(m, par, t, v_up)

  # solve the ODE
  println("Solving ODE...")
  prob = ODEProblem(rhs, moments_S_init, tspan, ODE_parameters)
  alg = Tsit5()
  sol = solve(prob, alg, dt=1e-6, callback=cb)

If you’d like to see the full code and related modules, check it out here: https://github.com/CliMA/Cloudy.jl/blob/edj/aerosol-activation/aerosols-emily/runtests.jl

My current assumption is that the callback condition is not checked until after a potential time step is identified and tested. Could it be possible that all time steps attempted (as identified by print statements within my derivative function) are rejected before the solver crashes such that the callback is never triggered?

7

I grabbed that file but GammaPrimitiveParticleDistribution isn’t defined so I can’t run it. Still looking things over.

8

No interaction with the callback or anything, I mean that instead of

f(u,p,t) = sqrt(u)

you’d do

f(u,p,t) = sqrt(abs(u))

so it doesn’t error, and then you’d tell it to be positive. If it errors then we can’t recover.

You probably meant to use a ContinuousCallback with

  function condition(m,t,integrator) 
    println("Condition checked")
    t-1e-9
  end
9

Yes, things are still in progress on many ends… you’ll have to run it from your local Cloudy.jl directory, accessing the REPL as julia --project in order to have access to all the modules.

10

Thank you Chris, I see what you mean now and went ahead and did so with a try-catch (and setting the derivative to Inf in the catch case as you suggested).

I’ve tried both continuous and discrete callbacks. You’re right, the bouncing ball example uses a Continuous Callback. In neither case is the callback condition ever checked for my code though.

11

Interesting, I tried running the code but couldn’t because

ERROR: LoadError: LoadError: UndefVarError: MP_n_0 not defined
Stacktrace:
 [1] include(::Module, ::String) at .\Base.jl:377
 [2] include(::String) at C:\Users\accou\.julia\packages\ClimateMachine\CCjXq\src\ClimateMachine.jl:1
 [3] top-level scope at C:\Users\accou\.julia\packages\ClimateMachine\CCjXq\src\ClimateMachine.jl:15
 [4] include(::Module, ::String) at .\Base.jl:377
 [5] top-level scope at none:2
 [6] eval at .\boot.jl:331 [inlined]
 [7] eval(::Expr) at .\client.jl:449
 [8] top-level scope at .\none:3
in expression starting at C:\Users\accou\.julia\packages\ClimateMachine\CCjXq\src\Atmos\Parameterizations\CloudPhysics\Microphysics.jl:16
in expression starting at C:\Users\accou\.julia\packages\ClimateMachine\CCjXq\src\ClimateMachine.jl:15

Does the release not compile or is that just me? Anyways, I can help debug this from a distance. You’re saying it never hits the condition because it never does that print? How big is the ODE? Are the function calls completing? If you stick @show t in the ODE, what do you see?

12

I couldn’t say. I’ve had no issues with Cloudy.jl when using Julia 1.3.

Anyways, I can help debug this from a distance. You’re saying it never hits the condition because it never does that print?

Right.

How big is the ODE? Are the function calls completing? If you stick @show t in the ODE, what do you see?

It’s a system of 4 coupled ODEs. The function calls complete. With @show here is the progress: (the DiscreteCallback condition for simplicity here is t>=1e-9)

t = 0.0
mom_p = [340.0, 36805.301004788, 5.017347797795586e6, 0.05]
t = 1.61e-7
mom_p = [340.0, 36833.31686818363, 5.018875368095533e6, 0.05000000008410558]
t = 3.27e-7
mom_p = [340.0, 36862.59962920981, 5.020471989475298e6, 0.05000000017082313]
t = 9.0e-7
mom_p = [340.0, 36966.89436184011, 5.026032595097983e6, 0.050000000470155403]
t = 9.800255409045097e-7
mom_p = [340.0, 36982.12897954516, 5.026813866263827e6, 0.05000000051196034]
t = 1.0e-6
mom_p = [340.0, 36985.94483332327, 5.027009314388916e6, 0.050000000522394894]
t = 1.0e-6
mom_p = [340.0, 36985.6308188531, 5.027020303629402e6, 0.050000000522394894]
Condition checked
  2.760753 seconds (6.22 M allocations: 340.891 MiB, 7.25% gc time)

So the condition is only checked on the very last attempted time step (and hence the integration terminates). Does that mean that none of the previous attempted steps were actually accepted by the solver, and hence there was no reason to reference the callback in those steps? If so, this could ultimately explain my confusion.

13

If you use a DiscreteCallback, you have to use a tstop so it stops exactly at the time you wanted. Event Handling and Callback Functions · DifferentialEquations.jl is a good example of that.

14

I see that now, thank you Chris.