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, and developing fast parallel algorithms. The Poisson equation is a second-order elliptic equation and can be represented as
where \(\Delta x\) and \(\Delta y\) is the grid spacing in the \(x\) and \(y\) directions, and \(f_{i,j}\) is the source term at discrete grid
locations. If we write the above Equation at each grid point, we get a system of linear equations. For the Dirichlet boundary condition, we assume that the values of
\(u_{i,j}\) are available when \((x_{i},y_{j})\) is a boundary point. These eauations can be written in the standard matrix notation as
The boundary points of the domain are incorporated in the source term vector \(\mathbf{b}\) in the above equation.
We can solve the equation by standard methods for systems of linear equations, such as Gaussian elimination.
However, the matrix \(\mathbf{A}\) is very sparse and the standard methods are computationally expensive for large size of \(\mathbf{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.
We will further present a multigrid framework which scales linearly with a number of discrete grid points in the domain.
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,N_{x}]\times[1,N_{y}]\).
The solution at right \((i=N_{x}+1)\) and top \((j=N_{y}+1)\) boundary is obtained from left \((i=1)\) and bottom \((j=1)\) boundary, respectively.
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 \(\mathbf{u}\) and compute the source term \(\mathbf{f}\) at each grid location.
We then solve the Poisson equation for this source term \(\mathbf{f}\) and compare the numerically calculated field \(\mathbf{u}\) with the exact solution field \(\mathbf{u}\). The exact field \(\mathbf{u}\) and the corresponding source term \(\mathbf{f}\) used for direct Poisson solver are given below
This problem can have both periodic and Dirichlet boundary conditions. The computational domain is square in shape and is divided into \(512\times 512\) grid in \(x\) and \(y\) directions.
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 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\text{log}N)\), where \(N\) is the data size. In two-dimensions, the DFT for a square domain discretized equally in both directions is defined as
where \(u_{i,j}\) is the function in the spatial domain and the exponential term is the basis function corresponding to each point
\(\widetilde{u}_{m,n}\) 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
\(\widetilde{u}_{m,n}\) 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, \(\widetilde{u}_{0,0}\) represents the component of the function with the lowest wavenumber and \(\widetilde{u}_{N_{x}-1,N_{y}-1}\) 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
The three step procedure to develop a Poisson solver using FFT is given below
Apply forward FFT to find the Fourier coefficients of the source term in the Poisson equation (\(\widetilde{f}_{m,n}\)
from the grid values \({f}_{i,j}\))
Get the Fourier coefficients for the solution \(\widetilde{u}_{m,n}\) using below relation
\[\begin{split}\displaystyle u_{i-1,j}+u_{i+1,j}=\cfrac{2}{M} \cfrac{2}{N}\sum_{m=0}^{M-1}\sum_{n=0}^{N-1} \hat{u}_{m,n}\bigg\{\text{ sin }\cfrac{\pi m{(i-1)}}{M}\text{ sin }\cfrac{\pi n{j}}{N}+\text{ sin }\cfrac{\pi m{(i+1)}}{M}\text{ sin }\cfrac{\pi n{j}}{N}\bigg\}\\\end{split}\]
\[\begin{split}\begin{array}{l}
\alpha=\cfrac{\pi m{(i-1)}}{M},\quad \beta=\cfrac{\pi m{(i+1)}}{M}\\
\cfrac{\alpha+\beta}{2}= \cfrac{\pi mi}{M},\quad\cfrac{\alpha-\beta}{2}= -\cfrac{\pi m}{M}\\
\text{ sin }\cfrac{\pi m{(i-1)}}{M}+\text{ sin }\cfrac{\pi m{(i+1)}}{M}=2\text{ sin }\cfrac{\pi mi}{M}\text{cos}\cfrac{\pi m}{M}\\
\bigg\{\text{ sin }\cfrac{\pi m{(i-1)}}{M}+\text{ sin }\cfrac{\pi m{(i+1)}}{M}\bigg\}\text{ sin }\cfrac{\pi nj}{N}=2\text{ sin }\cfrac{\pi mi}{M}\text{ sin }\cfrac{\pi nj}{N}\text{cos}\cfrac{\pi m}{M}\\
\end{array}\end{split}\]
\[\begin{split}\begin{array}{l}
\alpha=\cfrac{\pi n{(j-1)}}{M},\quad \beta=\cfrac{\pi n{(j+1)}}{N}\\
\cfrac{\alpha+\beta}{2}= \cfrac{\pi nj}{N},\quad\cfrac{\alpha-\beta}{2}= -\cfrac{\pi n}{N}\\
\text{ sin }\cfrac{\pi n{(j-1)}}{N}+\text{ sin }\cfrac{\pi n{(j+1)}}{N}=2\text{ sin }\cfrac{\pi nj}{N}\text{cos}\cfrac{\pi n}{N}\\
\bigg\{\text{ sin }\cfrac{\pi n{(j-1)}}{N}+\text{ sin }\cfrac{\pi n{(j+1)}}{N}\bigg\}\text{ sin }\cfrac{\pi mi}{M}=2\text{ sin }\cfrac{\pi nj}{N}\text{ sin }\cfrac{\pi mi}{M}\text{cos}\cfrac{\pi n}{N}\\
\end{array}\end{split}\]
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.
