Solving a 1D diffusion equation with linear and nonlinear source terms

Question

I would like to numerically solve the following equation: $$\frac{\partial \rho (z,t)}{\partial t} = B(N_D \rho (z,t) + \rho(z,t)^2) + D \frac{\partial^2 \rho (z,t)}{\partial z^2}$$ with the boundary conditions being: $$D \frac{\partial \rho}{\partial z}|_{z=0} = S_f \rho (z=0,t)$$
$$D \frac{\partial \rho}{\partial z}|_{z=L} = -S_b \rho (z=L,t)$$
Here, the constants $D$, $G$, $B$, $S_f$, $S_b$, $N_d$ are given, and are assumed not to vary with $z$, $t$, or $\rho$

In attempting to numerically solve this, I used central differencing to convert the boundary conditions to a set of ghost points in the $z$ mesh as follows: $$\rho (-\Delta z, t) = \rho(\Delta z, t) - \frac{S_f \Delta z \rho(0, t)}{D}$$ and
$$\rho (L+\Delta z, t)=\rho(L - \Delta z,t)-(\frac{S_b \Delta z \rho(L,t)}{D})$$.

I then implemented the following code in python:

import scipy as sp
from scipy.integrate import odeint
import numpy as np

def ODEEq(rho, t, D, B, Nd):
    return (D*np.gradient(np.gradient(rho))-B*(Nd*rho+(rho**2)))
def Solver(delz, delt, L, T, G0, D, B, Nd, Sf, Sb, unitMaker=1e-9):

    tspace = np.linspace(0,T,int(T/delt))

    xm1Ghost = G0[1] - (Sf*delz*G0[0]/D)
    xLp1Ghost = G0[-2] - (Sb*delz*G0[-1]/D)
    rho = G0
    np.insert(rho, 0, xm1Ghost)
    np.append(rho, xLp1Ghost)

    sol = odeint(ODEEq, rho, tspace, args=(D, B, Nd), full_output=1)
    return sol

L = 4000
timeRange = 2000
D = 1.639e18
B = 2e13
Nd = 0
delz = 0.1
delt = 0.1
G0 = 640*np.exp((-6.4e-2)*np.linspace(0,L,int(L/delz)))
Sf = 0.64
Sb = 0.64

sol = Solver(delz, delt, L, timeRange, G0, D, B, Nd, Sf, Sb)

The initial conditions (G0) are visualized as:

However, there seems to be some issue with this implementation (aside from the fact that it is extremely slow).

I'm hoping to get some help in resolving the issues with this approach, and any ideas on how to make it faster.

Answer 1

I'm not gonna edit your code to make it work but based on your code, I think you have some misunderstanding about finite difference method. Your PDE equation is:

$$\frac{\partial \rho}{\partial t} = D \frac{\partial^{2} \rho}{\partial z^{2}} + G - B(N_{D} \rho + \rho^{2})$$

Your boundary conditions are:

$$-D \frac{\partial \rho}{\partial z} \Bigg |_{z=0} = -S_{f} \rho |_{z=0}$$

and

$$-D \frac{\partial \rho}{\partial z} \Bigg |_{z=L} = S_{b} \rho |_{z=L}$$

Remember you can't or at least you can't easily implement a central difference at the boundaries. I believe it's better to use just first order at the boundaries for the moment. Now, your main PDE could be discretized by using first order difference in time and second order difference in spatial as:

$$\frac{\rho^{t+\Delta t}_{z}-\rho^{t}_{z}}{\Delta t} = D \frac{\rho^{t}_{z+\Delta z}+\rho^{t}_{z-\Delta z}-2\rho^{t}_{z}}{\Delta z^{2}} + G - B (N_{D} \rho^{t}_{z}+(\rho^{t}_{z})^{2})$$

I did not understand what you want to use as initial condition but let's say it's whatever that is shown as:

$$\rho^{0}_{z} = f(z)$$

You use that discretized equation by using above initial condition to update the points in your mid domain. But for boundary conditions:

$$-D \frac{\rho^{t}_{\Delta z}-\rho^{t}_{0}}{\Delta z} = -S_{f} \rho^{t}_{0}$$

or:

$$\rho^{t}_{0} = (S_{f}-\frac{D}{\Delta z})^{-1} \frac{D}{\Delta z} \rho^{t}_{\Delta z}$$

Note that $\rho^{t}_{\Delta z}$ is already known from your discretized PDE equation.

For the second boundary condition you have similar procedure:

$$\rho_{L}^{t} = (S_{b}+\frac{D}{\Delta z})^{-1} \frac{D}{\Delta z} \rho^{t}_{L-\Delta z}$$

Again $\rho^{t}_{L-\Delta z}$ is already known from your discretized PDE equation. I think you could easily follow this procedure and implement it to make it work by using Python.

Answer 2

Let's try some alternative avoiding the forward and backward differentiation of the np.gradient function. You need to apply the boundary conditions in each computation of the time derivative, as that determines the derivative at the boundary points. Thus you have to use there either $$ \rho_{zz}(x_0,t)\approx 2\rho[x_{-1},x_0,x_1]=\frac{ρ(x_{-1},t)-2ρ(x_0,t)+ρ(x_1,t)}{\Delta z^2} $$ or without the ghost point $$ \rho_{zz}(x_0,t)\approx 2\rho[x_0,x_0,x_1] =2\frac{ρ(x_1,t)-ρ(x_0,t)-Δz ρ_z(x_0,t)}{\Delta z^2} $$ which both in the end give the same numerical formula.

def ODEEq(rho, t, D, B, Nd):
    # ghost cells implementing the boundary conditions
    rho0m1 = rho[1] - (Sf*delz*rho[0]/D)
    rhoLp1 = rho[-2] - (Sb*delz*rho[-1]/D)
    # construct Laplace/second derivative
    D2rho = -2*rho; 
    D2rho[1:-1] += rho[:-2]+rho[2:];
    D2rho[0] += rho0m1+rho[1];
    D2rho[-1] += rho[-2]+rhoLp1;
    return  D*D2rho/delz**2+G-B*(Nd*rho+rho**2);

Then you need to be consistent in the construction of the space discretization. The interval is $[0,L]$, there is a number of segments in that interval, this segmentation then gives the space step $Δz$.

Nz=int(L/delz+0.1);
z = np.linspace(0,L,Nz+1);
delz = z[1]-z[0];

G0 = 640*np.exp(-6.4e-2 * z)

Note that this initial condition does not satisfy the boundary conditions. This means that there will be some dynamical artifacts at the boundary that depend on the space step size.

Changing the constants to get a result that is computable and actually visually representable,

L = 40
delz = 0.1
D = 1.639e2
B = 2e1
G = B
Nd = 0

Sf = 0.64
Sb = 0.64

t = np.arange(0,1.1,0.025)
rho = odeint(lambda u,t: ODEEq(u,t, D, B, Nd), rho0, t, atol=1e-4, rtol=1e-9)

gives a solution that rapidly (depending on your view of scale) converges to a function that is close to the stable equilibrium $1$ over the most part and has the required logarithmic slopes at the end points.