relax v2 times
gauss_seidel_mg(nxf, nyf, dx, dy, f, un, v2)
end
end
✝ ✆
```
Figure 23 displays the comparison of exact solution and numerical solution for full V-cycle multigrid framework. We use 512 × 512 grid resolution and hence we use nine levels from fine mesh to coarsest mesh. When the coarsest mesh has 2 × 2 grid resolution (i.e., since boundaries are set zero, only one grid point value is unknown and that can be solved algebraically), then it is usually referred to as the full multigrid framework. Figure 24 shows that the residual drops below tolerance in just nine iterations. This does not consist of a number of iterations for relaxation operations done for each intermediate grid levels. Although we report the outer iteration counter for multigrid results, we highlight that the workload in multigrid is more than the outer iteration counts. Equation (121) gives the total iteration count for a typical V-cycle multigrid method:

$$\text{V-cycle multiplier cost} = V \left( (v\_1 + v\_2) \sum\_{d=0}^{L-2} \frac{1}{2^d} + v\_3 \frac{1}{2^{L-1}} \right), \tag{121}$$

where *V* is the outer iteration count, *v*<sup>1</sup> and *v*<sup>2</sup> are the number of iterations during prolongation and restriction, *v*<sup>3</sup> is the number of iterations for coarsest grid, and *L* is total number of levels in full multigrid cycle. In a typical application, the values of *v*<sup>1</sup> = 2, *v*<sup>2</sup> = 2, and *v*<sup>3</sup> = 1 can work effectively (i.e., a higher value of *v*<sup>3</sup> can be used if the coarsest resolution is larger than 2 × 2 and intended to be solved by the same relaxation method). Table 2 provides the comparison of different iterative solvers for the Poisson equation.

**Table 2.** Comparison of different iterative methods in terms of iteration count, residual and CPU time for solving the Poisson equation on a resolution of *Nx* = 512 and *Ny* = 512.


**Figure 23.** Comparison of exact solution and numerical solution computed using multigrid framework for the Poisson equation.

**Figure 24.** Residual history comparison for different iterative methods for Poisson equation with Dirichlet boundary condition. The stopping criteria for all iterative methods is that the residual should reduce below 10−10. The grid resolution for all cases is *Nx* = 512 and *Ny* = 512.

#### **7. Incompressible Two-Dimensional Navier-Stokes Equation**

The Navier-Stokes equations for incompressible flow can be written as

$$
\nabla \cdot \mathbf{u} = 0,\tag{122}
$$

$$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u}\_\prime \tag{123}$$

where **u** = [u, v] is the velocity vector, *t* is the time, *ρ* is the density, *p* is pressure, and *ν* is the kinematic viscosity. By taking curl of Equation (123) and using ∇ × **u** = *ω*, we can derive the vorticity equation

$$
\nabla \times \frac{\partial \mathbf{u}}{\partial t} + \nabla \times (\mathbf{u} \cdot \nabla \mathbf{u}) = -\nabla \times \left(\frac{1}{\rho} \nabla p\right) + \nabla \times \nu \nabla^2 \mathbf{u}.\tag{124}
$$

The first term on the left side and last term on the right side becomes

$$
\nabla \times \frac{\partial \mathbf{u}}{\partial t} = \frac{\partial \omega}{\partial t}; \quad \nabla \times \nu \nabla^2 \mathbf{u} = \nu \nabla^2 \omega. \tag{125}
$$

Applying the identity ∇×∇ *f* = 0 for any scalar function *f* , the pressure term vanishes for the incompressible flow since the density is constant

$$
\nabla \times \nabla \left(\frac{p}{\rho}\right) = 0.\tag{126}
$$

The second term in Equation (123) can be writtern as

$$\mathbf{u} \cdot \nabla \mathbf{u} = \frac{1}{2} \nabla (\mathbf{u} \cdot \mathbf{u}) - \mathbf{u} \times (\nabla \times \mathbf{u}) = \nabla \left(\frac{\mathbf{u}^2}{2}\right) - \mathbf{u} \times \boldsymbol{\omega} \text{ where } \mathbf{u}^2 = \mathbf{u} \cdot \mathbf{u}. \tag{127}$$

The second term in Equation (124) becomes

$$\begin{split} \nabla \times (\mathbf{u} \cdot \nabla \mathbf{u}) &= \nabla \times \nabla \left(\frac{\mathbf{u}^2}{2}\right) - \nabla \times \mathbf{u} \times \boldsymbol{\omega} = \nabla \times (\boldsymbol{\omega} \times \mathbf{u}) \\ &= (\mathbf{u} \cdot \nabla) \boldsymbol{\omega} - \underbrace{(\boldsymbol{\omega} \cdot \nabla) \mathbf{u}}\_{\text{(vortex strength)}} + \underbrace{\boldsymbol{\omega} \cdot (\nabla \cdot \mathbf{u})}\_{\nabla \cdot \mathbf{u} = 0 \text{ (incompressible flow)}} + \underbrace{\mathbf{u} \cdot (\nabla \cdot \boldsymbol{\omega})}\_{\nabla \cdot (\nabla \times \mathbf{u}) = 0}, \end{split} \tag{128}$$

and the vortex stretching term vanishes in two-dimensional flows (i.e., *ω<sup>x</sup>* = 0, *ω<sup>y</sup>* = 0, *ω<sup>z</sup>* = *ω*)

$$(\boldsymbol{\omega} \cdot \nabla)\mathbf{u} = \left(\underbrace{\omega\_x}\_{0}\frac{\partial}{\partial x} + \underbrace{\omega\_y}\_{0}\frac{\partial}{\partial y} + \omega\_z\underbrace{\frac{\partial}{\partial z}}\_{0}\right) = 0. \tag{129}$$

We can further use below relations for two-dimensional flows

$$(\mathbf{u} \cdot \nabla)\omega = \mathbf{u}\frac{\partial \omega}{\partial \mathbf{x}} + \mathbf{v}\frac{\partial \omega}{\partial y},\tag{130}$$

and introducing a streamfunction with the following definitions:

$$\mathbf{u} = \frac{\partial \boldsymbol{\psi}}{\partial \boldsymbol{y}} \text{ and } \mathbf{v} = -\frac{\partial \boldsymbol{\psi}}{\partial \mathbf{x}} \text{,} \tag{131}$$

then the advection term becomes (i.e., sometimes called nonlinear Jacobian)

$$
\partial(\mathbf{u}\cdot\nabla)\omega = \frac{\partial\psi}{\partial y}\frac{\partial\omega}{\partial x} - \frac{\partial\psi}{\partial x}\frac{\partial\omega}{\partial y}.\tag{132}
$$

The vorticity equation for two-dimensional incompressible flow becomes

$$
\frac{\partial\omega}{\partial t} + \frac{\partial\psi}{\partial y}\frac{\partial\omega}{\partial x} - \frac{\partial\psi}{\partial x}\frac{\partial\omega}{\partial y} = \nu \left(\frac{\partial^2\omega}{\partial x^2} + \frac{\partial^2\omega}{\partial y^2}\right). \tag{133}
$$

The kinematic relationship between streamfunction and vorticity is given by a Poisson equation

$$
\frac{
\partial^2
\psi
}{
\partial x^2
} + \frac{
\partial^2
\psi
}{
\partial y^2
} = -\omega.
\tag{134}
$$

The vorticity-streamfunction formulation has several advantages over solving Equations (122) and (123). It eliminates the pressure term from the momentum equation and hence, there is no odd-even coupling between the pressure and velocity. We can use vorticity-streamfunction formulation directly on the collocated grid instead of using staggered grid. The number of equations to be solved in the vorticity-streamfunction formulation is also less than primitive variable formulation as it satisfies the divergence-free condition.

We use third-order Runge-Kutta numerical scheme (as discussed in Section 2.2) for the time integration. The right hand side terms in Equation (133) is discretized using the second-order central difference scheme similar to the diffusion term in heat equation (refer to Section 2.1). The nonlinear terms in Equation (133) is defined as the Jacobian

$$J(\omega, \psi) = \frac{\partial \psi}{\partial y} \frac{\partial \omega}{\partial x} - \frac{\partial \psi}{\partial x} \frac{\partial \omega}{\partial y}. \tag{135}$$

We can use any discretization method like explicit finite difference scheme (refer to Section 2.1) or compact finite difference scheme (refer to Section 2.4) for each term in Equation (135). Here, we use a numerical scheme introduced by Arakawa [44] for the discretization of Equation (135). This numerical scheme has conservation of energy, enstrophy and skew symmetry property and avoids computational instabilities arising from nonlinear interactions. The second-order Arakawa scheme for Equation (135) is given below

$$J(\omega, \psi) = \frac{J\_1(\omega, \psi) + J\_2(\omega, \psi) + J\_3(\omega, \psi)}{3},\tag{136}$$

where the discrete parts of the Jacobian are

$$J\_1(\omega, \psi) = \frac{(\omega\_{i+1,j} - \omega\_{i-1,j})(\psi\_{i,j+1} - \psi\_{i,j-1}) - (\omega\_{i,j+1} - \omega\_{i,j-1})(\psi\_{i+1,j} - \psi\_{i-1,j})}{4\Lambda x \Lambda y},\tag{137}$$

$$J\_2(\omega, \psi) = \frac{1}{4\Lambda x \Delta y} [\omega\_{i+1,j}(\psi\_{i+1,j+1} - \psi\_{i+1,j-1}) - \omega\_{i-1,j}(\psi\_{i-1,j+1} - \psi\_{i-1,j-1})] \tag{138}$$

$$\begin{split}-\omega\_{i,j+1}(\Psi\_{i+1,j+1}-\Psi\_{i-1,j+1})+\omega\_{i,j-1}(\Psi\_{i+1,j-1}-\Psi\_{i-1,j-1}) \\ J\_3(\omega,\boldsymbol{\upmu}) = \frac{1}{4\Delta x\Delta y}[\omega\_{i+1,j+1}(\boldsymbol{\upmu}\_{i,j+1}-\boldsymbol{\upmu}\_{i+1,j})-\omega\_{i-1,j-1}(\boldsymbol{\upmu}\_{i-1,j}-\boldsymbol{\upmu}\_{i,j-1}) \\ -\omega\_{i-1,j+1}(\boldsymbol{\upmu}\_{i,j+1}-\boldsymbol{\upmu}\_{i-1,j})+\omega\_{i+1,j-1}(\boldsymbol{\upmu}\_{i+1,j}-\boldsymbol{\upmu}\_{i,j-1})].\end{split}\tag{139}$$

#### *7.1. Lid-Driven Cavity Problem*

We test our two-dimensional Navier-Stokes solver using the lid-driven cavity benchmark problem for viscous incompressible flow [45]. The problem deals with a square cavity consisting of three rigid walls with no-slip conditions and a lid moving with a tangential unit velocity. The density of the fluid is taken to be unity. Therefore, we get *ν* = 1/Re, where Re is the Reynolds number of flow. The vorticity equation for lid-driven cavity problem can be written as

$$\frac{\partial\omega}{\partial t} = -\left(\frac{\partial\psi}{\partial y}\frac{\partial\omega}{\partial x} - \frac{\partial\psi}{\partial x}\frac{\partial\omega}{\partial y}\right) + \frac{1}{\text{Re}}\left(\frac{\partial^2\omega}{\partial x^2} + \frac{\partial^2\omega}{\partial y^2}\right). \tag{140}$$

The computational domain is square in shape with (*x*, *y*) ∈ [0, 1]×[0, 1]. We divide the computational domain into 64×64 grid resolution. All the walls have Dirichlet boundary conditions. The time integration implementation in Julia for two-dimensional Navier-Stokes solver is similar to Listing 2. We perform time integration from time *t* = 0 to *t* = 10 to make sure that steady state is reached and the residual reaches below 10−6. The residual is defined as the *L*<sup>2</sup> norm of the difference between two consecutive solutions. At each step of the Runge-Kutta numerical scheme, we also update the boundary condition for vorticity (i.e., further details can be found in [46]) and solve Equation (134) to update streamfunction. Any of the Poisson solvers mentioned in Section 6 can be used to solve this equation. We use fast sine transform Poisson solver discussed in Section 6.1.2 to get streamfunction from vorticity field as we have Dirichlet boundary conditions for all four walls. The implementation of FST Poisson solver in Julia is given in Listing 20. The functions to compute the right-hand side term for Runge-Kutta numerical scheme and to update boundary conditions are given in Listings 27 and 28 respectively. Figure 25 shows the vorticity and streamfunction field for the lid driven cavity problem.

**Listing 27.** Computation of right-hand side term in Equation (140) in Julia for two-dimensional incompressible Navier-Stokes equations. ✞ ☎

```