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