Not getting correct numerical solution for Advection-Diffusion-Reaction eqn

Question

Objective: I am trying to numerically solve $C(x,y,t)$ from the following advection-diffusion-reaction equation in 2D space (x,y) and time. I will be testing my numerical solution with an approximate 1D analytical solution (given on 2nd page in this file). However, the issue is that I am not able to get correct numerical solution.

$$\begin{align} \text{ADR Equation: }\frac{\partial C}{\partial t} + \nabla\left(v C - D\nabla{C} \right)= \alpha C \end{align}$$

Method: I discretized the above ADR equation in 2D using finite-difference implicit scheme and I get the following equation in difference form.

$$\begin{align} p_1C^{n+1}_{i,j-1}+p_2C^{n+1}_{i-1,j}+p_3C^{n+1}_{i,j}+p_4C^{n+1}_{i+1,j}+p_5C^{n+1}_{i,j+1} = C^{n}_{i,j} \end{align}$$

I want to solve this as a system of equations using $A^{n+1} C^{n+1}=C^{n}$

Variables: $A$ is a matrix formed from variables $p_1, p_2, p_3, p_4, p_5$, which are constants in time but vary along $(x,y)$ grid. Superscript $n$ represents time index and subscripts $i,j$ represent space along $x,y$.

I.C: $$\textrm{$C(x,y,0) = C_i$, assume $C_i$ as 1}$$

B.C: $$\frac{\partial C}{\partial x} = \frac{\partial C}{\partial y} = 0 \text{ at boundaries}$$. $$C(0,y=1,t) = \begin{cases} C_0\quad & \text{if $0<t\leq t_0$} \\ 0 \quad & \text{if $t>t_0$} \end{cases}$$ $\textrm{assume $C_0$ and $t_0$ as 5000 and 2 days, respectively}$

Numerical Solution: I estimated $C(x,y,t)$ for $t=0$ to $t=100$ days, where dimensions of $x$ and $y$ are shown below in results.

$$\textrm{Fig-1: Snapshots of $C(x,y)$ at various $t$}$$

In Fig-2, I am showing $C(x,t)$ for $0<t\leq t_0$ and $t=0$ to $t=100$ days:

$$\textrm{Fig-2: Snapshots of $C(x,t)$ obtained after averaging $C(x,y,t)$ over Y-axis}$$

Analytical Solution: Here I am showing $C(x,t)$ for $0<t\leq t_0$ and $t=0$ to $t=100$ days.

$$\textrm{Fig-3: Snapshots of $C(x,t)$. I am trying to match this with Fig-2}$$

Issue: I have been trying to figure out if there is a problem in my numerical solution developed as per the implicit scheme discussed above; however, I have spent several hours and I couldn't figure it out. Please let me know if I missed providing anything.

Answer 1

One of the best ways to test a PDE solver is to use the method of manufactured solutions. Essentially, you modify the PDE (and discretization) by adding a source term that yields an exact solution known in advance. You can then compare your numerical solution against the exact solution for debugging purposes. You should test your numerical solution against a known exact solution, and you should also test to make sure that the order of accuracy of your discretization is roughly what you expect (e.g., if you have a second-order accurate discretization in space, you should expect the error to be proportional to the square of the mesh spacing when performing a refinement study). A good reference is this Sandia National Laboratory technical report.

If you want to be more thorough, you could also set up and execute unit tests on various terms in your discretization (e.g., the diffusion term), comparing the numerical value on a known input to known (exact) output.