adam_poisson_mg_object

ADAM, multigrid-preconditioned conjugate-gradient Poisson solver on AMR grids (issue #81, M8-P5).

Solves L_h phi = factor (f - mean(f)) on the leaf cells of a 2:1-balanced octree (quadtree in 2-D, null_z), periodic in every active direction, with L_h a composite second-order Laplacian that is symmetric and consistent (issue #81, P5a):

• between leaves of one level, the 7-point (5-point) stencil;
• across a 2:1 face, the mimetic mean flux: every fine cell f_k (k = 1..n, n = 2^(d-1)) on a coarse cell C's face gets the same flux (2/3) w (mean_k phi_fk - phi_C), w the fine face weight (area over the fine spacing), and C gets their sum. It is the hybrid mimetic finite-difference flux of the coarse cell seen as a polygon with hanging nodes, its face unknowns eliminated: consistency (every sub-face flux exact for linear phi) forces the coarse cell's sub-face block to rank one (it sees only the mean of the face unknowns), and the elimination gives the mean of the fine cells. Exact for linear phi face by face, so the truncation error stays bounded at the seams (phi second order, the seam face gradient first order); symmetric in the V-weighted form (rows times the cell volume), so the solution map phi = K f has V K symmetric, which exact energy conservation of self-gravity needs. Rejected on the way (#81, P5-0 and P5a): Tomida & Stone's flux matching (consistent, not symmetric: the fallback); its symmetrisation (B + B^T) / 2 (symmetric, inconsistent: O(1/h) truncation, a seam gravity that does not converge); the Galerkin operator P^T A_fine P (symmetric, first order).

The solve is flexible preconditioned conjugate gradients (Notay's FCG(1)) on the composite operator: correctness lives in its matrix-vector product only. The preconditioner is one V-cycle over the level hierarchy: level l holds the leaves at l and a shadow block for every tree node at l that has children (ADAM keeps leaves only); red-black Gauss-Seidel sweeps down and up on the plain level Laplacian, a face to a cell outside the level reading the coarser correction (piecewise constant) on the way up; restriction the sum of the children's residuals, prolongation piecewise constant; the coarsest tree level (every block at the minimum leaf level: the whole domain) coarsened further inside the blocks (issue #84, E5c): each level below it halves the cells per block along every active direction, the same blocks on the same ranks, the plain level Laplacian with the face kinds, the same smoother and transfers, while the cells per block are even and the domain has more than root_cells cells at the level (default: the root solver's root_cells_max, told the growth of the root's transform: on the host 128^3 transform cells, so a periodic root of 128^3 cells and a ghost-Dirichlet one of 64^3; on cuFFT 32^3 cells, #84 E6 and E6b); the last of them, the root (lbot), is solved exactly by the FFT solver (adam_poisson_fft_object%solve_blocks; of the caller's backend type with initialize(coarse=)): on every rank (redundant), from the source gathered by MPI_ALLGATHERV, when it has at most REDUNDANT_ROOT_MAX (32^3) cells, so no rank-0 gather and scatter per V-cycle; a larger root (root_cells above the default) on the FFT solver's distributed path, since every rank transforming it contends for the node's memory (FNL, 4 GPUs, the isolated 128^3 coarsest level: 4.7 -> 6.4 s per solve). A coarsest tree level of at most root_cells cells is its own root, solved as before. Measured (test_poisson_mg --sub, every level down to one cell per block): 1 -> 14 FCG iterations on a uniform 1:2:3-cell box (point smoothing on anisotropic cells), +1 to +13 on AMR (24-47 -> 25-53). Medooza CPU, 3.9 M AMR cells per solve: isolated 50.7 -> 31.2 s at np 1 and 39.9 -> 6.5 s at np 16 (its 258^3 coarsest level and 42 exact solves were 70% of it); periodic 9.2 -> 12.7 s at np 1, even at np 8, 2.6 -> 2.3 s at np 16 (the generic CSR smoother on the base level costs more than one FFT of it on few ranks; a stencil smoother is #84 E6). It stops at tolerance or at the round-off floor of the source, 100 eps max|factor f| / max|factor (f - mean f)|, whichever is larger (a small perturbation of a uniform density is known only to that floor, #81 P5b).

Backend hooks (issue #84, E6a). The FCG iteration and the V-cycle work on vectors allocated once (make_work) through overridable hooks: apply, smooth, the vector operations (vec_*, nopass), vec_update (a vector to or from where the hooks compute) and work_place (the operators and vectors put there after initialize). Here they are the host loops, bitwise those of before; adam_fnl_poisson_mg_object runs them on the device. The right-hand side, the root, the James step and the ghost fill's placement stay on the host.

Distributed and assembled. Every rank gathers the leaf list once (block coordinates, owner, local index), derives the shadow blocks and their owners (the owner of the first leaf below, in the gathered order: the same on every rank), and assembles, from the geometry only, the rows of its own cells: the operator, the ghost fill, the level operators and transfers. Columns are global cell keys; every operator owns a halo-exchange plan (the remote columns it reads, fetched from their owners by MPI_ALLTOALLV). Correct for any rank count. Cell keys carry the level offset by KEY_LEVEL_OFFSET (the levels inside the blocks may be negative).

Domain boundaries (issue #87, P1; per face, issue #84, E1b). Every active axis is periodic (both faces) or bounded, and each face of a bounded axis has a kind (adam_flail_parameters; initialize(..., bc=), by default ghost-Dirichlet on every bounded axis, the #87 behaviour). A boundary face of a leaf cell C (weight w, the face area over the cell spacing h) adds to C's row, and maps its boundary value g onto the right-hand side:

• FLAIL_BC_DIRICHLET_GHOST, g the value at the centre of the first ghost cell: -w phi_C, w g;
• FLAIL_BC_DIRICHLET, g the value at the face (phi_ghost = 2 g - phi_C): -2 w phi_C, 2 w g;
• FLAIL_BC_NEUMANN, g the outward normal derivative at the face (phi_ghost = phi_C + h g): nothing, w h g. So the operator A keeps the homogeneous rows, symmetric, and solve(..., ghost=g) solves A phi = factor V f - B g, B the map of the boundary values onto the boundary rows (apply_b; its transpose apply_bt, the boundary rows onto the faces). The faces are listed per rank (nbnd, bnd_row, bnd_w the weight in B, bnd_kind, bnd_dir the outward axis signed, bnd_x where g lives: the ghost centre or the face centre). Without a Dirichlet face the operator is singular, as on a periodic grid (singular): the solve removes the mean of the right-hand side factor V f - B g (the compatibility condition, exact for compatible data) and of the solution. The V-cycle's levels take the same face kinds, homogeneous, and its coarsest level the FFT solver with them (odd or even extensions; issue #84, E1a); on a uniform grid the two operators coincide and one V-cycle is the exact solve. The two faces of an axis are both periodic or neither, and both ghost-Dirichlet or neither (the FFT solver's scope).

Isolated boundaries (issue #87, P2), a symmetric James step. initialize(..., isolated=.true.): each solve is the zero-ghost solve phi0 = A^-1 factor V f; the screening charge s = B^T phi0 on the ghost cells (extended by zero, phi0 satisfies the infinite operator with the extra source -s there); the ghost values of the free-space solution g = -G s, G the inverse of the V-weighted operator on the infinite coarsest lattice from ghost cell to ghost cell (the lattice Green's function over the coarsest spacing; ghost_green, a zero-padded FFT of the coarsest box widened by the ghost layer, gfft); and the solve with them. The map (A^-1 + A^-1 B G B^T A^-1) factor V is symmetric in the V-weighted form (A and G are), so self-gravity keeps its exact E + W exchange; on a uniform grid it is the discrete free-space solution (the FFT solver's isolated mode), on AMR second order (G misses the refined interior, a truncation-order perturbation). G is exact on the boundary only if every boundary leaf is at the coarsest level (each face then has its own ghost cell of the coarse lattice): fatal otherwise (D-87-2). The textbook James method (the screening charge on the boundary surface, the potential at the ghost centres) is not symmetric (selfgravity/james_proto.py, the P0 prototype).

Scope (fatal otherwise). The face kinds above: without bc every active direction periodic or none (a mixed grid would get a silent zero on its bounded faces), with it the grid's periodicity and bc agree; isolated: 3-D, no periodic direction, no bc, cubic cells, every boundary leaf at the coarsest level; 2:1 balance across faces; fewer than 2^19 cells per direction at the finest level (the cell keys). On output phi holds the leaf solution and its first ghost layer per face: same level, the neighbour; at a domain face, the ghost of its kind, g, 2 g - phi_C or phi_C + h g (g = 0 when none is given); across a 2:1 face, the value that makes the one-sided face difference equal the operator's own face gradient, (mean_f - phi_C) / (1.5 h) (fine side: phi_f - (mean_f - phi_C) / 1.5; coarse side: phi_C + (mean_f - phi_C) / 0.75); deeper ghost layers copy the first.

Source: src/lib/common/adam_poisson_mg_object.F90

Dependencies

Contents

• hash_table
• dcsr
• rvec
• ivec
• triplets
• poisson_mg_object
• destroy
• initialize
• solve
• solve_packed
• solve_fcg
• apply_a
• apply_b
• apply_bt
• apply
• vec_fill
• vec_copy
• vec_lincomb
• vec_axpy
• vec_permute
• vec_update
• work_place
• assemble
• a_row
• face_cells
• ghost_green
• comp_value
• make_work
• locate
• smooth
• vcycle
• kdecode
• hash_init
• hash_free
• hash_put
• trip_add
• symmetry_defect
• vec_dot
• vec_relmax
• cell_volume
• classify
• dx_of
• face_weight
• lin
• lin_of
• ncb_of
• ncell
• nijk_of
• parent3
• outside
• wrap
• gdot
• bkey
• ckey
• hash_get

Variables

Name	Type	Attributes	Description
MG_SEAM_GHOST_SCALE	real(kind=R8P)	parameter	composite operator's face gradient (clients that invert it, as
CLS_SAME	integer(kind=I4P)	parameter	Neighbour position: a leaf of the same level.
CLS_FINER	integer(kind=I4P)	parameter	Neighbour position: refined (finer leaves).
CLS_COARSER	integer(kind=I4P)	parameter	Neighbour position: inside a leaf one level coarser.
SPACE_COMP	integer(kind=I4P)	parameter	Vector space: the composite leaf cells (levels are the other spaces;
KEY_LEVEL_OFFSET	integer(kind=I4P)	parameter	Level offset in the cell keys (the levels below lmin, E5c).
REDUNDANT_ROOT_MAX	integer(kind=I4P)	parameter	A root of at most this many cells is solved on every rank;

Derived Types

hash_table

Open-addressing hash, int64 key (>= 0) -> int32 value (0: absent).

Components

Name	Type	Attributes	Description
keys	integer(kind=I8P)	allocatable	Keys (-1: empty).
vals	integer(kind=I4P)	allocatable	Values.
cap	integer(kind=I4P)		Capacity (a prime).

dcsr

Distributed CSR operator: local rows; columns local (1..nloc) or halo (nloc+1..nloc+nhalo).

Components

Name	Type	Attributes	Description
nrow	integer(kind=I4P)		Rows, local columns, halo columns.
nloc	integer(kind=I4P)		Rows, local columns, halo columns.
nhalo	integer(kind=I4P)		Rows, local columns, halo columns.
rowptr	integer(kind=I4P)	allocatable	Row pointers [nrow+1].
col	integer(kind=I4P)	allocatable	Column slots.
val	real(kind=R8P)	allocatable	Values.
scount	integer(kind=I4P)	allocatable	Send counts and displacements per rank.
sdispl	integer(kind=I4P)	allocatable	Send counts and displacements per rank.
sidx	integer(kind=I4P)	allocatable	Local column indices to send.
rcount	integer(kind=I4P)	allocatable	Receive counts and displacements per rank (halo order).
rdispl	integer(kind=I4P)	allocatable	Receive counts and displacements per rank (halo order).

rvec

A real vector (one per level).

Components

Name	Type	Attributes	Description
v	real(kind=R8P)	allocatable	Values.

ivec

An integer vector (one per level).

Components

Name	Type	Attributes	Description
v	integer(kind=I4P)	allocatable	Values.

triplets

Assembly buffer: (row, column key, value).

Components

Name	Type	Attributes	Description
n	integer(kind=I4P)		Entries.
r	integer(kind=I4P)	allocatable	Rows.
k	integer(kind=I8P)	allocatable	Column keys.
v	real(kind=R8P)	allocatable	Values.

poisson_mg_object

Multigrid-preconditioned CG Poisson solver on AMR grids.

Inheritance

Components

Name	Type	Attributes	Description
is_initialized	logical		Initialization status.
nijk	integer(kind=I4P)		Cells per block per direction.
ncb	integer(kind=I4P)		Cells per block.
ngc	integer(kind=I4P)		Ghost cells per side.
nd	integer(kind=I4P)		Active directions.
act	integer(kind=I4P)		The active directions.
null	logical		Null directions.
nblk0	integer(kind=I4P)		Blocks per direction at level 0.
dx0	real(kind=R8P)		Cell spacing at level 0.
lmin	integer(kind=I4P)		Leaf levels range.
lmax	integer(kind=I4P)		Leaf levels range.
lbot	integer(kind=I4P)		The root level: lmin less the levels inside the blocks (E5c).
root_max	integer(kind=I4P)		Most cells of the root (E5c; the root solver's default).
nr	integer(kind=I4P)		Ranks, this rank.
me	integer(kind=I4P)		Ranks, this rank.
nleaf	integer(kind=I4P)		Leaves over the ranks.
leaf_c	integer(kind=I4P)	allocatable	Leaf block coordinates (i,j,k,level) [4,nleaf].
leaf_owner	integer(kind=I4P)	allocatable	Owner rank.
leaf_loc	integer(kind=I4P)	allocatable	Local block index (the field's) on its owner.
leaf_hash	type(hash_table)		Block key -> leaf id.
nnode	integer(kind=I4P)		Nodes.
node_c	integer(kind=I4P)	allocatable	Node block coordinates (i,j,k,level) [4,nnode].
node_owner	integer(kind=I4P)	allocatable	Owner rank.
node_loc	integer(kind=I4P)	allocatable	Local block index in its level space on its owner.
node_leaf	logical	allocatable	The node is a leaf.
node_hash	type(hash_table)		Block key -> node id.
nloc_lev	integer(kind=I4P)	allocatable	This rank's blocks per level [lbot:lmax].
lev_node	integer(kind=I4P)	allocatable	This rank's level blocks: node id [max blocks, lbot:lmax]
nleaf_loc	integer(kind=I4P)		This rank's leaves.
leaf_lev	integer(kind=I4P)	allocatable	This rank's leaf b: its level.
leaf_lev_loc	integer(kind=I4P)	allocatable	This rank's leaf b: its local index in its level space.
a	type(dcsr)		The composite operator on the leaves.
g	type(dcsr)		Ghost fill: the first ghost layer of the leaves.
g_dst	integer(kind=I4P)	allocatable	Ghost fill destinations (i,j,k,b) [4,rows].
al	type(dcsr)	allocatable	Level operators [lbot+1:lmax].
el	type(dcsr)	allocatable	Level ghost terms from level l-1 [lbot+1:lmax].
rl	type(dcsr)	allocatable	Restriction l -> l-1 [lbot+1:lmax].
pl	type(dcsr)	allocatable	Prolongation l-1 -> l [lbot+1:lmax].
diag	type(rvec)	allocatable	Level operator diagonals [lbot+1:lmax].
colour	type(rvec)	allocatable	Red-black colour (0, 1) of the level cells [lbot+1:lmax].
vol	real(kind=R8P)	allocatable	Leaf cell volumes.
bc	integer(kind=I4P)		Kind per face [lo/hi, axis] (adam_flail_parameters).
periodic	logical		Periodic axes (a null axis counts as periodic).
singular	logical		No Dirichlet face: the operator's kernel is the constants.
emin	real(kind=R8P)		Domain lower corner.
nbnd	integer(kind=I4P)		This rank's boundary faces.
bnd_row	integer(kind=I4P)	allocatable	Boundary face: the local leaf cell (row) it bounds.
bnd_w	real(kind=R8P)	allocatable	Boundary face: its weight in B (see the module).
bnd_x	real(kind=R8P)	allocatable	Boundary face: where its value lives, ghost or face centre.
bnd_kind	integer(kind=I4P)	allocatable	Boundary face: its kind.
bnd_dir	integer(kind=I4P)	allocatable	Boundary face: the outward axis, signed (-3..3).
bnd_gs	real(kind=R8P)	allocatable	Boundary face: its value's weight in the ghost fill.
g_bnd	integer(kind=I4P)	allocatable	Ghost-fill row: its boundary face, 0 for an interior one.
isolated	logical		Free-space boundaries: the ghost values by the James step.
nc	integer(kind=I4P)		Cells per direction at the coarsest level.
bnd_gi	integer(kind=I4P)	allocatable	Boundary face: its ghost cell in the widened coarsest box.
u0	real(kind=R8P)	allocatable	The last zero-ghost solution (warm start of the first solve).
iterations0	integer(kind=I4P)		Iterations of the last zero-ghost solve.
fft	class(poisson_fft_object)	allocatable	The root solver, on every rank (FFTW, or the caller's backend).
gfft	class(poisson_fft_object)	allocatable
tolerance	real(kind=R8P)		Relative max-norm residual to stop.
max_iterations	integer(kind=I4P)		FCG iterations at most.
sweeps	integer(kind=I4P)		Red-black sweeps per level, down and up.
iterations	integer(kind=I4P)		Iterations of the last solve.
residual	real(kind=R8P)		Relative residual of the last solve.
tolerance_used	real(kind=R8P)		Stop tolerance of the last solve (the round-off floor of the
factor_mean	real(kind=R8P)		Mean convergence factor per iteration of the last solve.
u	real(kind=R8P)	allocatable	The last solution (warm start).
has_guess	logical		u holds a solution of the current grid.
wf	real(kind=R8P)	allocatable	previous p and A p (leaf space).
wr	real(kind=R8P)	allocatable	previous p and A p (leaf space).
wz	real(kind=R8P)	allocatable	previous p and A p (leaf space).
wp	real(kind=R8P)	allocatable	previous p and A p (leaf space).
wq	real(kind=R8P)	allocatable	previous p and A p (leaf space).
wpp	real(kind=R8P)	allocatable	previous p and A p (leaf space).
wqp	real(kind=R8P)	allocatable	previous p and A p (leaf space).
w1	real(kind=R8P)	allocatable	Ones (leaf space).
wg	real(kind=R8P)	allocatable	Ghost-fill values [g%nrow].
lr	type(rvec)	allocatable	V-cycle: level residual and correction [lbot:lmax].
le	type(rvec)	allocatable	V-cycle: level residual and correction [lbot:lmax].
lt1	type(rvec)	allocatable	V-cycle: level temporaries [lbot:lmax].
lt2	type(rvec)	allocatable	V-cycle: level temporaries [lbot:lmax].
mlev	type(ivec)	allocatable	Per level: the level-space and leaf-space cells of the
mleaf	type(ivec)	allocatable	Per level: the level-space and leaf-space cells of the
root_coords	integer(kind=I4P)	allocatable	V-cycle: the root blocks [3,blocks].
root_buf	real(kind=R8P)	allocatable	V-cycle: the root's source and solution (host).
root_out	real(kind=R8P)	allocatable	V-cycle: the root's source and solution (host).

Type-Bound Procedures

Name	Attributes	Description
destroy	pass(self)	Destroy the solver.
initialize	pass(self)	Gather the leaves, assemble the operators.
solve	pass(self)	Solve for phi from one variable of q.
solve_packed	pass(self)	Solve for phi from the packed interior source.
symmetry_defect	pass(self)	|x . A y - y . A x| / (|x| |A y|) (verification).
apply_a	pass(self)	y = A x on the leaves (verification).
apply_b	pass(self)	y = B g: the boundary values onto the boundary rows.
apply_bt	pass(self)	s = B^T x: the boundary rows onto the boundary faces.
apply	pass(self)	y = M x for a distributed operator.
smooth	pass(self)	Red-black Gauss-Seidel on a level.
vec_fill	nopass	x(1:n) = a.
vec_copy	nopass	y(1:n) = x(1:n).
vec_lincomb	nopass	w(1:n) = x(1:n) + a y(1:n).
vec_axpy	nopass	y(1:n) = y(1:n) + a x(1:n).
vec_dot	nopass	sum x y over the local n entries and the ranks.
vec_relmax	nopass	max |r / vol| over the ranks, over fmax.
vec_permute	nopass	y(iy(i)) = x(ix(i)).
vec_update	nopass	A work vector between the host and where the hooks compute.
work_place	pass(self)	Put (or release) the operators and work vectors where the hooks compute.
assemble	pass(self)	Finish a distributed operator from triplets.
a_row	pass(self)	The row of the composite operator at a cell.
cell_volume	pass(self)	Cell volume at a level.
classify	pass(self)	Classify a neighbour position.
comp_value	pass(self)	The level-l value of a cell as leaf-cell weights.
dx_of	pass(self)	Cell spacing at a level.
face_cells	pass(self)	The fine cells of a coarse face.
face_weight	pass(self)	V-weighted face weight (area over distance) at a level.
ghost_green	pass(self)	Isolated: g = -G s, the ghost-to-ghost Green's function.
solve_fcg	pass(self)	The FCG solve with given ghost values (or none).
lin	pass(self)	Linear index of a cell in a leaf block.
make_work	pass(self)	Allocate the work vectors and the level maps.
lin_of	pass(self)	Linear index of a cell in a block at a level.
ncb_of	pass(self)	Cells per block at a level.
ncell	pass(self)	Cells along a direction over the domain at a level.
nijk_of	pass(self)	Cells per block per direction at a level.
outside	pass(self)	A cell position beyond a domain face.
locate	pass(self)	Owner and local index of a cell key in a space.
parent3	pass(self)	Parent cell coordinates.
vcycle	pass(self)	The preconditioner.
wrap	pass(self)	Periodic wrap of cell coordinates at a level.

Subroutines

destroy

Destroy the solver.

subroutine destroy(self)

Arguments

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

Call graph

initialize

Gather the leaves, derive the hierarchy, assemble the operator, the ghost fill and the preconditioner, set up the coarsest-level FFT solver: of the dynamic type of coarse when given (the FNL backend passes its cuFFT solver), the FFTW one otherwise. With isolated (issue #87), also the ghost-to-ghost Green's function of the James step. bc (issue #84, E1b): the kind of every face (see the module); a periodic axis of the grid must be periodic in bc and the other way round. Without it, every bounded axis takes ghost-Dirichlet faces. root_cells (#84 E5c): the levels inside the blocks below lmin stop once the root has at most this many cells (default: the root solver's root_cells_max for the root's transform growth; tests force the levels with a small one, the bench measures others).

subroutine initialize(self, grid, field, tolerance, max_iterations, coarse, isolated, bc, root_cells)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
grid	type(grid_object)	in		The grid.
field	type(field_object)	in		The field (block coordinates).
tolerance	real(kind=R8P)	in	optional	Relative residual to stop.
max_iterations	integer(kind=I4P)	in	optional	FCG iterations at most.
coarse	class(poisson_fft_object)	in	optional	Mold of the coarsest-level solver.
isolated	logical	in	optional	Free-space boundaries (James), else as the grid.
bc	integer(kind=I4P)	in	optional	Kinds per face [lo/hi, axis].
root_cells	integer(kind=I4P)	in	optional	Most cells of the root (default: the backend's).

Call graph

solve

Solve A phi = factor V q(iv) - B ghost on the leaves, the mean of the right-hand side removed when the operator is singular (periodic, or no Dirichlet face); phi's ghosts filled (see the module).

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

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_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 V f phi (all ranks).
ghost	real(kind=R8P)	in	optional	The boundary value per face.

Call graph

solve_packed

Solve from the packed interior source of this rank's leaves (block-major, i fastest); then the ghost fill. Isolated (issue #87, P2), the James step: the zero-ghost solve phi0 = A^-1 factor V f, the screening charge s = B^T phi0 on the ghost cells, the ghost values of the free-space solution g = -G s (ghost_green), and the solve with them, phi = A^-1 (factor V f - B g): the map is (A^-1 + A^-1 B G B^T A^-1) factor V, symmetric in the V-weighted form since A and G are. Each of the two solves warm-starts from its own last solution.

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

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
field	type(field_object)	in		The field.
buf	real(kind=R8P)	in		Packed interior source.
factor	real(kind=R8P)	in		Source factor.
phi	real(kind=R8P)	inout		Potential [i,j,k,nb].
energy	real(kind=R8P)	out	optional	sum_cells V f phi (all ranks).
ghost	real(kind=R8P)	in	optional	The boundary value per boundary face.

Call graph

solve_fcg

Solve from the packed interior source of this rank's leaves (block-major, i fastest): FCG on the operator, preconditioned by the V-cycle, from the previous solution when there is one; then the ghost fill.

subroutine solve_fcg(self, field, buf, factor, phi, energy, ghost)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
field	type(field_object)	in		The field.
buf	real(kind=R8P)	in		Packed interior source.
factor	real(kind=R8P)	in		Source factor.
phi	real(kind=R8P)	inout		Potential [i,j,k,nb].
energy	real(kind=R8P)	out	optional	sum_cells V f phi (all ranks).
ghost	real(kind=R8P)	in	optional	The boundary value per boundary face.

Call graph

apply_a

y = A x on this rank's leaf cells (V-weighted), through the work vectors (where the hooks compute).

subroutine apply_a(self, x, y)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
x	real(kind=R8P)	in		Leaf vector.
y	real(kind=R8P)	inout		Result.

Call graph

apply_b

y = B g on this rank's leaf cells: every boundary face adds its weight times its boundary value to its row (A phi + B g is the operator with the boundary values g).

subroutine apply_b(self, g, y)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
g	real(kind=R8P)	in		Boundary value per boundary face.
y	real(kind=R8P)	inout		Result, a leaf vector.

apply_bt

s = B^T x: every boundary face takes its weight times the value of its row.

subroutine apply_bt(self, x, s)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
x	real(kind=R8P)	in		Leaf vector.
s	real(kind=R8P)	inout		Result, one value per boundary face.

Call graph

apply

y = M x: fetch the halo of x from the owners, then the local rows.

subroutine apply(self, m, x, y)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
m	type(dcsr)	in		Operator.
x	real(kind=R8P)	in		Local column-space vector.
y	real(kind=R8P)	inout		Result (rows).

Call graph

vec_fill

x(1:n) = a.

subroutine vec_fill(x, a, n)

Arguments

Name	Type	Intent	Attributes	Description
x	real(kind=R8P)	inout		Vector.
a	real(kind=R8P)	in		Value.
n	integer(kind=I4P)	in		Entries.

Call graph

vec_copy

y(1:n) = x(1:n).

subroutine vec_copy(x, y, n)

Arguments

Name	Type	Intent	Attributes	Description
x	real(kind=R8P)	in		Source.
y	real(kind=R8P)	inout		Destination.
n	integer(kind=I4P)	in		Entries.

Call graph

vec_lincomb

w(1:n) = x(1:n) + a y(1:n) (w distinct from x and y).

subroutine vec_lincomb(x, a, y, w, n)

Arguments

Name	Type	Intent	Attributes	Description
x	real(kind=R8P)	in		First vector.
a	real(kind=R8P)	in		Coefficient.
y	real(kind=R8P)	in		Second vector.
w	real(kind=R8P)	inout		Result.
n	integer(kind=I4P)	in		Entries.

Call graph

vec_axpy

y(1:n) = y(1:n) + a x(1:n).

subroutine vec_axpy(a, x, y, n)

Arguments

Name	Type	Intent	Attributes	Description
a	real(kind=R8P)	in		Coefficient.
x	real(kind=R8P)	in		Vector.
y	real(kind=R8P)	inout		Updated vector.
n	integer(kind=I4P)	in		Entries.

Call graph

vec_permute

y(iy(i)) = x(ix(i)), i = 1..n.

subroutine vec_permute(x, ix, y, iy, n)

Arguments

Name	Type	Intent	Attributes	Description
x	real(kind=R8P)	in		Source.
ix	integer(kind=I4P)	in		Source entries.
y	real(kind=R8P)	inout		Destination.
iy	integer(kind=I4P)	in		Destination entries.
n	integer(kind=I4P)	in		Entries.

Call graph

vec_update

A work vector to (to_device) or from where the hooks compute: nothing on the host.

subroutine vec_update(x, to_device)

Arguments

Name	Type	Intent	Attributes	Description
x	real(kind=R8P)	inout		Work vector.
to_device	logical	in		Direction.

Call graph

work_place

Put (attach) the operators and the work vectors where the hooks compute, or release them: nothing on the host.

subroutine work_place(self, attach)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
attach	logical	in		Put, or release.

Call graph

assemble

Finish a distributed operator: rows sorted; column keys resolved to local slots, or to halo slots grouped by owner; the exchange plan agreed with the owners (they receive the keys and return their local indices).

subroutine assemble(self, t, nrow, space, m)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
t	type(triplets)	inout		Triplets (emptied).
nrow	integer(kind=I4P)	in		Rows.
space	integer(kind=I4P)	in		Column space.
m	type(dcsr)	inout		Operator.

Call graph

a_row

The row of the composite operator at cell cc of level lv (V-weighted): same-level faces, and the fine and the coarse side of 2:1 faces (see the module).

subroutine a_row(self, lv, cc, keys, vals)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.
cc	integer(kind=I4P)	in		Cell coordinates at lv.
keys	integer(kind=I8P)	inout	allocatable	Column keys.
vals	real(kind=R8P)	inout	allocatable	Values.

Call graph

face_cells

The keys of the 2^(nd-1) fine cells (level lf) on the coarse face that contains the fine cell fc's face normal to d: fc's normal index, every pair of tangential indices of its coarse parent.

subroutine face_cells(self, lf, fc, d, keys_f)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lf	integer(kind=I4P)	in		Fine level.
fc	integer(kind=I4P)	in		A fine cell of the face.
d	integer(kind=I4P)	in		Normal direction.
keys_f	integer(kind=I8P)	out	allocatable	Keys.

Call graph

ghost_green

the infinite lattice, from ghost cell to ghost cell (G(n) / H, the unit lattice's Green's function over the coarsest spacing). Rank 0 gathers the charges into the coarsest box widened by the ghost layer (every face has its own ghost cell: a boundary leaf is at the coarsest level), convolves with the padded FFT (gfft, which solves L_H phi = factor f by the kernel H^2 G, so factor = -1 / H^3), and scatters the values back.

subroutine ghost_green(self, s, g)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
s	real(kind=R8P)	in		Screening charge per boundary face.
g	real(kind=R8P)	inout		Ghost value per boundary face.

Call graph

comp_value

Append w times the level-lv value of cell cc as leaf-cell weights: the leaf itself; the mean of the leaves inside, if refined; the containing coarser leaf, if inside one.

Attributes: recursive

subroutine comp_value(self, lv, cc, w, keys, wts)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.
cc	integer(kind=I4P)	in		Cell coordinates.
w	real(kind=R8P)	in		Weight.
keys	integer(kind=I8P)	inout	allocatable	Keys.
wts	real(kind=R8P)	inout	allocatable	Weights.

Call graph

make_work

Allocate the work vectors (FCG, ghost fill, the V-cycle's levels and root) and the maps of the leaves into the level spaces (#84 E6a).

subroutine make_work(self)

Arguments

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

Call graph

locate

Owner rank and local index of a cell key in a space (the composite leaves, or a level).

subroutine locate(self, space, key, owner, lidx)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
space	integer(kind=I4P)	in		Space.
key	integer(kind=I8P)	in		Cell key.
owner	integer(kind=I4P)	out		Owner rank.
lidx	integer(kind=I4P)	out		Local index on the owner.

Call graph

smooth

Red-black Gauss-Seidel sweeps on level lv: A_l e = rhs (the ghost term already in rhs).

subroutine smooth(self, lv, e, rhs, order)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
lv	integer(kind=I4P)	in		Level.
e	real(kind=R8P)	inout		Correction.
rhs	real(kind=R8P)	in		Right-hand side.
order	integer(kind=I4P)	in		0: colour 0 first, 1: colour 1 first.

Call graph

vcycle

The preconditioner: z = M r, one V-cycle over the levels (see the module), down to the root (lbot), solved exactly by FFT on every rank (#84 E5c: below lmin the levels coarsen inside the blocks). On the work vectors, through the hooks (#84 E6a); the root goes through the host.

subroutine vcycle(self, r_leaf, z_leaf)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
r_leaf	real(kind=R8P)	in		Leaf residual (V-weighted).
z_leaf	real(kind=R8P)	inout		Leaf correction.

Call graph

kdecode

Decode a cell key.

Attributes: pure

subroutine kdecode(key, lv, cc)

Arguments

Name	Type	Intent	Attributes	Description
key	integer(kind=I8P)	in		Key.
lv	integer(kind=I4P)	out		Level.
cc	integer(kind=I4P)	out		Cell coordinates.

Call graph

hash_init

A table for about n keys (capacity: the first prime above 2 n).

subroutine hash_init(h, n)

Arguments

Name	Type	Intent	Attributes	Description
h	type(hash_table)	inout		Table.
n	integer(kind=I4P)	in		Expected keys.

Call graph

hash_free

Free a table.

subroutine hash_free(h)

Arguments

Name	Type	Intent	Attributes	Description
h	type(hash_table)	inout		Table.

Call graph

hash_put

Insert or overwrite (linear probing).

subroutine hash_put(h, key, val)

Arguments

Name	Type	Intent	Attributes	Description
h	type(hash_table)	inout		Table.
key	integer(kind=I8P)	in		Key (>= 0).
val	integer(kind=I4P)	in		Value (> 0).

Call graph

trip_add

Append a triplet.

subroutine trip_add(t, row, key, val)

Arguments

Name	Type	Intent	Attributes	Description
t	type(triplets)	inout		Buffer.
row	integer(kind=I4P)	in		Row.
key	integer(kind=I8P)	in		Column key.
val	real(kind=R8P)	in		Value.

Functions

symmetry_defect

|x . A y - y . A x| / (|x| |A y|) over the ranks, for two leaf vectors: the operator's symmetry, verified.

Returns: real(kind=R8P)

function symmetry_defect(self, x, y) result(defect)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	inout		The solver.
x	real(kind=R8P)	in		First vector.
y	real(kind=R8P)	in		Second vector.

Call graph

vec_dot

sum_i x_i y_i over the local n entries and the ranks.

Returns: real(kind=R8P)

function vec_dot(x, y, n) result(d)

Arguments

Name	Type	Intent	Attributes	Description
x	real(kind=R8P)	in		First vector.
y	real(kind=R8P)	in		Second vector.
n	integer(kind=I4P)	in		Local entries.

Call graph

vec_relmax

max |r / vol| over the ranks, over fmax (when positive).

Returns: real(kind=R8P)

function vec_relmax(r, vol, n, fmax) result(rm)

Arguments

Name	Type	Intent	Attributes	Description
r	real(kind=R8P)	in		Residual.
vol	real(kind=R8P)	in		Cell volumes.
n	integer(kind=I4P)	in		Local entries.
fmax	real(kind=R8P)	in		Normalisation.

Call graph

cell_volume

Cell volume at level lv (the product over the active directions).

Attributes: pure

Returns: real(kind=R8P)

function cell_volume(self, lv) result(v)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.

Call graph

classify

Classify the level-lv cell position nb3: a leaf of lv, refined, or inside a leaf of lv-1 (fatal otherwise: the grid is not 2:1 balanced across faces).

Returns: integer(kind=I4P)

function classify(self, lv, nb3) result(cls)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.
nb3	integer(kind=I4P)	in		Cell coordinates.

Call graph

dx_of

Cell spacing at level lv (null directions keep the level-0 one).

Attributes: pure

Returns: real(kind=R8P)

function dx_of(self, lv) result(dx)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.

Call graph

face_weight

The V-weighted face weight at level lv along d: the face area over the centre distance, V / dx_d^2.

Attributes: pure

Returns: real(kind=R8P)

function face_weight(self, lv, d) result(w)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.
d	integer(kind=I4P)	in		Direction.

Call graph

lin

Linear index of cell (i,j,k) in a block, i fastest.

Attributes: pure

Returns: integer(kind=I4P)

function lin(self, ijk) result(l)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
ijk	integer(kind=I4P)	in		Cell (1-based).

Call graph

lin_of

Linear index of a cell (1-based) in a block at level lv (nijk_of(lv) cells).

Attributes: pure

Returns: integer(kind=I4P)

function lin_of(self, lv, ijk) result(l)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.
ijk	integer(kind=I4P)	in		Cell (1-based).

Call graph

ncb_of

Cells per block at level lv.

Attributes: pure

Returns: integer(kind=I4P)

function ncb_of(self, lv) result(n)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.

Call graph

ncell

Cells along direction d over the domain at level lv: the blocks at max(lv, lmin) times their cells.

Attributes: pure

Returns: integer(kind=I4P)

function ncell(self, lv, d) result(n)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.
d	integer(kind=I4P)	in		Direction.

Call graph

nijk_of

Cells per block per direction at level lv: nijk at lmin and above, halved per level below it along the active directions (the levels inside the blocks, #84 E5c).

Attributes: pure

Returns: integer(kind=I4P)

function nijk_of(self, lv) result(nl)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.

Call graph

parent3

Parent cell coordinates (one level coarser) along the active directions.

Attributes: pure

Returns: integer(kind=I4P)

function parent3(self, cc) result(pc)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
cc	integer(kind=I4P)	in		Cell coordinates.

Call graph

outside

A cell position at level lv beyond a domain face along a bounded active direction (never along a periodic one).

Attributes: pure

Returns: logical

function outside(self, lv, cc) result(beyond)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.
cc	integer(kind=I4P)	in		Cell coordinates.

Call graph

wrap

Periodic wrap of cell coordinates at level lv along the periodic active directions (bounded ones unchanged).

Attributes: pure

Returns: integer(kind=I4P)

function wrap(self, lv, cc) result(cw)

Arguments

Name	Type	Intent	Attributes	Description
self	class(poisson_mg_object)	in		The solver.
lv	integer(kind=I4P)	in		Level.
cc	integer(kind=I4P)	in		Cell coordinates.

Call graph

gdot

sum_i x_i y_i over the local n entries and the ranks.

Returns: real(kind=R8P)

function gdot(x, y, n) result(s)

Arguments

Name	Type	Intent	Attributes	Description
x	real(kind=R8P)	in		First vector.
y	real(kind=R8P)	in		Second vector.
n	integer(kind=I4P)	in		Local entries.

Call graph

bkey

Key of a block or cell: (i, j, k, level) packed in 19 + 19 + 19 + 7 bits.

Attributes: pure

Returns: integer(kind=I8P)

function bkey(c4) result(key)

Arguments

Name	Type	Intent	Attributes	Description
c4	integer(kind=I4P)	in		Coordinates (i, j, k, level).

Call graph

ckey

Key of a cell at level lv.

Attributes: pure

Returns: integer(kind=I8P)

function ckey(lv, cc) result(key)

Arguments

Name	Type	Intent	Attributes	Description
lv	integer(kind=I4P)	in		Level.
cc	integer(kind=I4P)	in		Cell coordinates.

Call graph

hash_get

The value of a key, 0 when absent.

Attributes: pure

Returns: integer(kind=I4P)

function hash_get(h, key) result(val)

Arguments

Name	Type	Intent	Attributes	Description
h	type(hash_table)	in		Table.
key	integer(kind=I8P)	in		Key.

Call graph