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.

enter image description here

$$\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: enter image description here

$$\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.

enter image description here

$$\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.

  • 1
    $\begingroup$ Have you checked if your code behaves correctly for simpler special cases first, e.g., constant coefficient heat equation? $\endgroup$
    – Kirill
    Jul 19, 2015 at 22:29
  • $\begingroup$ I've checked and it gives mixed indications on its output. The concentrations fall off very quickly to zero. $\endgroup$
    – user5510
    Jul 19, 2015 at 22:50
  • 1
    $\begingroup$ You should compare it to known exact solutions. Checking whether a solution approaches zero is not conclusive. $\endgroup$
    – Kirill
    Jul 19, 2015 at 22:57
  • $\begingroup$ @Pupil: Your first link has a permissions issue. $\endgroup$ Jul 20, 2015 at 7:22
  • $\begingroup$ The problem was the way I was forming the transmissibility matrix $A$ in $A\times X=B$ $\endgroup$
    – user5510
    Aug 10, 2015 at 1:44

1 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.

  • 3
    $\begingroup$ @Pupil: (1) If you have to read only one section, read Section 3. It's 7 pages, and describe the essentials of the method. The remaining pages are either background, or extended examples. (2) The extensive information you provided isn't helpful; it isn't diagnostic. Most of it is necessary for posing your question, but not for pinpointing what is wrong. (3) I am not going to test the code for you. Testing is an important skill you should learn so that you can debug problems like this one in the future, because they will occur again. $\endgroup$ Jul 20, 2015 at 8:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.