ADAM Documentation

Appearance

The library ​

The src/lib directory contains the ADAM SDK: the physics-agnostic, hardware-portable core upon which all ADAM solver applications are built. It is organised into four subdirectories — one portable common layer and three hardware backends — described in detail below. PRISM and FLUME build the common layer and the FNL backend; no maintained application builds NVF or GMP.

Design principles ​

The library is structured around two orthogonal axes of variability:

  • Physics axis — handled entirely by the application layer (src/app). The SDK is equation-agnostic: it knows nothing about Navier-Stokes, Maxwell, or any other PDE.
  • Hardware axis — handled by backend selection at compile time. _FNL selects the GPU backend: OpenACC with nvfortran (*-fnl-nvf modes), or the same kernels as OpenMP offloading with amdflang (*-fnl-omp-amd modes, DEV_OMP). _NVF and _GMP select the CUDA Fortran and OpenMP-offloading libraries, which no PRISM or FLUME mode builds.

The result is a matrix of N_equations × N_backends combinations served by a single shared source tree.

Directory overview ​

DirectoryRoleParallelism
common/Portable core objects: grid, tree, field, numerics, I/O, forestCPU · MPI · OpenMP
fnl/GPU backend (FUNDAL-based), built by PRISM and FLUMEOpenACC or OpenMP target · MPI
nvf/CUDA Fortran backend; no build mode builds itCUDA Fortran · MPI
gmp/OpenMP target offloading backend; no build mode builds itOpenMP target · MPI

The common library ​

The common library is the portable, physics-agnostic core of the ADAM SDK. The GPU objects (fnl; the unbuilt nvf, gmp) are separate types, initialised from these objects passed as arguments; every solver application (PRISM and FLUME are the maintained ones) composes them through realm_object. No GPU-specific code lives here — the entire directory compiles with any standard Fortran 2003+ compiler.

The aggregate entry point adam_common_library re-exports all modules below; a single use adam_common_library statement in application code exposes the complete API.

Contents ​

AMR data structures ​

adam_tree_object — octree/quadtree hierarchy ​

The tree is a CPU-resident hash-map keyed on 64-bit Morton indices that linearise the four octree coordinates (level, bx, by, bz) into a single integer. It provides O(1) block insertion, deletion, and neighbour lookup during AMR update phases. Morton ordering clusters geometrically adjacent blocks in index space, minimising ghost-cell communication volume across MPI ranks.

The refinement ratio is selectable: 4 (quadtree) or 8 (octree). Full 26-neighbourhood support covers all face, edge, and corner adjacencies needed by high-order stencils.

The tree is never accessed by numerical kernels during time integration — it is consulted only during AMR update cycles (marking → refinement/coarsening → load rebalancing).

adam_tree_node_object — individual tree node ​

Represents a single block in the octree. Key fields:

FieldTypePurpose
codeI8PMorton key
refinement_neededI4PAMR marker (TO_BE_REFINED, TO_BE_DEREFINED, TO_NOT_TOUCH)
myrankI4POwning MPI rank
block_indexI8PCompact index into the field array
neighbor(26)tree_node_neighbor_objectFace/edge/corner adjacency data
surface_stl_distanceR8PSigned distance to nearest IB surface

adam_tree_bucket_object — hash-map bucket ​

Double-linked list container underpinning the tree's hash-map. Stores nodes by Morton key and provides O(1) add_node, remove_node, and has_code operations. Supports iterator-based traversal via the traverse method and a user-supplied callback conforming to iterator_interface.

Grid and field ​

adam_grid_object — block-structured Cartesian grid ​

Manages the geometric description of a single structured block and the global domain layout. Provides cell/node coordinate computation, mesh metrics, ghost-cell layout, and boundary condition assignment.

Key configuration ([grid] INI section):

ParameterMeaning
ni, nj, nkCells per block per direction
ngcGhost-cell width (typically 3 for WENO5, 4 for WENO7)
emin_*, emax_*Domain extents
bc_type(6)Boundary condition type on each of the 6 domain faces

Notable methods: block_emin/emax, cell_xyz, node_xyz, compute_metrics, compute_weight_neighbor, load_from_ini_file.

adam_field_object — block geometry, ghost exchange, refinement ​

The field object holds the per-block geometry (emin, emax, dxyz, node and cell coordinates), the ghost exchange (update_ghost_local, update_ghost_mpi, and the reverse folds reduce_ghost_local, reduce_ghost_mpi), the refinement and derefinement of block data (adapt), and the restart blocks (save_blocks, load_blocks). The state arrays are owned by the application (e.g. flume_common_object%q), dense and contiguous over all blocks of an MPI rank:

q(nv, 1-ngc:ni+ngc, 1-ngc:nj+ngc, 1-ngc:nk+ngc, nb)
DimensionMeaning
nvPhysical variables (density, momenta, energy, …)
ni, nj, nkCell indices within a block, including ghost cells
nbCompact sequential block index (1 … current block count)

Fortran column-major order keeps nv contiguous in memory, enabling coalesced reads across CUDA threads or OpenACC/OpenMP SIMD lanes. The application allocates it once, to the block count nb that realm_object%initialize derives from the memory budget, and reuses it across all Runge-Kutta stages — no dynamic allocation during time integration.

Communication and maps ​

adam_maps_object — index translation and MPI communication patterns ​

Bridges the tree (Morton-key addressed) and the field (compact-index addressed). Rebuilt only when the AMR topology changes; invisible to numerical kernels during time integration.

Local maps (updated every AMR cycle):

MapPurpose
local_mapBlock index changes of the rank's nodes (old → new) for block reordering
local_map_ghostBlock-to-ghost-block adjacency for exchange
local_map_ghost_cellCell-level ghost stencil offsets
local_map_bc_face/edge/corner/crownBoundary condition application stencils

MPI communication lists (per-rank send/receive schedules):

  • comm_map_send/recv — block-level exchange for AMR rebalancing
  • comm_map_send/recv_ghost — ghost-cell exchange buffers before stencil sweeps
  • send_buffer_ghost, recv_buffer_ghost — pre-allocated MPI staging buffers

Key methods: initialize, make_comm_local_maps, blocks_reorder.

Numerical methods: spatial ​

adam_weno_object — WENO reconstruction ​

Weighted Essentially Non-Oscillatory finite difference reconstruction. Upwind and centered variants are available at orders 1, 3, 5, 7, and 9 (upwind) and 2, 4, 6, and 8 (centered); the nonlinear weights are Jiang–Shu (js, default) or scale-invariant (si). The stencil half-width S_MAX governs the ghost-cell requirement.

Scheme selectors:

ConstantScheme
WENO_U_1 … WENO_U_9Upwind WENO, orders 1/3/5/7/9
WENO_C_2 … WENO_C_8Centered WENO, orders 2/4/6/8

Reference: Chi-Wang Shu, NASA/CR-97-206253 (1997).

adam_fdv_operators_library — finite difference/volume operators ​

Spatial derivative operators at configurable order for use in application physics layers. All operators are available in three flavours:

FlavourDescription
Finite difference, centeredStandard FD stencils on cell centres
Finite volume, centeredCell-average-based FV operators
Finite volume, upwindWENO-reconstructed upwind FV fluxes

Operators: gradient, divergence, curl, Laplacian, and scalar derivatives of orders 1–6. Abstract compute_*_fdv_interface procedure pointers allow backend dispatch without branching.

adam_riemann_euler_library — Euler equation Riemann solvers ​

Convective interface flux computation for the compressible Euler equations. All solvers accept left/right primitive states and return the numerical flux; the optional lmax output gives the maximum wave speed for CFL estimation.

ProcedureSolver
compute_riemann_euler_exactExact solver (Newton iteration on Rankine-Hugoniot)
compute_riemann_euler_hllcHLLC (Harten–Lax–van Leer–Contact)
compute_riemann_euler_hllc_lmHLLC with low-Mach correction
compute_riemann_euler_hllemHLLEM (entropy-fix variant)
compute_riemann_euler_llfLocal Lax–Friedrichs
compute_riemann_euler_tsTrivial (first-order upwind)

Numerical methods: temporal ​

adam_rk_object — explicit Runge-Kutta integration ​

ConstantSchemeStagesOrder
RK_1Forward Euler11
RK_2Heun22
RK_3TVD-RK333
RK_SSP_11SSP(1,1), forward Euler11
RK_SSP_22SSP(2,2)22
RK_SSP_33SSP(3,3)33
RK_SSP_54SSP(5,4)54
RK_YOSHIDAYoshida symplectic44

Stage storage q_rk(nv, ni, nj, nk, nb, nrk) is pre-allocated; no allocation occurs during time-stepping. Configuration section: [runge_kutta].

The next three integrators are library objects that realm_object declares but does not initialize (the calls are commented out); no application can use them today, and PRISM refuses them at init (issue #55).

adam_leapfrog_object — leapfrog with RAW filter ​

Second-order leapfrog scheme with optional Robert–Asselin–Williams noise filter:

U^{n+2} = U^{n} + 2·Δt·R(t^{n+1}, U^{n+1})
Δ = ν/2·(U^{n} - 2·U^{n+1} + U^{n+2})    ! RAW correction

Filter coefficient nu (default 0.01) and Williams parameter alpha (default 0.53) are configurable. Configuration section: [leapfrog].

adam_blanes_moan_object — Blanes–Moan splitting integrator ​

Geometric integrator for separable Hamiltonian systems H = A + B. Two schemes available:

ConstantOrderStages
BM_SCHEME444
BM_SCHEME666

Coefficients a(:) and b(:) implement the alternating A-B operator splitting. Configuration section: [blanes_moan].

adam_cfm_object — Commutator-Free Magnus integrators ​

Exponential integrators for non-autonomous ODEs of the form dU/dt = A(t)·U. Computes U^{n+1} = exp(α·R)·U^n without evaluating the matrix exponential directly (commutator-free formulation).

ConstantOrderStagesExponentials
CFM_44——
CFM_66——
CFM_88——

Stage residuals dq and field buffer q are pre-allocated. Configuration section: [commutator_free_magnus].

Physics support ​

adam_eos_ic_object — ideal compressible gas EOS ​

Thermodynamic model for a calorically perfect ideal gas. Stores primary properties and derives all secondary quantities on initialisation.

PropertySymbolDescription
cp, cv—Specific heats at constant pressure/volume
gγHeat capacity ratio
gm1, gp1γ−1, γ+1Pre-computed combinations
delta, eta(γ−1)/2, 2γ/(γ−1)Isentropic flow coefficients
muμDynamic viscosity
kdkThermal conductivity

Key procedures: density, pressure, temperature, internal_energy, total_energy, speed_of_sound, primitive2conservative, conservative2primitive. Configuration section: [physics_specie_*].

adam_ib_object — immersed boundary method ​

Manages solid bodies immersed in the Cartesian grid. Computes the signed-distance field phi(nv, ni, nj, nk, nb) via an eikonal solver and applies wall boundary conditions at solid surfaces.

Supported solid definitions:

TypeConstantDescription
Sphere / circleIB_ANALYTICAL_SPHEREDefined by centre, radius, axis
Circle (XY plane)IB_ANALYTICAL_CIRCLE2D equivalent
RectangleIB_ANALYTICAL_RECTANGLEDefined by centre, edge lengths, axis
Surface mesh—file.off: no implementation, refused at load (issue #57)

Wall boundary condition types:

ConstantValueType
BCS_VISCOUS1Viscous (no-slip) wall
BCS_EULER2Inviscid (slip) wall

Configuration section: [solid]. Number of eikonal integration steps configurable via n_eikonal.

adam_amr_object — AMR refinement markers ​

Configures the criteria that drive block refinement and coarsening. Multiple markers can be combined; each marker specifies an independent criterion applied at every AMR update cycle.

Marker modes:

ConstantValueCriterion
AMR_GEO1Geometry: an IB solid (AMR_GEO_SOLID), a primitive box (AMR_GEO_PRIMITIVE_BOX, up to target_level) or an STL surface (AMR_GEO_STL, refused by both PRISM and FLUME). In PRISM the AMR_GEO_SOLID and AMR_GRAD branches do nothing (it implements the primitive box and AMR_TV)
AMR_GRAD2Field gradient magnitude
AMR_TV3Field total variation
AMR_LOHNER4Löhner normalised second derivative (refine_tol, derefine_tol; issue #74)
AMR_JEANS5Jeans length resolution (jeans_cells, jeans_cells_derefine; issue #91)
AMR_USER6A criterion the program registered (name; issue #110)

Registered criteria (issue #110, adam_amr_criterion_object). A refinement criterion can be code instead of a mode: a type extending the abstract amr_criterion_object (a name and a deferred vote(realm, b) returning TO_BE_REFINED, TO_NOT_TOUCH or TO_BE_DEREFINED for block b of realm), registered with the realm before the simulation starts and selected in the input by [amr_marker_N] mode = 6, name = <its name>. No edit of the library or of the app: the program, a copy of the app's driver, registers it. In FLUME, register_amr_criterion(criterion) (on every realm of a forest; the realm keeps a pointer, so the program owns the instance, which must have the TARGET attribute and outlive the run: a variable of the main program or a SAVEd one); a name nobody registered stops the run at initialisation, and mark_runtime combines its votes with the other markers' (refine if any marker votes refine, derefine if all do, never below iu_ref_levels). The vote runs on the host, for the blocks of the rank, at runtime regrids; realm%adam gives the grid, the tree and the block metrics, and a criterion that needs the app's fields selects the app's type. The built-in markers are criteria of the same kind (issue #110, P4.2): amr_criterion_object also has mark(realm, votes, initial), the votes of every block of the rank at once, whose default asks vote block by block, and runtime_only (skipped by the AMR of the initial condition; the default, for a registered criterion). FLUME builds one list, amr_marker_list, of the criteria of its [amr_marker_N] sections in input order, a registered criterion where mode = 6 and a FLUME type for each built-in (flume_box_marker, whose per-block vote is the box; the solid, gradient, Löhner and Jeans markers, which compare neighbours or dilate their marks, override mark and call the marking methods they wrap); amr_update (the initial AMR) and mark_runtime loop that list. A marker FLUME does not implement (total variation, an STL geometry, an unknown mode) stops the run when the list is built, at initialisation. Example and test: src/tests/flume/verification/regrid/test_flume_user_criterion.F90 (fobis mode test-flume-user-criterion-cpu-gnu), a criterion that votes as the built-in box marker, bitwise the same run on one realm and on two (multirealm/check.sh leg 10).

Each marker carries its own thresholds and a reference to the field variable. The markers refine the grid at initialization, iterated iters times to enforce 2:1 level balance. A realm that regrids during the run (FLUME, issue #74) re-applies them every frequency time steps through regrid_forest; PRISM refines at initialization only. Configuration section: [amr].

adam_flail_object — FLAIL, the elliptic solver library (issue #84) ​

One entry point for the elliptic solves of every app, the realm component flail (self%flail). It is built on the two solvers below, verified in #81/#87: the FFT solver on uniform grids and the multigrid-preconditioned CG solver on AMR grids. The legacy smoothers (next section) were deleted once PRISM moved here (#84, E4).

ProcedureRole
setup(grid, field, isolated, backend[, option_name][, bc][, dg_order])Choose the solver for the current grid (FLAIL_BACKEND_AUTO: the FFT solver on a uniform grid, the multigrid otherwise; _FFT, _MULTIGRID); call it on the initial grid and after every regrid. bc(2,3) (E1c): the kind of every face, forwarded (FLAIL_BC_PERIODIC, _DIRICHLET, _NEUMANN, _DIRICHLET_GHOST, re-exported here); without it, periodic or isolated as in E0; isolated with bc is fatal. dg_order (E2a): the wide D(G) of order 2, 4, 6, 8, on the FFT path only (fatal on AMR: not defined across 2:1 faces; #84 E7, closed with no consumer)
solve(field, q, iv, factor, phi, energy[, ghost]), solve_packed(field, buf, factor, phi, energy[, ghost])Forward to the chosen solver; energy = sum f phi dV; ghost the boundary value of each boundary face (mg%nbnd, mg%bnd_kind, mg%bnd_dir, mg%bnd_x), multigrid path only (fatal on the FFT path, which solves the homogeneous problem)
iterations(), is_multigrid(), destroy()The last FCG iteration count (0 for the FFT solver), the solver in use, cleanup

A vector whose components take different kinds (PRISM's vector potential at a PEC face: Neumann for the normal component, Dirichlet for the others) takes one FLAIL per set of kinds: the operators differ (the FFT extensions, the singularity), and the instances coexist (test_poisson_mg --flail, bitwise). The mapping from an app's boundaries to the kinds belongs to the app.

PRISM's t0 solves (#84, E3b). impose_pic_fields_time_zero (adam_prism_common_object.F90) solves D(G) phi = -rho/eps (PIC, then D <- D - eps G phi) and D(G) A_d = -mu J_d per component (then B <- B + curl A) through the realm's flail, one setup per component, dg_order = fdv_order, on the FFT path: periodic stays periodic, PRISM's Dirichlet is FLAIL_BC_ZERO_GHOSTS, PEC is FLAIL_BC_DIRICHLET (odd) for the potential and the tangential components of A, FLAIL_BC_NEUMANN (even) for the normal one, as the legacy ghost rules; an exact/open face (PRISM's Neumann, Silver_Muller) is FLAIL_BC_FREE_GHOSTS (E3c), a null axis periodic. A refined grid is refused (the wide operator across 2:1 faces is not defined; #84 E7 closed with no consumer). The FNL backend allocates the cuFFT solver; every PRISM mode links an FFT library (FFTW on the host, cuFFT on NVIDIA FNL). Regression cases pic-pec-t0, rmf-pec-t0, pic-open-t0 and rmf-open-t0 enter the two solves. Three limits: the CG solves (zero-ghost and free-space faces) gather the grid on rank 0 (and on its GPU on FNL), about 24 B per transform cell (#84 decision 3; the PEC solves are distributed over the z-planes since E5d, 8 N cells with reflected faces); and PRISM's centred div∘grad annihilates the checkerboard, so with PEC faces the component of rho along (-1)^(i+j+k) stays in div D - rho (no D can carry it; the legacy smoother could not either); the free-space ghost sum is a direct one on rank 0, O(ghost-slab cells x cells) per component (a lattice free-space solve is not implemented: #84 E7 closed with no consumer).

A backend chooses its FFT engine by allocating the fft component before setup (FLUME's FNL backend allocates the cuFFT one). The boundary kinds per face, FLAIL_BC_PERIODIC | _DIRICHLET | _NEUMANN | _DIRICHLET_GHOST | _ZERO_GHOSTS | _FREE_GHOSTS, live in adam_flail_parameters (a module of their own: the solvers use them, and the facade uses the solvers). Phases and status: issue #84.

Scaling (#84, E5-E6). Before E5 both solvers gathered to rank 0: the FFT solver the whole uniform grid, the multigrid its coarsest level (the whole base grid) in every V-cycle and the James step's ghost charges in every isolated solve. Since E5d the uniform FFT is distributed over z-slabs; the multigrid's root is redundant on every rank up to 32^3 cells and on the slab FFT above it. poisson_fft_object%wall(3) accumulates the gather, rank-0 solve and scatter times of solve_blocks on each rank (the caller resets it). exe/bench_poisson (modes bench-poisson-gnu, bench-poisson-fnl-nvf) times FLAIL's solves on the cube [-1,1]^3 (Gaussian source, uniform or with the central box refined), and src/tests/poisson/bench.sh [cpu|fnl] [--weak] [--strong] [--np-max N] runs the ladders and prints one row per run: solve time (min over the repetitions, max over the ranks) and rank 0's phases. Baseline (E5a, medooza CPU, strong scaling, 2.1 M uniform or 3.9 M AMR cells): the solve time stays flat from 1 to 16 ranks on the FFT paths (rank 0's part 83% periodic, 99% isolated at 16 ranks) and on the isolated multigrid (36.6 s of 39.9 s per solve in rank 0's coarsest level, a 258^3 odd extension solved 42 times); the periodic multigrid falls from 9.4 s to 2.6 s, 45% of it rank 0's. E5b sized the isolated transforms for speed (above); E5c coarsens the multigrid inside the blocks below the base level and solves the root on every rank (medooza CPU, 3.9 M AMR cells: isolated 39.9 -> 6.5 s per solve at 16 ranks, periodic 2.6 -> 2.3 s, but 9.2 -> 12.7 s at one rank, where one FFT of the 128^3 base level beats smoothing it with the generic CSR sweep; so the cuFFT root kept the exact coarsest level up to 128^3, until E6 measured 32^3 faster there too); E5d distributes the uniform FFT (medooza CPU, 2.1 M cells, 1 -> 16 ranks: periodic 0.115 -> 0.026 s, isolated 0.729 -> 0.108 s per solve, where they stayed at 0.12-0.14 s and 0.71-0.73 s; weak, np 8, 16.7 M cells: 0.81 -> 0.24 s and 5.7 -> 1.5 s). E6 runs the multigrid on the device and sets the cuFFT root to 32^3 (medooza FNL, 3.9 M AMR cells, 1 / 2 / 4 GPUs, E5d -> E6: periodic 7.94 / 4.21 / 2.46 -> 0.33 / 0.19 / 0.12 s, isolated 18.2 / 9.9 / 7.3 -> 0.85 / 0.52 / 0.43 s); E6b sizes the host root by its transform, 128^3 cells, now that the transform is distributed (medooza CPU, np 1 / 4 / 16: periodic 12.4 / 3.9 / 2.5 -> 8.6 / 2.7 / 1.6 s with a 128^3 root, isolated 30.4 / 9.6 / 6.3 -> 26.8 / 8.5 / 5.2 s with 64^3, its transform 130^3).

Legacy FLAIL smoothers (deleted in issue #84, E4) ​

adam_flail_legacy_object (the realm component flail_legacy, the [linear-algebra] section: multigrid, Gauss-Seidel and SOR smoothers of PRISM's D(G)) is deleted: PRISM's t0 solves run on FLAIL since #84 E3, [PIC] elliptic_correction was retired in E3d, and the only other caller was the dead divergence_correction = poisson code. FLAIL has no counterpart to the [linear-algebra] keys (the FFT solve is direct, the CG tolerance its own), so the section is gone from the committed inputs and the removed-keys lint rejects it; an input that still carries it is not read there.

adam_poisson_fft_object — FFT Poisson solver (uniform grids: periodic, Dirichlet/Neumann faces, or isolated) ​

Solves L_h phi = factor (f - mean f) (periodic) or L_h phi = factor f with phi -> 0 at infinity (isolated) for one variable f of a conservative array q(nv,i,j,k,b), with L_h the second-order 7-point Laplacian (5-point in 2-D; a null direction has no derivative) (#81).

  • Periodic (P2): each Fourier mode is divided 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 is at round-off.

  • Isolated (P4), Hockney–Eastwood: the source sits in a zero-padded array of m >= 2 (n + ngc) cells per direction (the smallest even size with no prime factor above 13, #84 E5b: FFTW took 3.3 times as long on 262 = 2 x 131 as on 264) and is convolved (one r2c, a product, one c2r) with the lattice Green's function of L_h, h^2 G(n), whose spectrum is computed once. The circular convolution equals the free-space one on every cell the scatter reads, ghosts included. G comes from adam_lattice_green_library (below): its residual is about 1e-13 per offset, so the residual of phi is that summed over the source (2e-11 of factor f at 64^3), not round-off.

  • Dirichlet/Neumann faces (issue #84, E1a): bc(2,3) gives each face a homogeneous kind; an axis is periodic, or each of its faces is Dirichlet (zero at the face) or Neumann (zero normal derivative), 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 faces agree and 4n when they differ, and the periodic transform of the extension solves the 7-point operator with those conditions exactly (round-off residual, face placement of the wall: second order). Without a Dirichlet face the source mean is dropped, as on a periodic grid. The energy sums the interior only.

  • The wide operator D(G) (issue #84, E2a): dg_order = p (2, 4, 6, 8; 0 the compact Laplacian, the default) solves PRISM's elliptic operator, the order-p centred first derivative composed with itself per axis (a stencil of 2p + 1 cells per axis), eigenvalue -sum_d s_p(theta_d)^2 / dx_d^2, s_p(theta) = sum_j 2 a_j sin(j theta). Every centred derivative annihilates the checkerboard, s_p(pi) = 0, set exactly (sin(pi) is 1e-16): the modes with each theta_d in {0, pi} are in the kernel and dropped, Dirichlet faces or not (G of them is zero). The extensions are reflections, which the wide stencil commutes with, so the kinds per face hold, the scatter filling p ghost layers; isolated boundaries take the compact operator only. Why a direct solve: on the prototype (src/tests/poisson/wide_dg_proto.py) a compact preconditioner, even inverted exactly, needs O(N) PCG iterations for D(G) (86 → 677 at p = 2, N = 32 → 256: the checkerboard branch it does not see), while the wide p = 2 operator (2^d compact Laplacians on the sub-lattices of spacing 2h) is spectrally equivalent (17–30 iterations, independent of N). PRISM's zero-ghost Dirichlet face is not a reflection: wide_dg_zero_ghost.py measures it (#84 E2b).

  • 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 (one layer read; the multigrid maps it the same way). With the wide D(G) the operator is the stencil truncated, symmetric and nonsingular, and no transform diagonalises it: rank 0 solves it on the gathered array by conjugate gradients (tolerance 1e-12 on the max-norm residual, iterations), 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 2^z sub-problems are FFT solves by sub-solvers of the solver's own type (cuFFT on FNL). Spectrally equivalent: 1, 17, 24, 30 iterations at p = 2, 4, 6, 8, independent of the grid and of the companion axes' kinds (wide_dg_zero_ghost.py, wide_dg_e2b.py). The ghost layers returned carry the faces' rules. A zero-ghost face against another kind on one axis is fatal (a reflection would couple the parity classes).

  • 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|), which rank 0 sums over the gathered source (free_ghosts: the values do not depend on the decomposition; PRISM's legacy sum ran over each rank's blocks, and its t0 fields at 1 and 2 ranks differed by up to 3.5 relative on rmf-open-t0). The operator is the zero-ghost one, the two kinds mixing on an axis; the ghosts move to the right-hand side, D(G)_0 u = factor f - D(G) u_g, and the returned ghost layers carry g. Edges and corners of the ghost region stay zero (D(G) reads along the axis lines only). Wide D(G) only (fatal with the compact operator; the multigrid refuses it). The PCG and the sub-lattice split skip a null axis.

All are symmetric linear maps of the source, which is what exact energy and momentum conservation of self-gravity needs (FLUME models).

ProcedureRole
initialize(grid, field, isolated, bc, dg_order)Check the scope, size the solver, make the plans (and, isolated, the kernel's spectrum) on rank 0; bc(2,3) optional, consistent with the grid's periodicity; dg_order the wide D(G) (E2a)
solve(field, q, iv, factor, phi, energy)Gather q(iv) to rank 0, transform, scatter phi(1-ngc:n+ngc, ..., b) with every ghost layer filled from the global solution; optionally energy = sum f phi
solve_packed(field, buf, factor, phi, energy)The same from the interior source already packed (block-major, i fastest: what a device kernel produces)
initialize_blocks(nblk, nijk, ngc, dxyz, level, null_xyz, isolated, dirichlet, bc, redundant)Size the solver for an explicit uniform block lattice, no grid object and no scope check (the caller owns it): the multigrid's root; redundant (#84 E5c) every rank holds the arrays and plans, and solve_blocks gathers the source to every rank and solves there, no scatter (compact operator, not isolated, no energy); dirichlet (#87) zero at the first ghost centres, by the odd extension to 2 (n + 1) cells per active direction (the type-I sine transform: FLAIL_BC_DIRICHLET_GHOST on every axis); bc the kinds per face (E1a)
solve_blocks(nbl, coords, buf, factor, gw, out, energy)Solve on an explicit list of this rank's blocks (integer coordinates, packed source), returning each block's potential with gw ghost layers; initialize/solve_packed are wrappers over these two
destroy()Free the arrays and the plans
  • Data path, the distributed slab transform (#84 E5d): each rank owns z-planes of the transform array (the grid with its extensions or padding); block planes go to the owners of the planes that map to them, the planes are transformed in 2-D, transposed to z-lines (MPI_ALLTOALLV) and transformed in 1-D, divided by the eigenvalues (or multiplied by the kernel's spectrum), and back; every block fetches its cells with their ghosts from the owners, so no ghost exchange follows; a rank's own planes and cells are copied directly, outside MPI. The stages are backend hooks (slab_update, slab_planes, slab_transpose, slab_lines, slab_spectral, slab_keep_kernel): the FNL twin keeps the slab arrays on the device. One plan per plane and per row of lines, each executed through a fixed buffer: the potential is bitwise the same on every rank count (poisson.sh leg slab: 26 of 26 potentials, CPU np 1-3, FNL np 1, 2, 4). At most m(3) ranks hold planes. The CG solves, the multigrid's redundant root and the CG's sub-problems keep the rank-0 path (MPI_Gatherv to rank 0, solve, MPI_Scatterv); global_potential returns the transform array on rank 0 in either path. Measured (medooza, 2.1 M uniform cells, per solve, E5c -> E5d): CPU periodic 0.138 -> 0.115 s on 1 rank, 0.124 -> 0.026 s on 16; isolated 0.734 -> 0.729 s and 0.733 -> 0.108 s; FNL periodic 0.108 -> 0.067 s on 1 GPU, 0.080 -> 0.030 s on 4; isolated 0.149 -> 0.138 s and 0.115 -> 0.112 s.
  • Scope (fatal otherwise): every block at one level and the blocks covering the domain. Without bc, every active direction periodic; one periodic face against a non-periodic one, or a ghost-Dirichlet face against a face kind, on one axis are fatal. Isolated: 3-D, every direction non-periodic, cubic cells (the kernel is the unit lattice's). AMR grids: adam_poisson_mg_object (below).
  • Build: -D_FFTW and -lfftw3 (FFTW 3, FFTW_ESTIMATE plans: deterministic). Without an FFT library, initialize stops with "built without an FFT library". The FNL twin is adam_fnl_poisson_fft_object (cuFFT). Every mode builds it since the FLAIL facade (#84, E0) is a realm component; FLUME's *-fft modes and every PRISM mode (#84, E3b) link an FFT library.
  • Verification: SV-2 (periodic) and SV-6 (isolated), src/tests/poisson/test_poisson_fft.F90 (--isolated) driven by src/tests/flume/verification/selfgravity/poisson.sh; the face kinds by src/tests/poisson/test_flail_bc.F90 (E1a: per-axis discrete eigenfunctions of every kind mix, the ghost layer included, at round-off on 1-3 ranks; the D/D order against the continuum 1.99; refusals).

adam_poisson_mg_object — Poisson solver on AMR grids (multigrid-preconditioned CG) ​

Solves A phi = factor V f - B g on the leaf cells of a 2:1-balanced octree (quadtree with null_z), periodic (#81, P5a) or with a boundary kind per face (#87, #84 E1b), the mean of the right-hand side removed when A is singular. A is the composite second-order Laplacian, V-weighted (finite-volume form): the 7-point stencil between leaves of one level and, across a 2:1 face, the mimetic mean flux — every fine cell on a coarse cell's face gets the flux (2/3) w (mean of the face's fine cells - phi_C), w its face area over its spacing, and the coarse cell their sum. That coupling is the hybrid mimetic finite-difference flux of the coarse cell with its face unknowns eliminated: exact for linear phi face by face (consistent: bounded truncation, phi second order, the 2:1 face gradient first order) and symmetric (phi = K f with V K symmetric: the exact energy exchange of self-gravity). Tomida & Stone's flux matching is consistent but not symmetric; its symmetrisation is symmetric but inconsistent (O(1/h) truncation, a non-converging seam gravity): both measured in selfgravity/seam_operators.py.

ProcedureRole
initialize(grid, field, tolerance, max_iterations, coarse, isolated, bc, root_cells)Gather the leaf list (every rank), derive the level hierarchy (the tree levels, then the levels inside the blocks down to the root, at most root_cells cells, default the root solver's root_cells_max(growth), growth the cells of the root's transform per cell (2 or 4 per bounded axis): on the host a transform of 128^3 cells, so a periodic root of 128^3 and a ghost-Dirichlet one of 64^3; on cuFFT 32^3 cells; #84 E5c, E6, E6b), assemble the operator, the ghost fill and the V-cycle operators from the geometry; size the root FFT solver (on every rank when at most 32^3 cells), of the dynamic type of coarse when given (FLUME's FNL backend passes its cuFFT solver), FFTW otherwise; isolated also sizes the ghost-to-ghost Green's function of the James step (an isolated FFT solver of the same type on the coarsest box widened by the ghost layer); bc the kind of every face (E1b; default ghost-Dirichlet on the bounded axes)
solve(field, q, iv, factor, phi, energy, ghost)FCG on A, preconditioned by one V-cycle, from the previous solution when the grid has not changed; phi returned with every ghost filled; optionally energy = sum V f phi; ghost the boundary value of each boundary face, of its kind (below); isolated: the James step (no ghost)
solve_packed(field, buf, factor, phi, energy, ghost)The same from the packed interior source
symmetry_defect(x, y), apply_a(x, y)Verification: `
apply_b(g, y), apply_bt(x, s)y = B g (the boundary values onto the boundary rows) and s = B^T x; the boundary faces are nbnd, bnd_row, bnd_w (the weight in B), bnd_kind, bnd_dir (the outward axis, signed), bnd_x (where the value lives: the ghost or the face centre)
destroy()Free everything

  • Solver: flexible CG (Notay's FCG(1)), stopping on the relative max-norm residual: tolerance (default 1e-10) or, when the source's mean dominates it, its round-off floor 100 eps max|factor f| / max|factor (f - mean f)| (a 1e-6 perturbation of a uniform density is known to 1e-10 of itself; without the floor every solve ran to the cap, #81 P5b); tolerance_used keeps the one applied, and reaching max_iterations (200) above it prints a warning. The V-cycle runs over levels holding the leaves and a shadow block per refined node (ADAM keeps leaves only); red-black Gauss-Seidel (2 sweeps down, 2 up), restriction by sum, piecewise-constant prolongation. Below the coarsest tree level (the minimum leaf level: the whole domain) the levels coarsen inside the blocks (#84 E5c): the same blocks on the same ranks, half the cells per block per level along the active directions, while those are even and the level has more than root_cells cells; the last, the root (lbot), is solved by adam_poisson_fft_object%solve_blocks, with redundant on every rank (MPI_ALLGATHERV, no scatter) when it has at most 32^3 cells (the default), on the FFT solver's distributed path above that (a larger root_cells). A coarsest tree level of at most that many cells is its own root, as before E5c. The design is Athena++'s (Tomida & Stone 2023, §3: each MeshBlock coarsened to one cell, the root solved redundantly after an MPI_Allgather); its cost is FCG iterations (test_poisson_mg --sub: 1 -> 14 on a uniform 1:2:3-cell box, point smoothing on anisotropic cells; +1 to +13 on AMR), its gain the rank-0 gather of the whole base grid per V-cycle. Diagnostics after a solve: iterations, residual, factor_mean.

  • Ghosts: the first layer per face from the solution (across a 2:1 face, the value whose one-sided difference is the operator's own face gradient), deeper layers copying the first.

  • Distributed, assembled: CSR rows over global cell keys, each operator with a halo-exchange plan (MPI_ALLTOALLV); every rank owns its leaves' rows and the shadows whose first leaf it owns. Correct for any rank count (results equal to round-off across rank counts); not tuned for speed.

  • Backend hooks (#84 E6a): the FCG iteration and the V-cycle work on vectors allocated once by initialize (wf..wqp, the level vectors lr, le, lt1, lt2, the maps mlev/mleaf of the leaves into the level spaces, the root buffers) through overridable hooks: apply (halo exchange and CSR rows), smooth, vec_fill, vec_copy, vec_lincomb, vec_axpy, vec_dot, vec_relmax, vec_permute, vec_update (a vector to or from where the hooks compute) and work_place (the operators and vectors put there after initialize, released by destroy). On the host they are the loops of before, bitwise. The right-hand side, the root (through solve_blocks), the James step and the ghost fill's placement stay on the host. apply_a and symmetry_defect stage host vectors through wp/wq. The FNL twin adam_fnl_poisson_mg_object keeps it all on the device (fnl); the facade takes it from flail%mg_mold.

  • Boundary faces (#87 P1; per face, #84 E1b): an active axis is periodic or bounded; each face of a bounded axis has a kind (adam_flail_parameters), its row homogeneous (symmetric) and its value g on the right-hand side (A phi = factor V f - B g), w the face area over the spacing:

    KindgRowBGhost
    FLAIL_BC_DIRICHLET_GHOST (#87, the default)value at the ghost centre-w phi_Cwg
    FLAIL_BC_DIRICHLETvalue at the face-2 w phi_C2 w2 g - phi_C
    FLAIL_BC_NEUMANNoutward derivative at the facenonew hphi_C + h g

    Without a Dirichlet face A is singular as on a periodic grid: the mean of factor V f - B g and of the solution are removed. The V-cycle's levels take the same kinds, homogeneous, and the root the FFT solver with them (E1a): on a uniform grid that is its own root one V-cycle is the exact solve.

  • Isolated boundaries (#87, P2), the James step: phi0 = A^-1 factor V f (zero ghosts), the screening charge s = B^T phi0 on the ghost cells, their free-space values g = -G s (G the lattice Green's function over the coarsest spacing, from ghost cell to ghost cell, by the padded FFT on the coarsest box widened by the ghost layer, gathered on rank 0), then phi = A^-1 (factor V f - B g): the map (A^-1 + A^-1 B G B^T A^-1) factor V is symmetric in the V-weighted form, the exact discrete free-space solution on a uniform grid, second order on AMR. Two solves per call, each warm-started from its own last solution (iterations0, iterations).

  • Scope (fatal otherwise): without bc every active direction periodic or none; with it, the grid's periodicity and bc agree, the two faces of an axis both ghost-Dirichlet or neither (the FFT solver's scope); isolated: 3-D, no periodic direction, no bc, cubic cells, every leaf on the domain boundary at the coarsest level; 2:1 balance across faces; fewer than 2^19 cells per direction at the finest level. Host code: FLUME's FNL backend runs it on the host with a cuFFT coarsest level (#81 P5b).

  • Build: as adam_poisson_fft_object (-D_FFTW, not re-exported by adam_common_library; mode test-poisson-mg-gnu).

  • Verification: SV-9, SV-12 (Dirichlet), SV-13 (isolated) and the boundary kinds (E1b, --bc: the FFT solver's result on uniform grids, a manufactured potential on AMR), src/tests/poisson/test_poisson_mg.F90, driven by src/tests/flume/verification/selfgravity/poisson.sh (legs sv9, sv13, bcmg).

adam_lattice_green_library — lattice Green's function of the 7-point Laplacian (generated) ​

lattice_green(n) returns G(n) on the unit cubic lattice, L_1 G = delta_(n,0), G -> -1/(4 pi |n|): the kernel of the isolated FFT solver. Generated, do not edit: scripts/lattice_green.py writes it and validates it (--check compares the committed file byte for byte). For max |n_d| <= 32 it is a table (the wedge n1 >= n2 >= n3 >= 0, 6545 values): the Dirichlet problem on the cube [-32, 32]^3, solved exactly by a type-I sine transform, with the asymptotic expansion on the faces. Beyond, the expansion itself: G0 + G1 + G2 + G3, G0 = -1/(4 pi r), G1 = -5 (n1^4 + n2^4 + n3^4 - 3 r^4/5)/(32 pi r^7), G2 = P8/r^13, G3 = P12/r^19, each order of the Taylor series of L_1 solved symbolically; the two coefficients the series leaves free (one cubic harmonic of degree 4 in G2, one of degree 6 in G3) are fitted to the table. Checked by the generator: the table's residual 2.2e-16; the residual of the composite kernel on [-64, 64]^3, across the faces where table and expansion meet, 1.0e-13; G(0) against Watson's closed form for the simple cubic lattice, 2.0e-13 relative.

Infrastructure ​

adam_mpih_object — MPI wrapper ​

Thin convenience layer over MPI. Provides rank/size query, tic-toc performance timing, rank-prefixed console output, barrier, and abort/error_stop helpers. The ghost exchange itself lives in adam_field_object (update_ghost_mpi), with its request flags (req_send_recv).

adam_io_object — parallel I/O ​

Holds the [IO] configuration (output basenames and cadences, auxiliary-field switches, the seam div(B) monitor keys), writes the history files (residuals, energy, divergence, charge conservation), and provides the XH5F field writer save_field (rank 4 and 5, kinds R8P, R4P, I8P, I4P, I2P, I1P). Field output is one HDF5 file per rank, <basename>-procNNNNNN.h5, indexed by one <basename>.xdmf (realm_object%open_file_xh5f, save_q_xh5f). Restart files are written by adam_object.

adam_slices_object — 2D diagnostic slices ​

Extracts planar 2D slices from the 3D field at configurable spatial bounds and resolution (emin, emax, nijk). Save frequency is controlled by n_save. Output: raw binary written with MPI-IO (adam_object%save_slice, extension .mat), converted by ASCOT or scripts/slices_converter.py. Configuration section: [slices].

adam_adam_object — top-level orchestrator ​

Composer that binds all SDK objects into a single adam_object and drives the simulation lifecycle:

  • Initialisation: initialize — sets up grid, tree, maps, field, and I/O from INI file
  • AMR cycle: adapt, amr_update, refine_uniform — marking, refinement/coarsening, load rebalancing
  • Communication: make_comm_local_maps_ghost_bc, mpi_redistribute, mpi_gather_refinement_needed, blocks_reorder
  • I/O: save_restart_files, load_restart_files (Fortran stream files: blocks per rank -procNNNNNN.fbd, tree nodes .tnd, time .time), save_slice, save_vtk
  • Geometry: interpolate_at_point, prune

adam_parameters — global constants ​

Named integer constants shared across all modules: boundary condition type flags, AMR marker values (TO_BE_REFINED, TO_BE_DEREFINED, TO_NOT_TOUCH), face enumeration codes (FEC_1_6_ARRAY), and coordinate-direction mappings (FEC_TO_DELTA, DELTA_TO_FEC).

adam_common_library — aggregate entry point ​

Re-exports every module listed above, including adam_mpih_global. Applications need only:

fortran
use adam_common_library

adam_globals re-exports only mpih, without the type-definition modules.

Program-scope singletons ​

The only program-scope singletons are the MPI handlers, one per process: mpih (adam_mpih_global, mpih_object) and, on the FNL backend, mpih_fnl (adam_fnl_mpih_global), which also owns the device context. Every other object is a component of a realm (self%adam%grid, self%adam%field, self%adam%maps, self%rk, self%weno, self%ib, …) and reaches library routines as an argument:

fortran
call self%adam%field%update_ghost_local(q=q, grid=self%adam%grid, maps=self%adam%maps)

The former adam_{grid,field,maps,tree,weno,ib,rk}_global modules no longer exist; do not add new *_global modules.

Forest orchestration (multi-realm) ​

The forest machinery coordinates more than one independent simulation domain — a realm — through a synchronized time loop. The user-facing conceptual overview, the manifest schema, and the α/β cadence trade-offs are documented in the forest concepts page; this section is the library-developer reference for the contract surface.

adam_realm_object — the per-realm contract ​

realm_object is the (concrete) parent every solver application extends. It owns every library object as a component (adam with grid, tree, field, maps; io, amr, slices, weno, ib, rk, flail, …) and declares the _forest-suffixed TBP family the forest dispatches against. The defaults of the core and staged TBPs (apply_reflux_to_stage_forest included) and of the seam pack/unpack/fill TBPs error_stop with the TBP's name; finalize_mpi_forest (mpih%finalize), after_topology_build_forest (nothing), coupling_descriptor_forest (sentinel), runtime_amr_forest (false) and regrid_forest (never) have usable defaults.

TBPRole
initialize_forest(filename, realms_number, memory_avail, nv, verbose)Per-realm prologue: load INI, allocate state, place initial conditions, open IO.
finalize_forest()Per-realm epilogue: close IO, release resources.
finalize_mpi_forest()Process-global MPI finalize, called once by the forest after every realm's finalize_forest (FNL overrides it).
compute_local_dt_forest(dt_local)Report this realm's CFL-limited dt; the forest min-reduces across realms.
advance_one_step_forest(dt)Fast path (one realm, no 2:1 face): realm owns the entire step internally.
open_step_forest(dt)Staged path: per-step prologue (external fields, RK stage buffer init, time bookkeeping).
begin_stage_forest(k, K_total, dt, realm)Staged path: open integrator stage k; sets self%stage_active = k.
end_stage_forest(k, K_total, dt, realm, flux_register)Staged path: residuals + stage assignment for stage k.
close_step_forest(dt)Staged path: per-step epilogue (state assembly, BC, div-clean, time advance); clears self%stage_active = 0.
fill_seam_from_peer_forest(peer, p_idx)Receive-side roundtrip: copy peer's interior into self's seam ghosts. Source/target buffer chosen by self%stage_active/peer%stage_active (>0 ⇒ rk%q_rk(:,...,stage), else committed q).
pack_seam_cells_forest(p_idx, buf), unpack_seam_cells_forest(p_idx, buf)Pack own cells of the cross-rank seam send rows; unpack the receive rows into own ghosts.
after_topology_build_forest()Backend hook fired once after the seam maps are built — FNL uses it to refresh the device-resident seam maps (field_fnl%maps%copy_cpu_gpu).
apply_reflux_to_stage_forest(stage, dt, flux_register)Apply Berger-Colella reflux at the realm's end-of-step (α.r1 gate).
stages_per_step_forest()Report K (the realm's integrator stride; rk%nrk for SSP-RK).
coupling_descriptor_forest(scheme_time, rk_scheme, nv)Report (scheme_time, rk_scheme, nv) for β admissibility (issue #18). Default returns ''/''/-1 — a realm can't join β unless it overrides.
is_done_forest(done)Report whether the simulation has reached its termination predicate.
post_step_forest(dt, t, it, do_save_state, do_save_residuals, do_save_restart, do_amr, realm)Per-step diagnostics / IO block (savers, residual prints).
runtime_amr_forest(), regrid_forest(regridded)Whether the realm regrids during the run, and the regrid itself, called by forest%regrid after post_step_forest (issue #74).

Each _forest TBP is pass(self)-dispatched. The forest receives type(realm_pointer), intent(inout) :: realm(:) as a parameter and never owns realm state: the driver allocates the realms as a concrete app type (e.g. type(prism_cpu_object), allocatable, target :: realm(:)) and points one realm_pointer at each, so element is is realm(is)%ptr. The forest never takes a polymorphic array class(realm_object) :: realm(:): nvfortran (25.11, 26.9) can address its elements with the declared type's stride and dispatch on the wrong realm (#96).

The fast path bypasses every per-stage TBP, but it is taken only by a single realm without 2:1 AMR faces (flux_register%nfaces == 0, evolve_one_step); a single realm with a 2:1 face runs the staged path. An app that refines its grid must therefore implement the staged TBPs; only the seam TBPs can stay at their defaults in a single-realm app.

adam_forest_object — the orchestrator ​

forest_object owns no realm state: its components are the realm count n, the flux register (flux_register, the coarse-fine faces of every realm and seam) and the timing report. All forest operations are TBP dispatches against the realm array passed in.

Public methods:

MethodPurpose
initialize(realm, filename)Shared-INI path: every realm reads the same INI. Used by legacy single-INI drivers (e.g. PRISM's no-manifest form).
initialize_from_manifest(realm, manifest)Manifest path: each realm reads its own INI from the parsed forest_manifest_t; topology is wired afterwards.
simulate(realm, filename)Top-level entry: initialize → time loop → finalize.
simulate_from_manifest(realm, manifest)Top-level entry: manifest variant.
compute_global_dt(realm, dt)Min-reduce per-realm CFL across realms (intra-process) and ranks (MPI_ALLREDUCE).
evolve_one_step(realm)One global timestep: Phase 0 (α seam fill, open) → K-loop (Phases 1–3) → Phase 4; per-seam cadence gating inside.
post_step(realm)Per-step diagnostics block across realms.
regrid(realm)Dispatch each realm's regrid_forest; rebuild the flux register when a tree changed (issue #74).
is_done(realm, done)AND-reduce per-realm termination across realms and ranks.
finalize(realm)Sequence each realm's finalize_forest.

Private helpers:

HelperPurpose
populate_inter_realm_topology(realm, manifest)Translate manifest face-pairs into per-realm maps%inter_realm_neighbors and the derived seam_local_* arrays.
check_beta_admissibility(realm, manifest)Enforce β admissibility on every stage_coincident seam (PRD #18 M2).
apply_reflux_corrections(realm, stage, dt)Dispatch each realm's apply_reflux_to_stage_forest (realm-side body gates on end-of-step under α.r1).

Manifest data flow ​

forest.ini                          (text manifest)
  │
  ▼  adam_forest_manifest.F90 ─ read_forest_manifest
forest_manifest_t                   (transient struct: realm_ini paths + face_pairs)
  │
  ▼  forest_object%populate_inter_realm_topology
realm(is)%adam%maps%inter_realm_neighbors(:)              (per-realm symmetric topology)
  │
  ▼  build_inter_realm_ghost_cell_map + build_seam_local_map
realm(is)%adam%maps%seam_local_map_ghost_cell(:,:)        (sorted per-cell ghost map)
realm(is)%adam%maps%seam_local_peer_realm(:)              (one entry per distinct peer)
realm(is)%adam%maps%seam_local_peer_row_start(:)
realm(is)%adam%maps%seam_local_peer_row_count(:)
realm(is)%adam%maps%seam_local_cadence(:)                 (CADENCE_END_OF_STEP | CADENCE_STAGE_COINCIDENT)
realm(is)%adam%maps%seam_local_send_buf / recv_buf        (per-peer pack/unpack buffers, allocated once)

The orchestrator's per-step Phase 2 (β) and Phase 0 (α) loops consume the seam_local_* arrays directly; no manifest reads happen at step time.

Per-seam coupling cadence ​

Two integer parameters on adam_maps_object tag each peer slot:

  • CADENCE_END_OF_STEP = 0 — α default: seam fill once per global step, at its top before open_step_forest (Phase 0); FLUME's open_step_forest copies the ghosts into every stage buffer (seam_ghosts_to_stages, #113).
  • CADENCE_STAGE_COINCIDENT = 1 — β opt-in: seam fill at every RK substage, before end_stage_forest (Phase 2 inside the K loop).

The orchestrator gates per-seam: Phase 2 fires only on STAGE_COINCIDENT seams (once per substage); Phase 0 fires only on END_OF_STEP seams (once per step).

β admissibility — same scheme_time, rk_scheme, nv, and K between both endpoint realms — is verified at initialize_from_manifest time via check_beta_admissibility. Mismatch → immediate error_stop; no silent downgrade.

Intra-realm AMR seams (a second, independent seam family) ​

The inter-realm machinery above couples two realms. A separate family — the intra-realm AMR seam — handles a 2:1 resolution jump inside a single realm's Morton octree (a level-ℓ block abutting a level-ℓ+1 block, produced by AMR markers). It is not manifest-declared and uses different code:

  • Registration: forest_object%register_intra_realm_amr_seams(realm) (private) does a two-pass registration of the coarse-fine faces into the flux register with SEAM_KIND_INTRA_REALM_AMR (adam_flux_register_object.F90), and allocates/fills maps%amr_seam_quadrant. The array is allocated for every realm with blocks, so the seam div(B) guard-rail (below) arms on its tolerance alone.
  • Ghost fill: the fine block's coarse-overlapping ghosts are filled through field_object%update_ghost_local's flag-4 map-row path → interp_seam_ghost. The interpolation regime is maps%seam_ghost_fill, one of SEAM_FILL_INJECTION (0, copy), SEAM_FILL_COMPATIBLE (1, q=2 restriction-compatible), SEAM_FILL_TRICUBIC (2, q=4) — declared in adam_seam_interpolation_library.F90, selected by [amr] seam_ghost_fill. Effective default is tricubic: adam_object%initialize overwrites the maps struct default (SEAM_FILL_INJECTION) with SEAM_FILL_TRICUBIC when the key is absent. The FNL device path carries the regime via update_ghost_local_gpu(..., seam_ghost_fill=..., ...). Seam fill guard (#101): when an app sets maps%seam_density_var to the index of its density, with the three momenta at the next indices (FLUME does; the default 0 turns it off, as in PRISM), the fill is clamped in the manner of Suresh and Huynh by seam_mp_bound (adam_seam_interpolation_library): the density to the range of the coarse cells that bracket the ghost, widened by the data's curvature, and the velocity (momentum over that density) to the same bound of the coarse cells' velocities; a ghost whose clamped density is below SEAM_DENSITY_GUARD = 0.1 times its anchor's (the coarse cell holding it) takes, in every variable, the value of that anchor (the injection value). Each variable recomputes the density of its ghost, so the decision is the same in every row and in any order: on the CPU, in the parallel FNL kernels, in the MPI pack and in the inter-realm seam pack (pack_seam_rows). The guard belongs to the state: update_ghost_local / update_ghost_mpi (and their FNL twins) take seam_guard = .false. for any other field filled with the same maps (FLUME's positivity-limiter factors, whose zero is legitimate). Cost: the ghost's part (seam_guard_prepare) is computed once per ghost on the CPU (once per row on FNL), and both clamps first test the range of the bracketing cells, inside the bound, so the full bound is computed only when a value leaves that range; on shock-cylinder-ib, the regression case with the most seam ghosts, the run is 20 % slower than without the guard (the others within noise).
  • Reflux: Berger-Colella reflux at the 2:1 face, dispatched by apply_reflux_corrections at each realm's own end-of-step (α.r1 gate), identical machinery to the inter-realm register.

The two families are structurally disjoint — a single-realm run never touches fill_seam_from_peer_forest; a same-resolution forest never touches interp_seam_ghost. They can coexist (a forest of internally-AMR-refined realms).

Seam div(B) guard-rail ([IO], issue #29) ​

The 2:1 AMR seam breaks the mimetic div_h(curl_h) ≡ 0 identity (coarse Δx vs fine Δx/2 stencils no longer commute) and injects an O(h^p) div(B) source that is refinement-convergent but unbounded in t at fixed h — a constraint violation, not an energy instability. No local collocated-FD fix exists (five classes ruled out); the resolution is accept-truncation. To make the otherwise-silent growth visible, save_divergence_history (adam_prism_common_object.F90) carries an opt-in monitor: [IO] seam_divB_tol (default -1.0, off) and seam_divB_error (default .false., warn-only). It arms when seam_divB_tol > 0: its second condition, allocated(maps%amr_seam_quadrant), holds for every realm with blocks, so it covers intra-realm 2:1 faces and inter-realm refined seams alike (without 2:1 seams div(B) stays at round-off and it stays silent); on exceedance it error_stops (seam_divB_error = .true.) or warns. seam_divB_tol/seam_divB_error live on the base adam_io_object.

Load-bearing invariants ​

  • Phase 2 → Phase 3 ordering (β): Phase 2 must complete on ALL realms before Phase 3 starts on ANY realm; otherwise the read-after-overwrite race between fill_seam_from_peer_forest writes and compute_residuals reads returns. The serial inner loops within a rank give this for free under the Phase-A replicated-forest layout.
  • Inter-realm seam diagnostic refill (issue #31): any div(B) diagnostic that runs self%update_ghost before reading the divergence stencil must ALSO refill the inter-realm seam (fill_seam_from_peer_forest per peer), because update_ghost fills intra-realm ghosts + physical BCs but NOT the inter-realm seam (the forest owns it). Omitting it reports a spurious seam-skin div(B) on a div-free field. present(realm)-guarded so single-realm runs are untouched. A 1:1 same-resolution inter-realm mirror seam must use stage_coincident (β) to stay div-free — α leaves the seam ghosts unfilled during RK substages.
  • stage_active discipline: begin_stage_forest(k) MUST set self%stage_active = k; close_step_forest(dt) MUST clear it back to 0. The fill_seam_from_peer_forest buffer-selection logic reads self%stage_active and peer%stage_active to pick between q (committed) and rk%q_rk(:,...,stage_active). Mid-step writes to q outside this discipline corrupt peer reads.
  • α.r1 reflux gate: apply_reflux_to_stage_forest and accumulate_seam_fluxes_fv MUST gate on stage == rk%nrk (the realm's own final substage). flux_register's third axis is collapsed to 1 under α.r1; per-stage RK-weighted reflux (Wang 2018) is deferred. Both α and β share this gate — β does not undo α.r1.
  • Fast path: when size(realm) == 1 and the flux register has no faces, evolve_one_step routes through advance_one_step_forest directly with zero exposure to the per-stage TBPs. A single realm with a 2:1 AMR face takes the staged path, so an app that refines must implement the staged TBPs.
  • No library singletons: apart from mpih (and mpih_fnl), every library object is a realm component, passed as an argument. See "Program-scope singletons" above.

Tracked through ​

PRDScopeStatus
#10Phase D inception — forest orchestrator designshipped
#13Coarse-fine interface machinery (A/B/C)A + C shipped; B characterized (→ #29)
#16α end-of-step barrier + asymmetric Kshipped
#18β stage-coincident recovery + cross-config oracleshipped
#21tricubic coarse→fine seam ghost fill (seam_ghost_fill)shipped
#22/#23FNL seam fill + FV residual/reflux device twinsshipped
#28np>1 reflux register ownership-gatingshipped
#292:1 seam div(B) — accept-truncation + guard-railshipped (resolution)
#31inter-realm 1:1 seam div-free (β + diagnostic refill); FNL fWLayer HtoDshipped

Module summary ​

ModuleCategoryFile
adam_parametersInfrastructureadam_parameters.f90
adam_tree_node_objectAMRadam_tree_node_object.f90
adam_tree_bucket_objectAMRadam_tree_bucket_object.f90
adam_tree_objectAMRadam_tree_object.F90
adam_grid_objectGrid/Fieldadam_grid_object.F90
adam_field_objectGrid/Fieldadam_field_object.F90
adam_maps_objectCommunicationadam_maps_object.F90
adam_weno_objectNumerics — spatialadam_weno_object.F90
adam_fdv_operators_libraryNumerics — spatialadam_fdv_operators_library.F90
adam_riemann_euler_libraryNumerics — spatialadam_riemann_euler_library.F90
adam_rk_objectNumerics — temporaladam_rk_object.F90
adam_leapfrog_objectNumerics — temporaladam_leapfrog_object.F90
adam_blanes_moan_objectNumerics — temporaladam_blanes_moan_object.F90
adam_cfm_objectNumerics — temporaladam_cfm_object.F90
adam_eos_ic_objectPhysicsadam_eos_ic_object.F90
adam_ib_objectPhysicsadam_ib_object.F90
adam_amr_objectPhysicsadam_amr_object.F90
adam_flail_objectPhysicsadam_flail_object.F90
adam_flail_parametersPhysicsadam_flail_parameters.F90
adam_poisson_fft_objectPhysicsadam_poisson_fft_object.F90
adam_poisson_mg_objectPhysicsadam_poisson_mg_object.F90
adam_lattice_green_libraryPhysicsadam_lattice_green_library.F90
adam_mpih_objectInfrastructureadam_mpih_object.F90
adam_io_objectInfrastructureadam_io_object.F90
adam_slices_objectInfrastructureadam_slices_object.F90
adam_adam_objectInfrastructureadam_adam_object.F90
adam_realm_objectForestadam_realm_object.F90
adam_forest_objectForestadam_forest_object.F90
adam_forest_manifestForestadam_forest_manifest.F90
adam_flux_register_objectForestadam_flux_register_object.F90
adam_common_libraryEntry pointadam_common_library.F90
adam_mpih_globalSingletonadam_mpih_global.F90
adam_globalsSingleton aggregatoradam_globals.F90