adam_poisson_fft_object

ADAM, FFT Poisson solver on uniform grids.

Solves L_h phi = factor (f - mean(f)) for a cell field f (one variable of a conservative array q) on a grid whose blocks all sit at one level and whose active directions are all periodic. L_h is the second-order 7-point Laplacian (5-point in 2-D; a null direction contributes nothing). The solve divides each Fourier mode by the eigenvalue of L_h itself, lambda_k = sum_d (2 cos(2 pi k_d / N_d) - 2) / dx_d^2, so the discrete residual of the returned phi is at round-off; the modes with lambda_k = 0 (the mean, and every mode along a null direction) are set to zero.

Data path: the distributed slab transform (issue #84, E5d). Rank r owns the z-planes [zr(r), zr(r+1)) of the transform array (m cells per direction: the grid, its extensions, its padding). Every block plane goes to the owners of the target planes that map to it, and each rank builds its planes by the per-axis maps (tmap, tsgn: a target cell's source cell and sign, the images of extend_global and the zero padding). The planes are transformed in 2-D (r2c), transposed to z-lines (MPI_ALLTOALLV, the y-rows [yr(r), yr(r+1)) on rank r) and transformed in 1-D, divided by the eigenvalues (or multiplied by the kernel's spectrum, kept in the pencil layout), and back. Every block then fetches its potential with all its ghost layers from the owners of its planes (wrapped modulo m), so no ghost exchange follows; the energy sums on the owners. One plan for one plane and one for one row of lines, executed on every plane and row in turn through fixed buffers (slab_planes, slab_lines, the backend hooks): the arithmetic of a value does not depend on how many planes or rows a rank holds, so the potential is bitwise the same on every rank count (poisson.sh leg slab). At most m(3) ranks hold planes (a pencil layout would lift it). PRISM's CG solves (pcg, below), the multigrid's redundant root and the CG's sub-problems keep the rank-0 path of the first version (issue #81, D-81-2): rank 0 gathers the source (MPI_Gatherv), places it, transforms (transform) and scatters every block's phi (MPI_Scatterv); global_potential returns the transform array on rank 0 in either path (diagnostics).

Backends. This object transforms with FFTW 3 (-D_FFTW, -lfftw3; FFTW_ESTIMATE plans, deterministic). The FNL twin adam_fnl_poisson_fft_object extends it and overrides the slab hooks (slab_make_plans, slab_planes, slab_lines, slab_destroy_plans: each rank's GPU) and the rank-0 ones (make_plans, transform, destroy_plans) with cuFFT. Built without either, initialize stops: "built without an FFT library".

Entry points. solve packs one variable of a host conservative array; solve_packed takes the interior source already packed (block-major, i fastest: what a device kernel produces on FNL). Both can return energy = sum_cells f Phi (summed on the slab owners, every rank gets it): W = energy dV / 2 is the gravitational energy of self-gravity on a uniform grid.

Isolated boundaries (issue #81, M8-P4), Hockney and Eastwood. initialize(..., isolated=.true.) solves L_h phi = factor f with phi -> 0 at infinity (no mean removed): the source sits in a zero-padded array of m >= 2 (n + ngc) cells per direction, and the circular convolution with the lattice Green's function of L_h (adam_lattice_green_library, h^2 G(n)) equals the free-space one on every cell the scatter reads, the ghost layers included, since m exceeds the largest offset, n - 1 + ngc, twice. m is the smallest even size at or above 2 (n + ngc) with no prime factor above 13 (#84 E5b): the transform cost no longer depends on the factors of n + ngc (medooza, n = 128: 262 = 2 x 131 cells cost FFTW 3.3 times as long as 264; the kernel entries beyond the offsets the read cells need meet only the zero padding). The kernel's spectrum (khat, real: the kernel is even) is computed once. Inside the domain the discrete residual is that of the kernel, about 1e-11 of the source (scripts/lattice_green.py).

Homogeneous Dirichlet boundaries (issue #87, P1) are for the multigrid's coarsest level only (initialize_blocks(..., dirichlet=.true.)): phi = 0 at the centres of the first ghost layer. The source is extended odd about those centres to m = 2 (n + 1) cells per active direction, and the periodic transform of m cells solves it: the odd extension is the type-I sine transform, exact for the 7-point operator with zero ghosts.

Boundary kinds per face (issue #84, E1a). bc(2,3) (adam_flail_parameters) gives each face a homogeneous kind: an axis is periodic, or each of its faces is Dirichlet (phi = 0 at the face) or Neumann (zero normal derivative at the face), possibly different at the two ends. The source is extended odd about a Dirichlet face and even about a Neumann one, over 2n cells when the two faces agree and 4n when they differ (a quarter-wave extension), and the periodic transform of the extension solves the 7-point operator with those boundary conditions exactly; the scatter reads the ghost layers from the extension, so they carry the conditions. The modes of zero eigenvalue are dropped as on a periodic grid: without a Dirichlet face the source mean is removed (the operator is singular). The energy sums the interior only. Inhomogeneous values are the clients' (FLAIL's) work.

The wide operator D(G) (issue #84, E2a). dg_order = p (2, 4, 6, 8; 0, the default, is the compact 7-point Laplacian) solves PRISM's elliptic operator instead: the order-p centred first derivative D1 composed with itself per axis (the stencil reaches p cells, as PRISM's D and G read it), eigenvalue -sum_d s_p(theta_d)^2 / dx_d^2 with s_p(theta) = sum_j 2 a_j sin(j theta) (the a_j of D1). Every centred D1 annihilates the checkerboard (s_p(pi) = 0, set exactly: sin(pi) is 1e-16, not 0), so the modes with every theta_d in {0, pi} are in the kernel and dropped, Dirichlet faces or not: G of them is zero, so the gradient a client takes is unique. The odd and even extensions are reflections, which the wide stencil commutes with, so the kinds per face hold as for the compact operator (the scatter fills p ghost layers from the extension); PRISM's zero-ghost Dirichlet face is not a reflection (below). Measured on the prototype (wide_dg_proto.py): a compact preconditioner needs O(N) iterations for D(G) (the checkerboard branch), so the direct solve here is the method on uniform grids.

Zero-ghost faces (issue #84, E2b). FLAIL_BC_ZERO_GHOSTS (both faces of an axis) is PRISM's elliptic Dirichlet face: every ghost layer zero. With the compact operator it is the ghost-Dirichlet kind (the 7-point stencil reads one layer). With the wide D(G) it is no reflection: the operator is the stencil truncated, symmetric and nonsingular, and no transform diagonalises it. Rank 0 then solves it, on the gathered array, by conjugate gradients preconditioned by the wide p = 2 operator with the same faces, inverted exactly: along a zero-ghost axis that operator splits into the two parity sub-lattices, each the compact Laplacian at spacing 2h with zero at the ghost centres (the ghost-Dirichlet kind), along the other axes it is the p = 2 operator with their own kinds; the 2^z sub-problems (z zero-ghost axes) are FFT solves by sub-solvers of this solver's own type (sub). Spectrally equivalent (kappa <= 2.8 at p = 4, 7.1 at p = 8): 1, 17, 24, 30 iterations at p = 2, 4, 6, 8, independent of the grid (wide_dg_zero_ghost.py, wide_dg_e2b.py, with periodic and with reflected companion axes). The ghost layers returned carry the faces' rules (zero, reflection, periodic wrap).

Free-space ghosts (issue #84, E3c). FLAIL_BC_FREE_GHOSTS is PRISM's exact/open face: every ghost layer holds the continuum free-space potential of the source, g(x) = -sum_y factor f(y) dV / (4 pi |x - y|). The operator is the zero-ghost one (the two kinds mix on an axis, and share the CG path and its preconditioner); the ghost values move to the right-hand side, D(G)_0 u = factor f - D(G) u_g (u_g: g in the ghost slabs, zero inside), and the returned ghost layers carry g. Rank 0 sums over the gathered source (free_ghosts), so the values do not depend on the decomposition; the cost is O(ghost-slab cells x cells), a direct sum (a lattice free-space solve is #84 E7). The wide D(G) only. Edges and corners of the ghost region stay zero: D(G) reads along the axis lines only.

Scope (fatal otherwise). One level for every block, the blocks covering the domain. Without bc, every active direction periodic. Isolated: 3-D, every direction non-periodic, cubic cells (the kernel is the unit lattice's), the compact operator only. AMR is the multigrid's (adam_poisson_mg_object), which calls this solver on its coarsest level.

Source: src/lib/common/adam_poisson_fft_object.F90

Dependencies

Contents

• fft_slot
• poisson_fft_object
• destroy
• initialize
• initialize_blocks
• solve
• solve_packed
• solve_blocks
• solve_everywhere
• slab_setup
• solve_slabs
• slab_transform
• global_potential
• make_kernel
• destroy_plans
• make_plans
• extend_global
• set_plane
• fill_ghosts
• apply_wide
• precondition
• solve_pcg
• free_ghosts
• subtract_free_ghosts
• put_free_ghosts
• transform
• transform_kernel
• slab_make_plans
• slab_destroy_plans
• slab_planes
• slab_lines
• slab_transpose
• slab_spectral
• slab_keep_kernel
• slab_update
• root_cells_max
• smooth_size
• axis_eigenvalue
• dg_stencil

Variables

Name	Type	Attributes	Description
PI	real(kind=R8P)	parameter	Pi.
D1_COEF	real(kind=R8P)	parameter	a_j [j, p/2].

Derived Types

fft_slot

One sub-solver of the E2b preconditioner (issue #84). A scalar polymorphic pointer in a slot, not a polymorphic array component: nvfortran 26.9 rejects a recursive allocatable component (NVFORTRAN-S-0155) and crashed deallocating a polymorphic pointer array allocated with mold (pgf90_dealloc_mbr03a, on the FNL type). Owned: destroy frees it.

Components

Name	Type	Attributes	Description
s	class(poisson_fft_object)	pointer	The sub-solver.

poisson_fft_object

FFT Poisson solver on uniform periodic grids.

Inheritance

Components

Name	Type	Attributes	Description
is_initialized	logical		Initialization status.
n	integer(kind=I4P)		Global cells per direction.
m	integer(kind=I4P)		Transform cells per direction: n, or >= 2 (n + ngc) isolated.
isolated	logical		Isolated boundaries (zero padding), else periodic.
dirichlet	logical		Homogeneous Dirichlet boundaries (odd extension, #87).
bc	integer(kind=I4P)		Kinds per face [lo/hi, axis] (issue #84, E1a).
extended	logical		Some axis is extended (Dirichlet or Neumann faces).
dg_order	integer(kind=I4P)		0: the compact Laplacian; p: the wide D(G) of order p (E2a).
axis_order	integer(kind=I4P)		The operator per axis (dg_order, but in E2b's sub-problems).
pcg	logical		Some axis has zero-ghost faces and the operator is wide.
free	logical		Some face has free-space ghosts (E3c).
g_free	real(kind=R8P)	allocatable	Rank 0, E3c: the free-space ghost slabs, zero elsewhere.
null_xyz	logical		Null directions.
stencil	real(kind=R8P)		The 1-D D1 D1 stencil (times h^2), offsets -p..p.
tolerance	real(kind=R8P)		Relative max-norm residual to stop.
max_iterations	integer(kind=I4P)		CG iterations at most.
iterations	integer(kind=I4P)		Iterations of the last solve.
residual	real(kind=R8P)		Relative residual of the last solve.
redundant	logical		Every rank holds the global arrays and solves (#84 E5c).
wall	real(kind=R8P)		Wall time (s) of solve_blocks, accumulated over calls on
nsub	integer(kind=I4P)		Preconditioner sub-problems, 2^z.
sub_par	integer(kind=I4P)	allocatable	Parity per axis of each sub-problem [3,nsub] (-1: no split).
sub	type(fft_slot)	allocatable	The preconditioner's sub-problems [nsub] (see fft_slot).
nt	integer(kind=I4P)		Cells of the periodic transform the eigenvalues refer to: n,
nijk	integer(kind=I4P)		Cells per block per direction.
ngc	integer(kind=I4P)		Ghost cells per side.
level	integer(kind=I4P)		Level of every block.
dxyz	real(kind=R8P)		Cell spacing per direction.
lam1	real(kind=R8P)	allocatable	Eigenvalue part along x [0:n(1)/2].
lam2	real(kind=R8P)	allocatable	Eigenvalue part along y [0:n(2)-1].
lam3	real(kind=R8P)	allocatable	Eigenvalue part along z [0:n(3)-1].
f_global	real(kind=R8P)	allocatable	Rank 0: the source, then the potential [m(1),m(2),m(3)].
f_source	real(kind=R8P)	allocatable	Rank 0: a copy of the source, for the energy.
f_hat	complex(kind=R8P)	allocatable	Rank 0: the transform [m(1)/2+1,m(2),m(3)].
khat	real(kind=R8P)	allocatable	Rank 0, isolated: the kernel's spectrum over prod(m).
plan_r2c	type(c_ptr)		FFTW forward plan (rank 0).
plan_c2r	type(c_ptr)		FFTW backward plan (rank 0).
slab	logical		Every rank transforms its z-planes (see the module).
h1	integer(kind=I4P)		Complex cells along x after r2c, m(1)/2 + 1.
zr	integer(kind=I4P)	allocatable	First z-plane (0-based) of each rank's slab [0:ranks].
yr	integer(kind=I4P)	allocatable	First y-row (0-based) of each rank's pencils [0:ranks].
zown	integer(kind=I4P)	allocatable	Owner rank of each z-plane [0:m(3)-1].
tmap	integer(kind=I4P)	allocatable	Per axis, target cell -> source cell, -1 for zero [0:,3].
tsgn	real(kind=R8P)	allocatable	Per axis, the sign of the image [0:,3].
slab_r	real(kind=R8P)	allocatable	This rank's z-planes [m(1),m(2),planes].
slab_c	complex(kind=R8P)	allocatable	Their 2-D transforms [h1,m(2),planes].
pen	complex(kind=R8P)	allocatable	This rank's z-lines [m(3),h1,rows].
khat_pen	real(kind=R8P)	allocatable	Isolated: the kernel's spectrum in pencils [m(3),h1,rows].
w_r	real(kind=R8P)	allocatable	FFTW: the plane buffer [m(1),m(2)].
w_c	complex(kind=R8P)	allocatable	FFTW: the transformed plane buffer [h1,m(2)].
w_l	complex(kind=R8P)	allocatable	FFTW: the row-of-lines buffer [m(3),h1].
plan_pr2c	type(c_ptr)		FFTW: one plane, r2c.
plan_pc2r	type(c_ptr)		FFTW: one plane, c2r.
plan_lfwd	type(c_ptr)		FFTW: one row of z-lines, forward.
plan_lbwd	type(c_ptr)		FFTW: one row of z-lines, backward.

Type-Bound Procedures

Name	Attributes	Description
destroy	pass(self)	Destroy the solver.
initialize	pass(self)	Check the grid, size the solver, make the plans.
solve	pass(self)	Solve for phi from one variable of q.
solve_packed	pass(self)	Solve for phi from the packed interior source.
initialize_blocks	pass(self)	Size the solver from block counts and spacing (no grid, no field).
solve_blocks	pass(self)	Solve from an explicit list of blocks (coordinates and source).
make_kernel	pass(self)	Rank 0, isolated: the kernel's spectrum.
extend_global	pass(self)	Rank 0: the source extended over the transform cells.
apply_wide	pass(self)	Rank 0: y = D(G) x with the faces' ghost rules (E2b).
fill_ghosts	pass(self)	Rank 0: the ghost layers of a padded array by the faces.
precondition	pass(self)	Rank 0: z = M^-1 r, the sub-lattice solves (E2b).
solve_pcg	pass(self)	Rank 0: f_global <- factor D(G)^-1 f_global by PCG (E2b).
free_ghosts	pass(self)	Rank 0: g_free, the free-space ghost slabs (E3c).
subtract_free_ghosts	pass(self)	Rank 0: f <- f - D(G) u_g (E3c).
put_free_ghosts	pass(self)	Rank 0: the free-space ghost slabs into a padded array.
solve_everywhere	pass(self)	solve_blocks with redundant: every rank solves.
slab_setup	pass(self)	The slab layout, maps, arrays and plans (E5d).
solve_slabs	pass(self)	solve_blocks with slab: the distributed solve.
slab_transform	pass(self)	slab_r <- factor * L^-1 slab_r, distributed.
global_potential	pass(self)	Gather the slab solution to rank 0 (diagnostics).
destroy_plans	pass(self)	Destroy the transform plans.
make_plans	pass(self)	Make the transform plans (rank 0).
transform	pass(self)	f_global <- factor * L_h^-1 f_global (rank 0).
transform_kernel	pass(self)	khat <- the spectrum of f_global over prod(m) (rank 0, isolated).
root_cells_max	nopass	Most cells of a multigrid root solved by this backend (#84 E5c, E6b).
slab_make_plans	pass(self)	Make the slab plans: one plane, one row of lines (E5d).
slab_destroy_plans	pass(self)	Destroy the slab plans.
slab_planes	pass(self)	2-D transform of every local plane, slab_r slab_c.
slab_lines	pass(self)	1-D transform of every local z-line, in place in pen.
slab_transpose	pass(self)	slab_c pen (MPI_ALLTOALLV).
slab_spectral	pass(self)	pen <- the spectral solve (eigenvalues or kernel) of pen.
slab_keep_kernel	pass(self)	khat_pen <- pen over prod(m) (isolated, at setup).
slab_update	pass(self)	slab_r between the host and where the stages run.

Subroutines

destroy

Destroy the solver (and its sub-solvers).

Attributes: recursive

subroutine destroy(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

initialize

Check the grid (one level, the domain covered, the boundaries), size the solver, make the plans and, isolated, the kernel. bc (issue #84, E1a): the kind of every face; a periodic axis of the grid must be periodic in bc and the other way round. Without it every active axis must be periodic (or the solve isolated).

subroutine initialize(self, grid, field, isolated, bc, dg_order)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
grid	type(grid_object)	in		The grid.
field	type(field_object)	in		The field (block coordinates).
isolated	logical	in	optional	Isolated boundaries (default: periodic).
bc	integer(kind=I4P)	in	optional	Kinds per face [lo/hi, axis] (default periodic).
dg_order	integer(kind=I4P)	in	optional	0 compact (default); 2, 4, 6, 8: the wide D(G) (E2a).

Call graph

initialize_blocks

Size the solver for a uniform grid of nblk blocks per direction, nijk cells each, at level with spacing dxyz, and make the plans (and, isolated, the kernel) on rank 0. No check of a grid: the caller owns the scope (the multigrid's coarsest level calls it directly, issue #81 P5; with dirichlet, issue #87).

bc (issue #84, E1a) gives the boundary kind of every face (adam_flail_parameters), homogeneous: an axis is periodic (both faces), or its faces are Dirichlet/Neumann at the face (the source extended odd/even about each face: 2n cells when the two faces agree, 4n when they differ), or ghost-Dirichlet (both faces, odd about the ghost centres, 2(n+1) cells: what dirichlet sets on every active axis). Periodic against non-periodic, and ghost-Dirichlet against a face kind, on one axis are fatal; a null axis takes no extension.

redundant (issue #84, E5c): every rank holds the global arrays and the plans, and solve_blocks gathers the source to every rank (MPI_ALLGATHERV), transforms on every rank and keeps each rank's own blocks: no scatter, and no rank waits on another's transform. For a small grid (the multigrid's root); the compact operator only, not isolated, no energy.

Attributes: recursive

subroutine initialize_blocks(self, nblk, nijk, ngc, dxyz, level, null_xyz, isolated, dirichlet, bc, dg_order, axis_dg_order, redundant)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
nblk	integer(kind=I4P)	in		Blocks per direction.
nijk	integer(kind=I4P)	in		Cells per block per direction.
ngc	integer(kind=I4P)	in		Ghost cells per side (isolated: the padding margin).
dxyz	real(kind=R8P)	in		Cell spacing.
level	integer(kind=I4P)	in		Level of every block.
null_xyz	logical	in		Null directions.
isolated	logical	in		Isolated boundaries, else periodic.
dirichlet	logical	in	optional	Homogeneous ghost-Dirichlet everywhere (default false).
bc	integer(kind=I4P)	in	optional	Kinds per face [lo/hi, axis] (default periodic).
dg_order	integer(kind=I4P)	in	optional	0 compact (default); 2, 4, 6, 8: the wide D(G).
axis_dg_order	integer(kind=I4P)	in	optional	The operator per axis (E2b's sub-problems).
redundant	logical	in	optional	Every rank solves (default false).

Call graph

solve

Solve L_h phi = factor (q(iv) - mean(q(iv))) (periodic) or L_h phi = factor q(iv) with phi -> 0 at infinity (isolated); every ghost layer of phi is filled from the global solution.

subroutine solve(self, field, q, iv, factor, phi, energy)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
field	type(field_object)	in		The field (coordinates).
q	real(kind=R8P)	in		Conservative [nv,i,j,k,nb].
iv	integer(kind=I4P)	in		Variable of q, the source.
factor	real(kind=R8P)	in		Source factor (4 pi G).
phi	real(kind=R8P)	inout		Potential [i,j,k,nb].
energy	real(kind=R8P)	out	optional	sum_cells f phi (all ranks).

Call graph

solve_packed

Solve L_h phi = factor (f - mean(f)) (periodic) or L_h phi = factor f (isolated) from the packed interior source of this rank's blocks (block-major, i fastest, field%blocks_number * ni * nj * nk values); every ghost layer of phi is filled from the global solution (the periodic wrap, or the padded free-space solution).

subroutine solve_packed(self, field, buf, factor, phi, energy)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
field	type(field_object)	in		The field (block coordinates).
buf	real(kind=R8P)	in		Packed interior source of this rank.
factor	real(kind=R8P)	in		Source factor (4 pi G for gravity).
phi	real(kind=R8P)	inout		Potential [i,j,k,nb].
energy	real(kind=R8P)	out	optional	sum_cells f phi (all ranks).

Call graph

solve_blocks

The solve on an explicit list of this rank's blocks: their integer coordinates at the solver's level and their packed interior source (block-major, i fastest); out returns each block's potential with gw ghost layers (block-major, i fastest, from 1-gw to n+gw), filled from the global solution. Optionally the energy sum_cells f phi (rank 0 computes it, every rank gets it).

subroutine solve_blocks(self, nbl, coords, buf, factor, gw, out, energy)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
nbl	integer(kind=I4P)	in		Blocks of this rank.
coords	integer(kind=I4P)	in		Block coordinates [3,nbl].
buf	real(kind=R8P)	in		Packed interior source of this rank.
factor	real(kind=R8P)	in		Source factor (4 pi G for gravity).
gw	integer(kind=I4P)	in		Ghost layers of out.
out	real(kind=R8P)	inout		Packed potential, ghosts included.
energy	real(kind=R8P)	out	optional	sum_cells f phi (all ranks).

Call graph

solve_everywhere

solve_blocks with redundant (issue #84, E5c): the source of every rank gathered to every rank, placed and transformed on every rank (the same arrays, the same plans: the same values as rank 0 computes in solve_blocks), and each rank's blocks packed from its own global solution. wall: gather, solve, 0.

subroutine solve_everywhere(self, nbl, coords, buf, factor, gw, out)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
nbl	integer(kind=I4P)	in		Blocks of this rank.
coords	integer(kind=I4P)	in		Block coordinates [3,nbl].
buf	real(kind=R8P)	in		Packed interior source of this rank.
factor	real(kind=R8P)	in		Source factor.
gw	integer(kind=I4P)	in		Ghost layers of out.
out	real(kind=R8P)	inout		Packed potential, ghosts included.

Call graph

slab_setup

The distributed slab transform (issue #84, E5d): rank r owns the z-planes [zr(r), zr(r+1)) of the transform array and, after the transpose, the y-rows [yr(r), yr(r+1)) of its z-lines; the per-axis maps from a target cell of the transform array to its source cell and sign (the extensions of extend_global, the isolated padding); the arrays, the plans and, isolated, the kernel's spectrum in the pencil layout.

subroutine slab_setup(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

solve_slabs

solve_blocks with slab (issue #84, E5d): every block plane sent to the owners of the target planes that map to it, each rank's planes built (the extensions and the padding by the maps), transformed (slab_transform), and every block's potential with its ghost layers fetched back from the owners of its planes (wrapped modulo m, as solve_blocks reads f_global). The arithmetic of every value does not depend on the rank count (one plan per plane and per row of lines). wall: redistribution, transform, return.

subroutine solve_slabs(self, nbl, coords, buf, factor, gw, out, energy)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
nbl	integer(kind=I4P)	in		Blocks of this rank.
coords	integer(kind=I4P)	in		Block coordinates [3,nbl].
buf	real(kind=R8P)	in		Packed interior source of this rank.
factor	real(kind=R8P)	in		Source factor.
gw	integer(kind=I4P)	in		Ghost layers of out.
out	real(kind=R8P)	inout		Packed potential, ghosts included.
energy	real(kind=R8P)	out	optional	sum_cells f phi (all ranks).

Call graph

slab_transform

slab_r <- factor L^-1 slab_r on the distributed slabs: the planes' 2-D transforms, the transpose to z-lines, their 1-D transforms, the spectral solve, and back; the stages run where the backend keeps the slab arrays (slab_update moves slab_r there and back).

subroutine slab_transform(self, factor)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
factor	real(kind=R8P)	in		Source factor.

Call graph

global_potential

Rank 0: the transform array after the last solve, [m(1),m(2),m(3)] (the slabs gathered; f_global otherwise). For diagnostics (test_poisson_fft --dump); other ranks get an empty array.

subroutine global_potential(self, g)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	in		The solver.
g	real(kind=R8P)	inout	allocatable	The array.

make_kernel

Rank 0, isolated: place the kernel h^2 G(o) at every offset o, wrapped into (-m/2, m/2] per direction, and keep its spectrum over prod(m) (transform_kernel).

subroutine make_kernel(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

destroy_plans

Destroy the FFTW plans.

subroutine destroy_plans(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

make_plans

Make the FFTW plans on rank 0's global arrays (FFTW_ESTIMATE: the arrays are not touched, the plans deterministic).

subroutine make_plans(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

extend_global

Rank 0: the source extended over the m cells of every extended axis (cell index t = 0..m-1, array index t + 1), one axis after the other (the extensions are separable):

• ghost-Dirichlet (m = 2 (n + 1)): f(t) for t < n, 0 at t = n and t = m - 1 (the ghost centres), -f(2 n - t) for n < t < m - 1;
• face kinds: f(t) for t < n; the reflection about the high face, s_hi f(2 n - 1 - t), for n <= t < 2 n; with m = 4 n (the two faces differ) the reflection of [0, 2 n) about the low face, s_lo g(4 n - 1 - t), for 2 n <= t < 4 n. The sign s is -1 at a Dirichlet face, +1 at a Neumann one; the periodic extension then has both reflection symmetries, and so does its periodic solution: the boundary conditions on [0, n).

subroutine extend_global(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

set_plane

The plane t of axis d_ set to sg times the plane src (cell indices; sg = 0 zeroes it).

Attributes: pure

subroutine set_plane(f, d_, t, src, sg)

Arguments

Name	Type	Intent	Attributes	Description
f	real(kind=R8P)	inout		The array.
d_	integer(kind=I4P)	in		Axis, target and source cell indices.
t	integer(kind=I4P)	in		Axis, target and source cell indices.
src	integer(kind=I4P)	in		Axis, target and source cell indices.
sg	real(kind=R8P)	in		Sign (0: zero).

Call graph

fill_ghosts

Rank 0: the ng ghost layers of the padded global array u by the faces' rules, one axis after the other: periodic (and null) the wrap, Dirichlet and Neumann at the face the odd and even reflections about it, ghost-Dirichlet the odd reflection about the ghost centre (zero there), zero-ghost zero.

subroutine fill_ghosts(self, u, ng)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	in		The solver.
u	real(kind=R8P)	inout		Padded array [1-ng:n+ng per axis].
ng	integer(kind=I4P)	in		Ghost layers.

Call graph

apply_wide

Rank 0: y = D(G) x on the global interior, the ghosts by the faces' rules (fill_ghosts).

subroutine apply_wide(self, x, y)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	in		The solver.
x	real(kind=R8P)	in		Interior values [n].
y	real(kind=R8P)	out		Result [n].

Call graph

precondition

Rank 0: z = M^-1 r, M the wide p = 2 operator with the same faces: the 2^z parity sub-problems, each solved by its sub-solver's transform (the extension of its kinds, then the division by its eigenvalues).

subroutine precondition(self, r, z)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
r	real(kind=R8P)	in		Residual [n].
z	real(kind=R8P)	inout		Correction [n].

Call graph

solve_pcg

Rank 0: f_global <- u, D(G) u = factor f_global with zero-ghost faces (E2b), by conjugate gradients preconditioned by precondition, from zero, to tolerance (the max-norm residual relative to the source). D(G) and M are negative definite here (the zero ghosts remove every null mode), which the iteration takes as is.

subroutine solve_pcg(self, factor)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
factor	real(kind=R8P)	in		Source factor.

Call graph

free_ghosts

Rank 0, E3c: g_free, to depth gd, the ghost slab of every free-space face (the cells outside along that face's axis, inside along the others): the continuum free-space potential of the gathered source, g(x) = -sum_y factor f(y) dV / (4 pi |x - y|) over the interior cell centres, as PRISM's exact/open face.

subroutine free_ghosts(self, factor, gd)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
factor	real(kind=R8P)	in		Source factor.
gd	integer(kind=I4P)	in		Ghost layers.

Call graph

subtract_free_ghosts

Rank 0, E3c: f <- f - D(G) u_g on the interior, u_g the free-space ghost slabs (g_free) and zero inside.

subroutine subtract_free_ghosts(self, f)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	in		The solver.
f	real(kind=R8P)	inout		Right-hand side [n].

Call graph

put_free_ghosts

Rank 0, E3c: the ng ghost layers of every free-space face of the padded array u set to g_free.

subroutine put_free_ghosts(self, u, ng)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	in		The solver.
u	real(kind=R8P)	inout		Padded array [1-ng:n+ng per axis].
ng	integer(kind=I4P)	in		Ghost layers.

Call graph

transform

Rank 0: f_global <- factor * L_h^-1 f_global. Periodic: the modes of zero eigenvalue (the mean) set to zero; isolated: the convolution with the kernel (khat).

subroutine transform(self, factor)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
factor	real(kind=R8P)	in		Source factor.

Call graph

transform_kernel

Rank 0, isolated: khat <- FFT(f_global) / prod(m), real (the kernel is even in every direction).

subroutine transform_kernel(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

slab_make_plans

The slab plans (#84 E5d): one 2-D plane (r2c, c2r) and one row of z-lines (h1 lines of m(3) cells, c2c), each made once on its own buffer and executed on every plane and row in turn (FFTW_ESTIMATE: deterministic, and the arithmetic of a value does not depend on how many planes or rows a rank holds).

subroutine slab_make_plans(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

slab_destroy_plans

Destroy the slab plans.

subroutine slab_destroy_plans(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

Call graph

slab_planes

The 2-D transform of every local plane: forward slab_r -> slab_c (r2c), backward slab_c -> slab_r (c2r), one plane at a time through the plane buffers.

subroutine slab_planes(self, forward)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
forward	logical	in		Direction.

Call graph

slab_lines

The 1-D transform of every local z-line, in place in pen, one row of h1 lines at a time through the row buffer.

subroutine slab_lines(self, forward)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
forward	logical	in		Direction.

Call graph

slab_transpose

own block is copied directly.

subroutine slab_transpose(self, forward)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
forward	logical	in		Direction.

Call graph

slab_spectral

The spectral solve on this rank's z-lines: the division by the eigenvalues (the zero modes dropped) or, isolated, the product with the kernel's spectrum.

subroutine slab_spectral(self, factor)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
factor	real(kind=R8P)	in		Source factor.

Call graph

slab_keep_kernel

Isolated, at setup: khat_pen <- real(pen) / prod(m), the kernel's spectrum (real: the kernel is even).

subroutine slab_keep_kernel(self)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.

slab_update

Move slab_r to (to_device) or from where the stages run: nothing on the host.

subroutine slab_update(self, to_device)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_fft_object)	inout		The solver.
to_device	logical	in		Direction.

Call graph

Functions

root_cells_max

The most cells of a multigrid root on this backend (#84 E5c, E6b), for a root whose transform (the FFT solver's extension of its bounded axes, 2 or 4 times per axis) has growth cells per cell: on the host a transform of at most 128^3 cells, where the distributed transform (E5d) of a large root beats smoothing more levels (medooza CPU, 3.9 M AMR cells, np 1 / 4 / 16): periodic 12.4 / 3.9 / 2.5 s with a 32^3 root, 8.6 / 2.7 / 1.6 s with 128^3; isolated (the root ghost-Dirichlet, its transform 8 times its cells) 30.4 / 9.6 / 6.3 s with 32^3, 26.8 / 8.5 / 5.2 s with 64^3, 49.2 / 15.0 / 6.2 s with 128^3. cuFFT overrides it.

Attributes: pure

Returns: integer(kind=I4P)

function root_cells_max(growth) result(n)

Arguments

Name	Type	Intent	Attributes	Description
growth	integer(kind=I4P)	in		Transform cells per root cell.

Call graph

smooth_size

The smallest even size m >= k with no prime factor above 13: FFTW has hard-coded transforms for these factors (#84 E5b, medooza: a 7-smooth rule moved the James step's 260 = 2^2 5 13 to 270 and cost it 15%; cuFFT took the same time on 262 and 270).

Attributes: pure

Returns: integer(kind=I4P)

function smooth_size(k) result(m)

Arguments

Name	Type	Intent	Attributes	Description
k	integer(kind=I4P)	in		Least size.

Call graph

axis_eigenvalue

The eigenvalue part of mode k of an axis on the nt-cell periodic transform: the compact operator, (2 cos(theta) - 2) / dx^2, or the wide D(G) of order p, -s_p(theta)^2 / dx^2, exactly 0 at theta = 0 and pi.

Attributes: pure

Returns: real(kind=R8P)

function axis_eigenvalue(k, nt, dx, dg_order) result(lam)

Arguments

Name	Type	Intent	Attributes	Description
k	integer(kind=I4P)	in		Mode, transform cells.
nt	integer(kind=I4P)	in		Mode, transform cells.
dx	real(kind=R8P)	in		Spacing.
dg_order	integer(kind=I4P)	in		0 compact, else the order of D(G).

Call graph

dg_stencil

The 1-D stencil of D1 D1 (times h^2), offsets -p..p, the centred D1 of order p (D1_COEF) composed with itself.

Attributes: pure

Returns: real(kind=R8P)

function dg_stencil(p) result(c)

Arguments

Name	Type	Intent	Attributes	Description
p	integer(kind=I4P)	in		Order.

Call graph