ADR equation implicit solution: Penta-diagonal matrix for a 2D $N\times N$ system

Question

Objective: I am trying to simulate the following advection-diffusion-reaction equation in 2D space (x,y) and time. $$\begin{align} \text{ADR Equation: }\frac{\partial C}{\partial t} + \nabla\left(v.C - D\nabla{C} \right)= \alpha.C \end{align}$$ I discretized the above ADR equation in 2D using finite-difference implicit scheme and as a result I get the following discretized equation.

$$\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}$$ where, $p_1, p_2, p_3, p_4, p_5$ are constants in time.

I want to solve this as a system of equations using $A^{n+1}.C^{n+1}=C^{n}$, with no-flow i.e. $C=0$ outside the boundary domain. Here, $A^{n+1}$ would be a penta-diagonal, symmetric (not sure about this) and a diagonally dominant matrix. I have derived matrix $A$ for $2\times2$, $3\times3$ and $4\times4$ systems. For example, below you can see matrix $A$ for $3\times3$ and $4\times4$ systems, respectively.

Issue:

1. I am not sure if the form of matrix $A$ I have derived is correct because as per my understanding it should be symmetric, however, it's not as per my derivations.
2. Owing to an unsymmetric form of the matrix $A$ I need help to efficiently form $A$ for $N\times N$ system.

Would appreciate if someone could use their awesome numerical skills to answer these issues.

Answer 1

This is mostly just a matter of assigning an index to each vertex. If you have $n\times n$ unknowns $C_{i,j}$, $0\leq i,j\leq n-1$, then you would typically map the variable corresponding to the vertex $(i,j)$ to something like an index $I(i,j) = i + n j$. So each variable $C_{i,j}$ that is indexed by a pair of numbers $(i,j)$ is equivalent to the variable $C_{I(i,j)}$ indexed by a single number.

With this mapping, you rewrite the equations as $$ p_1 C_{I(i,j)} + p_{2:3} C_{I(i\pm1,j)} + p_{4:5} C_{I(i,j\pm1)} = \cdots, $$

This means that the $I(i,j)$-th row (the $(i+jn)$-th row) of the matrix $A^{n+1}$ has the entries $p_{1,2,3,4,5}$ in columns $I(i,j)$, $I(i\pm1,j)$, $I(i,j\pm1)$.

With this indexing method the nonzero entries will all fall into the diagonals $\pm n$ (coefficients $p_4$, $p_5$), $\pm1$ ($p_2$, $p_3$), and $0$ ($p_1$). The main thing to bear in mind, when constructing such a matrix knowing its diagonal elements, is that the diagonal elements for non-existent boundary variables such as $I(i,n)$, $I(i,-1)$, $I(-1,j)$, $I(n,j)$ should be zeroed out.

So for a $3\times 3$ grid, you would get the matrix (I dropped the minus signs here) $$ \begin{pmatrix} p_1 & p_2 & 0 & p_4 & 0 & 0 & 0 & 0 & 0 \\ p_3 & p_1 & p_2 & 0 & p_4 & 0 & 0 & 0 & 0 \\ 0 & p_3 & p_1 & 0 & 0 & p_4 & 0 & 0 & 0 \\ p_5 & 0 & 0 & p_1 & p_2 & 0 & p_4 & 0 & 0 \\ 0 & p_5 & 0 & p_3 & p_1 & p_2 & 0 & p_4 & 0 \\ 0 & 0 & p_5 & 0 & p_3 & p_1 & 0 & 0 & p_4 \\ 0 & 0 & 0 & p_5 & 0 & 0 & p_1 & p_2 & 0 \\ 0 & 0 & 0 & 0 & p_5 & 0 & p_3 & p_1 & p_2 \\ 0 & 0 & 0 & 0 & 0 & p_5 & 0 & p_3 & p_1 \end{pmatrix} $$