Heat / diffusion equation with co-ordinate dependent diffusion coefficient

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:
2020-06-14 15_54_27-HeatEqEx1D-05.jl — C___MyOwnDocs_SoftwareDevelopment_MebraneFlowPDEs_MebraneFlow

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:
2020-06-14 16_15_53-Plots — C___MyOwnDocs_SoftwareDevelopment_MebraneFlowPDEs_MebraneFlowPDEs.jl — A

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

Could you open an issue at https://github.com/SciML/DiffEqOperators.jl ? This looks like it could be a bug.

This is because the cache that you had setup was not compatible with autodiff because there wasn’t enough storage space for numbers larger than Float64. So either do KenCarp4(autodiff=false) or

https://diffeq.sciml.ai/latest/basics/faq/#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods-1

That describes what to do.

What do you mean?

Chris, thank you!

For the first question: the bug was in my code:

function stp(c,p,t)
   r = Δ2*bc_c*(c.*D) + Δ1*bc_c*c*Δ1*bc_D*D
   return r
end

in the first term taking (D*c)" is wrong

function stp(c,p,t)
   D.*(Δ2*bc_c*c) + Δ1*bc_c*c*Δ1*bc_D*D
end

works as expected.

Another question: Is it the optimal way to use DiffEqOperators?

Let’s set D[30:34] .= 0.1

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase C:\Users\eben\.julia\packages\DiffEqBase\KnYSY\src\integrator_interface.jl:329

2020-06-14 19_04_05-Plots — C___MyOwnDocs_SoftwareDevelopment_MebraneFlowPDEs_MebraneFlowPDEs.jl — A

That’s not a stable way to discretize that portion of the equation.

The last plot was actually produced by the code using Δ1*bc_c*c*Δ1*bc_D*D , but that’s probably equivalent - all three codes share the problem

What would you suggest?

You do forward derivative in one of them and backward derivative in the other.

Like (u[i+1]-u[i])*(D[i]-D[i-1])/h^2 ? - That produced stable but wrong results with different slopes on the left and right side of the curve.

However (u[i+1]-u[i])*(D[i+1]-D[i])/h^2 or (u[i]-u[i-1])*(D[i]-D[i-1])/h^2 looked OK at any D ratios, so I took the mean of both of them:
((u[i+1]-u[i])*(D[i+1]-D[i])+(u[i]-u[i-1])*(D[i]-D[i-1]))/(2h^2)
and now it integrates stably and correctly down to D-ratios of at least 1 to 0.001

Thank you Chris