The isolated Poisson solver

ID	Group	Milestone	Backends	Figure	Regression
SV-6	Self-gravity	M8 (#81), P4	CPU, FNL	yes	selfgravity-star (through the star it solves for)

Case and scope

On a periodic domain a self-gravitating object feels its own images, and the mean density must be removed. Isolated boundaries give the free-space potential, $\mathrm{\Phi }\to 0$ at infinity, of the mass on the grid: what a star, a collapsing cloud or a binary needs. P4 adds them to the FFT solver and verifies the solver (SV-6, this page), a star in equilibrium (SV-7) and a collapsing one (SV-8, both on the n = 1 polytrope).

The method (poisson_bc = isolated). Hockney and Eastwood's zero padding: the density sits in an array of $m=2(N+n_{gc})$ cells per direction, zero outside the domain, and is convolved, by one r2c transform, a product and one c2r, with the lattice Green's function of $L_h$, $h^{2}G(\mathbf{n})$. The circular convolution equals the free-space one on every cell the scatter reads, ghosts included, because $m$ exceeds twice the largest offset. Using the lattice kernel of $L_h$ rather than the sampled continuum $-1/(4\pi r)$ makes $\mathrm{\Phi }$ satisfy the discrete Poisson equation, not just approximate it.

The kernel. $G$ is generated once, offline, by scripts/lattice_green.py into adam_lattice_green_library (6545 values, the wedge $n_1\ge n_2\ge n_3\ge 0$):

• near field, $max|n_d|\le 32$: the Dirichlet problem $L_{1}G=\delta$ on the cube $[-32,32{]}^3$, solved exactly by a type-I sine transform, with the far-field series on its faces;
• far field: the series $G_0+G_1+G_2+G_3$, $G_0=-1/(4\pi r)$, $G_1=-5(n_1^4+n_2^4+n_3^4-\frac{3}{5}{r}^4)/(32\pi r^7)$, $G_2=P_8/r^{13}$, $G_3=P_{12}/r^{19}$, each term solving its order of the Taylor series of $L_1$ (checked symbolically). The series leaves one coefficient free in $G_2$ and one in $G_3$ (the cubic harmonics of degree 4 and 6); the lattice fixes them, and they are fitted to the table (least squares over $12\le r\le 26$).

What the case catches:

• a kernel that is not the lattice Green's function of $L_h$ (a sampled continuum kernel, a bad table/series seam): the residual of $L_h$ leaves the kernel's level;
• a padding too small for the offsets the scatter reads: the error against the exact potential on the ghosts and the padded far field;
• a wrong mass normalisation or sign: the error against the exact potentials outside the source ($-GM/r$);
• a solver that accepts a configuration it cannot solve: the scope is fatal otherwise (3-D, every direction non-periodic, cubic cells: the kernel is the unit lattice's).

Models and numerics

Item	Choice	Code
Equation	$L_h\mathrm{\Phi }=4\pi G\rho$, $\mathrm{\Phi }\to 0$ at infinity	models
Solver, CPU	FFTW, zero padding to $m=2(N+n_{gc})$, convolution with the kernel's spectrum	initialize_blocks, make_kernel, transform, transform_kernel, adam_poisson_fft_object.F90
Solver, FNL	cuFFT on rank 0's GPU	transform_kernel_fnl, spectral_convolve_dev, kernel_spectrum_dev, adam_fnl_poisson_fft_object.F90
Kernel	the lattice Green's function: table for $max|n_d|\le 32$ (LGF_RADIUS), series $G_0+G_1+G_2+G_3$ beyond	adam_lattice_green_library.F90, generated by scripts/lattice_green.py
Selection	[gravity] poisson_bc = isolated	load_from_file, adam_flume_gravity_object.F90

Initial and boundary conditions

Item	Value
Domain	the cube $[-1,1{]}^3$, non-periodic
Grid	octree levels 1–3 (${16}^3$ to ${64}^3$ cells), 8 cells per block side, 3 ghosts
Ranks	CPU: 1, 2, 3; FNL: 1, 2
Factor	$4\pi G$ with $G=1$
Source 1: the n = 1 polytrope	$\rho =\sin (kr)/(kr)$, $k=\pi /a$, $a=1/2$ (smooth; its derivative jumps at $r=a$); exact $\mathrm{\Phi }=-(4\pi /k^2)\sin (kr)/(kr)-M/a$ inside, $-M/r$ outside, $M=4a^3/\pi$
Source 2: the uniform sphere	$\rho =1$, $r<a$, as cell averages ($8^3$ samples per cell); exact $-2\pi (a^2-r^2/3)$ inside, $-M/r$ outside, $M=4\pi a^3/3$
Refusals	--refuse-isolated-periodic (a periodic grid), --refuse-isolated-anisotropic (cells $1/8\times 2/8\times 3/8$), --refuse-isolated-2d (null $z$): each must stop with exit 1 and its message

Generator checks (deterministic, scripts/lattice_green.py):

Check	Bound	Measured
residual of $L_1$ on the table	${10}^{-13}$	$2.2\cdot {10}^{-16}$
residual of $L_1$ on the composite kernel (table and series, as the Fortran evaluates it) on $[-64,64{]}^3$, across the faces where they meet	${10}^{-12}$	$1.0\cdot {10}^{-13}$
$G(\mathbf{0})$ against Watson's closed form for the simple cubic lattice, relative	${10}^{-12}$	$2.0\cdot {10}^{-13}$
$G(1,0,0)-G(\mathbf{0})-1/6$	${10}^{-14}$	$2.8\cdot {10}^{-17}$

Input files

File	Role
src/tests/poisson/test_poisson_fft.F90	the unit test, --isolated (and the three --refuse-isolated-* flags)
src/tests/flume/verification/selfgravity/poisson.sh	leg sv6: the test on every rank count of --np-list, and the refusals; --dump DIR covers SV-6 too
src/tests/flume/verification/selfgravity/poisson_oracle.py	compare DIR_CPU DIR_FNL on the padded arrays
scripts/lattice_green.py	generates adam_lattice_green_library.F90 and runs the generator checks
src/tests/flume/verification/make_doc_runs.sh, group isolated	sv6-poisson: the test with its log and dumps
src/tests/flume/verification/make_doc_figures.py, fig_selfgravity_isolated	draws the figure (--only selfgravity-isolated)

python3 scripts/lattice_green.py --check            # generator checks; the committed module against a fresh emission
cd src/tests/flume/verification/selfgravity
./poisson.sh --legs sv6                                                                   # CPU (FFTW)
POISSON_EXE=$PWD/../../../../../exe/test_poisson_fft_fnl ./poisson.sh --legs sv6 --np-list "1 2"  # FNL (cuFFT)

Results

Left: the error against the exact potentials falls as $\mathrm{\Delta }{x}^2$ (the sphere less regularly); the residual sits at the kernel's level. Middle: the polytrope's potential along the line through the centre over the whole padded array, on the exact curve, $-GM/r$ outside the star. Right: the error on the plane through the centre, $6.1\cdot {10}^{-4}$ at most, at the centre; smooth, with no trace of the domain's edge or of the padding.

Asserted	Bound	CPU (FFTW)	FNL (cuFFT, A30 and the WSL GPU)
residual of $L_h$, relative to $max|4\pi G\rho |$	${10}^{-10}$	$\le 2.1\cdot {10}^{-11}$	$\le 2.1\cdot {10}^{-11}$
least-squares order against the exact $\mathrm{\Phi }$, every cell and ghost, polytrope	$\ge 1.9$	2.03	2.03
the same, uniform sphere	$\ge 1.2$	1.35	1.35
ghost cells left unfilled	0	0	0
isolated on a periodic grid, with non-cubic cells, in 2-D	refused, exit 1	refused	refused
CPU against FNL (poisson_oracle.py compare, the padded arrays)	${10}^{-13}$	$\le 1.6\cdot {10}^{-15}$	$\le 1.6\cdot {10}^{-15}$

Related cases: the periodic solver is SV-2; the stars it holds and collapses are SV-7 and SV-8; isolated boundaries on AMR are the multigrid's James step (#87, SV-13, SV-14).

Findings

• The kernel series needs $G_2$ and $G_3$. With $G_0+G_1$ alone the composite residual was $7\cdot {10}^{-10}$, concentrated on the cube's faces, and summed over a ${64}^3$ source it showed as a residual of $3.4\cdot {10}^{-10}$ in $\mathrm{\Phi }$; with $G_2$ and $G_3$ they are $1.0\cdot {10}^{-13}$ and $2.1\cdot {10}^{-11}$.
• The residual is the kernel's, not round-off. It is ${10}^{-14}$ while every offset lies inside the table (${16}^3$, ${32}^3$) and $2\cdot {10}^{-11}$ at ${64}^3$, where the offsets reach the series: about ${10}^{-13}$ per offset, summed over the source cells. At ${128}^3$ (measured once, 6 ranks) it is $5\cdot {10}^{-11}$.
• The polytrope converges at second order (2.14, 1.92, 2.02 per pair from ${16}^3$ to ${128}^3$; the dip is the kink of its density at the surface). The uniform sphere converges irregularly (0.6, 2.1, 1.75): its surface cuts cells, and the bound only asserts better than first order.
• Conservation needs symmetry, not a round-off residual. The convolution with an even kernel is a symmetric linear map, which is what the exact exchange of $E$ and $W$ needs (models): SV-8 measures it.
• The error is largest at the centre, where the potential is deepest.