# Enforce positivity in simple PDE solution

**URL:** <https://discourse.julialang.org/t/enforce-positivity-in-simple-pde-solution/121738>\
**Category:** Modelling & Simulations\
**Tags:** package, differentialequation, nonlinearsolve\
**Created:** [October 25, 2024, 9:26am UTC](https://discourse.julialang.org/t/enforce-positivity-in-simple-pde-solution/121738 "2024-10-25T09:26:16Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![ify](https://avatars.discourse-cdn.com/v4/letter/i/a88e4f/32.png) [@ify](https://discourse.julialang.org/u/ify)\
**Post date:** [October 25, 2024, 9:26am UTC](https://discourse.julialang.org/t/enforce-positivity-in-simple-pde-solution/121738/1 "2024-10-25T09:26:16Z")

</div>

Hi all,  
Pretty new to julia. I am trying to solve a simplified PDE that I normally would do in R. In R, the numerical steady state solution can be constrained to obtain solution concentration. However, when I translate this code to julia and solve numerically with differentialEquations.jl and NonlinearSolve (NewtonRaphson()), the solutions are always negative (which it shouldn’t be as the concentrations are positive.

Can anyone help with this.

```julia
thick = 0.05
Intdepth = collect(0:thick:10)
Nint = length(Intdepth)
Depth = 0.5*(Intdepth[1:(Nint-1)] + Intdepth[2:Nint])
N = length(Depth)

porosity = 0.8 # -
minrate = 1000 # mmol/m3/d - oxygen consumption rate
O2BW = 300 # mmol/m3 - bottom water oxygen concentration
ks = 1 # mmol O2/m3 - Half saturation constant
DispO2 = 1.289186
w = 0.1/1000/365

function O2model(du, conc, p, t)
    O2 = @view conc[:, 1]
    du = @view du[:, 1]

    O2tran = (-DispO2 .* porosity .* diff([O2BW ; O2 ; O2[N] ])./thick) + (w .* porosity .* [O2BW; O2]) 

    # Respiration
    O2cons = minrate.*(O2./(O2.+ks))    
    
    dO2 = -diff(O2tran)./thick./porosity - O2cons

    # du[:,1] = max.(0.0, dO2)
    du[:,1] = dO2
end 

using DifferentialEquations, NonlinearSolve
using CairoMakie

conc = ones(N, 1)
prob = SteadyStateProblem(O2model, conc)
#sol = solve(prob, NewtonRaphson(), reltol = 1e-8, abstol = 1e-8)
#Test with isoutofdomain option
sol = solve(prob, NewtonRaphson(), reltol = 1e-4, abstol = 1e-4, isoutofdomain=(conc,parms,t) -> any(x -> x < 0, conc))

```

```julia
retcode: Success
u: 200×1 Matrix{Float64}:
    -87.89480788353349
   -473.8280950051232
   -857.818076779688
  -1239.866591215867
  -1619.9743362518211
  -1998.141679374796
      ⋮
 -38650.39626156616
 -38660.09255274791
 -38667.84958571944
 -38673.667360470194
 -38677.54587698711
 -38679.48513525461

```

---

<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:** [October 28, 2024, 1:38am UTC](https://discourse.julialang.org/t/enforce-positivity-in-simple-pde-solution/121738/2 "2024-10-28T01:38:40Z")

</div>

Change your initial guess to positive values?

---

<div class="post-metadata">

**Author:** ![Vitor\_Patricio\_Canta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vitor_patricio_canta/32/50861_2.png) [@Vitor\_Patricio\_Canta](https://discourse.julialang.org/u/Vitor_Patricio_Canta)\
**Post date:** [October 28, 2024, 9:14am UTC](https://discourse.julialang.org/t/enforce-positivity-in-simple-pde-solution/121738/3 "2024-10-28T09:14:13Z")

</div>

Hi,

I am also experimenting with reactive transport models. I would say your model looks fine. Can you please check your parameter values and see if they are feasible?  
For example playing around with your max consumption rate and and your velocity I get a plot like this:

```julia
porosity = 0.8 # -
minrate = 5 # mmol/m3/d - oxygen consumption rate
O2BW = 300 # mmol/m3 - bottom water oxygen concentration
ks = 1 # mmol O2/m3 - Half saturation constant
DispO2 = 1.289186
w = 20/365

```

 ![test_values](https://global.discourse-cdn.com/julialang/original/3X/f/3/f372fd0bdaae8b3b74554370062dbf0ac3d082c0.png)

---

<div class="post-metadata">

**Author:** ![ify](https://avatars.discourse-cdn.com/v4/letter/i/a88e4f/32.png) [@ify](https://discourse.julialang.org/u/ify)\
**Post date:** [October 28, 2024, 9:54am UTC](https://discourse.julialang.org/t/enforce-positivity-in-simple-pde-solution/121738/4 "2024-10-28T09:54:06Z")

</div>

Hi @ChrisRackauckas I tried changing the initial guess value and nothing change (The solution is still negative). Moreover, `conc(N, 1)` or `fill(10, N)` will always produce a positive value right?

---

<div class="post-metadata">

**Author:** ![ify](https://avatars.discourse-cdn.com/v4/letter/i/a88e4f/32.png) [@ify](https://discourse.julialang.org/u/ify)\
**Post date:** [October 28, 2024, 10:01am UTC](https://discourse.julialang.org/t/enforce-positivity-in-simple-pde-solution/121738/5 "2024-10-28T10:01:49Z")

</div>

Hi @Vitor_Patricio_Canta Yes, you are right. Playing with the oxygen consumption rate and/or the velocity makes the solution positivity. What I don’t understand is why the R version for the same previous parameter values give positive result but the julia doesn’t.

By the way, great to know other folks within the julia community working on reactive transport models.

---

<div class="post-metadata">

**Author:** ![Vitor\_Patricio\_Canta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vitor_patricio_canta/32/50861_2.png) [@Vitor\_Patricio\_Canta](https://discourse.julialang.org/u/Vitor_Patricio_Canta)\
**Post date:** [October 28, 2024, 10:35am UTC](https://discourse.julialang.org/t/enforce-positivity-in-simple-pde-solution/121738/6 "2024-10-28T10:35:57Z")

</div>

Same here @ify ! Its hard to tell you what could have gotten wrong. Could be a typo or a wrong multiplication somewhere. Your dispersion coefficient looks high also, but maybe it makes sense with your choice of parameters.