Appearance
The FNL library (GPU)
The FNL backend is ADAM's GPU acceleration layer, built on top of the FUNDAL GPU memory management library. The same sources compile as OpenACC with nvfortran (*-fnl-nvf modes) or as OpenMP target offloading with amdflang (*-fnl-omp-amd modes, DEV_OMP): every kernel carries an !$acc directive and its !$omp twin. It follows a consistent two-tier pattern for every subsystem: a wrapper object that holds the device-resident arrays, paired with a kernel module that contains the device subroutines implementing the actual computation. The wrappers are separate types, not extensions of the common classes (the Poisson solvers excepted): each is a component of the application's FNL type, initialised from the realm's common objects passed as arguments.
No physics or algorithmic logic is duplicated — all equations, coefficients, and data structures are defined once in src/lib/common and mirrored to the GPU by the FNL layer.
The aggregate entry point adam_fnl_library re-exports the entire FNL API together with adam_common_library, including adam_fnl_mpih_global; a single use adam_fnl_library statement in application code exposes both layers.
Contents
- Design and memory model
- Program-scope singletons
- Field
- WENO reconstruction
- Runge-Kutta integration
- Immersed boundary method
- Communication maps
- Finite difference/volume operators
- Poisson solver
- MPI handler
- Module summary
Program-scope singletons
The only FNL singleton is mpih_fnl (adam_fnl_mpih_global, type mpih_fnl_object): the GPU-aware MPI handler, which also owns the device context and must be initialised exactly once per process. The other FNL objects are components of the application's FNL type, one set per realm, initialised after the common initialisation from the realm's own objects (PRISM, adam_prism_fnl_object.F90):
fortran
if (.not.mpih_fnl_is_initialized) then ! once per process
call mpih_fnl%initialize(do_mpi_init=..., do_device_init=.true., verbose=.true.)
mpih_fnl_is_initialized = .true.
endif
call self%prism_common_object%initialize(filename=filename, memory_avail=memory_avail_, verbose=.true.)
call self%field_fnl%initialize(grid=self%adam%grid, field=self%adam%field, maps=self%adam%maps, verbose=.true.)
call self%ib_fnl%initialize(grid=self%adam%grid, field=self%adam%field, ib=self%ib)
call self%rk_fnl%initialize(grid=self%adam%grid, field=self%adam%field, rk=self%rk)
call self%weno_fnl%initialize(weno=self%weno)App-specific FNL objects follow the same pattern (PRISM's coil_fnl, pic_fnl, pml_fnl). The former adam_fnl_{field,ib,rk,weno}_global modules no longer exist.
Design and memory model
CPU vs GPU array layout
The common library stores field data in Fortran column-major order with variables as the fastest-varying index:
CPU: q_cpu(nv, ni, nj, nk, nb) — stride-1 on nvOn the GPU the block index is moved to the front so that threads mapped to adjacent blocks access adjacent memory locations, improving warp-level coalescing:
GPU: q_gpu(nb, ni, nj, nk, nv) — blocks contiguousThe device state arrays (q_gpu, dq_gpu, …) belong to the application's FNL type (e.g. flume_fnl_object%q_gpu); field_fnl%initialize(..., q_gpu=) can allocate one in this layout. The application transposes between layouts with FUNDAL's buffered transposed copies (dev_memcpy_to_device(..., ij=[1,5], buf=), with a staging buffer of the destination's extent, issue #63) or with the device kernels pack_blocks_transposed_dev / unpack_blocks_transposed_dev of adam_fnl_field_kernels; the copy_transpose_cpu_gpu / copy_transpose_gpu_cpu routines in adam_fnl_field_object are not bound as TBPs. This copy happens only at I/O and AMR update boundaries; during time integration the GPU layout is used exclusively.
FUNDAL integration
GPU memory allocation and host-device transfers route through FUNDAL utilities:
| Utility | Purpose |
|---|---|
dev_alloc | Allocate device array |
dev_assign_to_device | Copy host array to device |
dev_memcpy_to_device | Raw device-to-device copy |
DEVICEVAR(array) macro | Mark array as already device-resident (suppresses implicit copy) |
The DEVICEVAR macro is defined in fundal.H and appears at the top of every kernel file. It is the primary mechanism by which the compiler is told that a given array lives on the device, preventing spurious host-copy insertion by the OpenACC runtime.
IB solid masking
Every kernel that modifies field variables carries an implicit guard: cells inside a solid body (phi_gpu(b,i,j,k,all_solids) < 0) are skipped. This masking is applied at the innermost loop level so that it incurs no branch divergence in fluid-only regions.
Field
adam_fnl_field_object — GPU field wrapper
A plain type (not an extension of field_object), initialised from the realm's grid, field and maps. Holds the device-side geometry and maps needed for field operations (its maps component is a maps_fnl_object) and provides the ghost exchange on the device.
Device arrays:
| Array | Shape | Purpose |
|---|---|---|
x_cell_gpu, y_cell_gpu, z_cell_gpu | (nb, 1-ngc:ni+ngc) (per axis) | Cell centroid coordinates |
dxyz_gpu | (nb, 3) | Block mesh spacing (dx, dy, dz) |
fec_1_6_array_gpu | (26) | Maps the face/edge/corner codes 1-26 to the boundary faces 1-6 |
Key methods:
| Method | Purpose |
|---|---|
initialize(grid, field, maps[, nv_aux][, q_gpu][, verbose]) | Allocate the device arrays via dev_alloc (and q_gpu(nb,…,nv) when passed) |
copy_cpu_gpu | Transfer the coordinate arrays and maps to device |
reduce_ghost_local_gpu, reduce_ghost_mpi_gpu | Fold deposited ghosts back onto their owner cells, locally and across ranks |
update_ghost_local_gpu | Apply intra-rank ghost-cell updates entirely on device |
update_ghost_mpi_gpu | Pack send buffer on device, perform MPI exchange, unpack on device |
compute_q_gradient(b, ivar, dx, dy, dz, q_gpu, gradient) | AMR refinement criterion: `max |
update_ghost_mpi_gpu, compute_q_gradient, copy_cpu_gpu and initialize take the realm-local CPU grid/field/maps (and the per-rank comm_map_*_ptr_ghost) as dummy arguments, not from program-scope singletons — so a multi-realm caller always exchanges with its own decomposition. Only mpih_fnl remains a singleton (MPI is genuinely program-global).
GPU-direct ghost exchange and the WSL2 rendezvous abort (issue #12)
update_ghost_mpi_gpu posts MPI_Isend/MPI_Irecv directly on the device-resident ghost buffers (send_buffer_ghost_gpu / recv_buffer_ghost_gpu) — the GPU-direct path that lets CUDA-aware MPI use GPUDirect RDMA. On a healthy InfiniBand stack this is the fast default and the only path: there is no host-staging fallback in the code.
On WSL2 with >=2 ranks this path aborts: UCX moves the (large) device ghost buffers via its rendezvous protocol, whose device-memory transports (cuda_copy / gdr) cannot get the GPU primary context through the /dev/dxg shim — SIGABRT in ucp_proto_rndv_send_start at the first exchange. The fault is in the UCX rendezvous path on a broken WSL stack, not in ADAM or the request handles.
The WSL workaround is export UCX_RNDV_THRESH=inf, set by src/tests/prism/regression/run-fnl-local.sh. It forces every message eager, so rendezvous is never entered and no device pointer reaches the broken transport. It is a blunt, WSL-only crutch — a pure performance regression on real IB + GPUDirect RDMA, where rendezvous is the fast path — so it is confined to the local run wrapper and must never reach a cluster job script or the application. Verified: rmf regression passes -np 1 and -np 2 (single realm) with the field digest rank-count-invariant within tolerance.
adam_fnl_field_kernels — field device kernels
All routines carry !$acc parallel loop independent and use DEVICEVAR on every device pointer argument.
| Kernel | Purpose |
|---|---|
compute_q_gradient_dev | Centred-difference gradient magnitude with reduction(max:) |
compute_normL2_residuals_dev | L2 norm √(Σ dq²) per variable with reduction(+:) |
copy_transpose_gpu_cpu_dev | Transpose (nb,ni,nj,nk,nv) → (nv,ni,nj,nk,nb) on device |
populate_send_buffer_ghost_gpu_dev | Pack ghost-cell values into MPI send buffer; supports 1-cell and 8-cell AMR averaging |
receive_recv_buffer_ghost_gpu_dev | Unpack MPI receive buffer into ghost cells |
update_ghost_local_gpu_dev | Apply intra-rank block-to-block ghost updates; supports AMR coarse↔fine averaging |
WENO reconstruction
adam_fnl_weno_object — GPU WENO coefficient wrapper
A plain type, initialised from the realm's weno (initialize(weno)). The CPU object computes all WENO coefficients once during initialisation; the FNL object mirrors them to device memory and holds the ROR (Reduced-Order Reconstruction) tables used near solid boundaries.
Device arrays:
| Array | Shape | Purpose |
|---|---|---|
a_gpu | (2, 0:S-1, S) | Optimal WENO weights per sub-stencil and face |
p_gpu | (2, 0:S-1, 0:S-1, S) | Polynomial reconstruction coefficients |
d_gpu | (0:S-1, 0:S-1, 0:S-1, S) | Smoothness indicator coefficients |
ror_schemes_gpu | (:) | ROR fallback scheme orders near solid walls |
ror_ivar_gpu | (:) | Variable indices checked by ROR |
cell_scheme_gpu | (nb, ni, nj, nk) | Per-cell effective reconstruction order |
adam_fnl_weno_kernels — WENO device kernels
| Kernel | Directive | Purpose |
|---|---|---|
weno_reconstruct_upwind_dev(S, a, p, d, zeps, V, VR) | !$acc routine seq | Reconstruct left (VR(1)) and right (VR(2)) interface values from stencil V |
The reconstruction follows the standard three-step algorithm:
- Compute
Spolynomial reconstructions from overlapping sub-stencils - Compute smoothness indicators from second-derivative sums
- Weight and convolve:
VR(f) = Σ_k w(f,k) · VP(f,k)
!$acc routine seq marks the procedure as callable from within a parallel region without launching a new kernel — one thread per stencil, invoked inside the outer loop over cells.
Runge-Kutta integration
adam_fnl_rk_object — GPU RK stage manager
A plain type, initialised from the realm's grid, field and rk (initialize(grid, field, rk)). Manages stage storage on device and drives the per-stage updates.
Device arrays:
| Array | Shape | Purpose |
|---|---|---|
q_rk_gpu | (nb, ni, nj, nk, nv, nrk) | Stage values (1 stage for low-storage, nrk for SSP) |
alph_gpu | (nrk, nrk) | SSP alpha coefficients |
beta_gpu | (nrk) | SSP beta coefficients |
gamm_gpu | (nrk) | SSP gamma coefficients |
Supported schemes:
| Scheme | Storage mode | Stages |
|---|---|---|
RK_1, RK_2, RK_3 | Low-storage | 1 — overwrites q_gpu in place |
RK_SSP_11, RK_SSP_22, RK_SSP_33 | Multi-stage | 1 / 2 / 3 |
RK_SSP_54 | Multi-stage | 5 |
Key methods:
| Method | Purpose |
|---|---|
initialize(grid, field, rk) | Allocate q_rk_gpu sized to scheme requirements |
initialize_stages(q_gpu) | Broadcast q_gpu into all stage slots |
assign_stage(s, q_gpu, phi_gpu) | Copy q_gpu into stage s, skipping solid cells |
compute_stage(s, dt, phi_gpu) | Accumulate stages 1…s−1 into stage s (SSP) |
compute_stage_ls(s, dt, phi_gpu, dq_gpu, q_gpu) | Low-storage update: q = ark·q_n + brk·q + dt·crk·dq |
update_q(s, dt, phi_gpu) | Final assembly: q += dt·beta(s)·q_rk(:,:,:,:,:,s) |
adam_fnl_rk_kernels — RK device kernels
All kernels carry !$acc parallel loop independent and mask solid cells via phi_gpu.
| Kernel | Purpose |
|---|---|
rk_assign_stage_dev | q_rk(:,s) ← q_gpu (fluid cells only) |
rk_initialize_stages_dev | q_rk(:,all_s) ← q_gpu |
rk_compute_stage_dev | q_rk(:,s) += dt·α(s,ss)·q_rk(:,ss) for ss = 1…s−1 |
rk_compute_stage_ls_dev | q = ark·q_n + brk·q + dt·crk·dq (low-storage) |
rk_update_q_dev | q += dt·β(s)·q_rk(:,s) for s = 1…nrk |
Immersed boundary method
adam_fnl_ib_object — GPU IB wrapper
A plain type, initialised from the realm's grid, field and ib (initialize(grid, field, ib)). Manages the signed-distance field phi_gpu on device and drives the eikonal enforcement cycle.
Device arrays:
| Array | Shape | Purpose |
|---|---|---|
phi_gpu | (nb, ni, nj, nk, n_solids+1) | Signed-distance function; last slice holds max over all solids |
q_bcs_vars_gpu | (:,:) | Boundary condition state variables per solid |
The sign convention: phi < 0 inside the solid (ghost region), phi > 0 in the fluid.
Key methods:
| Method | Purpose |
|---|---|
initialize(grid, field, ib) | Allocate phi_gpu and q_bcs_vars_gpu (dimensions from the grid and field arguments); copy BCS data to device |
evolve_eikonal(dxyz_gpu, dq_gpu, q_gpu) | Advance eikonal equation inside solid: q -= ∇φ·(q_bc − q); dxyz_gpu is passed explicitly by the caller (typically field_fnl%dxyz_gpu) |
invert_eikonal(q_gpu) | Enforce wall BC at solid surface (φ > 0): reflect momentum |
Wall BC modes applied by invert_eikonal:
BCS_VISCOUS(no-slip):(u, v, w) → (−u, −v, −w)BCS_EULER(inviscid):u → u − 2(u·n̂)n̂
adam_fnl_ib_kernels — IB device kernels
| Kernel | Purpose |
|---|---|
compute_phi_analytical_sphere_dev | φ = −(‖x − xc‖ − R) — negative inside sphere |
compute_phi_all_solids_dev | φ_all = max(φ₁, φ₂, …, φ_ns) — union of all solids |
compute_eikonal_dq_phi_dev | Gradient-weighted residual: `dq = |
evolve_eikonal_q_phi_dev | q -= dq inside solid (φ > 0) |
invert_eikonal_q_phi_dev | Momentum reflection at surface (BCS_VISCOUS or BCS_EULER) |
move_phi_dev | Level-set advection: ∂φ/∂t = −v·∇φ for moving bodies |
reduce_cell_order_phi_dev | Lower reconstruction order in cells adjacent to solid surface |
Communication maps
adam_fnl_maps_object — GPU maps wrapper
A plain type, the maps component of field_fnl_object, initialised from the realm's maps. Mirrors all communication index tables to device memory so that ghost-cell packing and unpacking happen entirely on the GPU, eliminating CPU staging for MPI buffers.
Device arrays:
| Array | Columns | Content |
|---|---|---|
local_map_ghost_cell_gpu | 9 | (b_src, b_dst, i_src, j_src, k_src, i_dst, j_dst, k_dst, mode) |
comm_map_send_ghost_cell_gpu | 7 | (b_src, i, j, k, v_offset, buf_idx, mode) |
comm_map_recv_ghost_cell_gpu | 6 | (buf_idx, b_dst, i, j, k, v_offset) |
send_buffer_ghost_gpu | — | 1D packed MPI send staging buffer |
recv_buffer_ghost_gpu | — | 1D packed MPI receive staging buffer |
local_map_bc_crown_gpu | — | Boundary condition crown ghost-cell map |
The mode column distinguishes two cases:
mode = 1— one-to-one cell correspondence (same refinement level)mode = 8— eight-cell average (fine block → coarse block at AMR interface)
Key methods:
| Method | Purpose |
|---|---|
initialize(maps) | Initialise and call copy_cpu_gpu |
copy_cpu_gpu(verbose) | Transfer all map arrays from CPU to device via dev_assign_to_device |
Finite difference/volume operators
adam_fnl_fdv_operators_library — device-callable spatial operators
Provides the same operators as adam_fdv_operators_library (common) in a form callable from within OpenACC parallel regions. All routines carry !$acc routine seq — no internal parallelism, one thread per cell.
Available operators:
| Operator | FD centred | FV centred |
|---|---|---|
Gradient ∇q | compute_gradient_fd_centered_dev | compute_gradient_fv_centered_dev |
Divergence ∇·q | compute_divergence_fd_centered_dev | compute_divergence_fv_centered_dev |
Curl ∇×q | compute_curl_fd_centered_dev | compute_curl_fv_centered_dev |
Laplacian ∇²q | compute_laplacian_fd_centered_dev | compute_laplacian_fv_centered_dev |
Each routine accepts a stencil half-width s and the local mesh spacing dxyz, allowing the accuracy order to be selected at call time without recompilation.
Poisson solver
adam_fnl_poisson_mg_object — multigrid Poisson solver on the GPU (#84 E6a)
Extends poisson_mg_object (common): the leaves, the hierarchy and every CSR operator are assembled on the host; work_place_fnl then enters the operators (composite operator, ghost fill, level operators, ghost terms, transfers), the level diagonals, colours and leaf maps, the cell volumes and the work vectors into device memory, with halo buffers sized for the largest operator, until destroy. The FCG iteration and the V-cycle run on the device through the common hooks: apply_fnl packs the halo on the device, exchanges it by MPI_ALLTOALLV on the device buffers (GPU-aware MPI; none on one rank) and computes one row per thread in the host's order; smooth_fnl runs one kernel per colour after its halo exchange. That parallel update equals the host's sequential Gauss-Seidel because no row couples to a cell of its own colour: the colour is the global cell parity and the level operators are 7-point (measured: no such entry in 5704 smoother calls over the AMR, isolated, boundary-kind and sub-block suites). The vector hooks are device kernels, the dots and the max-norm device reductions. The root goes through the host to the FFT solver: cuFFT by default (initialize_fnl passes poisson_fft_fnl_object when no coarse is given), at most 32^3 cells whatever its transform's growth (root_cells_max_fnl, #84 E6: a 128^3 cuFFT root, E5c's choice while the multigrid smoothed on the host, cost 2.4 to 6.8 times more per solve on 3.9 M AMR cells; the host sizes its root by the transform instead, E6b). So do the right-hand side, the James step and the ghost fill's placement. initialize_fnl releases and places the device copies itself, since the parent's initialize, called through the parent component, dispatches to the host's no-op hooks. Results equal the host's to round-off (another summation order in the dots, FMA contraction). FLUME's FNL backend and bench_poisson_fnl set flail%mg_mold to it; test-poisson-mg-fnl-nvf builds test_poisson_mg on it.
adam_fnl_poisson_fft_object — FFT Poisson solver on the GPU (cuFFT)
Extends poisson_fft_object (common). The distributed slab transform (#84 E5d) keeps the CPU redistribution and return of the block planes on the host; the transform runs on each rank's GPU with its arrays resident there (slab_r, slab_c, pen, khat_pen and the transpose buffers, created with the plans): slab_r goes up (slab_update_fnl), the planes (slab_planes_fnl: one D2Z/Z2D plan for one plane, executed on every plane), the transpose (slab_transpose_fnl: packed and unpacked on the device, MPI_ALLTOALLV on the device buffers, so GPU-aware MPI as the ghost exchange; none on one rank), the lines (slab_lines_fnl: one Z2Z plan for one row of z-lines, executed on every row), the spectral solve (slab_spectral_fnl), back, and slab_r comes down. The result does not depend on the rank count. The rank-0 path below remains for PRISM's CG solves and the multigrid's small root. The global source goes to the device, a cuFFT D2Z transform, a kernel divides by the eigenvalues of L_h (the host-computed tables of the CPU solver), a Z2D transform, and the potential comes back. cuFFT is bound to the stream of the synchronous OpenACC queue (acc_get_cuda_stream(acc_async_sync)), so the transforms and the kernel between them are ordered.
- Build:
-D_CUFFTand-cudalib=cufft(modetest-poisson-fft-fnl-nvf);cufftandopenaccare declared intrinsic modules infobos.d/templates.fobos. cuFFT is NVIDIA only: without_CUFFT(an AMD build) making the plans stops the run. - Not FUNDAL-agnostic, by necessity: the device arrays live in an OpenACC data region around the cuFFT calls (
host_data use_device), not in FUNDAL allocations. - Isolated (P4): the padded transforms (
m >= 2 (n + ngc)per direction, the smallest 13-smooth even size, #84 E5b); the kernel's spectrumkhatis made once on the device (a D2Z of the kernel) and stays there (enter data, released bydestroy_plans); each solve multiplies by it instead of dividing by the eigenvalues. - Verification: SV-2 asserts it against the discrete solution at round-off, SV-6 against the exact potentials of two spheres; both against the FFTW solver (
selfgravity/poisson_oracle.py compare: isolated,1.6e-15).
MPI handler
adam_fnl_mpih_object — FUNDAL MPI alias
fortran
use :: fundal_mpih_object, only : mpih_fnl_object => mpih_objectA direct re-export of FUNDAL's MPI handler under the FNL-namespaced type alias mpih_fnl_object. Provides rank/size queries, rank-prefixed console output, and timing utilities. No FNL-specific extensions are needed because FUNDAL's handler already covers GPU-aware MPI requirements.
Module summary
| Module | Role | Extends |
|---|---|---|
adam_fnl_library | Aggregate entry point | — |
adam_fnl_field_object | GPU field wrapper: device geometry, maps, ghost exchange | — |
adam_fnl_field_kernels | Gradient, L2 norm, ghost-cell pack/unpack | — |
adam_fnl_weno_object | GPU WENO coefficient mirror + ROR tables | — |
adam_fnl_weno_kernels | Upwind WENO reconstruction (!$acc routine seq) | — |
adam_fnl_rk_object | GPU RK stage storage and update dispatch | — |
adam_fnl_rk_kernels | Stage assign, accumulate, low-storage, final update | — |
adam_fnl_ib_object | GPU distance field + eikonal BC wrapper | — |
adam_fnl_ib_kernels | Eikonal evolution, sphere distance, momentum inversion | — |
adam_fnl_maps_object | GPU communication index tables + MPI buffer staging | — |
adam_fnl_mpih_object | FUNDAL MPI handler alias | — |
adam_fnl_poisson_fft_object | FFT Poisson solver, cuFFT on each rank's GPU (the slab transform; not re-exported by adam_fnl_library) | poisson_fft_object |
adam_fnl_poisson_mg_object | Multigrid Poisson solver on AMR, the FCG iteration and the V-cycle on each rank's GPU (not re-exported by adam_fnl_library) | poisson_mg_object |
adam_fnl_fdv_operators_library | Device-callable FD/FV spatial operators | — |
adam_fnl_mpih_global | Singleton — mpih_fnl | — |