The FFT Poisson solver

ID	Group	Milestone	Backends	Figure	Regression
SV-2	Self-gravity	M8 (#81), P2	CPU, FNL	yes	—

Case and scope

Self-gravity needs a Poisson solver that is exact for the discrete equation $L_h\mathrm{\Phi }=4\pi G(\rho -\overline{\rho })$, with $L_h$ the 7-point Laplacian the gravity stencils telescope against, so that the scheme is the one P0 analysed. P2 wrote this as the condition of exact conservation; P4 found the condition weaker: any linear, symmetric solve conserves $E+\frac{1}{2}\rho \mathrm{\Phi }$ and momentum exactly in the semi-discrete form (see SV-6). The round-off residual remains a property of this solver, verified here.

What the case catches (the unit test test_poisson_fft):

• a wrong eigenvalue of $L_h$ (a missing factor, a swapped axis, a wrong spacing): the error against the discrete solution leaves round-off. The domain is anisotropic, $[0,1]\times [0,2]\times [0,3]$, so an axis mix-up cannot cancel;
• a misplaced block in the gather, or a wrong ghost wrap in the scatter: the error over every cell and ghost, and the residual from the scattered ghosts;
• a ghost layer the scatter leaves unfilled (a sentinel value);
• a solver that accepts a grid it cannot solve: a non-periodic grid must stop the run.

Models and numerics

Item	Choice	Code
Equation	$L_h\mathrm{\Phi }=4\pi G(\rho -\overline{\rho })$, $L_h$ the 7-point Laplacian (5-point in 2-D)	models
Solver, CPU	FFTW. Rank 0 gathers the block interiors and transforms; each Fourier mode is divided by the eigenvalue of $L_h$ itself, ${\lambda }_k=\sum _d(2\cos (2\pi k_d/N_d)-2)/\mathrm{\Delta }{x}_d^2$; the mean (${\lambda }_k=0$) is set to zero; the inverse transform follows, then the scatter fills every block with all its ghost layers by the periodic wrap	initialize, solve, transform, adam_poisson_fft_object.F90
Solver, FNL	cuFFT on rank 0's GPU, the same gather and scatter	transform_fnl, spectral_solve_dev, adam_fnl_poisson_fft_object.F90
Scope	uniform, periodic grids; anything else stops the run. The gather makes the solve serial on rank 0: correct for any rank count, not scalable. Isolated boundaries come in P4 (SV-6), AMR in P5 (SV-9)	initialize, adam_poisson_fft_object.F90

Initial and boundary conditions

Item	Value
Domain	$[0,1]\times [0,2]\times [0,3]$, periodic on every face
Source	$\rho =1.5+\sin (k_{1}x)\cos (k_{1}y)\cos (k_{1}z)+0.5\cos ({\mathbf{k}}_2\cdot \mathbf{x})$, ${\mathbf{k}}_1=2\pi (1,1,\frac{1}{3})$, ${\mathbf{k}}_2=2\pi (2,\frac{1}{2},\frac{1}{3})$; in 2-D the z wave numbers are 0
Factor	$4\pi G$ with $G=1$
Grids	3-D octree at levels 1–3 (${16}^3$ to ${64}^3$ cells); 2-D quadtree (null_z) at levels 1–3 (${16}^2$ to ${64}^2$); a single $8^3$ block (level 0: the wrap inside one block; on 2+ ranks, ranks without blocks). 8 cells per block side, 3 ghosts
Ranks	CPU: 1, 2, 3 (the poisson.sh default list), and 4; FNL: 1, 2
Exact solutions	each mode is an eigenfunction of $L_h$: the discrete solution divides it by its ${\lambda }_k$, the continuum one by $-|\mathbf{k}{|}^2$
Refusal	--no-periodic: one 3-D configuration without periodic boundaries; the run must stop with exit 1, naming the direction (direction 1 is not periodic)

Input files

File	Role
src/tests/poisson/test_poisson_fft.F90	the unit test (modes test-poisson-fft-gnu, test-poisson-fft-fnl-nvf); --dump DIR writes rank 0's potential of every level
src/tests/flume/verification/selfgravity/poisson.sh	runs the test on every rank count of --np-list (leg sv2), and the refusal; POISSON_EXE selects the build
src/tests/flume/verification/selfgravity/poisson_oracle.py	compare DIR_CPU DIR_FNL: the CPU = FNL check on the dumped potentials, and the figure's reader
src/tests/flume/verification/make_doc_runs.sh, group selfgravity	sv2-poisson: the test with its log and dumps
src/tests/flume/verification/make_doc_figures.py, fig_selfgravity_poisson	draws the figure (--only selfgravity-poisson)

cd src/tests/flume/verification/selfgravity
./poisson.sh --legs sv2 --dump <dir-cpu>                                    # CPU (FFTW)
POISSON_EXE=$PWD/../../../../../exe/test_poisson_fft_fnl ./poisson.sh --legs sv2 --np-list "1 2" --dump <dir-fnl>
python3 poisson_oracle.py compare <dir-cpu> <dir-fnl>                       # CPU = FNL

Results

Left: against N, the error against the continuum solution falls as $\mathrm{\Delta }{x}^2$, while the error against the discrete solution and the residual stay at round-off. Middle and right: the 2-D potential at $N=64$ and its error against the continuum solution, $4.3\cdot {10}^{-4}$, largest where $|\mathrm{\Phi }|$ is.

Asserted	Bound	CPU (FFTW)	FNL (cuFFT, A30 and the WSL GPU)
$max|\mathrm{\Phi }-{\mathrm{\Phi }}_h|/max|{\mathrm{\Phi }}_h|$, every cell and ghost	${10}^{-12}$	$\le 2.8\cdot {10}^{-15}$	$\le 2.8\cdot {10}^{-15}$
residual of $L_h$ from the scattered ghosts, relative	${10}^{-10}$	$\le 3.7\cdot {10}^{-14}$	$\le 4.3\cdot {10}^{-14}$
order against the continuum solution, $N=32\to 64$	$\ge 1.9$	1.98 (3-D), 1.98 (2-D)	1.98, 1.98
ghost cells the scatter left unfilled	0	0	0
a non-periodic grid	refused, exit 1	refused	refused
CPU against FNL, $max|{\mathrm{\Phi }}_{\text{CPU}}-{\mathrm{\Phi }}_{\text{FNL}}|/max|\mathrm{\Phi }|$ (poisson_oracle.py compare)	${10}^{-13}$	$\le 9.8\cdot {10}^{-15}$	$\le 9.8\cdot {10}^{-15}$

Related cases: the isolated variant of the same solver is SV-6; the multigrid solver on AMR grids, checked against this one on uniform grids, is SV-9; the coupling to the scheme is SV-3 to SV-5.

Findings

• The printed numbers do not depend on the rank count (np 1 and np 4 agree to every digit), and the A30 potential is bitwise that of the WSL GPU.
• The solver is exact for $L_h$, and the scheme is second order through $L_h$.
• The residual grows from ${10}^{-15}$ to $4\cdot {10}^{-14}$ between $N=16$ and 64: $L_h$ amplifies the round-off of $\mathrm{\Phi }$ by $(k\mathrm{\Delta }x{)}^{-2}$. That is the floor the conservation of P3 inherits.
• The solve is serial on rank 0 (the gather): correct on any rank count, not scalable.