I have trying for quite some time to write a finite difference solver for the following advection-diffusion-reaction differential equation:
$$ \frac{\partial C}{\partial x} = \frac{1}{u(z) + u_e(x,z)}\left( \frac{\partial}{\partial z} \left( K(z)\frac{\partial C}{\partial z} \right) - w(z) \frac{\partial C}{\partial z} - \lambda C \right) + S(x,z) $$ one the domain $(x,z) \in [0,X]\times [0,H]$, where \begin{align*} K(z) = K_r &\left( \frac{z}{z_r} \right)^{1-\alpha_p}\\ u(z) = u_r &\left( \frac{z}{z_r} \right)^{\alpha_p}\\ w(z) = \frac{az}{\alpha_p+1}&\left( \frac{z}{z_r}\right)^{\alpha_p}\\ u_e(x,z) = -ax&\left( \frac{z}{z_r} \right)^{\alpha_p} \end{align*} Additionally, $\lambda \sim 10^{-5}$, $a\sim 10^{-3}$, $z_r,u_r,K_r \sim 1$, and $\alpha_p \in (0,1)$. My first step in trying to create this solver was to rewrite the PDE in a conservative form, I can first say that $u_0(x)/s(z) = u_e+u$, where $u_0(x) = (u_r-ax)$ and $s(z) = \left( \frac{z}{z_r} \right)^{-\alpha_p}$. From there I can write $$ \frac{\partial C}{\partial x} = \frac{s(z)}{u_0(x)}\left( \frac{\partial}{\partial z} \left( K(z)\frac{\partial C}{\partial z} \right) - w(z) \frac{\partial C}{\partial z} - \lambda C \right) + S(x,z) $$ Furthermore, $$ s(z)\left( K(z) \frac{\partial C}{\partial z} \right) = \frac{\partial}{\partial z} \left(s(z) K(z) \frac{\partial C}{\partial z} \right) - s'(z)\cdot K(z) \frac{\partial C}{\partial z} $$ We define $\mathcal{D}(z) = s(z)\cdot K(z)$. Next we can write $$ - \left( s'(z) \cdot K(z) + s(z) w(z) \right) \frac{\partial C}{\partial z} = -\frac{\partial}{\partial z} \left( (s'(z) \cdot K(z) + s(z) w(z)) C \right) + \frac{\partial}{\partial z} \left( s'(z) \cdot K(z) + s(z) w(z) \right) \cdot C $$ where we define $\mathcal{A}(z) = -(s'(z) \cdot K(z) + s(z) w(z))$, and $\mathcal{R}(z) = \frac{\partial}{\partial z} \left( s'(z) \cdot K(z) + s(z) w(z) \right) - s(z)\lambda$. Allowing for the equation to become $$ \frac{\partial C}{\partial x} = \frac{1}{u_0} \left( \frac{\partial}{\partial z} \left(\mathcal{D} \frac{\partial C}{\partial z} + \mathcal{A} C \right) + \mathcal{R}C \right) + S $$ There are only three boundary conditions, one is a Dirichlet, one is Robin, and one is Neumann. \begin{align} C(0,z) = C_0\\ \lim_{z\to0^+} K(z) \frac{\partial C}{\partial z} = v_d C(x,0) - Q(x)\\ K(H) \cdot \frac{\partial C}{\partial z}\Big|_{z=H} = 0 \end{align} where $Q(x)$ is a boundary source function and $v_d \sim 10^{-3}$. Due to the fact that the second condition is very difficult to deal with in this context it is replaced with the condition \begin{gather} K(z^*) \frac{\partial C}{\partial z}\Big|_{z=z^*} = v_d C(x,z^*) - Q(x) \end{gather} where $z^*>0$ and is reasonably close to $0$. Moreover, the domain of the problem then becomes $(x,z) \in [0,X] \times [z^*,H]$. Since the PDE lacks a condition at the $x=X$ boundary, the finite difference solver must be "time-stepped" along $x$. To make the scheme even further conservative, the boundary conditions are written in a conservative form as \begin{align*} \left[\frac{\partial}{\partial z} (K\cdot C) - K'\cdot C\right]_{z=z^*} = v_d C(x,z^*) - Q(x)\\ \left[\frac{\partial}{\partial z} (K\cdot C) - K'\cdot C\right]_{z=H} = 0 \end{align*} The next step I took was to perform nondimensionalization, where $x=x_s\widetilde{x}$, $z = z_s\widetilde{z}+z^*$, and $C = C_s \widetilde{C}$. The PDE can be written $$ \frac{C_s}{x_s} \frac{\partial \widetilde{C}}{\partial \widetilde{x}} = \frac{1}{u_0} \left( \frac{1}{z_s}\frac{\partial}{\partial \widetilde{z}} \left( \mathcal{D} \frac{C_s}{z_s} \frac{\partial \widetilde{C}}{\partial \widetilde{z}} + \mathcal{A}C_s \widetilde{C} \right) + \mathcal{R} C_s \widetilde{C} \right) + S $$ This can simplify to $$ \frac{\partial \widetilde{C}}{\partial \widetilde{x}} = \frac{x_s}{z_s^2 u_0} \left( \frac{\partial}{\partial \widetilde{z}} \left( \mathcal{D} \frac{\partial \widetilde{C}}{\partial \widetilde{z}} + z_s\mathcal{A} \widetilde{C} \right) + z_s^2\mathcal{R} \widetilde{C} \right) + \frac{x_sS}{C_s} $$ With this, I define the functions to be the following \begin{align*} \widetilde{u_0}(\widetilde{x}) = \frac{z_s^2 u_0(x_s\widetilde{x})}{x_s}\\ \widetilde{\mathcal{D}}(\widetilde{z}) = \mathcal{D}(z_s\widetilde{z}+z^*)\\ \widetilde{\mathcal{A}}(\widetilde{z}) = z_s\mathcal{A}(z_s\widetilde{z}+z^*)\\ \widetilde{\mathcal{R}}(\widetilde{z}) = z_s^2\mathcal{R}(z_s\widetilde{z}+z^*)\\ \widetilde{\mathcal{S}}(\widetilde{x},\widetilde{z}) = \frac{x_sS(x_s\widetilde{x},z_s\widetilde{z}+z^*)}{C_s} \end{align*} The shifting now allows for the domain to become $(\widetilde{x},\widetilde{z})\in[0,X/x_s]\times [0,(H-z^*)/z_s]$. The boundary conditions need to be nondimensionalized too, so \begin{align*} \frac{\partial}{\partial \widetilde{z}}\left( \widetilde{K} \widetilde{C} \right)\Big|_{\widetilde{z} = 0} - z_s \left[\widetilde{K}' \widetilde{C}\right]_{\widetilde{z}=0} = z_sv_d \widetilde{C}(\widetilde{x},0) - \widetilde{Q}(\widetilde{x})\\ \frac{\partial}{\partial \widetilde{z}}\left( \widetilde{K} \widetilde{C} \right)\Big|_{\widetilde{z} = 0} - z_s \left[\widetilde{K}' \widetilde{C}\right]_{\widetilde{z}=0} \end{align*} Where $\widetilde{K}(\widetilde{z}) = K(z_s\widetilde{z}+z^*)$, $\widetilde{K}'(\widetilde{z}) = K'(z_s\widetilde{z}+z^*)$, and $\widetilde{Q}(\widetilde{x}) = z_sQ(x_s\widetilde{x})/C_s$. For the scheme's discretization, I chose to do a Crank-Nicolson scheme. Giving the form $$ 2\frac{\widetilde{C}_{i+1,j} - \widetilde{C}_{ij}}{\Delta \widetilde{x}} = \frac{1}{[u_0]_{i+1}} \left[ \frac{\partial}{\partial \widetilde{z}} \left( \widetilde{\mathcal{D}} \frac{\partial \widetilde{C}}{\partial\widetilde{z}} + \widetilde{\mathcal{A}} \widetilde{C} \right) + \widetilde{\mathcal{R}} \widetilde{C} \right]_{i+1,j} + \frac{1}{[u_0]_{i}} \left[ \frac{\partial}{\partial \widetilde{z}} \left( \widetilde{\mathcal{D}} \frac{\partial \widetilde{C}}{\partial\widetilde{z}} + \widetilde{\mathcal{A}} \widetilde{C} \right) + \widetilde{\mathcal{R}} \widetilde{C} \right]_{ij} + \widetilde{S}_{i+1,j} + \widetilde{S}_{ij} $$ Due to the length of this post, I will not show the scheme's derivation, but I have the following $$ E_{ij}\widetilde{C}_{i+1,j+1} + F_{ij}\widetilde{C}_{i+1,j} + G_{ij}\widetilde{C}_{i+1,j-1} = A_{ij}\widetilde{C}_{i,j+1} + B_{ij}\widetilde{C}_{ij} + D_{ij}\widetilde{C}_{i,j-1} + \mathcal{H}_{ij} $$ where \begin{align*} E_{ij} = -\Delta \widetilde{x}[u_0]_i\left(2\widetilde{\mathcal{D}}_{j+1/2} + \Delta \widetilde{z} \widetilde{\mathcal{A}}_{j+1/2} + \Delta \widetilde{z}^2 \widetilde{\mathcal{R}}_j\right)\\ F_{ij} = [u_0]_i\left(4 \Delta \widetilde{z}^2 [u_0]_{i+1} - \Delta\widetilde{x}(\Delta \widetilde{z}(\widetilde{\mathcal{A}}_{j+1/2} - \widetilde{\mathcal{A}}_{j-1/2}) - 2(\widetilde{\mathcal{D}}_{j+1/2} + \widetilde{\mathcal{D}}_{j-1/2}))\right)\\ G_{ij} = -\Delta \widetilde{x}[u_0]_i\left(2\widetilde{\mathcal{D}}_{j-1/2} - \Delta \widetilde{z} \widetilde{\mathcal{A}}_{j-1/2} + \Delta \widetilde{z}^2 \widetilde{\mathcal{R}}_j\right)\\ A_{ij} = \Delta \widetilde{x}[u_0]_{i+1}\left(2\widetilde{\mathcal{D}}_{j+1/2} + \Delta \widetilde{z} \widetilde{\mathcal{A}}_{j+1/2} + \Delta \widetilde{z}^2 \widetilde{\mathcal{R}}_j\right)\\ B_{ij} = [u_0]_{i+1}\left(4 \Delta \widetilde{z}^2 [u_0]_{i} + \Delta\widetilde{x}(\Delta \widetilde{z}(\widetilde{\mathcal{A}}_{j+1/2} - \widetilde{\mathcal{A}}_{j-1/2}) - 2(\widetilde{\mathcal{D}}_{j+1/2} + \widetilde{\mathcal{D}}_{j-1/2}))\right)\\ D_{ij} = \Delta \widetilde{x}[u_0]_{i+1}\left(2\widetilde{\mathcal{D}}_{j-1/2} - \Delta \widetilde{z} \widetilde{\mathcal{A}}_{j-1/2} + \Delta \widetilde{z}^2 \widetilde{\mathcal{R}}_j\right)\\ \mathcal{H}_{ij} = 2[u_0]_{i+1}[u_0]_i(\widetilde{S}_{ij} + \widetilde{S}_{i+1,j} ) \end{align*} The discretization of the boundary conditions yields the following terms referencing ghost nodes \begin{align*} \widetilde{C}_{i,-1} = T\widetilde{C}_{i,1} + [S_b]_i\\ \widetilde{C}_{i,n+1} = V \widetilde{C}_{i,n-1} \end{align*} with \begin{align*} T = \frac{\widetilde{K}_1 - z_s \Delta \widetilde{z} \left( \widetilde{K}'_0 +v_d\right)}{\widetilde{K}_{-1} + z_s \Delta \widetilde{z} \left( \widetilde{K}'_0 +v_d\right)}\\ [S_b]_i = \frac{2\Delta \widetilde{z}\widetilde{Q}_i}{\widetilde{K}_{-1} + z_s \Delta \widetilde{z} \left( \widetilde{K}'_0 +v_d\right)}\\ V = \frac{\widetilde{K}_{n-1} + z_s \Delta \widetilde{z} \left( \widetilde{K}'_n +v_d\right)}{\widetilde{K}_{n+1} - z_s \Delta \widetilde{z} \left( \widetilde{K}'_n +v_d\right)} \end{align*} Since $\widetilde{K}$ is ill-defined for values $\widetilde{z}<-z^*/z_s$, it re-expressed as $$ \widetilde{K}(\widetilde{z}) = K_r\cdot \textrm{sign}{(z_s\widetilde{z}+z^*)} \left| \frac{z_s\widetilde{z}+z^*}{z_r} \right|^{1-\alpha_p} $$ Furthermore, the the functions $\widetilde{\mathcal{D}}$, $\widetilde{\mathcal{A}}$, and $\widetilde{\mathcal{R}}$ are problematic for values $\widetilde{z}<0$, due to their asymptotic behavior (barring $\widetilde{D}$ for $\alpha_p\geq 0.5$). To circumvent this issue, linear extrapolation is performed using the output and derivatives of the functions at $\widetilde{z}=z^*/z_s$. Doing this allows for the evaluation of the three functions at the staggered node $j=-1/2$.
At each time-step we have the matrix equation \begin{multline*} \begin{bmatrix} F_{i,0} & E_{i,0} + TG_{i,0} & 0 & \cdots\\ G_{i,1} & F_{i,1} & E_{i,1} & 0 & \cdots\\ & \ddots & \ddots& \ddots\\ & & G_{i,n-1} & F_{i,n-1} & E_{i,n-1}\\ & & & G_{i,n} + VE_{i,n} & F_{i,n} \end{bmatrix} \begin{bmatrix} \widetilde{C}_{i+1,0}\\ \widetilde{C}_{i+1,1}\\ \widetilde{C}_{i+1,2}\\ \vdots\\ \widetilde{C}_{i+1,n}\\ \end{bmatrix} \\= \begin{bmatrix} B_{i,0} & A_{i,0} + TD_{i,0} & 0 & \cdots\\ D_{i,1} & B_{i,1} & A_{i,1} & 0 & \cdots\\ & \ddots & \ddots& \ddots\\ & & D_{i,n-1} & B_{i,n-1} & A_{i,n-1}\\ & & & D_{i,n} + VA_{i,n} & B_{i,n} \end{bmatrix} \begin{bmatrix} \widetilde{C}_{i,0}\\ \widetilde{C}_{i,1}\\ \widetilde{C}_{i,2}\\ \vdots\\ \widetilde{C}_{i,n}\\ \end{bmatrix} + \begin{bmatrix} \mathcal{H}_{i,0} + D_{i,0}[S_b]_i - G_{i,0}[S_b]_{i+1}\\ \mathcal{H}_{i,1}\\ \mathcal{H}_{i,2}\\ \vdots\\ \mathcal{H}_{i,n}\\ \end{bmatrix} \end{multline*} I am implementing this scheme using MATLAB and despite this implementation, I continue to get terrible relative errors when testing the method using the method of manufactured solutions, for a function $$ C(x,z) = C_0 + x\cdot (H-z)^3 $$ The convergence should be both order 2 in $x$ and $z$, but I can't even get the scheme to consistently converge. I've picked the following values: $x_s=10^6,z_s=10^6,C_s=10^3,\lambda=4\cdot 10^{-5},a=2\cdot 10^{-3},v_d=6\cdot 10^{-3},z_r=10,u_r=2,H=200, X=1000$. I've played around with values of $\alpha_p$ in it's range and $0.5<z^*<20$. I've been working on this problem for quite a while and have tried different approaches to the scheme, all with no luck. If you happen to see anything off in my derivation, please point it out to me, or if you have any insight as to why this problem would be so unruly, please let me know. I'm fairly certain that I've shown my derivation used correctly, but there is always a possibility that I've made some typographical error.
Here's my MATLAB code for the scheme:
function Ctilde = myCDR(X,H,m,n,xs,zs,Cs,vd,alpha_p,Kr,C0,Q,zanchor,Stilde)
K = @(zvar) K_func(zvar,alpha_p,Kr);
dK = @(zvar) dK_func(zvar,alpha_p,Kr);
Qtilde = @(xvar) zs*Q(xs*xvar)/Cs;
dKtilde = @(zvar) dK(zs*zvar+zanchor);
u0tilde = @(xvar) zs^2/xs*u0_func(xs*xvar);
Ktilde = @(zvar) sign(zs*zvar+zanchor).*K(abs(zs*zvar+zanchor));
sDtilde = @(zvar) sD(zs*zvar+zanchor);
sAtilde = @(zvar) zs*sA(zs*zvar+zanchor);
sRtilde = @(zvar) zs^2*sR(zs*zvar+zanchor);
sDbelow = @(zvar) sD_below(zs*zvar+zanchor);
sAbelow = @(zvar) zs*sA_below(zs*zvar+zanchor);
Htilde = (H-zanchor)/zs;
dztilde = Htilde/n;
ztilde = 0:dztilde:Htilde;
Xtilde = X/xs;
dxtilde = Xtilde/m;
xtilde = 0:dxtilde:Xtilde;
gammaTerm = zeros(1,n+1);
betaTerm = zeros(1,n+1);
alphaTerm = 2*sDtilde(ztilde+0.5*dztilde) ...
+ dztilde*sAtilde(ztilde+0.5*dztilde) ...
+ dztilde^2*sRtilde(ztilde);
betaTerm(2:end) = dztilde*(sAtilde(ztilde(2:end)+0.5*dztilde) - sAtilde(ztilde(2:end)-0.5*dztilde)) ...
- 2*(sDtilde(ztilde(2:end)+0.5*dztilde) + sDtilde(ztilde(2:end)-0.5*dztilde));
betaTerm(1) = dztilde*(sAtilde(0.5*dztilde) - sAbelow(-0.5*dztilde)) ...
- 2*(sDtilde(0.5*dztilde) + sDbelow(-0.5*dztilde));
gammaTerm(2:end) = 2*sDtilde(ztilde(2:end)-0.5*dztilde) ...
- dztilde*sAtilde(ztilde(2:end)-0.5*dztilde) ...
+ dztilde^2*sRtilde(ztilde(2:end));
gammaTerm(1) = 2*sDbelow(-0.5*dztilde) ...
- dztilde*sAbelow(-0.5*dztilde) ...
+ dztilde^2*sRtilde(0);
T = (Ktilde(dztilde) - zs*dztilde*(dKtilde(0) + vd)) ...
/(Ktilde(-dztilde) + zs*dztilde*(dKtilde(0) + vd));
Sb = 2*dztilde*Qtilde(xtilde) ...
/(Ktilde(-dztilde) + zs*dztilde*(dKtilde(0) + vd));
V = (Ktilde(ztilde(n-1)) + zs*dztilde*dKtilde(ztilde(n)))...
/(Ktilde(ztilde(n)+dztilde) - zs*dztilde*dKtilde(ztilde(n)));
Ctilde = zeros(n+1,m+1);
Ctilde(:,1) = C0/Cs;
for i=1:m
A = dxtilde*u0tilde(xtilde(i+1))*alphaTerm;
E = -dxtilde*u0tilde( xtilde(i) )*alphaTerm;
B = u0tilde(xtilde(i+1))*(4*dztilde^2*u0tilde( xtilde(i) ) + dxtilde*betaTerm);
F = u0tilde( xtilde(i) )*(4*dztilde^2*u0tilde(xtilde(i+1)) - dxtilde*betaTerm);
D = dxtilde*u0tilde(xtilde(i+1))*gammaTerm;
G = -dxtilde*u0tilde( xtilde(i) )*gammaTerm;
Adiag = [nan,A(1) + T*D(1),A(2:end-1)];
Ediag = [nan,E(1) + T*G(1),E(2:end-1)];
Bdiag = B;
Fdiag = F;
Ddiag = [D(2:end-1),D(end) + V*A(end),nan];
Gdiag = [G(2:end-1),G(end) + V*E(end),nan];
sH = 2*dztilde^2*dxtilde*u0tilde(xtilde(i))*u0tilde(xtilde(i+1)) ...
*(Stilde(xtilde(i),ztilde)+Stilde(xtilde(i+1),ztilde));
sH(1) = sH(1) + D(1)*Sb(i)-G(1)*Sb(i+1);
M = spdiags([Ddiag.', Bdiag.', Adiag.'],-1:1,n+1,n+1);
N = spdiags([Gdiag.', Fdiag.', Ediag.'],-1:1,n+1,n+1);
if n==20 && i==1
full(M)
full(N)
full(sH.')
end
Ctilde(:,i+1) = N\(M*Ctilde(:,i) + sH.');
end
end
Here's the code I'm using to run it:
clear variables
n = 20*2.^(0:7);
m = 5*2.^(0:7);
%%% Nondimensionalization Factors %%%
Cs = 1e6;
xs = 1e6;
zs = 1e6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Problem's Constants %%%
X = 10;
H = 200;
vd = 6e-3;
alpha_p = 0.61;
L = 1/0.0385;
zanchor = 100;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
Kr = Kr_constant(alpha_p,L);
generate_conservative_form_functions(alpha_p,L,zanchor);
K = @(zvar) K_func(zvar,alpha_p,Kr);
schemeFunc = @myCDR;
myError = zeros(1,length(n));
C0 = 1;
Ctrue_func = @(xvar,zvar) C0 ...
+ xvar.*(H-zvar).^3;
syms xx zz
K0 = K_func(zanchor,alpha_p,Kr);
Cweight = sqrt(double(int(int((Ctrue_func(xs*xx,zs*zz+zanchor)/Cs).^2,zz,0,(H-zanchor)/zs),xx,0,X/xs)));
C_check_BC = K0*diff(Ctrue_func(xx,zz),zz);
C_check_BC = matlabFunction(C_check_BC);
if any(Ctrue_func(0,0:H) ~= C0)
Ctrue_func(0,0:H/10:H)
error("Recheck your true solution to make sure it equals 0 at x=0")
elseif any(C_check_BC(0:X/10:X,H) ~= 0) % This technically shouldn't have the K(zanchor) and instead have the K(H) but it won't matter since the derivative needs to be zero
error("Recheck your true solution to make sure the flux is 0 through z=H")
else
Ctrue_func(0,0:H/10:H)
C_check_BC(0:X/10:X,H)
Q = @(xvar) vd*Ctrue_func(xvar,zanchor) - C_check_BC(xvar,zanchor);
end
mySource(xx,zz) = diff(Ctrue_func(xx,zz),xx) ...
- (diff(sD(zz)*diff(Ctrue_func(xx,zz),zz) ...
+ sA(zz)*Ctrue_func(xx,zz),zz) + sR(zz)*Ctrue_func(xx,zz))./u0_func(xx);
S = matlabFunction(mySource);
Stilde = @(xvar,zvar) S(xs*xvar,zs*zvar+zanchor)*xs/Cs;
for i=1:length(n)
Ctilde = schemeFunc(X,H,m(i),n(i),xs,zs,Cs,vd,alpha_p,Kr,C0,Q,zanchor,Stilde);
dz = H/n(i);
dx = X/m(i);
[XX,ZZ] = meshgrid(0:X/m(i):X,0:H/n(i):H);
Ctilde_true = Ctrue_func(xs*XX,zs*ZZ+zanchor);
myError(i) = ...
norm(Ctilde_true-Ctilde,"fro")*sqrt(dz/zs*dx/xs)/Cweight;
end
fig = figure;
loglog(1./n(1:end-0),myError(1:end-0),'-*b')
grid on
ylabel("Rel. Error","Interpreter","latex")
xlabel("$\Delta z$","Interpreter","latex")
myTitle = sprintf("$|| \\widetilde{C} - \\widetilde{C}_{true} ||_2\\:/\\:|| \\widetilde{C}_{true} ||_2$ with $\\alpha=%0.2f$",alpha_p);
tt = title(myTitle,"Interpreter","latex");
Running this gave the following errors
Please note that the zanchor value which is the same as $z^*$ is $z^*=100$, but I've tried this with all sorts of values and none produce anything even remotely correct.
I listed all the codes used here on a nearly identical MathWorks post, with the intention of seeing if anyone can spot a bug in my code.