Interface fluxes (Rusanov scheme)
for i = 1:nx+1 for m = 1:3
f[i,m] = 0.5*(fR[i,m]+fL[i,m]) - 0.5*ps[i]*(qR[i,m]-qL[i,m])
end end
end
function wavespeed(nx,gamma,uL,uR,ps)
```
*Fluids* **2019**, *4*, 159

```
gm = gamma-1.0
for i = 1:nx+1
#Left and right states:
rhLL = uL[i,1]
uuLL = uL[i,2]/rhLL
eeLL = uL[i,3]/rhLL
ppLL = gm*(eeLL*rhLL - 0.5*rhLL*(uuLL*uuLL))
hhLL = eeLL + ppLL/rhLL
rhRR = uR[i,1]
uuRR = uR[i,2]/rhRR
eeRR = uR[i,3]/rhRR
ppRR = gm*(eeRR*rhRR - 0.5*rhRR*(uuRR*uuRR))
hhRR = eeRR + ppRR/rhRR
alpha = 1.0/(sqrt(abs(rhLL)) + sqrt(abs(rhRR)))
uu = (sqrt(abs(rhLL))*uuLL + sqrt(abs(rhRR))*uuRR)*alpha
hh = (sqrt(abs(rhLL))*hhLL + sqrt(abs(rhRR))*hhRR)*alpha
aa = sqrt(abs(gm*(hh-0.5*uu*uu)))
ps[i] = abs(aa + uu)
end
end
✝ ✆
```
#### **6. Two-Dimensional Poisson Equation**

In this section, we explain different methods to solve the Poisson equation which is encountered in solution to the incompressible flows. The Poisson equation is solved at every iteration step in the solution to the incompressible Navier-Stokes equation due to the continuity constraint. Therefore, many studies have been done to accelerate the solution to the Poisson equation using higher-order numerical methods [37], and developing fast parallel algorithms [38]. The Poisson equation is a second-order elliptic equation and can be represented as

$$
\frac{
\partial^2 
\mu
}{
\partial x^2
} + \frac{
\partial^2 
\mu
}{
\partial y^2
} = f.
\tag{82}
$$

Using the second-order central difference formula for discretization of the Poisson equation, we get

$$\frac{u\_{i+1,j} - 2u\_{i,j} + u\_{i-1,j}}{\Delta x^2} + \frac{u\_{i,j+1} - 2u\_{i,j} + u\_{i,j-1}}{\Delta y^2} = f\_{i,j},\tag{83}$$

where Δ*x* and Δ*y* is the grid spacing in the *x* and *y* directions, and *fi*,*<sup>j</sup>* is the source term at discrete grid locations. If we write Equation (83) at each grid point, we get a system of linear equations. For the Dirichlet boundary condition, we assume that the values of *ui*,*<sup>j</sup>* are available when (*xi*, *yj*) is a boundary point. These equations can be written in the standard matrix notation as

$$
\begin{bmatrix} D\_{xy} & D\_y & 0 & \cdots & D\_x & \cdots & \cdots & 0 \\ D\_y & D\_{xy} & D\_y & \ddots & & \ddots & & \vdots \\ 0 & D\_y & D\_{xy} & D\_y & \ddots & & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ D\_x & & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & & \ddots & & \ddots & D\_y & D\_{xy} & D\_y & 0 \\ \vdots & & & & & & & \\ \vdots & & & & & & \ddots & D\_y & D\_{xy} & D\_{xy} \\ \vdots & & & & & & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & D\_x & \cdots & 0 & D\_y & D\_{xy} \end{bmatrix} \begin{bmatrix} u\_{1,1} \\ u\_{1,2} \\ \vdots \\ \vdots \\ u\_{2,1} \\ \vdots \\ \vdots \\ u\_{N\_x, N\_y} \\ \vdots \\ \vdots \\ \vdots \\ \vdots \end{bmatrix}\_{N\_{\mathcal{X}}, N\_{\mathcal{Y}}} = \begin{bmatrix} f\_{1,1} \\ f\_{1,2} \\ \vdots \\ \vdots \\ f\_{2,1} \\ \vdots \\ f\_{N\_{\mathcal{X}}, N\_{\mathcal{Y}}} \end{bmatrix},\tag{84}$$

where

$$D\_{xy} = \frac{-2}{\Delta x^2} + \frac{-2}{\Delta y^2}; \quad D\_x = \frac{1}{\Delta x^2}; \quad D\_y = \frac{1}{\Delta y^2}.$$

The boundary points of the domain are incorporated in the source term vector *b* in Equation (84). We can solve Equation (84) by standard methods for systems of linear equations, such as Gaussian elimination [39]. However, the matrix *A* is very sparse and the standard methods are computationally expensive for large size of *A*. In this paper, we discuss direct methods based on fast Fourier transform (FFT) and fast sine transform (FST) for periodic and Dirichlet boundary condition. We also explain iterative methods to solve Equation (84). We will further present a multigrid framework which scales linearly with a number of discrete grid points in the domain.

#### *6.1. Direct Solver*

Figure 15 shows the finite difference grid for two-dimensional CFD problems. We use two different boundary conditions for direct Poisson solver: periodic and Dirichlet boundary condition. For the Dirichlet boundary condition, the nodal values of a solution are already known and do not change, and hence, we do not do any calculation for the boundary points. In case of periodic boundary condition, we do calculation for grid points (*i*, *j*) between [1, *Nx*] × [1, *Ny*]. The solution at right (*i* = *Nx* + 1) and top (*j* = *Ny* + 1) boundary is obtained from left (*i* = 1) and bottom (*j* = 1) boundary, respectively.

**Figure 15.** Finite difference grid for two-dimensional problems. The calculation is done only for points shown by blue color. The solution at black points is extended from left and bottom boundary for periodic boundary condition. The solution is already available at all four boundaries for Dirichlet boundary condition.

We perform the assessment of direct solver using the method of manufactured solution. We assume certain field *u* and compute the source term *f* at each grid location. We then solve the Poisson equation for this source term *f* and compare the numerically calculated field *u* with the exact solution field *u*. The exact field *u* and the corresponding source term *f* used for direct Poisson solver are given below

$$u(x,y) = \sin(2\pi x)\sin(2\pi y) + \frac{1}{16^2}\sin(32\pi x)\sin(32\pi y),\tag{85}$$

$$f(\mathbf{x}, y) = -8\pi^2 \sin(2\pi \mathbf{x}) \sin(2\pi y) - 8\pi^2 \sin(32\pi \mathbf{x}) \sin(32\pi y). \tag{86}$$

This problem can have both periodic and Dirichlet boundary conditions. The computational domain is square in shape and is divided into 512 × 512 grid in *x* and *y* directions.

#### 6.1.1. Fast Fourier Transform

There are two different ways to implement fast Poisson solver for the periodic domain. One way is to perform FFTs directly on the Poisson equation, which will give us the spectral accuracy. The second method is to discretize the Poisson equation first and then apply FFTs on the discretized equation. The second approach will give us the same spatial order of accuracy as the numerical scheme used for discretization. We use the second-order central difference scheme given in Equation (83) for developing a direct Poisson solver.

The Fourier transform decomposes a spatial function into its sine and cosine components. The output of the Fourier transform is the function in its frequency domain. We can recover the function from its frequency domain using the inverse Fourier transform. We use both the function and its Fourier transform in the discretized domain which is called the discrete Fourier transform (DFT). Cooley-Tukey proposed an FFT algorithm that reduces the complexity of computing DFT from *O*(*N*2) to *O*(*N*log*N*), where *N* is the data size [17]. In two-dimensions, the DFT for a square domain discretized equally in both directions is defined as

$$\tilde{\mu}\_{m,n} = \sum\_{i=0}^{N\_x - 1} \sum\_{j=0}^{N\_y - 1} \mu\_{i,j} e^{-i2\pi \left(\frac{\mu i}{N\_x} + \frac{nj}{N\_y}\right)},\tag{87}$$

where *ui*,*<sup>j</sup>* is the function in the spatial domain and the exponential term is the basis function corresponding to each point *u*˜*m*,*<sup>n</sup>* in the Fourier space, *m* and *n* are the wavenumbers in Fourier space in *x* and *y* directions, respectively. This equation can be interpreted as the value in the frequency domain *u*˜*m*,*<sup>n</sup>* can be obtained by multiplying the spatial function with the corresponding basis function and summing the result. The basis functions are sine and cosine waves with increasing frequencies. Here, *u*˜0,0 represents the component of the function with the lowest wavenumber and *u*˜*Nx*−1,*Ny*−<sup>1</sup> represents the component with the highest wavenumber. Similarly, the function in Fourier space can be transformed to the spatial domain. The inverse discrete Fourier transform (IDFT) is given by

$$\mu\_{i,j} = \frac{1}{N\_x} \frac{1}{N\_y} \sum\_{m=-N\_x/2}^{N\_x/2-1} \sum\_{n=-N\_y/2}^{N\_y/2-1} i\mathbb{1}\_{m,n} \mathfrak{e}^{\left(\frac{\pi i}{N\_x} + \frac{nj}{N\_y}\right)}\,,\tag{88}$$

where 1/(*NxNy*) is the normalization term. The normalization can also be applied to forward transform, but it should be used only once. If we use Equation (88) in Equation (83), we get

$$\eta\_{m,n} \left( \frac{e^{l\frac{2\pi n}{N\_x}} - 2 + e^{l\frac{-2\pi n}{N\_x}}}{\Delta x^2} + \frac{e^{l\frac{2\pi n}{N\_y}} - 2 + e^{l\frac{-2\pi n}{N\_y}}}{\Delta y^2} \right) = \tilde{f}\_{m,n} \tag{89}$$

in which we use the definition of cosine to yield

$$\ddot{a}\_{m,n}\left(\frac{2\cos\left(\frac{2\pi n}{N\_x}\right)-2}{\Delta x^2} + \frac{2\cos\left(\frac{2\pi n}{N\_y}\right)-2}{\Delta y^2}\right) = \tilde{f}\_{m,n}.\tag{90}$$

If we take the forward DFT of Equation (90), we get a spatial function *ui*,*j*. The three step procedure to develop a Poisson solver using FFT is given below


$$\Pi\_{m,n} = \frac{f\_{m,n}}{\left(\frac{2}{\Delta x^2}\cos\left(\frac{2\pi m}{N\_x}\right) + \frac{2}{\Delta y^2}\cos\left(\frac{2\pi n}{N\_y}\right) - \frac{2}{\Delta x^2} - \frac{2}{\Delta y^2}\right)}.\tag{91}$$

• Apply inverse FFT to get the grid values *ui*,*<sup>j</sup>* from the Fourier coefficients *u*˜*m*,*n*.

Implementation of Poisson solver using FFT is given in Listing 18. Figure 16 shows the exact and numerical solution for for the test problem.

**Listing 18.** Implementation of Julia function for the Poisson solver using fast Fourier transform (FFT) for the domain with periodic boundary condition. ✞ ☎

```