if the convergence critera (rms<tolerance) is satidfied, stop the iteration
if (rms/initial_rms) <= 1e-12
break
end
end
✝ ✆
```


1: Given *b* - Given source term in space, i.e., *b* = *f*(*x*, *y*) in Equation (82) 2: Given matrix operator A - Given discretization by Equation (107) 3: *u*(0) = *b* - Initialize the solution 4: *<sup>r</sup>*(0) <sup>=</sup> *<sup>b</sup>* − A*u*(0) - Initialize the residual 5: *p*(0) = *r*(0) - Initialize the conjugate 6: *k* = 0 7: *<sup>ρ</sup>*<sup>0</sup> = *<sup>r</sup>*(0)*<sup>T</sup> r*(0) 8: **while** tolerance met (or *k* < N ) **do** 9: *<sup>q</sup>*(*k*) <sup>=</sup> <sup>A</sup>*p*(*k*) 10: *<sup>α</sup><sup>k</sup>* <sup>=</sup> *<sup>ρ</sup>k*/*p*(*k*)*<sup>T</sup> q*(*k*) 11: *u*(*k*+1) = *u*(*k*) + *αkp*(*k*) - Update the solution 12: *<sup>r</sup>*(*k*+1) <sup>=</sup> *<sup>r</sup>*(*k*) <sup>−</sup> *<sup>α</sup>kq*(*k*) 13: *<sup>ρ</sup>k*+<sup>1</sup> <sup>=</sup> *<sup>r</sup>*(*k*+1)*<sup>T</sup> r*(*k*+1) 14: *β<sup>k</sup>* = *ρk*<sup>+</sup>1/*ρ<sup>k</sup>* 15: *p*(*k*+1) = *r*(*k*+1) + *βkp*(*k*) - Update the residual 16: check convergence; continue if necessary 17: *k* ← *k* + 1 18: **end while**

#### *6.3. Multigrid Framework*

We saw in Section 6.2 that the rate of convergence for iterative methods depends upon the iteration matrix *M*−1*N*. Point iterative methods like Jacobi and Gauss-Seidel methods have large eigenvalue and hence the slow convergence. As the grid becomes finer, the maximum eigenvalue of the iteration matrix becomes close to 1. Therefore, for very high-resolution simulation, these iterative methods are not feasible due to the large computational time required for residuals to go below some specified tolerance.

The multigrid framework is one of the most efficient iterative algorithm to solve the linear system of equations arising due to the discretization of the Poisson equation. The multigrid framework works on the principle that low wavenumber errors on fine grid behave like a high wavenumber error on a coarse grid. In the multigrid framework, we restrict the residuals on the fine grid to the coarser grid. The restricted residual is then relaxed to resolve the low wavenumber errors and the correction to the solution is prolongated back to the fine grid. We can use any of the iterative methods like Jacobi,

Gauss-Seidel method for relaxation. The algorithm can be implemented recursively on the hierarchy of grids to get faster convergence.

Let *u*˜Δ*<sup>x</sup>* denote the approximate solution after applying few steps of iterations of a relaxation method. The fine grid solution will be

$$\mathcal{A}(\vec{u}^{\Lambda x} + e^{\Lambda x}) = b^{\Lambda x},\tag{111}$$

$$\mathcal{A}\mathcal{e}^{\Lambda x} = b^{\Lambda x} - \mathcal{A}\vec{u}^{\Lambda x},\tag{112}$$

where *e*Δ*<sup>x</sup>* is the error at fine grid and the operator *A* is defined in Equation (84) for the Cartesian grid. We define the residual vector *<sup>r</sup>*Δ*<sup>x</sup>* <sup>=</sup> <sup>A</sup>*e*Δ*x*. This residual vector contains mostly low wavenumber errors. The correction at the fine grid is obtained by restricting the residual vector to a coarse grid (2Δ*x* grid spacing) and solving an equivalent system to obtain the correction. The equivalent system on a coarse grid can be written as

$$\mathcal{A}^{2\Delta x}\overline{\varepsilon}^{2\Delta x} = \mathcal{R}(r^{\Delta x}),\tag{113}$$

where *<sup>A</sup>*2Δ*<sup>x</sup>* is the elliptic operator on coarse grid, and <sup>R</sup> is the restriction operator. The discretized form of Equation (113) can be written as

$$\frac{\vec{\mathbb{E}}\_{i+1,j}^{2\Delta x} - 2\vec{\mathbb{E}}\_{i,j}^{2\Delta x} + \vec{\mathbb{E}}\_{i-1,j}^{2\Delta x}}{2\Delta x^2} + \frac{\vec{\mathbb{E}}\_{i,j+1}^{2\Delta x} - 2\vec{\mathbb{E}}\_{i,j}^{2\Delta x} + \vec{\mathbb{E}}\_{i,j-1}^{2\Delta x}}{2\Delta y^2} = \mathcal{R}(r\_{i,j}^{\Lambda x}).\tag{114}$$

We compute the approximate error at an intermediate grid using some relaxation operation. Once the approximation to *e*˜ <sup>2</sup>Δ*<sup>x</sup>* is obtained, it is prolongated back to the fine grid using

$$
\hat{\epsilon}^{\Delta x} = \mathcal{P}(\hat{\epsilon}^{2\Delta x}),
\tag{115}
$$

where <sup>P</sup> is the prolongation operator. The approximate solution at fine grid *<sup>u</sup>*˜Δ*<sup>x</sup>* is corrected with *<sup>e</sup>*<sup>ˆ</sup> <sup>Δ</sup>*x*. The approximate solution is relaxed for a specific number of iterations and convergence criteria is checked. If the residuals are below specified tolerance we stop or we repeat the steps. The procedure can be easily extended to more number of levels. The illustration of the V-cycle multigrid framework for two levels is shown in Figure 22. In this study, we use two iterations of the Gauss-Seidel method for relaxation at every grid level. The number of iterations can be set to a different number or we can also set the residual tolerance instead of a fixed number of iterations.

The residuals corresponding to the fine grid can be projected on to the coarse grid either using full-weighting. The full-weighting restriction operator is given by

$$r\_{\hat{i,j}}^{2\Delta x} = \frac{4r\_{2i-1,2j-1}^{\Delta x} + 2(r\_{2i-1,2j}^{\Delta x} + r\_{2i-1,2j-2}^{\Delta x} + r\_{2i-2,2j-1}^{\Delta x} + r\_{2i-2,2j-1}^{\Delta x}) + r\_{2i,2j}^{\Delta x} + r\_{2i-2,2j}^{\Delta x} + r\_{2i-2,2j-2}^{\Delta x}}{16}.\tag{116}$$

For boundary points, the residuals are directly injected from fine grid to coarse grid. The implementation of restriction operator in Julia is given in Listing 23.

The prolongation operator transfers the error from the coarse grid to the fine grid. We use direct injection for points which are common to both coarse and fine grid and bilinear interpolation for points which are on the fine grid but not on the coarse grid. The prolongation operations are given by below equations

$$\varepsilon\_{2i-1,2j-1}^{\text{Ax}} = \varepsilon\_{i,j}^{2\text{\AA}x} \, , \tag{117}$$

$$
\epsilon\_{2i-1,2j-1+1}^{\Delta x} = \frac{\epsilon\_{i,j}^{2\Delta x} + \epsilon\_{i,j+1}^{2\Delta x}}{2},
\tag{118}
$$

$$
\sigma\_{2i-1+1,2j-1}^{\Delta x} = \frac{c\_{i,j}^{2\Delta x} + c\_{i+1,j}^{2\Delta x}}{2},
\tag{119}
$$

$$e\_{2i-1+1,2j-1+1}^{\Delta x} = \frac{e\_{i,j}^{2\Delta x} + e\_{i,j+1}^{2\Delta x} + e\_{i+1,j}^{2\Delta x} + e\_{i+1,j+1}^{2\Delta x}}{2}.\tag{120}$$

The implementation of prolongation operation in Julia is given in Listing 24.

**Listing 23.** Implementation of restriction operation in multigrid framework. ✞ ☎

**Figure 22.** Illustration of multigrid V-cycle framework with three levels between fine and coarse grid. For full multigrid cycle, the coarsest mesh has 2 × 2 grid resolution and is solved directly.

**Listing 24.** Implementation of prolongation operation in multigrid framework. ✞ ☎

```