adam_fnl_poisson_fft_object

ADAM, FNL FFT Poisson solver on uniform periodic grids (cuFFT).

Extends poisson_fft_object (issue #81, P2): the gather to rank 0, the placement and the scatter of phi with its ghosts are the CPU ones; the spectral solve runs on rank 0's GPU through cuFFT (-D_CUFFT, -cudalib=cufft). The global source goes to the device, a D2Z transform, the division by the eigenvalues of L_h (the same host-computed lam1..3 as the CPU solver), a Z2D transform, and the potential comes back. cuFFT runs on the stream of the synchronous OpenACC queue, so the transforms and the kernel between them are ordered.

Isolated boundaries (issue #81, M8-P4): the transforms are the padded ones (m), the kernel's spectrum khat is computed once on the device (transform_kernel) and stays there (enter data, released by destroy_plans); the solve multiplies by it instead of dividing by the eigenvalues.

The distributed slab transform (issue #84, E5d) keeps the CPU redistribution and return of the block planes (adam_poisson_fft_object, on the host); the transform runs on each rank's GPU, its arrays resident there (slab_r, slab_c, pen, khat_pen, the transpose buffers: enter data with the plans): slab_r goes up, the local planes (one D2Z or Z2D plan for one plane, executed on every plane in turn), the transpose to z-lines (packed and unpacked on the device, MPI_ALLTOALLV on the device buffers: GPU-aware MPI, as the ghost exchange), the local z-lines (one Z2Z plan for one row of h1 lines, on every row in turn), the spectral solve, back, and slab_r comes down. The arithmetic of a value does not depend on the rank count.

cuFFT is NVIDIA only: built without _CUFFT (an AMD build), making the plans stops the run.

Source: src/lib/fnl/adam_fnl_poisson_fft_object.F90

Dependencies

Contents

• poisson_fft_fnl_object
• destroy_plans_fnl
• make_plans_fnl
• transform_fnl
• transform_kernel_fnl
• slab_make_plans_fnl
• slab_destroy_plans_fnl
• slab_planes_fnl
• slab_lines_fnl
• slab_transpose_fnl
• slab_spectral_fnl
• slab_keep_kernel_fnl
• slab_update_fnl
• root_cells_max_fnl

Derived Types

poisson_fft_fnl_object

FNL FFT Poisson solver on uniform periodic grids (cuFFT).

Inheritance

Extends: poisson_fft_object

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.
plan_d2z	integer(kind=c_int)		cuFFT forward plan (rank 0).
plan_z2d	integer(kind=c_int)		cuFFT backward plan (rank 0).
has_plans	logical		The cuFFT plans exist.
khat_on_device	logical		khat is on the device (isolated).
plan_sd2z	integer(kind=c_int)		cuFFT slab plan: one plane, D2Z (#84 E5d).
plan_sz2d	integer(kind=c_int)		cuFFT slab plan: one plane, Z2D.
plan_sz2z	integer(kind=c_int)		cuFFT slab plan: one row of z-lines, Z2Z.
has_slab_plans	logical		The cuFFT slab plans exist.
slab_on_device	logical		The slab arrays are on the device.
yown	integer(kind=I4P)	allocatable	Owner rank of each y-row [0:m(2)-1].
t_sb	complex(kind=R8P)	allocatable	Device: the transpose's send buffer.
t_rb	complex(kind=R8P)	allocatable	Device: the transpose's receive buffer.

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).
global_potential	pass(self)	Gather the slab solution to rank 0 (diagnostics).
destroy_plans	pass(self)	Destroy the cuFFT plans.
make_plans	pass(self)	Make the cuFFT plans (rank 0).
transform	pass(self)	The spectral solve on the GPU (rank 0).
transform_kernel	pass(self)	The kernel's spectrum on the GPU (isolated).
root_cells_max	nopass	32^3: the device smooths, the host transforms.
slab_make_plans	pass(self)	The cuFFT slab plans (#84 E5d).
slab_destroy_plans	pass(self)	Destroy them.
slab_planes	pass(self)	Every local plane on the GPU.
slab_lines	pass(self)	Every local z-line on the GPU.
slab_transpose	pass(self)	The transpose on the GPU.
slab_spectral	pass(self)	The spectral solve on the GPU.
slab_keep_kernel	pass(self)	khat_pen on the GPU.
slab_update	pass(self)	slab_r to and from the GPU.

Subroutines

destroy_plans_fnl

Destroy the cuFFT plans.

subroutine destroy_plans_fnl(self)

Arguments

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

make_plans_fnl

Make the cuFFT plans (D2Z and Z2D, cuFFT is row-major: the slowest dimension first) on the synchronous queue.

subroutine make_plans_fnl(self)

Arguments

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

Call graph

transform_fnl

Rank 0: f_global <- factor * L_h^-1 f_global on the GPU, the modes of zero eigenvalue (the mean) set to zero.

subroutine transform_fnl(self, factor)

Arguments

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

Call graph

transform_kernel_fnl

Rank 0, isolated: khat <- FFT(f_global) / prod(m) on the GPU (real: the kernel is even); khat stays on the device for every solve.

subroutine transform_kernel_fnl(self)

Arguments

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

Call graph

slab_make_plans_fnl

The cuFFT slab plans (#84 E5d): one plane D2Z and Z2D, one row of z-lines Z2Z (h1 lines of m(3) cells, contiguous), on the synchronous OpenACC queue's stream; the slab arrays and the transpose buffers on the device.

subroutine slab_make_plans_fnl(self)

Arguments

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

Call graph

slab_destroy_plans_fnl

Destroy the cuFFT slab plans.

subroutine slab_destroy_plans_fnl(self)

Arguments

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

slab_planes_fnl

Every local plane on the GPU: forward slab_r -> slab_c (D2Z), backward slab_c -> slab_r (Z2D), on the device.

subroutine slab_planes_fnl(self, forward)

Arguments

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

Call graph

slab_lines_fnl

Every local z-line on the GPU, in place in pen (Z2Z, forward or inverse), on the device.

subroutine slab_lines_fnl(self, forward)

Arguments

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

Call graph

slab_transpose_fnl

The transpose on the device: forward slab_c(x, y, z local) -> pen(z, x, y local), backward the reverse; the blocks packed by destination rank (the CPU layout), MPI_ALLTOALLV on the device buffers (skipped on one rank).

subroutine slab_transpose_fnl(self, forward)

Arguments

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

Call graph

slab_spectral_fnl

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

subroutine slab_spectral_fnl(self, factor)

Arguments

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

Call graph

slab_keep_kernel_fnl

Isolated, at setup: khat_pen <- real(pen) / prod(m) on the device, where it stays.

subroutine slab_keep_kernel_fnl(self)

Arguments

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

Call graph

slab_update_fnl

slab_r to (to_device) or from the device.

subroutine slab_update_fnl(self, to_device)

Arguments

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

Functions

root_cells_max_fnl

The most cells of a multigrid root on cuFFT (#84 E6, E6b): 32^3 whatever the transform's growth. The multigrid smooths on the device and the root's transform goes through the host's redistribution, so a small root wins (medooza FNL, 3.9 M AMR cells, 1 / 2 / 4 GPUs: periodic 0.33 / 0.19 / 0.12 s with a 32^3 root, 0.80 / 0.45 / 0.37 s with 128^3); for the isolated root 16^3 gains only 3-10% over 32^3 (0.83 / 0.51 / 0.39 against 0.85 / 0.52 / 0.43 s), not worth a uniform 32^3 grid losing its exact one-V-cycle solve.

Attributes: pure

Returns: integer(kind=I4P)

function root_cells_max_fnl(growth) result(n)

Arguments

Name	Type	Intent	Attributes	Description
growth	integer(kind=I4P)	in		Transform cells per root cell (unused: the cap is in cells).