11 Poisson Solvers – IPPL Manual

11.1 Main classes

Class	Role
`Poisson`	Common Poisson solver interface (base class).
`FFTPeriodicPoissonSolver`	FFT-based periodic solver for Poisson.
`FFTOpenPoissonSolver`	FFT-based open-boundary solver based on Hockney-Eastwood and a new Vico-Greengard spectral solver [1].
`FFTTruncatedGreenPeriodicPoissonSolver`	An FFT-based solver for long-range interaction only, using a truncated Green function. This is needed in the P3M approach to include collisions.
`PoissonCG`	(Preconditioned) Conjugate Gradient solver based on a finite differences discretization of the Poisson equation.
`FEMPoissonSolver`	Matrix-free FEM-based Poisson solver (using Conjugate Gradient as the iterative solver).
`PreconditionedFEMPoissonSolver`	FEM Poisson solver with preconditioning (using PCG).

For more information on Poisson solvers in IPPL, see [2].

11.2 Solver hierarchy

classDiagram
  class Poisson {
    +setLhs(lhs)
    +setRhs(rhs)
    +setGradient(grad)
    +updateParameter(key, value)
    +mergeParameters(params)
    +solve()
  }
  class FFTPeriodicPoissonSolver {
    +setRhs(rhs)
    +solve()
  }
  class FFTOpenPoissonSolver {
    +setRhs(rhs)
    +solve()
    +greensFunction()
  }
  class FFTTruncatedGreenPeriodicPoissonSolver {
    +setRhs(rhs)
    +solve()
    +greensFunction()
  }
  class PoissonCG {
    +setRhs(rhs)
    +setSolver(lhs)
    +solve()
  }
  class FEMPoissonSolver {
    +setRhs(rhs)
    +solve()
  }
  class PreconditionedFEMPoissonSolver {
    +setRhs(rhs)
    +solve()
  }

  Poisson <|-- FFTPeriodicPoissonSolver
  Poisson <|-- FFTOpenPoissonSolver
  Poisson <|-- FFTTruncatedGreenPeriodicPoissonSolver
  Poisson <|-- PoissonCG
  Poisson <|-- FEMPoissonSolver
  Poisson <|-- PreconditionedFEMPoissonSolver

The base Poisson<FieldLHS, FieldRHS> class stores LHS, RHS, and solver parameters. One can also pass a solver output type as a parameter, which is an enum:

Output type	Meaning
`SOL`	Get the solution of the Poisson equation, in-place (replaces the rhs with the solution).
`GRAD`	Get the gradient of the potential only i.e. the electric field (in the electrostatic case). Note that this returns the negative gradient ($\mathbf{E} = -\nabla\phi$).
`SOL_AND_GRAD`	Return both solution and gradient.

11.3 FFT solver parameters

FFT-based Poisson solvers inherit the FFT parameter set (see Chapter 10):

params.add("use_heffte_defaults", false);
params.add("use_pencils", true);
params.add("use_reorder", false);
params.add("use_gpu_aware", true);
params.add("comm", ippl::a2av);
params.add("r2c_direction", 0);

For the FFTOpenPoissonSolver, one can also choose the algorithm type:

Solver	Parameter	Values
`FFTOpenPoissonSolver`	`algorithm`	`HOCKNEY`, `VICO`, `BIHARMONIC`, or `DCT_VICO`.
`FFTOpenPoissonSolver`	`greens_function`	`STANDARD` or `INTEGRATED`; only used by `HOCKNEY`.

The HOCKNEY algorithm is the standard Hockney-Eastwood method for solving the Poisson equation with open boundary conditions [3]. The VICO algorithm is a novel method with spectral accuracy properties, described in [1], and the DCT_VICO is its memory-optimized counterpart. We suggest the user to use DCT_VICO as the algorithm of choice, as it has similar performance and memory footprint as the Hockney-Eastwood method, with the added benefit of higher-accuracy.

11.3.1 Hockney Green-function models

The HOCKNEY open-boundary solver can build either a point-sampled or a cell-integrated real-space Green function before transforming the doubled-grid convolution kernel. The selection is made with the integer greens_function parameter:

Value	Meaning
`FFTOpenPoissonSolver::STANDARD`	Point-sampled free-space kernel. This is the default and preserves historical IPPL behavior.
`FFTOpenPoissonSolver::INTEGRATED`	Cell-averaged kernel following Qiang et al. [4], [5]. This is currently implemented for 3D `HOCKNEY` only.

The standard kernel samples \[ G(\mathbf{r}) = -\frac{1}{4\pi |\mathbf{r}|}. \]

The integrated kernel assumes the charge density is uniform over each source cell and uses the cell average \[ \overline{G}(\mathbf{r}) = -\frac{1}{4\pi h_x h_y h_z} \int_{-h_x/2}^{h_x/2} \int_{-h_y/2}^{h_y/2} \int_{-h_z/2}^{h_z/2} \frac{dx' dy' dz'} {\sqrt{(x-x')^2 + (y-y')^2 + (z-z')^2}}. \]

IPPL evaluates this integral by the closed-form eight-corner difference of the antiderivative \[ \begin{aligned} F(x,y,z) ={}& -\frac{z^2}{2}\tan^{-1}\!\left(\frac{xy}{zr}\right) -\frac{y^2}{2}\tan^{-1}\!\left(\frac{xz}{yr}\right) -\frac{x^2}{2}\tan^{-1}\!\left(\frac{yz}{xr}\right) \\ &+ yz\ln(x+r) + xz\ln(y+r) + xy\ln(z+r), \end{aligned} \] where $r=\sqrt{x^2+y^2+z^2}$. This removes the point singularity from the cell containing the source and is useful for high-aspect-ratio beams where point sampling of $1/r$ can be inefficient near narrow dimensions.

A minimal setup is:

using Solver_t = ippl::FFTOpenPoissonSolver<VField_t, Field_t>;

ippl::ParameterList params;
params.add("algorithm", Solver_t::HOCKNEY);
params.add("greens_function", Solver_t::INTEGRATED);

Solver_t solver;
solver.mergeParameters(params);
solver.setRhs(rho);
solver.solve();

Calling shiftedGreensFunction(shift) replaces the cached Hockney kernel by the shifted kernel $G(\mathbf{r}-\mathbf{s})$. With greens_function == INTEGRATED, the same shift is applied before the cell average is evaluated. This is intended for workflows such as image-charge corrections, where the caller constructs the image contribution around an IPPL solve with a shifted kernel. The shifted kernel remains active until greensFunction() is called again or the mesh spacing changes.

11.4 (Preconditioned) Conjugate Gradient solver

PoissonCG solves the Poisson problem with a finite-difference Laplace operator. The Poisson equation is discretized using a second-order centered finite differences scheme, which in the 1D case looks like:

\[ \frac{\phi_{i+1} - 2\phi_i + \phi_{i-1}}{\Delta x^2} = - \rho_i \hspace{5mm} \forall i \in \{0,..., n-1\}. \] This is also known as the 3-point stencil. In the 2D case, this becomes a 5 point stencil, while in the 3D case, it becomes a 7D stencil.

When applying this to every grid point of the field, we obtain a linear system of equations of type $A \mathbf{\phi} = \mathbf{\rho}$, where $A$ is an $n\times n$ sized matrix: \[ \underbrace{\frac{1}{\Delta x^2}\begin{bmatrix} 2 & -1 & \\ -1 & 2 & -1 & \\ & -1 & \ddots & \ddots &\\ & & \ddots & \ddots & -1 \\ & & & -1 & 2 \\ \end{bmatrix}}_\text{$A$} \cdot \underbrace{\begin{bmatrix} \phi_0 \\ \phi_1 \\ \vdots \\ \phi_{n-2}\\ \phi_{n-1} \\ \end{bmatrix}}_\text{$\vec{\phi}$} = \underbrace{\begin{bmatrix} \rho_0 \\ \rho_1 \\ \vdots \\ \rho_{n-2}\\ \rho_{n-1} \\ \end{bmatrix}}_\text{$\vec{\rho}$}. \]

Instead of building the matrix, one can solve this linear system of equations in a manner by defining an operator for the action of the matrix $A$ on a vector $\vec{x}$: \[ [F_A(\vec{x})]_j = \frac{\vec{x}_{j+1} - 2\vec{x}_j + \vec{x}_{j-1}}{h_x^2}. \]

This is what is done in PoissonCG, as building the matrix is an expensive task, especially when the number ofgrid points in the simulation is large. By defining a Laplace operator, one can simply pass the action of the matrix on a vector to the CG iterative solver, which uses this operator to simply compute any needed matrix-vector products on the fly.

This solver can either be run using plain CG, or with a preconditioner to reduce the condition number of the matrix and therefore make the iterative solver converge faster (PCG).

The preconditioners which are currently available are: Jacobi, Richardson, Gauss-Seidel, and SSOR (see [6], [7], [8] for more information on these preconditioners and their implementation).

The preconditioner can be set using the solver ParameterList. A list of the input parameters which one can pass are:

Key	Meaning
`solver`	`non-preconditioned` for CG or `preconditioned` for PCG.
`preconditioner_type`	`jacobi`, `richardson`, `gauss-seidel`, `ssor`
`richardson_iterations`	Number of Richardson preconditioner iterations.
`gauss_seidel_inner_iterations`	Inner Gauss-Seidel or SSOR iterations.
`gauss_seidel_outer_iterations`	Outer Gauss-Seidel or SSOR iterations.
`ssor_omega`	SSOR relaxation parameter.
`communication`	This is a bool (`0` or `1`) which when enabled turns on halo communication in the preconditioners.
`max_iterations`	Iteration cap for the linear solve.
`tolerance`	Residual tolerance.

TestCGSolver verifies convergence against a manufactured periodic solution. Its command-line parser accepts a grid-size exponent and optional scaling/preconditioner selectors; for example j selects Jacobi, r selects Richardson, g selects Gauss-Seidel, and s selects SSOR. The test prints relative error, residual, and iteration count.

11.5 FEMPoissonSolver

The FEM-based Poisson solver is based on the FEM framework in IPPL (src/FEM). It employs a first-order Lagrange space to discretize the Poisson equation in its variational form:

\[ \int_{\Omega} \nabla \phi \cdot \nabla v = \int_{\Omega} \rho \cdot v, \hspace{4mm} \forall v\in V. \] where $v$ is known as the test function, and $\phi$ the trial function.

When applying a finite element discretization to this variational form, one obtains a linear system of equations $A\mathbf{x} = \mathbf{b}$, where $A$ is known as the stiffness matrix and $\mathbf{b}$ as the load vector. Assembling the stifness matrix is very expensive, therefore we choose to employ a matrix-free approach, where the product $A\mathbf{x}$ is computed on-the-fly. More information on this approach can be found in [2].

11.6 Minimal working example of a solver

FFTPeriodicPoissonSolver solves periodic Poisson problems using real-to-complex FFTs. test/solver/TestFFTPeriodicPoissonSolver.cpp verifies the solver against a smooth manufactured solution on [-1,1]^3 over multiple grid sizes.

Minimal pattern:

using Field_t  = ippl::Field<double, Dim, Mesh_t, Centering_t>;
using VField_t = ippl::Field<ippl::Vector<double, Dim>, Dim, Mesh_t, Centering_t>;
using Solver_t = ippl::FFTPeriodicPoissonSolver<VField_t, Field_t>;

Field_t rhs(mesh, layout);

ippl::ParameterList params;
params.add("output_type", Solver_t::SOL);
params.add("use_heffte_defaults", false);
params.add("use_pencils", true);
params.add("use_gpu_aware", true);
params.add("comm", ippl::a2av);
params.add("r2c_direction", 0);

Solver_t solver;
solver.mergeParameters(params);
solver.setRhs(rhs);
solver.solve();
# rhs should now contain the solution of the Poisson equation.

With output_type == SOL, the RHS field is overwritten with the potential in the current test pattern.

11.7 Verification material

The correctness and convergence of all solvers are tested in test/solver and test/solver/fem(for the FEM-based solver). All test files have instructions on how to run them at the beginning of the file.