I am starting with the example at DiffEqOperators . I am interested in the steady state only. Let us have Neumann border condition (physically, a constant source of e.g. heat) on the left, and Dirichlet BC on the right. Here the code, which works as expected and plots a straight line:
using DiffEqOperators, OrdinaryDiffEq, SteadyStateDiffEq
using Plots
gr()
nknots = 64
h = 1.0/(nknots+1)
knots = range(h, step=h, length=nknots)
ord_deriv = 2
ord_approx = 2
Dc = 1.0
const Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)
const bc = RobinBC((0.0, 1.0, -1.0), (1.0, 0.0, 0.0), h, 1)
# Neumann BC left & Dirichlet BC right
u0 = fill(0.0, nknots)
function stp(u,p,t)
r = Δ*bc*(u.*Dc)
return r
end
alg=SSPRK932()
prob = SteadyStateProblem(stp, u0)
@time sol = solve(prob,DynamicSS(alg ;abstol=1e-4,reltol=1e-3), maxiters=1e5)
plot(sol)
Now, I am interested in a case where Dc is not a constant - D(x). If J is flow, u temperature (or concentration) J=D u' and du/dt = J' = 0 for steady state. J' = D u" + D' u' . Let’s try to implement it:
using DiffEqOperators, OrdinaryDiffEq, SteadyStateDiffEq
using Plots
gr()
nknots = 64
h = 1.0/(nknots+1)
knots = range(h, step=h, length=nknots)
u0 = fill(0.0, nknots)
D_0 = 1.0
D = fill(D_0, nknots)
D[30:34] .= 0.5 # "heat insulator" or diffusion barrier in the middle
const Δ2 = CenteredDifference(2, 2, h, nknots)
const Δ1 = CenteredDifference(1, 2, h, nknots)
const bc_c = RobinBC((0.0, 1.0, -1.0), (1.0, 0.0, 0.0), h, 1)
const bc_D = DirichletBC(D_0, D_0)
function stp(c,p,t)
r = Δ2*bc_c*(c.*D) + Δ1*bc_c*c*Δ1*bc_D*D
return r
end
alg=SSPRK932()
prob = SteadyStateProblem(stp, u0)
@time sol = solve(prob,DynamicSS(alg ;abstol=1e-4,reltol=1e-3), maxiters=1e6)
plot(sol)
What goes wrong here? - the result is not as expected:
Lets try an implementation without DiffEqOperators:
using OrdinaryDiffEq, SteadyStateDiffEq
using Plots
gr()
const nknots = 64
const h = 1.0/(nknots+1)
const Dc = 1.0
const D = fill(Dc, nknots+2)
D[30:34] .= 0.5
const ugp = zeros(nknots+2) # preallocation: extend the u array by boundary points left and right
u0 = fill(0.0, nknots)
function lcrviews(v)
ln = size(v)[1]
c = view(v, 2:ln-1)
l = view(v, 1:ln-2)
r = view(v, 3:ln)
return (l=l, c=c, r=r, ln=ln)
end
function x1der(x, h)
(x.r-x.l)/(2h)
end
function x2der(x, h)
(x.l+x.r-2x.c)./h^2
end
function stp(u,p,t)
uv = lcrviews(ugp)
Dv = lcrviews(D)
uv.c .= u
ugp[end] = 0.0 # Dirichlet
ugp[1] = ugp[2]+h # Neumann
r = x2der(uv, h).*Dv.c + x1der(uv, h).*x1der(Dv, h)
end
# alg = KenCarp4()
alg=SSPRK932()
prob = SteadyStateProblem(stp, u0)
@time sol = solve(prob,DynamicSS(alg ;abstol=1e-4,reltol=1e-3), maxiters=1e5)
plot(sol)
The result looks correct:
However trying to use KenCarp4 algorithm throws an error: LoadError: MethodError: no method matching...
So let’s rewrite the code in a more explicit way:
using OrdinaryDiffEq, SteadyStateDiffEq
using Plots
gr()
const nknots = 64
const h = 1.0/(nknots+1)
const Dc = 1.0
const D = fill(Dc, nknots)
D[30:34] .= 0.5
function stp!(du, u,p,t)
u_end = 0.0 # Dirichlet
u_0 = u[1]+h # Neumann
l = length(du)
for i in 2:l-1
du[i] = D[i]*(u[i-1]+u[i+1]-2u[i])/h^2 + ((u[i+1]-u[i-1])/2h)*((D[i+1]-D[i-1])/2h)
end
du[1] = D[1]*(u_0+u[2]-2u[1])/h^2
du[end]=D[end]*(u[end-1]+u_end-2u[end])/h^2
return nothing
end
alg = KenCarp4()
# alg=SSPRK932()
prob = SteadyStateProblem(stp!, u0)
@time sol = solve(prob,DynamicSS(alg ;abstol=1e-4,reltol=1e-3), maxiters=1e5)
plot(sol)
Now it works with KenCarp4 as well.
However in both cases it works only for D[30:34] down to about 0.25, then diverges.
So I have three questions:
- What was wrong with the first code (DiffEqOperators)?
- Why the second code doesn’t work with KenCarp4?
- What happens if the D-step is too large and what can be done to improve (e.g. slope instead of step?).
