adam_poisson_fft_object

ADAM, FFT Poisson solver on uniform periodic 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 (issue #81, D-81-2, first version). Every rank packs the interior f of its blocks with their integer block coordinates; MPI_Gatherv brings them to rank 0, which places them by the global cell index, transforms (transform, real-to-complex and back), and MPI_Scattervs every block's phi with all its ghost layers already filled by the periodic wrap of the global solution, so no ghost exchange follows. Correct for any rank count, ranks without blocks included; not scalable (rank 0 holds the whole grid).

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 make_plans, transform and destroy_plans with cuFFT on the GPU of rank 0. 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 (computed on rank 0 from the global arrays, broadcast): 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. 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.

Scope (fatal otherwise). One level for every block, the blocks covering the domain. Periodic: every active direction periodic. Isolated: 3-D, every direction non-periodic, cubic cells (the kernel is the unit lattice's). 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

• poisson_fft_object
• destroy
• initialize
• initialize_blocks
• solve
• solve_packed
• solve_blocks
• make_kernel
• destroy_plans
• make_plans
• transform
• transform_kernel

Variables

Name	Type	Attributes	Description
PI	real(kind=R8P)	parameter	Pi.

Derived Types

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).
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).

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.
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).

Subroutines

destroy

Destroy the solver.

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.

subroutine initialize(self, grid, field, isolated)

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).

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).

subroutine initialize_blocks(self, nblk, nijk, ngc, dxyz, level, null_xyz, isolated, dirichlet)

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 Dirichlet boundaries (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

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

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