f: source term of the Poisson equation
function ps_fft(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)
aa = -2.0/(dx*dx) - 2.0/(dy*dy)
bb = 2.0/(dx*dx)
cc = 2.0/(dy*dy)
hx = 2.0*pi/nx #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) # FFT of the data
```
*Fluids* **2019**, *4*, 159

```
e[1,1] = 0.0
for i = 1:nx
for j = 1:ny
data1[i,j] = e[i,j]/(aa + bb*cos(kx[i]) + cc*cos(ky[j]))
end
end
u = real(ifft(data1)) # inverse FFT of the data1
return u
end
```
✝ ✆

**Figure 16.** Comparison of exact solution and numerical solution computed using fast Fourier transform method for the Poisson equation with periodic boundary condition.

We develop the Poisson solver using a second-order central difference scheme for discretization of the Poisson equation. Instead of using the discretized equation, we can directly use Fourier analysis for the Poisson equation. This will compute the exact solution and the solution will not depend upon the grid size. The only change will be the change in Equation (91). We can easily derive the equation for Fourier coefficients in spectral solver by performing forward FFT of Equation (82). The Fourier coefficients for spectral solver can be written as

$$\vec{u}\_{m,n} = -\frac{\vec{f}\_{m,n}}{m^2 + n^2}.\tag{92}$$

We can get the solution by performing inverse FFT of the above equation. The implementation of spectral Poisson solver in Julia is outlined in Listing 19.

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

```