Enforce positivity in simple PDE solution

Specific Domains Modelling & Simulations
, ,
1

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.

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

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
2

Change your initial guess to positive values?

3

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:

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
test_values1600×1200 48.2 KB

4

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?

5

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.

6

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.