## 3 Numerical Algorithm
### The Physics-Based Preconditioners
#### 3.4.2. Based on the Crank-Nicolson Method  
To obtain the preconditioning equation, we recast the discretization (8) as  
\[\frac{c_{i,j}^{k+1}-c_{i,j}^{n}}{\Delta t}-\frac{[((D(x,y)c_{x}))_ {x}]_{i,j}^{k+1}+[((D(x,y)c_{x}))_{x}]_{i,j}^{n}}{2}\] \[-\frac{[((D(x,y)c_{y}))_{y}]_{i,j}^{k+1}+[((D(x,y)c_{y}))_{y}]_{i,j}^{n}}{2}\approx\frac{f_{i,j}^{k}+f_{i,j}^{n}}{2}. \tag{22}\]  
If we follow the similar steps as in Section 3.4.1, we obtain the following update formula  
\[-c\delta c_{i-1,j}\,-\,b\delta c_{i+1,j}+a\delta c_{i,j}-d\delta c_{i,j+1}-e \delta c_{i,j-1}\approx-\text{res}c_{i,j}, \tag{23}\]where  
\[\text{res}c_{i,j} \approx c_{i,j}^{k}-c_{i,j}^{n}-\frac{\Delta t}{2}([((D(x,y)c_{x}))_{x}]_{i, j}^{k}+[((D(x,y)c_{x}))_{x}]_{i,j}^{n}) \tag{24}\] \[-\frac{\Delta t}{2}([((D(x,y)c_{y}))_{y}]_{i,j}^{k}+[((D(x,y)c_{y} ))_{y}]_{i,j}^{n})-\frac{\Delta t}{2}(f_{i,j}^{k}+f_{ij}^{n}),\]  
and  
\[a\,=\,1\,+\,\frac{\Delta t}{2\Delta x^{2}}D_{i+\frac{1}{2},j}\,+\,\frac{ \Delta t}{2\Delta x^{2}}D_{i-\frac{1}{2},j}+\frac{\Delta t}{2\Delta y^{2}}D_{ i,j+\frac{1}{2}}+\frac{\Delta t}{2\Delta y^{2}}D_{i,j-\frac{1}{2}},\]  
\[b=\frac{\Delta t}{2\Delta x^{2}}D_{i+\frac{1}{2},j^{\prime}}\,\,\,\,c=\frac{ \Delta t}{2\Delta x^{2}}D_{i-\frac{1}{2},j^{\prime}}\,\,\,\,d=\frac{\Delta t }{2\Delta y^{2}}D_{i,j+\frac{1}{2}},\,\,\,\,e=\frac{\Delta t}{2\Delta y^{2}}D_ {i,j-\frac{1}{2}}.\]  
for \(i\,=\,1,\cdots,M_{x}-1\) and \(j\,=\,1,\cdots,M_{y}-1\). Notice that we need boundary conditions when solving Equations (21) and (23) for \(\delta c\). Essentially, we define few ghost cell values similarly to those described in Section 3.1 as \(\delta c_{-1,j}\approx\delta c_{0,j}\) and \(\delta c_{M_{x},j}\,\approx\,\delta c_{M_{x}-1,j}\) for \(j\,=\,1,\cdots,M_{y}-1\) and \(\delta c_{i,-1}\approx\delta c_{i,0}\) and \(\delta c_{i,M_{y}}\approx\delta c_{i,M_{y}-1}\) for \(i\,=\,1,\cdots,M_{x}-1\). The update Equations (21) and (23) are solved iteratively to produce the preconditioned Krylov vector \(z\). The iterative procedure is based on a multi-grid V-cycle method specifically designed for our applications. A typical multi-grid V-cycle method consists of few recursively applied Gauss-Seidel (GS) or other basic such as Jacobi iterations from finest to coarsest and from coarsest to finest grid levels [41; 43]. It is well known that the GS or other basic iterative methods (sometimes referred to as smoothers) are very good at smoothing the high frequency errors in just a few iterations. However, they are extremely slow to produce fully converged solutions. Often, the iteration count for convergence is related to the problem size (e.g., \(M_{x}\times M_{y}\) in our 2D cases) which is impractical. We take advantage of the GS smoother to design our multi-grid V-cycle in the following way. At the beginning of each V-cycle, we perform _two_ Gauss-Seidel iterations to remove or smooth out the high-frequency errors occurring at the finest grid. Then, we restrict ourselves to a coarser grid at which we form coarse grid solutions and the related coarse residuals to apply _two_ more GS sweeps. We note that, at this coarse level, the previously low-frequency errors act like high-frequency ones and the application of a small number of GS sweeps are sufficient to smooth them out. This continues recursively until we reach the coarsest level (two grid cells in each direction) at which we perform _hundred_ GS iterations. At the end of the coarsest level iterations, we recursively apply the prolongation step to interpolate the coarse level solutions to the finer levels at which we carry out _two_ more GS sweeps. The completion of these steps defines _one_ multi-grid V-cycle. In our test problems, we carry out _five_ V-cycles that we found sufficient to decrease the GMRES iteration count to very small numbers. However, one can always increase or decrease the number of multi-grid V-cycles depending on the physical model and its numerical discretization. **Remark 2**.