# Can't get proper numerical convergence for complicated Advection-​Diffusion-​Reaction PDE

2 views (last 30 days)
David Gillcrist on 29 Nov 2023
Edited: David Gillcrist on 5 Dec 2023
I have written an advection-diffusion-reaction PDE using a Crank-Nicolson finite difference scheme. The detail of which and the derivation can be found here: My Stack Exchange Question. I'm writing this to see if anyone see's an issue with the code I've written to solve this problem. The functions I have are
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;
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);
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
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)
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");
This produces the error plot:
Additionally, I use the following functions attached to run the setup.
Torsten on 30 Nov 2023
Edited: Torsten on 30 Nov 2023
For use of pdepe, you have to multiply your equation by u_0~ to get the derivative term "free" of a multiplying factor. And I'd define the flux f in "pdepe" as D*dC/dz, not as D*dC/dz + A*C. Instead, put the d/dz(A*C) together with S into the source term s in "pdepe". Otherwise, you will come into trouble when defining the boundary conditions.
How close 0 do you think $u_0$ can be before issues?
I usually set x~ and/or t~ to a small number to avoid division by 0 - it should be no problem to try which value(s) is/are appropriate and see whether it influences the result when the coding is finished.
David Gillcrist on 1 Dec 2023
Edited: David Gillcrist on 5 Dec 2023
@Torsten This ends up doing great! Thank you. Though I'm still concerned about why my scheme approach doesn't work.

### Categories

Find more on Boundary Conditions in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!