# SailFFish: A Lightweight, Parallelised Fast Poisson Solver Library

Joseph Saverin
Hermann Föttinger Institute
Technische Universität Berlin
Berlin, Germany
<EMAIL_ADDRESS>

###### Abstract

A solver for the Poisson equation for 1D, 2D and 3D regular grids is
presented. The solver applies the convolution theorem in order to efficiently
solve the Poisson equation in spectral space over a rectangular computational
domain. Conversion to and from the spectral space is achieved through the use
of discrete Fourier transforms, allowing for the application of highly
optimised $\mathcal{O}(N\log N)$ algorithms. The data structure is configured
to be modular such that the underlying interface for operations to, from and
within the spectral space may be interchanged. For computationally demanding
tasks, the library is optimised by making use of parallel processing
architectures. A range of boundary conditions can be applied to the domain
including periodic, Dirichlet, Neumann and fully unbounded. In the case of
Neumann and Dirichlet boundary conditions, arbitrary inhomogeneous boundary
conditions may be specified. The desired solution may be found either on
regular (cell-boundary) or staggered (cell-centre) grid configurations. For
problems with periodic, Dirichlet or Neumann boundary conditions either a
pseudo-spectral or a second-order finite difference operator may be applied.
For unbounded boundary conditions a range of Green’s functions are available.
In addition to this, a range of differential operators may be applied in the
spectral space in order to treat different forms of the Poisson equation or to
extract highly accurate gradients of the input fields. The underlying
framework of the solver is first detailed, followed by a range of validations
for each of the available boundary condition types. Finally, the performance
of the library is investigated. The code is free and publicly available under
a GNU v3.0 license at github.com/ZeppSav/SailFFish.

_K_ eywords Poisson Solver $\cdot$ Fast Fourier Transform $\cdot$ Spectral
Solver, Numerical Integration

##  1 Introduction

For any given vector field $\vec{v}$ which demonstrates $C^{2}$ smoothness–
implying that the second derivative (or higher) is continuous– the field may
be expressed using the Helmholtz decomposition:

$\vec{v}=\nabla\times\vec{\psi}-\nabla\phi\text{ ,}\hskip 14.22636pt\text{
where }\hskip 14.22636pt\nabla\cdot\vec{\psi}=0\text{ .}$ (1)

Here $\vec{\psi}$ is a vector potential which represents the solenoidal
(divergence-free) part of the field and $\phi$ is a scalar potential which
represents the irrotational part of the field. The divergence free constraint
on the vector potential is necessary to ensure a unique decomposition. By
taking the curl of (1), one retrieves the Poisson equation: a linear, second
order partial differential equation (PDE):

$\nabla^{2}\vec{\psi}=\vec{f}\text{ , }$ (2)

where $\vec{f}=-\nabla\times\vec{v}$ is the curl of the vector field. By
taking the divergence of (1), one again retrieves a Poisson equation relating
the Laplacian of the scalar potential and the divergence of the vector field
$\nabla\cdot\vec{v}$. This direct connection to continuous vector fields
illustrates precisely why the Poisson equation is so ubiquitous in physical
modelling. The conservation of physical field quantities leads to vector and
scalar potentials which are amenable to solution with the Poisson equation or
simplifications thereof. Applications include fluid dynamics, electrostatics,
electromagnetism, gravitation, diffusion and heat transfer, amongst many
others. A range of methods exist for the solution of the Poisson equation
depending on the desired accuracy of the solution and boundary conditions
(BCs) of the domain. These may be broken down broadly into three categories:
multigrid methods, fast multipole methods and spectral methods [1]. For the
work here spectral methods have been applied, however the other solution
methods shall briefly be described for comparison.

###  1.1 Multigrid Methods

Multigrid (MG) methods can be applied for the solution of a range of both
nonlinear and linear PDEs on arbitrary grid configurations over a range of
domain topologies [2, 3, 4, 5]. The main concept of the MG method is to solve
the desired PDE on a range of grid resolutions $\mathcal{H}_{i=0,1,\cdots n}$.
This shall be described first for a two grid levels $\mathcal{H}_{0}$,
$\mathcal{H}_{1}$. An initial solution is found on a disjoint portion of the
finer grid using e.g. finite differences. A residual error is calculated by
comparing source and preliminary solution. This residual is anterpolated
(adjoint interpolation) from the base grid $\mathcal{H}_{0}$ to
$\mathcal{H}_{1}$. The solution is then found on the coarser grid and
interpolated back down to the finer grid to improve the initial guess. This
gives rise to a iterative process of progressively smoothing the solution
until a convergence criteria is satisfied. This process is generally extended
to higher grid levels to further improve efficiency. Numerous schemes are
available for optimally traversing between grid levels e.g. $V$ and $W$-cycles
depending on the problem type [6]. It can be demonstrated that for a grid with
$N$ grid points, the computational expense scales as $\mathcal{O}(N)$ [7]. An
additional advantage of such methods is their amenability to irregular domains
or domains with variable resolution.

###  1.2 Fast Multipole Methods

The Fast Multipole Method (FMM) was initially developed by Greengard & Rokhlin
[8, 9]. This method makes use of the Green’s function solution to the Poisson
equation:

$\psi(r)=\mathcal{G}*f=\int_{\mathcal{R}^{n}}f(r^{\prime})\,\mathcal{G}(r-r^{\prime})\,d^{n}r\text{
,}$ (3)

where ($*$) indicates the convolution operator. The Green’s function
$\mathcal{G}$ represents the fundamental solution to the unbounded Poisson
equation: $\nabla^{2}\mathcal{G}(r)=\delta(r)$ where $\delta(r)$ is the Dirac
delta function [10]. FMM exploits the nature of $\mathcal{G}$ to create a low-
rank approximation to the far-field influence. $\mathcal{G}(r)$ may be
expressed as a multipole expansion which, in the far-field requires a smaller
expansion for a desired accuracy $\epsilon$. The use of a hierarchical grid
allows coarser levels of approximation with larger interaction distances. The
difficulty in formulating the multipole expansion of the desired kernel
motivated development of the kernel-independent FMM (KIFMM) [11] or similar
methods such as the multi-level multi-interaction cluster (MLMIC) [12, 13]. An
advantage of the FMM is that it naturally handles unbounded BCs, although
periodic or homogeneous Neumann or Dirichlet BCs may also be treated by using
the method of images [14, 15]. A significant advantage of the FMM is the
ability to handle sparse or even disjoint source distributions. Typically the
FMM scales either as $\mathcal{O}(N\log N)$ or $\mathcal{O}(N)$ depending on
problem and setup [8].

###  1.3 Spectral Methods

The third class of methods exploit the linearity of the Poisson equation by
representing the solution as a sum of basis functions. It shall be assumed
without loss of generality that in the case of interest here, these are
trigonometric functions. The forward and inverse Fourier transforms of a
function are introduced:

$\text{Forward: }\hskip
14.22636pt\hat{f}(k)=\int_{-\infty}^{\infty}f(x)e^{-2\pi
ikx}\,dx\text{,}\hskip 28.45274pt\text{Inverse: }\hskip
14.22636ptf(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{f}(k)e^{2\pi
ikx}\,dk\text{ , }$ (4)

This transforms the representation of the function from real ($x$) space to
spectral ($k$) space. The Poisson equation may be written in the spectral
space as:

$-k^{2}\hat{\psi}(k)=\hat{f}\text{ ,}\hskip 14.22636pt\text{ where }\hskip
14.22636ptk^{2}=\begin{cases}k_{x}^{2}&\text{in 1D}\\\
k_{x}^{2}+k_{y}^{2}&\text{in 2D}\\\ k_{x}^{2}+k_{y}^{2}+k_{z}^{2}&\text{in
3D}\end{cases}\text{ ,}$ (5)

and $k_{i}$ are the angular wave numbers in the spectral space. Applying the
convolution theorem, convolution in the real space is represented in the
spectral space by point-wise multiplication. Herein lies the enormous
advantage of this approach. The appropriate Green’s function (3) need only be
calculated once, converted to the spectral representation and stored.
Thereafter the entire convolution of any given source field is equivalent to a
multiplication in the spectral space. In practice, the continuous transforms
of (4) are treated with a discrete Fourier transform (DFT) which represents
the Fourier transform in a discrete, truncated spectral space. Application of
the DFT allows one to exploit the fast Fourier transform (FFT) to achieve a
computational expense $\mathcal{O}(N\log N)$ [16]. Although this appears to
not be a significant improvement over MG or FMM approaches, the myriad
applications of the FFT in signal processing implies that highly optimised
libraries exist for the calculation of FFTs. Unlike the MG or FMM approaches
whereby a desired accuracy of the solution $\epsilon$ is achieved, the use of
the spectral approach allows, under suitable conditions, spectral accuracy to
be achieved. This implies the solution is practically exact, with the error
saturating at machine precision. A range of BCs may be represented including
Neumann, Dirichlet, periodic and unbounded. These shall be detailed in the
proceeding section. The method has the disadvantage that the entire grid must
be generated and stored for a given resolution which leads to large memory
overheads. This makes the method less suitable for sparse source
distributions. In addition, extensions of the method to non-rectangular
domains or irregular grids can be challenging [17].

For the work carried out here, the spectral method has been implemented within
SailFFish, an open-source library written in C++. The library has been written
in such a way that emphasis is placed on simplicity, modularity and
adaptability. The code is available for free online under a GNU v3.0 license
at github.com/ZeppSav/SailFFish. The rest of the paper is devoted to an
overview and validation of the SailFFish solver. In Section 2 the numerical
method implemented in the solver is detailed. In Section 3 the framework of
the solver is described including details on the application of object-
oriented programming, data type modularity and input/output interfaces. In
Section 4 a range of validation cases have been investigated to demonstrate
the accuracy of the solver. In Section 5 library performance is detailed along
with scaling behaviour. Finally in Section 6 conclusions are drawn and a scope
for future features and applications is described.

##  2 Method of Solution

The spectral method has been applied in SailFFish for the solution of the
Poisson equation. The methodology and solution steps are described in this
section along with the options for different solver types. For illustrative
purposes a simple 1D periodic case will first be described. As detailed in
Section 1.3, it shall be assumed that the function is expressed as a
trigonometric series:

$f(x)=\sum_{k=0}^{N}\hat{f}_{k}e^{2\pi i\frac{kx}{L_{x}}}=\sum_{k=0}^{N}\
\hat{f}_{k}\left(\cos 2\pi k\frac{x}{L}+i\sin 2\pi k\frac{x}{L}\right)\text{
,}$ (6)

where $L$ is the length of the domain. The periodicity of the solution over
the domain can be deduced from this representation. The trigonometric
functions represent solution eigenvectors. By substituting this representation
into (2), one obtains the following representation of the solution on the
grid:

$-\left(\frac{2\pi k}{L}\right)^{2}\hat{\psi}_{k}=\hat{f}_{k}\text{ .}$ (7)

The factor multiplying the $\hat{\psi}_{k}$ represents the spectral
eigenvalues corresponding to the solution eigenvectors [18]. Regarding the
solution accuracy, the (backwards) transformed values $\psi_{n}$ on the grid
will correspond exactly to the solution of the Poisson equation, to machine
precision. This is referred to as demonstrating _spectral_ accuracy. It is
observed that the case $k=0$ represents a division by zero. The Fourier mode
$k=0$ represents the average value of the function $\psi$ and taking this to
be zero leads to a zero mean and avoids numerical instabilities.
An alternative approach to the solution of the Poisson equation is to directly
apply the representation (6) with finite differences (FD) [19]. Unlike the
accuracy achieved with the spectral eigenvalues described above, this approach
converges with the accuracy of the FD scheme applied.
A third approach is to calculate the Green’s function (3) in real space, and
then transform this to retrieve the spectral space representation. This
approach has the advantage that mollified Green’s kernels may be applied which
allow for smoothing of the input field [20, 21, 22]. This approach has been
applied to unbounded solver cases in SailFFish – see Section 2.5.

###  2.1 Transformation between Real and Spectral Space Using DFTs

As opposed to the continuous forward and inverse Fourier transforms as defined
by (4), in practice the grid values are specified on a discrete, equispaced
grid. The conversions to and from the discrete spectral space are achieved
with the forward and inverse DFT, defined as:

$\text{Forward: }\hskip
14.22636pt\hat{f}_{d}(k)=\sum_{n=0}^{N-1}f_{d}(x)e^{-2\pi
ik\frac{n}{N}}\text{,}\hskip 28.45274pt\text{Inverse: }\hskip
14.22636ptf(x)=\frac{1}{N}\sum_{n=0}^{N-1}\hat{f}_{d}(k)e^{2\pi
ik\frac{n}{N}}\text{ , }$ (8)

Depending on the BCs being applied, along with the type of the input data
(purely real or complex), different numerical forms of the DFT may be chosen
in order to improve performance. As an example, whenever the input data is an
array of $n$ purely real values $f=[f_{1},f_{2},\dots,f_{n}]$, the output of
the DFT can be shown to demonstrate Hermitian redundancy. A real-to-complex
(R2C) DFT may be applied which expresses the output vector as a complex array
with $\frac{n}{2}+1$ elements. This simplification achieves approximately a
factor of two improvement in speed and memory overhead. The inverse transform
is achieved in a similar manner.
Depending on the dimension of the problem being solved, 1D, 2D or 3D DFTs are
executed. The 1D DFT is applied to a single, contiguous array of data. Higher
dimensional DFTs are applied by carrying out the DFT sequentially in
1D-pencils. A 2D DFT carried out on a set of data points $f(x_{i},y_{j})$ is
therefore calculated as follows:

  1. 1.

X forward transform: Forward DFT taken in $x$ direction to transform
$f(x_{i},y_{j})$ into $\hat{f}(k_{x},y_{j})$,

  2. 2.

Y forward transform: Forward DFT taken in $y$ direction to transform
$\hat{f}(k_{x},y_{j})$ into $\hat{f}(k_{x},k_{y})$.

The inverse procedure follows precisely the same template, however applying
the inverse DFT and proceeding in the reserve order. In the case of the
aforementioned 2D example, the inverse transform is executed as follows:

  1. 1.

Y inverse transform: Backward DFT taken in $y$ direction to transform
$\hat{f}(k_{x},k_{y})$ into $\hat{f}(k_{x},y_{j})$,

  2. 2.

X inverse transform: Backward DFT taken in $x$ direction to transform
$\hat{f}(k_{x},y_{j})$ into $f(x_{i},y_{j})$.

There are two approaches for handling these procedures in 2D. In the first
approach, the partially transformed data after step (1) is re-aligned (or
stored in an intermediate array) such that it is contiguous in memory for
application of the second DFT. This has the advantage of allowing more
flexibility in choosing BCs in each spatial direction [20]. The second
approach makes use of devoted 2D DFT algorithms which automatically apply the
DFTs to the correctly ordered data based on the specified grid dimensions. The
latter approach is applied in SailFFish, as will be detailed in Section 3. The
extension to 3D cases is straightforward.

Figure 1: A 2D paraboloid distribution of the source term $f$ (left) along
with its representation $\hat{f}$ in the spectral space (right). The logarithm
of $\hat{f}$ is shown purely for visual purposes.

###  2.2 Calculating and Retrieving the Solution

Within SailFFish the spatial convolution is carried out in the spectral space
as described in Section 1.3. This is always achieved by carrying out four
steps:

  1. 1.

Forward Transform: The source distribution $f$ is transformed to a spectral
representation $\hat{f}$ via a forward DFT;

  2. 2.

Convolution: Point-wise multiplication of $\hat{f}$ with the spectral space
multiplier (Either the corresponding eigenvalue or the Fourier transform of
the Green’s function $\hat{\mathcal{G}}$);

  3. 3.

Spectral Operations: If necessary, additional operations within the spectral
space are carried out;

  4. 4.

Backward Transform: The spectral solution $\hat{\mathcal{G}}\hat{f}$ is
transformed to the real space $f$ via an inverse DFT.

The general procedure is therefore relatively straightforward. A visualisation
of the outputs of step 1) are shown in Figure 1. Despite the simplicity of the
above approach, there are a range of solver options. Type of DFT: This is
determined by the form of the data (purely real, complex) being transformed
and which type of BCs are being enforced (R2C, R2R, DFT). Spectral space
convolution kernel: These can either be derived from the spectral or finite
difference approaches as described in the previous section, or by applying
solutions of the Green’s kernel (3). Spectral space operator: A range of
operators may be applied in the spectral domain to modify the form of the
Poisson equation under investigation or to extract spectral gradients. These
options are detailed in Section 2.3.3.

###  2.3 Solver Options in SailFFish

The specification of the solver type and solver setup in SailFFish depends on
the desired solution field, BCs and grid properties. The range of options
available are described in this section.

####  2.3.1 Solver Types

Currently SailFFish has the ability to resolve either 1D/2D/3D scalar fields
or 3D vector fields. In the case of a 3D vector field, the three components
are resolved by individual component-wise application of the corresponding
solver methodology for a 3D scalar solver. Keeping this is mind, four solver
classes are available:

  * •

Poisson_Periodic_1D/2D/3D/3DV: Periodic BCs are applied.

  * •

Poisson_Dirichlet_1D/2D/3D/3DV: Dirichlet (odd function extension) BCs are
applied.

  * •

Poisson_Neumann_1D/2D/3D/3DV: Neumann (even function extension) BCs are
applied.

  * •

Poisson_Unbounded_1D/2D/3D/3DV: The problem is treated with free-space BCs.

In all cases, the 3DV option represents the case of a 3D vector field,
otherwise scalar fields are resolved. The first three options described above
represent _bounded_ solver types – see Section 2.4. The final solver type is
treated using a different solution method and is separately described in
Section 2.5. Currently only uniform BCs are available within SailFFish. For
details of an implementation which allows for arbitrary BCs, the reader is
referred to the FLUPS library [20]. The reason for this is elaborated upon in
Section 3.

####  2.3.2 Grid Type

SailFFish solves the Poisson equation (2) on regular grids in 1D, 2D and 3D
cases. The resolution of the grid is always specified by the number of grid
cells, $N_{c}$. During solver initialization the upper and lower spatial
coordinate for each grid dimension must be specified. Additionally, one has
the option to specify the position on the grid where the solution is being
found. If the REGULAR grid option is specified, the solution will be
calculated at the cell boundary positions. The number of resolved grid
position $N_{s}$ changes depending on the type of solver being applied – see
Section 2.4. If the STAGGERED grid option is specified, the solution is found
at the centre of the cells. In this case the solution is found at $N_{c}$
positions. Numerically, the grid positions are given by:

$X_{i}=\begin{cases}X_{l}+i(X_{u}-X_{l})/N_{c}&\text{for }i=0,1,\dots
N_{s}\hskip 28.45274pt\text{{REGULAR} grid option}\\\
X_{l}+(i+0.5)(X_{u}-X_{l})/N_{c}&\text{for }i=0,1,\dots N_{c}-1\hskip
11.38109pt\text{{STAGGERED} grid option}\end{cases}\text{ , }$ (9)

where $X_{l}$ and $X_{u}$ are the lower and upper values of the grid domain,
respectively. This scheme is applied for each of the resolved spatial
directions.

####  2.3.3 Differential Operators

In some cases it may be desirable to apply spectral differentiation to the
transformed source field $\hat{f}$. Inspecting (4), one observes:

$\frac{d}{dx}f(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{d}{dx}\left(\hat{f}(k)e^{2\pi
ikx}\right)\,dk=\frac{1}{2\pi}\int_{-\infty}^{\infty}ik\,\hat{f}(k)e^{2\pi
ikx}\,dk\text{ , }$ (10)

The process of multiplying by $ik$ in the frequency domain hence has the
result of extracting the (spectrally accurate) derivatives of the function
$f$. This may also be applied with a convolution kernel $\hat{\mathcal{G}}(k)$
in order to modify the form of the Poisson equation being solved. A range of
options exist in SailFFish to extract the spatial gradients:

  * •

DIV: Resolves the divergence of the input field

  * •

GRAD: Resolves the gradients of the input field

  * •

CURL: Resolves the curl of the input field

  * •

NABLA: Resolves the Laplacian of the input field

The application of spectral differentiation to the point-wise solution to PDEs
is commonly referred to as the Pseudo-spectral method [23].

###  2.4 Bounded Solvers

If either the nature or the numerical value of the BCs on the domain are
specified, the solution on the domain shall be referred to as a bounded
solution. Three forms of bounded solvers can be specified: Dirichlet (odd
symmetry), Neumann (even symmetry), or periodic. The choice of transform types
and their corresponding eigenvalues are described below.

####  2.4.1 Transform Types

The Dirichlet and Neumann type solvers allow the use of a real-to-real (R2R)
DFT. Purely real forms of (6) are utilised whereby only the $\sin$ (discrete
sine transform – DST) or $\cos$ (discrete cosine transform – DCT) expansions
are applied for symmetric odd or symmetric even BCs, respectively. Depending
on whether a regular or staggered grid option is chosen also dictates the
solution eigenvectors and therewith the type of DST or DCT being applied.
These are summarised for all boundary condition and grid type variations in
Table 2.4.1. Further references are provided here for a detailed description
of the corresponding eigenvector functions [18, 20, 24, 25]. By nature of the
fact that it is composed of periodic trigonometric functions, a DFT is itself
periodic. For periodic BCs the direct complex-to-complex (C2C) standard DFT is
applied. Although this may be simplified for purely real data by applying the
R2C transform described in Section 2.1, this is configured for complex inputs
for generality for future application cases. Real data inputs are stored as
complex data inputs with imaginary component zero.

Table 1: DFT Type and Eigenvalues for Bounded Solver Types

BC |  |  |  |  |
---|---|---|---|---|---
Type |  Grid |  |  |  |
|  DFT |  |  |  |
|  DFT |  |  |  |
|  $\mathcal{G}_{i}$ |  |  |  |
|  $\mathcal{G}_{i}$ |  |  |  |
| (FD2) |  |  |  |
Dirichlet | Regular | DST-I | DST-I | $-\left(\pi\frac{k+1}{L}\right)^{2}$ | $-\frac{4}{H^{2}}\sin^{2}\frac{\pi(k+1)}{2(N+1)}$
Dirichlet | Staggered | DST-II | DST-III | $-\left(\pi\frac{k+1}{L}\right)^{2}$ | $-\frac{4}{H^{2}}\sin^{2}\frac{\pi(k+1)}{2N}$
Neumann | Regular | DCT-I | DCT-I | $-\left(\pi\frac{k}{L}\right)^{2}$ | $-\frac{4}{H^{2}}\sin^{2}\frac{\pi k}{2(N-1)}$
Neumann | Staggered | DCT-II | DST-III | $-\left(\pi\frac{k}{L}\right)^{2}$ | $-\frac{4}{H^{2}}\sin^{2}\frac{\pi k}{2N}$
Periodic |  Both |  DFT |  Inverse DFT |  $-\left(2\pi\frac{k}{L}\right)^{2}$: k $<\frac{N}{2}$ |
|  $-\frac{4}{H^{2}}\sin^{2}\frac{\pi k}{N}$ |  |  |  $-\left(2\pi\frac{N-k}{L}\right)^{2}$: k $\geq\frac{N}{2}$ |

####  2.4.2 Choice of Eigenvalues

For bounded problems, two options for the eigenvalues are available. The first
are pseudo-spectral (PS) eigenvalues. As was illustrated above for the
periodic case, there are derived directly from the eigenvalue resulting from
differentiation in the real space of the trigonometric expansion. These vary
due to the corresponding zeros of the eigenvectors for particular BC-grid type
combinations. The second option is the use of the second-order FD (FD2)
eigenvalues. Use of this option has the advantage that it allows for the use
of inhomogeneous BCs for Dirichlet and Neumann solver types. This is a result
of the fact that the finite differences are built into the spectral solution.
In order to enforce the desired boundary values a simple weighting of the
boundary source terms may be applied [19]. This is automatically accounted for
in SailFFish based on the solver and grid type. An illustration of
inhomogeneous BCs is given in Fig. 2. The correct eigenvalues to all available
BC and grid type combinations are also provided in Table 2.4.1.

Figure 2: The solution to a 2D Dirichlet problem with homogeneous (left) and
inhomogeneous (right) boundary conditions using FD2 type eigenvalues. The
source field $f$ is equivalent in both cases.

In order to ensure the periodicity of the desired solution for spectral
eigenvalues the application of a cyclic Green’s function is required. This can
be observed for the periodic solvers PS eigenvalue, which is mirrored about
the point $k=\frac{N}{2}$.

###  2.5 Unbounded Solvers

If free-space BCs are specified, the solution on the domain shall be referred
to as an unbounded solution. There are numerous approaches to handling this
problem including the approach due to James & Lackner [26, 27], the Vico-
Greengard method [28], or the Hockney-Eastwood method [29, 30]. The latter has
been applied in SailFFish. In order to mimic the effect of free boundary
conditions, the domain of the Green’s function is doubled in the dimensions of
interest and mirrored about the point $N=k$, as was done for the periodic
spectral eigenvalues. Due to the periodic nature of the DFT, this has the
effect that the convolution of any point in the frequency domain extends to
the boundary as it would for the unbounded solution. In order to ensure that
the convolution in the frequency domain occurs without reflection, the source
field must also be doubled, and the additional solver regions are zero-padded
to avoid any artificial influence. A visualisation of this is given in Fig. 3.
After the periodic solution has been found on the doubled domain, the solution
within the desired region is extracted and the remaining region is discarded.
For an elegant explanation of this concept in 1D and for mixed BCs, the reader
is referred to the FLUPS documentation [20].

####  2.5.1 Transform Types

For unbounded cases, the direct DFT is applied. The source function and
Green’s function may therefore be complex values. Although not yet implemented
in SailFFish, this allows the solution of the Helmholtz equation [31, 20]. In
the case that purely real data is investigated, the R2C and C2R transforms are
applied to improve solver efficiency.

####  2.5.2 Greens Functions

Unlike the bounded solvers whereby the Green’s functions are available for
direct application in the frequency space, more care must be taken for
unbounded solutions. The problem of the undefined Green’s function at the
origin cannot by circumvented. In addition, the original approach of Hockney &
Eastwood only displays $\mathcal{O}(h^{2})$ convergence, where $h$ is the grid
size [30]. In order to overcome this, a range of options have been suggested.
Hejlesen et al. suggested a range of Gaussian smoothings which progressively
approximate the singular solution $\nabla^{2}\mathcal{G}=\delta(r)$ by
applying $P^{th}$ order approximations to the inverse Gaussian regularisation
[21]. The same author derived a spectrally compact kernel which achieves
spectral accuracy on finite grids [22]. Within SailFFish the spectrally
compact Hejlesen kernels are available for 1D, 2D and 3D cases. For 2D and 3D
solvers the Hejlesen Gaussian approximations of order $P=2,4,\dots 10$ may
also be selected.

Figure 3: Operations necessary for application of the Hockney-Eastwood method.
The doubling and mirroring of the Green’s function in real space is carried
out (left). The source domain is also doubled and the remaining space is zero-
padded (right). For the case here a bump function is being simulated –see
Section 4.3.

###  2.6 Semi-unbounded Solvers

It is possible to manually break 2D or 3D DFT operations down to a set of
successive 1D DFT processes. This is the approach taken in a range of other
Poisson solver packages [20, 32]. This approach has the advantage that the
Laplace operator may be split in the desired direction, yielding a 2D
Helmholtz equation which allows for the treatment of mixed boundary conditions
[31]. By applying combinations of DSTs, DCTs and DFTs this allows for the
implementation of fully arbitrary boundary conditions [20]. Within SailFFish
the optimised $N$-dimensional DFT algorithms are applied to reduce
communication times and allow simple development on a graphical processing
unit (GPU) – as explained in Section 3. This feature is therefore currently
not available in SailFFish, however it is planned to implement this in future
work.

##  3 Solver Framework

SailFFish has been prepared with the following objectives in mind:

  1. 1.

Simplicity, readability and abstraction;

  2. 2.

Optimise for shared-memory platforms;

  3. 3.

Minimise communication time and maximise application of existing libraries and
algorithms.

Abstraction has been exploited with the object-oriented structure of C++ to
avoid code repetition and improve readability. Class objects and inheritance
are described in the proceeding section. Preparing the software in this way
simplifies tailored application to specific problems for future applications.
SailFFish is designed to be run on shared-memory devices because, put simply,
not every user has access to cluster computing equipment and the overarching
concept should be efficiency rather than power. For all computationally
demanding tasks, the multi-platform shared-memory parallel programming
framework OpenMP has been utilised to parallelise calculations [33].

The final point refers to exploitation of existing optimised programs or APIs
in the interest of maintaining readable code. Take for example the execution
of 2D and 3D FFTs. If the software is configured to handle mixed BCs, then it
is necessary to carry out the transform into and out of spectral space
individually in 1D arrays. This requires reorganisation of data such that it
is contiguously aligned for the 1D FFT algorithms. This necessity both
complicates the code (particular for arbitrary boundary conditions) and
increases memory communication overhead as described in [20]. In SailFFish
this is avoided by directly exploiting 2D and 3D FFT algorithms which exist
within the chosen FFT engine. One advantage of this approach is that data may
be passed as very large monolithic blocks. An additional advantage is that
operations on the data (e.g. convolution) may be carried out along a single
index.

This approach is particularly relevant for the solution of unbounded problems
– see Section 2.5. The Hockney-Eastwood algorithm requires the doubling of the
domain in each dimension, leading to an increase in calculation expense and
memory overhead by a factor of approximately $2^{d}$ where $d$ is the
dimension of the problem [29]. This can be reduced to a factor of
approximately $2^{d-1}$ by carrying out the transforms in each direction
sequentially and discarding the additional information in the zero-padded
region [29, 21]. This necessarily greatly reduces the size of the FFTs. As
described above however it also requires reordering the data for the FFT
algorithms. This is avoided for now by simply treating the entire doubled
domain in order to adhere to the mantra of the simplification. It is the hope
of the author that in the near future the necessity to double the domain will
be avoided by making direct use of the James-Lackner method to treat the
unbounded problem [27, 26]. The proceeding sections are devoted to a
description of the architecture of the SailFFish library.

###  3.1 Architecture of SailFFish

There are two parent classes within SailFFish which contain the majority of
the functionality: the DataType class and the Solver class. These are
described here. A visualisation of the class data structure is given in Fig.
4.

####  3.1.1 DataType Class

The DataType parent class contains a range of virtual methods relevant to
accessing (input/output), manipulating (convolution) and carrying out
transforms of data within the spectral space. For a given choice of FFT tool,
a corresponding class is derived from the DataType parent class and the
corresponding virtual methods must be implemented. Currently two DataType
bases are available in SailFFish:

  * •

FFTW_DataType makes use of the highly optimised and popular FFTW library
(Fastest Fourier Transform in the West) [24, 25, 34]. This library has broad
functionality and displays consistently good performance;

  * •

cuFFT_DataType makes use of the cuda language to leverage highly accelerated
parallel execution on an NVIDIA graphical processing unit (GPU) [35]. Although
this library has far less functionality than FFTW, the speedups capable on a
GPU unit allow for extremely fast execution on a desktop.

Although these only represent two options, there are a range of other options
depending on target device, compilation environment and intended solver type
which may warrant the implementation of further DataType classes. A brief list
of possible options include the Intel MKL FFT framework [36], the OpenCL FFT
library [37], or implementation in other languages such as Python FFT [38].
Depending on the chosen data type, additional implementation may be necessary.

DataType ClassesDataType ParentFFTW DataTypecuFFT DataType–User-defined–Solver
ClassesSolver Parent1D_Scalar2D_Scalar3D_Scalar3D_Vector–User-defined– Figure
4: Architecture of SailFFish

####  3.1.2 Solver Class

The Solver parent class is derived directly from the user-defined DataType
parent class and contains all of the functionality and interfacing methods for
the FFT tool. By inheriting this in the solver implementations, an additional
layer of abstraction is achieved. This class furthermore contains all of the
virtual functions relevant to grid declaration and manipulation, solver option
specification, preparation of the Green’s functions and convolution execution.
Four types of solver are derived from the Solver parent class: 1D_Scalar,
2D_Scalar, 3D_Scalar and 3D_Vector, named for the type of the source (and
solution) data. The solvers relevant to specific BCs (see Section 2.3.1) are
then derived from these classes. Within each solver class the following
virtual methods must be declared and defined:

  * •

FFT_Data_Setup\- Allocates data and FFT Transform types;

  * •

Specify_Operator\- Specify desired spectral space operator types (Section
2.3.3);

  * •

Forward_Transform\- Form of the DFT from real to spectral space;

  * •

Backward_Transform\- Form of the DFT from spectral to real space;

  * •

Convolution\- Form of the convolution (purely real, complex);

  * •

Specify_Greens_Function\- Calculates the Green’s function on the grid for the
convolution step.

Although these are declared and defined for the aforementioned solver types in
SailFFish, the option exists for user defined solver classes with any
additional operations necessary for a desired application.

###  3.2 Generating and Executing a Solver

A short overview is given here of the routines which are carried out to
generate and execute a solver within SailFFish.

####  3.2.1 Solver Initialization

Numerous functions are called during the initialization of a solver. These are
listed here:

  1. 1.

X/Y/Z_Grid_Setup Grid positions and statistics are specified.

  2. 2.

Datatype_Setup Datatype-specific initializations are carried out.

  3. 3.

FFT_Data_Setup FFT Data arrays are allocated and FFT plans are specified.

  4. 4.

Specify_Greens_Function Specification of either eigenvalues (bounded BCs) or
Green’s functions (unbounded BCs) arrays.

  5. 5.

Prepare_Dif_Operators Spectral-space differential operators are prepared (if
appropriate).

The specification of the Green’s function is the most expensive process here,
as many expensive function evaluations must be carried out. This process is
parallelised with OpenMP. The computational expense of the initialization is
not of significant consequence as this process is only carried out once at the
beginning of a simulation. The solver may thereafter be executed arbitrarily
many times.

####  3.2.2 Solver Execution

The execution process is relatively straightforward for simplicity. The steps
are listed chronologically:

  1. 1.

Set_Input The values of $f$ on the grid are transferred into the necessary
arrays for transforming.

  2. 2.

Forward_Transform Forward FFTs are carried out on the data.

  3. 3.

Convolution Convolution in the spectral space is carried out. This includes
other procedures in the spectral space such as spectral differentiation.

  4. 4.

Backward_Transform Backward (inverse) FFTs are carried out on the data.

  5. 5.

Get_Output The values of $\vec{\psi}$ on the grid are transferred into the
necessary arrays for post-processing.

The most expensive steps here are the forward and backward FFTs. An overview
of solver timings is given in Section 5

##  4 Validation

To demonstrate the convergence of the solvers a set of validation calculations
have been carried out. Depending on the type of solver, a suitable source
function has been chosen for which an analytical solution to the Poisson
equation is known. It should be noted that the convergence achieved is a
function of the input field $f$. The accuracy of the numerical solution
depends on the number of continuous derivatives of the source function. It is
therefore helpful for a convergence study to select a source function which
displays $C^{\infty}$ smoothness.
For all presented cases SailFFish has been configured to have double floating
point accuracy. For all results shown, the convergence of the solver is
demonstrated by showing the infinite norm of the error:

$E_{\infty}=\text{sup}\\{\psi(x,y,z)-\psi_{a}(x,y,z)\\}\text{ , }$ (11)

where $\psi_{a}$ represents the known analytical solution. This error norm is
chosen as the $E_{\infty}$-norm bounds the $E_{2}$-norm. The absolute value of
error is scaled by the lowest value of the solution
$\psi_{a,min}=\text{inf}\\{\psi_{a}\\}$ on the grid. In order to make the
error results comparable, all values of $E_{\infty}$ have been scaled such
that they correspond to the grid solution with $\psi_{a,min}=1$.

###  4.1 Bounded Solvers: Homogeneous Boundary Conditions

The first class of solvers investigated are the bounded solvers. The following
1D test functions have been used:

$S(x,L_{x})=\sin\left(4\pi\frac{x}{L_{x}}\right)\text{, }\hskip
28.45274ptC(x,L_{x})=\cos\left(4\pi\frac{x}{L_{x}}\right)\text{, }$ (12)

where $L_{x}$ is the length of the domain in $x$ direction and the coordinate
$x$ is taken to be relative to the lower bound of the $x$-domain ($x_{l}$). In
higher dimensional cases, the product of this function in each direction is
taken. Derivation of the analytical solution is trivial. For homogeneous
Dirichlet BCs for example, the following function has been used:

$f_{a}(x,y,z)=S(x,L_{x})\,S(y,L_{y})\,S(z,L_{z})\hskip
8.53581pt\longrightarrow\hskip
8.53581pt\psi_{a}(x,y,z)=-(4\pi)^{2}\left(\frac{L_{x}+L_{y}+L_{z}}{L_{x}^{2}L_{y}^{2}L_{z}^{2}}\right)f_{a}(x,y,z)$
(13)

For homogeneous Neumann BCs, a similar expression to (12) is used however with
the cosine expression $C(x,L_{x})$. The analytical solution corresponding to
(13) then follows easily. The error for the pseudo-spectral (PS) and second
order finite difference (FD2) eigenvalues are shown in Fig. 2.4.1. Results are
shown for both regular and staggered grid configurations. As expected,
application of the PS eigenvalues to the spectral source distribution yields
spectral accuracy of the solution, as the error is seen to saturate at machine
precision. Application of the FD2 type eigenvalues results in the expected
quadratic convergence. For reference, the spectral 1D operator of Hejlesen is
shown for a 1D unbounded problem [21]. It is observed that the expected
behaviour is reproduced, independent of the grid configuration. The
corresponding plots for the 2D and 3D cases are also shown in Fig. 2.4.1.
Practically exactly the same behaviour is observed in these cases as with the
1D case. It is seen that the 3D unbounded case does not quite reach spectral
accuracy for the chosen grid size as the spatial resolution is lower than for
the 1D and 2D cases (due to increased memory demand). The 3D-vector type
solver results are not shown here as, by construction, these simply carry out
the same procedures as the 3D scalar solver type.

&
$2^{3}$$2^{4}$$2^{5}$$2^{6}$$2^{7}$$2^{8}$$2^{9}$$2^{10}$$2^{11}$$2^{12}$$N$Dirichlet
(SP)Neumann (SP)Periodic (SP)Unbounded (HEJ)Dirichlet (FD2)Neumann
(FD2)Periodic (FD2)
$2^{3}$$2^{4}$$2^{5}$$2^{6}$$2^{7}$$2^{8}$$2^{9}$$2^{10}$$2^{11}$$10^{-16}$$10^{-8}$$10^{0}$$N^{\frac{1}{2}}$$E_{\infty}$
$2^{3}$$2^{4}$$2^{5}$$2^{6}$$2^{7}$$2^{8}$$2^{9}$$2^{10}$$2^{11}$$N^{\frac{1}{2}}$
$2^{3}$$2^{4}$$2^{5}$$2^{6}$$2^{7}$$2^{8}$$2^{9}$$10^{-17}$$10^{-8}$$10^{1}$$N^{\frac{1}{3}}$$E_{\infty}$
$2^{3}$$2^{4}$$2^{5}$$2^{6}$$2^{7}$$2^{8}$$2^{9}$$N^{\frac{1}{3}}$

Figure 5: Error of 1D (top row), 2D (middle row) and 3D (bottom row) solver
types of staggered (left) and regular (right) grid configurations. The dashed
line indicates quadratic decay.

###  4.2 Bounded Solvers: Inhomogeneous Boundary Conditions

The option also exists within SailFFish to treat inhomogeneous BCs for the
Dirichlet and Neumann bounded solvers. In this case, the form of the FD2
eigenvalue can be exploited to account for the desired boundary condition by
simply modifying the interior boundary nodes values as described in Section
2.4.2. This however excludes the use of the PS eigenvalue type. The accuracy
of these solution shall be investigated here. Inhomogeneous BCs are specified
by applying the following input function:

$\psi_{a}(x,L_{x})=a\left(1-\frac{x}{L_{x}}\right)^{2}+b\left(\frac{x}{L_{x}}\right)^{2}\hskip
5.69054pt\longrightarrow\hskip
5.69054ptf_{a}(x,L_{x})=2.0\frac{a+b}{L_{x}^{2}}\text{.}$ (14)

This function represents a paraboloid with the values $\psi_{a}(0,L_{x})=a$
and $\psi_{b}(L_{x},L_{x})=b$ on the lower and upper limits of the domain,
respectively. In the 3D case, a function is simply added for each dimension.
This results in a parabolic boundary distribution on each domain face, with
corner values specified by the constants $a_{x},b_{x},a_{y},\cdots$. This
function furthermore has the advantage that the gradient function is easily
calculated as:

$\frac{\partial\psi_{a}}{\partial
x}(x,L_{x})=2a\left(\frac{x}{L_{x}}-1\right)+2b\frac{x}{L_{x}}\text{ . }$ (15)

For test cases with Dirichlet BCs, (14) is therefore applied to calculate the
numerical value on the boundary. For test cases with Neumann BCs, (15) is
applied. As higher derivatives of this function vanish, even the FD2 type
eigenvalue leads to spectral accuracy. For this reason, the test function is
superimposed with the corresponding distribution from (12). For Dirichlet BCs,
$S(\cdot)$ is applied and for Neumann BCs, $C(\cdot)$ is applied, such that
the total BC on each face remains unchanged by this addition. The results for
2D and 3D cases are shown for a range of grid sizes in Fig. 2.4.1. The 1D case
is omitted for brevity, as the results are indistinguishable from the 2D and
3D cases. The error can be seen to decay quadratically as is expected for the
FD2 eigenvalue type. The results here are seen to mirror those for homogeneous
BCs in Fig. 2.4.1. It should be noted that if all boundaries have Neumann-type
BCs, the solution can only be determined up to a constant. An optimum solution
which produces a solution with the lowest residual may be employed. This
however is not implemented in SailFFish for simplicity. The results displayed
in Fig. 2.4.1 are therefore achieved by shifting the solution values such that
$\text{sup}\\{\psi\\}=\text{sup}\\{\psi_{a}\\}$.

&
$2^{3}$$2^{4}$$2^{5}$$2^{6}$$2^{7}$$2^{8}$$2^{9}$$10^{-4}$$10^{-2}$$10^{0}$$N^{\frac{1}{3}}$Dirichlet
(staggered)Dirichlet (regular)Neumann (staggered)Neumann (regular)

Figure 6: Error of 2D (left) and 3D (right) solvers with inhomogeneous BCs for
both staggered and regular grid configurations with FD2 type eigenvalues. The
dashed line indicates quadratic decay.

###  4.3 Unbounded Solvers

As with the bounded solvers, it advantageous to apply a test function for
which the boundary values are finite, however display $C^{\infty}$ smoothness.
For the unbounded problems, a bump function has been applied. This is
described as:

$\psi_{a}(\zeta)=\begin{cases}e^{-cd}&\zeta<1\\\ 0&\zeta\geq
1\end{cases}\hskip 2.84526pt\rightarrow\hskip
2.84526ptf_{a}(\zeta)=\begin{cases}-\left(c^{2}\zeta^{2}d^{4}-2cd^{2}-8c\zeta^{2}d^{3}\right)e^{-cd}&\zeta<1\\\
0&\zeta\geq 1\end{cases}\text{ where: }d=\frac{1}{1-\zeta^{2}}$ (16)

and $c$ is a constant which was chosen for all cases to be 10. In the 1D case
$\zeta$ is chosen to be the absolute distance from the origin. In the 2D case
$\zeta$ is taken to be the radius from the origin in the $x$-$y$ plane with:
$\zeta=\sqrt{x^{2}+y^{2}}$ and the vorticity vector is aligned with the $z$
axis. In the 3D case a vortex ring of radius $r_{0}$ in the $y$-$z$ plane is
prescribed. The vorticity distribution in any plane $\theta=\tan^{-1}(z/y)=$
const. is described by (16) where $\zeta$ is the distance to the vortex ring
coordinate in that plane: $\zeta=\sqrt{(r-r_{0})^{2}+z^{2}}$ and the vorticity
vector is aligned with the tangent to the vortex ring at that position. Due to
the additional radial terms of the Laplacian in cylindrical coordinates, the
expressions for $f_{a}$ become more complex in 2D and 3D cases. These are
omitted here for brevity, however the expressions can be found in the
corresponding unbounded test function files in the library repository.
In the case of unbounded BCs, the Green’s function is first calculated in real
space, before being forward transformed to the spectral space. This procedure
is essentially independent of whether staggered or regular grid types are
applied (again, this is true for the currently implemented cases of uniform BC
types, the situation differs for mixed BCs [20]). For this reason, in the
following results only staggered results are shown for brevity. The results
are identical for regular grid types. As described in Section 2.5, a range of
kernels have been implemented for unbounded BCs including Gaussian
approximations from order 2-10 and the spectral kernel representation, all
attributable to Hejlesen [21]. For the 1D case only the spectral-type kernel
is employed. For brevity, 1D results are not shown here, however the spectral
convergence can be seen for the 1D cases in Fig. 2.4.1.

&
$2^{3}$$2^{4}$$2^{5}$$2^{6}$$2^{7}$$2^{8}$$2^{9}$$10^{-15}$$10^{-7}$$10^{1}$$N^{\frac{1}{3}}$SPECTG2G4G6G8G10

Figure 7: Error of 2D (left) and 3D (right) solvers with unbounded BCs for a
range of unbounded kernel types.

The results for the bump function in 2D and 3D are given in Fig. 2.4.1. It can
be seen that the convergence order of the different kernels follows precisely
the behaviour expected for the given kernel order. The grid size in the 3D
extends only up to $[512,512,512]$ due to the memory overhead of the Hockney-
Eastwood method and the available memory on the device employed. Despite this,
the convergence behaviour is precisely as expected.

&
$2^{3}$$2^{4}$$2^{5}$$2^{6}$$2^{7}$$2^{8}$$10^{-14}$$10^{-7}$$10^{0}$$N^{\frac{1}{3}}$SPECTG2G4G6G8G10

Figure 8: Error of 2D (left) and 3D (right) solvers with unbounded BCs with
the grad (2D) and curl (3D) differential operator for a range of unbounded
kernel types.

The correct functioning of the differential operators has been tested by
applying to the distributions described previously, however for the solution
of the curl Poisson equation:

$\nabla^{2}\vec{\psi}=\nabla\times\vec{f}\text{ .}$ (17)

In this case, the same solution procedure is carried out as for the standard
unbounded solver, however after convolution spatial differentiation is carried
out on the data in the spectral space via (10) in order to extract the
derivatives. This form of the equation may be applied, for example, when
extracting the velocity field from the vorticity field in fluid dynamics [31,
21]. In this case, analytical expressions also exist for the 2D and 3D bump
functions and are again omitted here for brevity. The 2D curl operator is
practically equivalent to the application of the gradient operator. As
compared to the previous unbounded tests cases, the application of these
gradient terms results in multiple output fields, corresponding to the
velocity components of the extracted field. In this case, $E_{\infty}$ is
applied to all velocity components to extract the maximum absolute error. The
results for the application of spectral differentiation to the bump function
in 2D and 3D are shown in Fig. 2.4.1. The error is again seen to scale with
the same behaviour as for the direct Poisson equation for the bump function.
This demonstrates the efficacy of the spectral differentiation implemented
with SailFFish.

##  5 Performance

The performance of the library is now investigated. An inspection of timings
is provided, scaling of the solver is discussed and finally a breakdown of the
solver timings is given for a range of simulation setups. As stated in Section
3.2.1, solver initialization only occurs once and thus the expense of this
step is of less consequence than the solver timings. For the following
discussion, steps 2-4 of the solver execution (see Section 3.2.2) are
considered for the timings, as the input and output steps are generally of low
computational expense and are also likely to be very application-specific.
These function serve only as an access point to the underlying data arrays.
The simulations carried out here were done so on a desktop computer with an
Intel(R) 12th Gen Core${}^{\text{TM}}$ i7-12700KF processor with a clock speed
of 3.60 GHz, 20 available threads and 64 GB of installed RAM. For simulations
utilizing the cuFFT_DataType, execution was carried out on an NVIDIA GeForce
RTW 3070 Ti graphics card with 8GB of on-device memory. All simulations were
carried out with single floating point accuracy. All data points shown were
gathered by averaging the execution times over ten simulation runs.

###  5.1 Execution Time

A plot of solver wall clock timings is given in Fig. 9. As mentioned
previously, these timings are the sum total of forward transformation,
convolution, and backward transformation. Results are shown for both
FFTW_DataType and cuFFT_DataType base classes. It is observed that for lower
$N$, activating multi-threading within FFTW has a negative effect due to the
additional overhead associated with the use OpenMP and the comparatively small
problem size. The advantages can however be seen for larger $N$. With the use
of eight threads, an improvement of an order of magnitude is observed as
compared to single-thread execution. It is also observed that the use of the
GPU for the calculation drastically improves performance. As compared to CPU
operation an improvement in calculation time of between one and two orders of
magnitude can be observed.

$10^{2}$$10^{3}$$10^{4}$$10^{5}$$10^{6}$$10^{7}$$10^{8}$$10^{2}$$10^{4}$$10^{6}$$10^{8}$$N$$t_{calc}$
[$\mu s$]FFTW (1 thread)FFTW (4 threads)FFTW (8 threads)CUDA Figure 9:
Calculation times in microseconds for a range of problem sizes using both FFTW
and CUDA datatypes. The solver type is a bounded 3D solver with periodic BCs.

As described in Section 3, it is desired to minimise solver communication
overheads. Ideally, memory allocation for the FFT data arrays should therefore
occur directly on the GPU and transfers between CPU and GPU should be avoided.
This limits directly the allowable problem size to the available memory on the
GPU. These limits may be extended with the use of half-floating point
precision, provided the GPU architecture allows for this.

###  5.2 Scalability

As SailFFish is configured for shared-memory execution, strong scaling shall
be discussed. The three components of the execution time are all easily multi-
threaded with OpenMP. The convolution procedure is nothing more than a vector-
vector multiplication, as are any of the spectral differentiation routines.
The scaling here is therefore device-relevant and constrained by the same
overhead limitation as described in the previous section. For the
FFTW_DataType, the forward and backward transform steps are nothing more than
a forward and backward execution within FFTW (e.g. fftw_execute_dft function)
and are therefore subject to the multi-threaded implementation of FFTW. It was
difficult to attain precisely linear scaling as there appeared to be a
significant overhead between the cases of a single and multiple threads. For
some problem sizes linear scaling was observed between 4-28 threads. For other
cases, a discontinuous jump in execution time resulted above a given number of
threads. Further investigations are required here to establish robust scaling
behaviour. In most cases however, as observed in Fig. 9, the use of multiple
threads greatly improves performance. Similar investigations are not possible
for the cuFFT_DataType, as specifying the number of active GPU threads cannot
be easily facilitated.

###  5.3 Breakdown of Execution Steps

In order to investigate the relative computational expense of the different
processes of the execution, these are broken down and relative execution times
are compared. A range of problem sizes and DataType were investigated, these
are summarised in Table 2. Two cases were considered with the curl operator,
as this is particularly relevant to the solution of e.g. (17). For the
(F-type) simulations carried out here, eight threads have been applied
concurrently for calculation. Similar behaviour was however observed that for
any choice of number of concurrent threads.

Table 2: Simulation Parameters for Fig. 10 Simulation | Datatype |  $N$ Nodes | BC Type | Operator
---|---|---|---|---
F1 | FFTW | $64^{3}$ | Periodic | —
F2 | FFTW | $128^{3}$ | Unbounded | CURL
F3 | FFTW | $512^{3}$ | Periodic | —
C1 | CUDA | $64^{3}$ | Periodic | —
C2 | CUDA | $128^{3}$ | Unbounded | CURL
C3 | CUDA | $512^{3}$ | Periodic | —

The results for these cases are shown in Fig. 10. A number of observations can
be made. Convolution procedures generally consume only a small portion of the
total execution time, including cases where a differential operator is being
applied in the spectral space. This portion increases for smaller problem
sizes. The convolution times of the cuFFT_DataType are generally smaller. This
is likely due to the (adapted) use of a highly optimised matrix-vector
multiplication algorithm for convolution procedures. It is also observed that
an asymmetry occurs for the cuFFT_DataType. The causes for this require
further investigation.

F1F2F3C1C2C3$0$$0.2$$0.4$$0.6$$0.8$ProportionForward FFTConvolution
OperationsBackward FFT Figure 10: Proportion of calculation time for a range
of problem sizes, datatype and differential operators. For solver parameters,
see Table 2.

##  6 Conclusions & Outlook

SailFFish is a lightweight, parallelised library for the solution of the
Poisson equation on a rectangular grid with regular grid spacing. The solver
is constructed with the goals of simplicity and adaptability and is optimised
for execution within shared-memory applications. The solver is capable of
handing bounded Poisson problems whereby Dirichlet, Neumann or periodic
boundary conditions are specified on the domain boundary. Spectral eigenvalue
specification allows for the solution of the Poisson equation with spectral
accuracy for homogeneous boundary conditions. The use of the second-order
finite difference eigenvalues additionally allows for the treatment of
arbitrary inhomogeneous boundary conditions of Dirichlet or Neumann type.
Unbounded boundary conditions can also be treated within SailFFish through
application of the Hockney-Eastwood method of domain doubling. In this case a
range of Green’s function are available. Currently only uniform boundary
conditions are available, however it is desired to implement mixed boundary
conditions in the near future. It is hoped that in the future an optimised
James-Lacker algorithm may be applied to reduce memory overhead for the case
of unbounded boundary conditions. All solver types can be executed with either
regular or staggered grid configurations. All solver configuration were
validated against analytical solutions to the Poisson equation and were found
to converge at the expected rate. A description of the solver framework was
given to provide insight for future adaptation of the solver. The solver is
currently implemented with two FFT frameworks for transforms and operations
within the spectral space:

  * •

FFTW_DataType, whereby the FFTW library is leveraged for optimised calculation
on a CPU; and

  * •

cuFFT_DataType, whereby the cuFFT library is leveraged for optimised
calculation on a GPU.

It is hope that within an open-source context, future DataType formats will be
implemented as new, optimised libraries become available. The development of
SailFFish is a first step towards the implementation of a GPU-accelerated
solver for the vorticity transport equation.

## Acknowledgments

The author would like to acknowledge Denis-Gabriel Caprace for many helpful,
interesting and enjoyable conversations on the topic of fast Poisson solvers
and their application.

## References

  * [1] A. Gholami, D. Malhotra, H. Sundar, and G. Biros.  FFT, FMM, or multigrid? a comparative study of state-of-the-art poisson solvers for uniform and nonuniform grids in the unit cube.  SIAM Journal on Scientific Computing, 38(3):C280–C306, jan 2016\.
  * [2] A. Brandt.  Guide to multigrid development.  In Multigrid Methods, pages 220–312, Köln-Porz, 11 1982.
  * [3] A Brandt and A.A Lubrecht.  Multilevel matrix multiplication and fast solution of integral equations.  Journal of Computational Physics, 90(2):348–370, 1990.
  * [4] Van Emden Henson.  Multigrid methods nonlinear problems: an overview.  In Charles A. Bouman and Robert L. Stevenson, editors, Computational Imaging, volume 5016, pages 36 – 48. International Society for Optics and Photonics, SPIE, 2003.
  * [5] C.H. Venner and A.A. Lubrecht.  Multilevel Methods in Lubrication, volume 37.  Elsevier Science, 1993.
  * [6] J.C. Adams.  Mudpack: Multigrid portable fortran software for the efficient solution of linear elliptic partial differential equations.  Appl. Math. & Comp., 34(2):113–146, 1989.
  * [7] Brandt. A.  Multiscale scientific computation; review 2001.  In T.J. Barth, T. Chan, and R. Haimes, editors, Multiscale and Multiresolution Methods: Theory and Applications, pages 3–95. Springer Berlin, Heidelberg, 1 edition, 2002.
  * [8] J. Carrier, L. Greengard, and V. Rokhlin.  A fast adaptive multipole algorithm for particle simulations.  SIAM Journal on Scientific and Statistical Computing, 9(4):669–686, 1988.
  * [9] L. Greengard and V. Rokhlin.  A new version of the fast multipole method for the laplace equation in three dimensions.  Acta Numerica, 6:229–269, 1997.
  * [10] E. Kreyszig.  Advanced Engineering Mathematics.  John Wiley & Sons, 10 edition, 7 2010.
  * [11] Y. Ying, G. Biros, and D. Zorin.  A kernel-independent adaptive fast multipole algorithm in two and three dimensions.  Journal of Computational Physics, 196(2):591–626, 2004.
  * [12] A. van Garrel, C.H. Venner, and H.W. Hoeijmakers.  Fast multilevel panel method for wind turbine rotor flow simulations.  In Proc. AIAA SciTech Forum, Grapevine, Texas, USA, 1 2017.
  * [13] J. Saverin, D. Marten, G. Pechlivanoglou, C.N. Nayeri, and C.O. Paschereit.  Multilevel simulation of aerodynamic singularity elements.  In Proc. AIAA SciTech Forum, Kissimee, Florida, USA, January 2018\.
  * [14] R. Cocle, G. Winckelmans, and G. Daeninck.  Combining the vortex-in-cell and parallel fast multipole methods for efficient domain decomposition simulations.  Journal of Computational Physics, 227(21):9091–9120, 2008.  Special Issue Celebrating Tony Leonard’s 70th Birthday.
  * [15] J. Saverin.  Multilevel vortex particle method for aerodynamic simulations.  PhD thesis, Technische Universität Berlin, Berlin, Germany, 11 2021.
  * [16] James W. Cooley and John W. Tukey.  An algorithm for the machine calculation of complex fourier series.  Mathematics of Computation, 19:297–301, 1965.
  * [17] G. Pöplau and Daniel Potts.  Fast poisson solvers on nonequispaced grids: Multigrid and fourier methods compared.  Proceedings of SPIE - The International Society for Optical Engineering, 12 2003.
  * [18] V. Fuka.  PoisFFT – a free parallel fast poisson solver.  Applied Mathematics and Computation, 267:356–364, 2015.
  * [19] A. Wiegmann.  Fast poisson, fast helmholtz and fast linear elastostatic solvers on rectangular parallelepipeds.  Technical report, Berkeley lab, 6 1999.  LBNL-43565.
  * [20] Denis-Gabriel Caprace, Thomas Gillis, and Philippe Chatelain.  Flups: A fourier-based library of unbounded poisson solvers.  SIAM Journal on Scientific Computing, 43(1):C31–C60, Jan 2021\.
  * [21] M.M. Hejlesen.  A high order regularisation method for solving the Poisson equation and selected applications using vortex method.  PhD thesis, DTU Mechanical Engineering, 2016.
  * [22] M.M. Hejlesen, G. Winckelmans, and J.H. Walther.  Non-singular green’s functions for the unbounded poisson equation in one, two and three dimensions.  Applied Mathematics Letters, 89:28–34, 2019.
  * [23] Steven A. Orszag.  Comparison of pseudospectral and spectral approximation.  Studies in Applied Mathematics, 51(3):253–259, 1972.
  * [24] M. Frigo and S. G. Johnson.  The design and implementation of FFTW3.  Proceedings of the IEEE, 93(2):216–231, 2005.  Special issue on “Program Generation, Optimization, and Platform Adaptation”.
  * [25] M. Frigo and S. G. Johnson.  The fastest Fourier transform in the west.  Technical Report MIT-LCS-TR-728, Massachusetts Institute of Technology, September 1997.
  * [26] K. Lackner.  Computation of ideal mhd equilibria.  Comp. Phys. Comm., 12:33–44, 1976.
  * [27] R.A James.  The solution of poisson’s equation for isolated source distributions.  Journal of Computational Physics, 25(2):71–93, 1977.
  * [28] F. Vico, Greengard L., and M. Ferrando.  Fast convolution with free-space green's functions.  Journal of Computational Physics, 323:191–203, oct 2016.
  * [29] R.W. Hockney and J.W. Eastwood.  Computer Simulation using Particles.  CRC Press, 1988.
  * [30] J.W. Eastwood and D.R.K. Brownrigg.  Remarks on the solution of poisson’s equation for isolated systems.  Journal of Computational Physics, 32(1):24–38, 1979.
  * [31] Chatelain C. and P. Koumoutsakos.  A fourier-based elliptic solver for vortical flows with periodic and unbounded directions.  Journal of Computational Physics, 229(7):2425–2431, 2010.
  * [32] M.M. Hejlesen.  poisson solver.  https://github.com/mmhej/poissonsolver, 2019.
  * [33] OpenMP Architecture Review Board.  Openmp.  https://www.openmp.org/, 2021.
  * [34] M. Frigo and S. G. Johnson.  Fftw, 2021.
  * [35] NVIDIA Corporation.  cufft, 2022.
  * [36] Intel Corporation.  Intel mkl.  https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/fourier-transform-functions.html, 2022\.
  * [37] B. Natarajan and C. Robeck.  clfft.  https://github.com/clMathLibraries/clFFT, 2022.
  * [38] NumPy.  Numpy fft.  https://numpy.org/doc/stable/reference/routines.fft.html, 2022.