# Using DiffEq isoutofdomain to enforce positivity

**URL:** <https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443>\
**Category:** General Usage\
**Created:** [July 2, 2020, 6:15pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443 "2020-07-02T18:15:03Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![emily\_dj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emily_dj/32/16280_2.png) [@emily\_dj](https://discourse.julialang.org/u/emily_dj)\
**Post date:** [July 2, 2020, 6:15pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/1 "2020-07-02T18:15:03Z")

</div>

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:

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 2, 2020, 6:26pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/2 "2020-07-02T18:26:48Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![emily\_dj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emily_dj/32/16280_2.png) [@emily\_dj](https://discourse.julialang.org/u/emily_dj)\
**Post date:** [July 9, 2020, 9:20pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/3 "2020-07-09T21:20:29Z")

</div>

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”?

---

<div class="post-metadata">

**Author:** ![emily\_dj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emily_dj/32/16280_2.png) [@emily\_dj](https://discourse.julialang.org/u/emily_dj)\
**Post date:** [July 9, 2020, 10:37pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/4 "2020-07-09T22:37:36Z")

</div>

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

```julia
`# 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?

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [July 9, 2020, 11:42pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/5 "2020-07-09T23:42:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![emily\_dj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emily_dj/32/16280_2.png) [@emily\_dj](https://discourse.julialang.org/u/emily_dj)\
**Post date:** [July 10, 2020, 1:01am UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/6 "2020-07-10T01:01:20Z")

</div>

Yes, of course.

```julia
  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](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?

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [July 10, 2020, 1:12am UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/7 "2020-07-10T01:12:17Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 10, 2020, 1:51am UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/8 "2020-07-10T01:51:43Z")

</div>

> [@emily\_dj](#):
>
> 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”?

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

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

```

you’d do

```julia
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.

> [@emily\_dj](#):
>
> 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?

You probably meant to use a `ContinuousCallback` with

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

```

---

<div class="post-metadata">

**Author:** ![emily\_dj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emily_dj/32/16280_2.png) [@emily\_dj](https://discourse.julialang.org/u/emily_dj)\
**Post date:** [July 10, 2020, 2:49am UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/9 "2020-07-10T02:49:16Z")

</div>

> [@klaff](#):
>
> I grabbed that file but `GammaPrimitiveParticleDistribution` isn’t defined so I can’t run it. Still looking things over.

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.

---

<div class="post-metadata">

**Author:** ![emily\_dj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emily_dj/32/16280_2.png) [@emily\_dj](https://discourse.julialang.org/u/emily_dj)\
**Post date:** [July 10, 2020, 2:53am UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/10 "2020-07-10T02:53:43Z")

</div>

> [@ChrisRackauckas](#):
>
> No interaction with the callback or anything, I mean that instead of
> 
> ```julia
> f(u,p,t) = sqrt(u)
> 
> ```
> 
> you’d do
> 
> ```julia
> 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.

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).

> [@ChrisRackauckas](#):
>
> You probably meant to use a `ContinuousCallback` with

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 10, 2020, 4:05am UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/11 "2020-07-10T04:05:32Z")

</div>

> [@emily\_dj](#):
>
> 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.

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

```julia
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?

---

<div class="post-metadata">

**Author:** ![emily\_dj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emily_dj/32/16280_2.png) [@emily\_dj](https://discourse.julialang.org/u/emily_dj)\
**Post date:** [July 10, 2020, 11:08pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/12 "2020-07-10T23:08:25Z")

</div>

> [@ChrisRackauckas](#):
>
> Does the release not compile or is that just me?

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`)

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 10, 2020, 11:26pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/13 "2020-07-10T23:26:34Z")

</div>

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](https://diffeq.sciml.ai/latest/features/callback_functions/#Example-1:-Interventions-at-Preset-Times-1) is a good example of that.

---

<div class="post-metadata">

**Author:** ![emily\_dj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/emily_dj/32/16280_2.png) [@emily\_dj](https://discourse.julialang.org/u/emily_dj)\
**Post date:** [July 12, 2020, 6:24pm UTC](https://discourse.julialang.org/t/using-diffeq-isoutofdomain-to-enforce-positivity/42443/14 "2020-07-12T18:24:11Z")

</div>

I see that now, thank you Chris.