f: source term of the Poisson equation
function ps_spectral(nx,ny,dx,dy,f)
eps = 1.0e-6
kx = Array{Float64}(undef,nx)
ky = Array{Float64}(undef,ny)
data = Array{Complex{Float64}}(undef,nx,ny)
data1 = Array{Complex{Float64}}(undef,nx,ny)
e = Array{Complex{Float64}}(undef,nx,ny)
```

```
u = Array{Complex{Float64}}(undef,nx,ny)
hx = 2.0*pi/(nx*dx) #wavenumber indexing
for i = 1:Int64(nx/2)
kx[i] = hx*(i-1.0)
kx[i+Int64(nx/2)] = hx*(i-Int64(nx/2)-1)
end
kx[1] = eps
ky = kx
for i = 1:nx
for j = 1:ny
data[i,j] = complex(f[i,j],0.0)
end
end
e = fft(data)
e[1,1] = 0.0
for i = 1:nx
for j = 1:ny
data1[i,j] = e[i,j]/(-kx[i]^2 -ky[j]^2)
end
end
u = real(ifft(data1))
return u
end
```
✝ ✆

**Figure 17.** Comparison of error for spectral method (Top) and second-order CDS (bottom) for three different grid resolutions. The error is defined as (*x*, *<sup>y</sup>*) = <sup>|</sup>*uexact* <sup>−</sup> *<sup>u</sup>numerical*|.

The spectral Poisson solver gives the exact solution and the error between exact and numerical solution is machine zero error. We demonstrate it by using spectral Poisson solver for three different grid resolutions. We also use the central-difference scheme (CDS) Poisson solver at these grid resolutions and compare the discretization errors for two methods. We systematically change the grid spacing by a factor of two in both *x* and *y* direction and compute the *L*<sup>2</sup> norm of the discretization error. Figure 17 shows the contour plot for the error at three grid resolutions for spectral and second-order CDS. We can notice that the error for the spectral method does not depend on the grid resolution and the value of error is of the same order as machine zero error. On the other hand, for CDS the error reduces as we increase the grid resolution. This is because the truncation error goes down with the decrease in grid spacing. Figure 18 plots the *L*<sup>2</sup> norm of the error for spectral method and CDS. The slope of the discretization error is two for CDS showing that the scheme has second-order of accuracy. This method is used in CFD to validate the code and numerical methods. For the spectral method, the *L*<sup>2</sup> norm of the error has the machine zero value for all three grid resolutions.

**Figure 18.** Order of accuracy for spectral method and second-order CDS. We calculate discretization error at five grid resolutions and change the grid size by a factor of 2 in both *x* and *y* direction. The CDS Poisson solver give second-order accuracy. Spectral method give exact solution and error is of the same order as machine zero.

#### 6.1.2. Fast Sine Transform

The procedure described in Section 6.1.1 is valid only when the problem has periodic boundary condition, i.e., the solution that satisfies *ui*,*<sup>j</sup>* = *ui*+*Nx*,*<sup>j</sup>* = *ui*,*j*+*Ny* . If we have homogeneous Dirichlet boundary condition for the square domain problem, we need to use the sine transform given by

$$\mathfrak{u}\_{m,n} = \sum\_{i=0}^{N\_x - 1} \sum\_{j=0}^{N\_x - 1} u\_{i,j} \sin \frac{\pi mi}{N\_x} \sin \frac{\pi nij}{N\_y} \tag{93}$$

where *m* and *n* are the wavenumbers. This satisfies the homogeneous Dirichlet boundary conditions that *u* = 0 at *i* = 0, *Nx* and *j* = 0, *Ny*. If we substitute Equation (93) in Equation (83) and use trigonometric identity (sin(*A* ± *B*) = sin(*A*)cos(*B*) ± sin(*B*)cos(*A*)) we get

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

We can compute the spatial field using inverse sine transform given by

$$
\mu\_{i,j} = \frac{2}{N\_x} \frac{2}{N\_y} \sum\_{m=0}^{N\_x - 1} \sum\_{n=0}^{N\_y - 1} \bar{u}\_{m,n} \sin \frac{\pi m i}{N\_x} \sin \frac{\pi n j}{N\_y}.\tag{95}
$$

The three-step procedure to develop Poisson solver using sine transform can be listed as


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

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

The implementation of Poisson solver for problems with homogeneous Dirichlet boundary condition using the sine transform is given in Listing 20. Figure 19 displays the comparison of exact and numerical solution for the test problem computed using sine transform.

**Listing 20.** Implementation of Julia function for the Poisson solver using fast sine transform (FST) for the domain with Dirichlet boundary condition. ✞ ☎

✝ ✆

```