Skip to content

Adaptive Mesh Refinement (AMR)#

This page describes SELF's approach to adaptive mesh refinement for discontinuous Galerkin spectral element models. Sections 1-5 document the 2-D implementation, which came first: the refinement trigger, and the staged design along which the invasive parts (dynamic mesh mutation, solution transfer, MPI re-partitioning, and GPU re-allocation) were reviewed and landed. Section 6 covers the 3-D stack - octrees and 2:1 face mortars - which is a deliberate transcription of the 2-D one and records only what differs.

AMR in SELF builds directly on the 2:1 nonconforming (mortar) interface support (see Nonconforming (Mortar) Interfaces). Refining a quadrilateral element into four children produces exactly the 2:1 hanging-node configurations that the mortar operators already couple correctly, so the flux-coupling layer that adaptive refinement needs is already in place and tested.

1. Scope and status#

Component Status
Refinement trigger (Legendre modal-decay indicator, CPU + GPU) Implemented
h-refinement primitives (isoparametric subdivision + refined connectivity) Implemented
Uniform h-refinement (conforming, serial) Implemented
Adaptive quad-forest (flagged refine / coarsen, level tracking) Implemented
Solution transfer (prolongation / restriction, conservative) Implemented
Forest face-neighbour queries + 2:1 balancing Implemented
Hanging-node / mortar-table + Mesh2D_t emission Implemented
Adaptation-epoch transfer plan (old-leaf → new-leaf mapping) Implemented
Model regrid (DGModel2D%Regrid) + AMR controller (serial, CPU/GPU) Implemented
Ultrasound point-source example + AMR visualization script Implemented
MPI dynamic re-partitioning / load balancing (v1: replicated forest, allgathered migration) Implemented
GPU device re-allocation for a changing element count Implemented (amortized high-water-mark storage, Stage 6b)
Device-side solution transfer (no host round trip on one GPU) Implemented (Stage 6a)
Geometry reuse for unchanged elements across an epoch Implemented (Stage 6c)

The full serial adaptive-refinement loop is wired into the library today: flag with the Stage-1 indicator, mutate the forest (Stage 2b), transfer the solution (Stage 3), balance and emit a runnable nonconforming Mesh2D_t (Stage 4). What remains is scaling that loop - dynamic MPI re-partitioning across ranks (Stage 5) and GPU device re-allocation for a changing element count (Stage 6) - each deferred behind this design so it lands as a self-contained, reviewable piece.

2. The refinement trigger (implemented)#

The trigger answers a single, element-local question: is the solution on this element well-resolved by the current polynomial degree, or does it need more resolution? It never mutates the mesh; it only produces a per-element flag. This makes it cheap, embarrassingly parallel, free of MPI communication, and identical in structure on CPU and GPU.

2.1 Legendre modal-decay (spectral) indicator#

For a nodal solution field \(u_{ij}\) on an element, we form its tensor-product Legendre modal expansion

\[ u(\xi,\eta) = \sum_{p=0}^{N}\sum_{q=0}^{N} \hat u_{pq}\, \tilde L_p(\xi)\,\tilde L_q(\eta), \]

where \(\tilde L_p\) are the \(L^2\)-normalized Legendre polynomials on \([-1,1]\) (\(\int_{-1}^1 \tilde L_p \tilde L_q\,d\xi = \delta_{pq}\)). The nodal-to-modal transform is the exact inverse of the Legendre Vandermonde \(V_{ip} = \tilde L_p(\xi_i)\) built from the interpolant control points, so the indicator is independent of the control-node type (Legendre–Gauss or Legendre–Gauss–Lobatto) and is exact for polynomial data. With the normalized basis the total modal energy equals the exact \(L^2\) energy on the reference element,

\[ E_{\text{tot}} = \sum_{p,q}\hat u_{pq}^2 = \lVert u\rVert_{L^2([-1,1]^2)}^2 . \]

The smoothness of the field is measured by how much of that energy sits in the highest modes. Defining the clipped energies that drop the highest one and two modes in each direction,

\[ E_{\text{clip1}} = \sum_{p,q\le N-1}\hat u_{pq}^2, \qquad E_{\text{clip2}} = \sum_{p,q\le N-2}\hat u_{pq}^2, \]

the smoothness ratio is the larger of the top-shell and next-shell energy fractions,

\[ S_e = \max\!\left(\frac{E_{\text{tot}}-E_{\text{clip1}}}{E_{\text{tot}}},\; \frac{E_{\text{clip1}}-E_{\text{clip2}}}{E_{\text{clip1}}}\right), \qquad \sigma_e = \log_{10} S_e . \]

A smooth, well-resolved field has energy concentrated in the low modes, so \(S_e\) is tiny and \(\sigma_e\) is very negative; an under-resolved field (a steep front or a discontinuity) leaves energy in the top modes, driving \(\sigma_e\) toward 0.

This is the modal-energy indicator of Persson & Peraire (2006). The second (next-shell) term is the robustification introduced by Hennemann, Rueda-Ramírez, Hindenlang & Gassner (2021) — the shock indicator used in Trixi.jl — which guards against the odd/even parity dropouts a single-mode measure suffers for symmetric data. It is only active for \(N\ge 2\).

2.1.1 The amplitude gate (relative energy floor)#

\(S_e\) is a ratio, so it is scale-free by construction: it reports how an element's energy is distributed across modes and says nothing about how much energy the element carries. Grid-scale residue at \(10^{-5}\) of the field's amplitude has exactly the same \(S_e\approx 1\) as a full-amplitude discontinuity.

Guarding the ratio with an absolute floor at machine epsilon does not fix this, because \(E_{\text{tot}}\) carries the square of the field amplitude while \(\varepsilon\) is a fixed \(2.2\times10^{-16}\): a wake at \(10^{-5}\) of the peak sits some ten orders of magnitude above the floor and is still judged on shape alone. Such elements never reach \(\sigma_{\text{coarsen}}\), so the mesh behind a passing front is never released and the refined region grows to the whole area the wave has swept (issue #162).

Each element therefore also carries a gate energy

\[ g_e = \sum_v w_v\,E_{\text{tot},e,v}, \]

a weighted sum of the per-variable modal energies. Since \(E_{\text{tot},e,v}\) is the exact \(L^2\) energy of variable \(v\) on the reference element, \(g_e\) is a convex quadratic functional of the state — a discrete entropy integral over the element when \(w\) is taken from the diagonal of an entropy Hessian. It is compared against a field-wide energy scale,

\[ g_e \;\le\; \max\!\left(\varepsilon,\; \texttt{relativeEnergyFloor}\cdot \texttt{energyScale}\right) \quad\Longrightarrow\quad S_e := 0 , \]

so a quiescent element is reported perfectly resolved (\(\sigma_e = \log_{10}\varepsilon\), hence COARSEN) whatever its modal shape.

Because energy is amplitude squared, relativeEnergyFloor is \(10^{dB/10}\) in amplitude terms. The default is \(10^{-12}\): amplitudes below \(10^{-6}\) (−120 dB) of the field scale count as quiescent. It sits below \(\varepsilon\) in real32, so in a single-precision build the absolute term dominates and the gate is inactive unless a caller raises it — the safe direction, since a real32 field cannot represent −120 dB structure anyway.

The energy hysteresis band#

As described so far the gate is a single hard cut, which reintroduces on the energy axis exactly the thrashing the two \(\sigma\) thresholds exist to prevent: an element whose energy drifts across that one value flips COARSEN ↔ REFINE on successive epochs, churning the mesh without the solution changing meaningfully. significantEnergyFloor opens a band instead:

gate energy \(g_e\) flag
\(\le\) quiescent floor COARSEN, whatever the spectrum says
between the floors KEEP — hold the mesh
\(>\) significant floor the spectrum decides, as before

The middle zone reads too weak to be worth spending levels on, too strong to declare resolved, and is stable under a small change in amplitude. Inside the band \(\sigma_e\) still reports the true spectrum — the band holds the flag, it does not falsify the diagnostic — so for those elements flag is deliberately not what thresholding \(\sigma_e\) would give.

It defaults to the quiescent floor, collapsing the band to the single cut, so behaviour is unchanged unless a caller opts in. test/refinement_indicator_2d_energy_hysteresis.f90 covers the three zones and asserts that an in-band energy drifting by a factor of ten either way does not move the flag, while the same drift across a degenerate band does.

Keep the band narrow. Its upper edge is functionally "do not spend levels below this energy" — the same knob as the floor, with the same cost in refinement depth. On the ultrasound benchmark's initial adaptation, with a quiescent floor of \(10^{-12}\):

significant floor 1e-12 3e-12 1e-11 3e-11 1e-10
elements / max level 328 / 2 328 / 2 268 / 1 268 / 1 268 / 1

A band up to ~3× the floor leaves depth untouched; at 10× it collapses a level, exactly as raising the floor to \(10^{-10}\) does. The band is a thrash damper, not a savings knob — widen it only as far as is needed to stop flags oscillating.

Raising the floor is not free. The gate cannot distinguish residue from the low-amplitude flank of a feature that is genuinely under-resolved, so an aggressive floor buys element count at the cost of refinement depth. Measured on the ultrasound point-source benchmark (examples/linear_euler2d_amr_ultrasound_pointsource, 30 epochs, maxLevel = 2):

relativeEnergyFloor initial adaptation mean elements final
0 (no gate) 328, level 2 1482 1732
1e-12 (default) 328, level 2 1057 (1.40x fewer) 1528
1e-8 268, level 1 2210 3304

At 1e-8 the gate also suppresses the source pulse's skirt, so the forest never reaches the level cap; RecommendedTimeStep is tied to the max level, so dt doubles, and the resulting time-integration error drives enough later refinement to roughly double the final mesh. Measure before raising this. A related trap: a 2-D cylindrical pulse leaves a genuine algebraically-decaying tail behind its front (Huygens' principle fails in even dimensions), which is physics rather than residue and should not be gated away.

energyScale defaults to the largest \(g_e\) over the elements, reduced with MPI_MAX when the indicator is given a communicator, so the gate follows a decaying front. Two consequences worth knowing:

  • the automatic scale makes the flags depend on a floating-point reduction, which can differ at round-off between rank counts — SetEnergyScale pins it and is the deterministic escape hatch;
  • once a wave has left the domain entirely the scale collapses onto the residue's own peak and the gate stops discriminating; the absolute \(\varepsilon\) floor is what bounds that case.

Note that raising \(\sigma_{\text{coarsen}}\) is not a substitute. In a low-amplitude wake \(\sigma_e\) is near \(0\), so any threshold high enough to coarsen the wake is also high enough to stop the front from refining; the two decisions cannot be separated by thresholds alone, which is what motivates an amplitude-based gate.

2.2 Refine / keep / coarsen semantics#

With user-supplied thresholds \(\sigma_{\text{refine}} > \sigma_{\text{coarsen}}\), each element is flagged

Condition Flag
\(\sigma_e > \sigma_{\text{refine}}\) SELF_AMR_REFINE (+1)
\(\sigma_e < \sigma_{\text{coarsen}}\) SELF_AMR_COARSEN (−1)
otherwise SELF_AMR_KEEP (0)

Two thresholds (rather than one) give a hysteresis band that prevents elements from thrashing between refine and coarsen on successive checks. Recommended starting values for double precision are \(\sigma_{\text{refine}}\approx-3\) and \(\sigma_{\text{coarsen}}\approx-8\); the best values are problem-dependent and are deliberately required arguments rather than hidden defaults.

2.3 API#

The indicator lives in SELF_RefinementIndicator_2D and follows the standard SELF backend pattern: a portable base type with a do concurrent implementation (SELF_RefinementIndicator_2D_t), a thin CPU extension, and a GPU extension that runs a device kernel (SELF_Refinement.cpp) and copies the resulting flags back to the host, where the mesh-adaptation logic runs.

use SELF_RefinementIndicator_2D

type(RefinementIndicator2D) :: amr

! interp : the model interpolant (Lagrange); nElem : rank-local element count
call amr%Init(interp, nElem, refineThreshold=-3.0_prec, coarsenThreshold=-8.0_prec)

! Amplitude gate (all optional; these are the defaults unless set).
call amr%SetRelativeEnergyFloor(1.0e-12_prec) ! 0 disables the gate
call amr%SetRelativeEnergyFloor(1.0e-12_prec, &  ! ... or open a hysteresis band
                                significantEnergyFloor=1.0e-10_prec)
call amr%SetEnergyScale(pRef**2)              ! pin the scale; ClearEnergyScale undoes it
call amr%SetEnergyWeights(w)                  ! per-variable gate weights, e.g. from the entropy

! solution : the model's MappedScalar2D (or any Scalar2D); ivar selects the driving
! variable, or SELF_AMR_ALLVARS (=0) reduces over all variables (most conservative).
! comm : optional, makes the automatic energy scale global (MPI_MAX).
! gate : optional, a caller-computed per-element gate energy replacing the weighted sum.
call amr%Estimate(solution, ivar=1)

! amr%indicator(1:nElem)  -- per-element sigma_e (host)
! amr%flag(1:nElem)       -- per-element SELF_AMR_REFINE / KEEP / COARSEN (host)
! amr%gate(1:nElem)       -- per-element gate energy actually used (host, diagnostic)
n = amr%CountFlagged(SELF_AMR_REFINE)

The estimate runs in two phases: an element-local, parallel (device) pass that forms the spectra, the raw \(S_e\) and \(g_e\); then a host pass that reduces for the energy scale, applies the gate, takes the \(\log_{10}\) and sets the flags. The second phase is host-side because the scale needs a reduction that no element-local pass can supply, and it is shared by both backends, so CPU and GPU produce identical flags. AMRController2D forwards relativeEnergyFloor and energyWeights from its own Init and re-applies them after each epoch's indicator resize; it also supplies the communicator so the scale is global.

The GPU device kernel uses per-thread scratch bounded to \(N \le 15\); higher degrees fail loudly at Init with a pointer to the bound (AMR2D_MAXNP in SELF_Refinement.cpp).

2.4 Validation#

test/refinement_indicator_2d_spectraldecay.f90 checks the indicator against fields with known spectra, on both Gauss and Gauss–Lobatto nodes:

  • a constant field (energy only in mode 0) floors out → COARSEN;
  • a pure highest tensor mode \(\tilde L_N(\xi)\tilde L_N(\eta)\) has \(S_e=1\) so \(\sigma_e = 0\) to roundoff → REFINE (this also verifies the transform recovers a single mode exactly);
  • a low-degree polynomial (degree \(\le N-2\)) has no top-shell energy → COARSEN;
  • a steep tanh front is under-resolved → REFINE.

2.5 h-refinement primitives and uniform refinement (implemented)#

The mesh-mutation core of Stage 2 is implemented as two additive modules that leave every existing mesh type and interface untouched:

  • SELF_RefinementPrimitives_2D — the element-local, dependency-light pieces:

    • SubdivideNodeCoords performs exact isoparametric subdivision of one element's geometry into its four children. The parent geometry is the degree-nGeo Lagrange interpolant through the element's geometry nodes; each child node coordinate is that interpolant evaluated at the corresponding point of the parent reference square. This is exact for straight-sided and curved (isoparametric) elements alike and for any control node type, so refinement never perturbs the represented domain.
    • RefineConnectivity builds the refined mesh's connectivity by pure deterministic integer bookkeeping. Sibling faces interior to a parent are same-orientation (flip 0); a child face on a parent boundary inherits the parent face's neighbor, side, and flip, with the sub-position pairing across the face taken from the parent flip. Because refinement is orientation preserving, this reproduces exactly the flips a corner-node matching pass would compute, without requiring globally consistent node ids (which structured meshes do not guarantee) and without any coordinate hashing.

  • SELF_MeshRefinement_2D — UniformRefineMesh(meshIn, meshOut) assembles a fully-formed Mesh2D_t with 4× the elements from those primitives: child geometry by isoparametric subdivision, connectivity and flips inherited from the base, boundary-condition and material metadata carried over, and a (serial) domain decomposition. The result is conforming (no hanging nodes), so it needs neither mortars nor 2:1 balancing and is immediately usable for geometry generation and time stepping. It is a genuine capability on its own (e.g. grid convergence studies) and it exercises all of the subdivision, connectivity, and array machinery that adaptive refinement will reuse.

Uniform refinement is currently serial (it replicates the input mesh's single-rank MPI state rather than re-decomposing); multi-rank refinement is Stage 5.

Validation. The primitives are unit-tested in isolation (exact child geometry for affine and curved elements; neighbor / side / flip / global-side-id reciprocity of the refined connectivity). The end-to-end UniformRefineMesh is covered by test/mesh2d_uniform_refine.f90: refining a structured mesh quadruples the element count, conserves the domain area (integral of the Jacobian) to roundoff, keeps every Jacobian strictly positive, produces reciprocal interior connectivity, and doubles the physical-boundary side count.

2.6 Adaptive quad-forest (implemented)#

SELF_QuadTreeMesh_2D provides the adaptive mesh-mutation data structure of Stage 2b. Each base element is the root of a quadtree; QuadTreeMesh2D stores a growable node pool (level, parent, per-node child pointers, originating root) plus the base geometry, and maintains the active leaf set.

  • Init(mesh) seeds one root per base element (all leaves at level 0).
  • AdaptFromFlags(flag) consumes a per-leaf flag array indexed exactly like the Stage-1 indicator's flag(:) (+1 refine, -1 coarsen, 0 keep): flagged leaves are subdivided into four children, and a family of four leaf siblings is merged back only when all four are flagged coarsen (the standard de-refinement rule). Refine and coarsen are resolved against the same pre-adaptation snapshot so they never interfere.
  • LeafCoords(i, geomInterp, coords) regenerates any leaf's physical geometry by repeated exact isoparametric subdivision of its root along the quadtree path (reusing the §2.5 primitive), so the forest never stores redundant coordinates.
  • Node storage grows by amortized doubling; the leaf set is always recovered by traversal from the roots, which makes nodes orphaned by coarsening invisible without an explicit free list.

This wires directly to the trigger: indicator%Estimate(solution, ivar) then forest%AdaptFromFlags(indicator%flag) performs one adaptation step. What the forest deliberately does not yet do is answer face-neighbour queries, enforce 2:1 balance, or emit a solver-ready Mesh2D_t with hanging-node mortars - an adaptively refined forest is generally nonconforming, and turning it into a runnable mesh is Stage 4.

The forest is unit-tested (refine/coarsen leaf and level bookkeeping, the four-sibling coarsening guard, amortized capacity growth, and leaf geometry matching direct subdivision at multiple levels).

2.7 Solution transfer (implemented)#

SELF_SolutionTransfer_2D moves the prognostic solution with the mesh when an element is refined or coarsened, reusing the mortar operators already built on Lagrange:

  • ProlongToChildren(interp, nVar, uParent, uChildren) samples the parent's degree-N nodal polynomial onto its four children - a tensor product of mortarR in each direction. Exact interpolation, no loss.
  • RestrictFromChildren(interp, nVar, uChildren, uParent) L2-projects the four children back onto the parent - a tensor product of mortarP (the adjoint of mortarR, already carrying the 1/2 per-direction sub-edge Jacobian for solution traces).

Both are element-local and portable (host do concurrent); the driver maps forest parent/child relations onto the element index ranges. From the 1-D mortar identities they inherit exactly:

  • Reversibility - RestrictFromChildren(ProlongToChildren(u)) = u to roundoff (refine then immediately coarsen leaves the solution unchanged), from sum_k P_k R_k = I per direction.
  • Conservation - the reference-cell integral of the restricted parent equals the sum of the children's, i.e. sum_ij w_i w_j u_parent = (1/4) sum_c sum_ij w_i w_j u_child_c; weighted by the geometry Jacobian this is conservation of the cell-integrated quantity.

Validated at two levels: a unit test drives the tensor operators with the exact mortar matrices (prolong reproduces a degree-N polynomial at the child nodes; prolong-then-restrict is the identity to roundoff; discrete conservation defect is zero), and test/solution_transfer_2d.f90 checks the end-to-end behaviour against Stage-2 geometry - prolonging a coarse field onto a UniformRefineMesh and back is reversible, and int u dA (with the geometry Jacobian) is identical on the coarse and refined meshes.

2.8 Forest face-neighbours and 2:1 balancing (implemented)#

SELF_QuadTreeMesh_2D also answers face-neighbour queries on the forest and enforces the 2:1 balance condition - the first part of Stage 4 and the prerequisite for mortar generation. The forest stores the base mesh's root face connectivity (rootNbr / rootNbrSide / rootFlip, from sideInfo; a conforming base is assumed).

  • FaceNeighbor(node, s, nbr, ns, nf) returns the equal-or-larger neighbour across local side s by the classic quadtree ascend/descend search: cross to a sibling when the face is interior to the parent, otherwise ascend to the parent's neighbour and descend one level, matching sub-positions across the face through the base flip. nbr is either a leaf at any level <= level(node) or an internal node at exactly level(node); ns / nf are the neighbour's facing side and the edge flip. Consequently a 2:1 hanging face is precisely "nbr is a leaf with level(nbr) = level(node)-1" and finer neighbours are precisely "nbr is internal" - exactly the classification Stage 4b needs to emit mortarInfo.
  • Balance2to1() iterates to a fixed point: any leaf whose equal-or-larger neighbour is a leaf two or more levels coarser refines that neighbour, and the (possibly rippling) refinement repeats until no face violates the condition.
  • MaxLevelJump() reports the largest level difference across any leaf face (0 conforming, 1 for a balanced adaptive forest) - a cheap invariant for tests and drivers.

Unit-tested standalone: level-0 neighbour queries on a structured base (including boundaries and directional reciprocity); a uniformly refined forest is conforming (MaxLevelJump = 0); an adaptive refinement that creates a two-level jump is reduced to one level by Balance2to1 (rippling into the coarse neighbour); and equal-level leaf-neighbour reciprocity holds across the balanced forest.

2.9 Mesh emission (implemented) — closing the loop#

SELF_AdaptiveMesh_2D's EmitMesh(forest, baseMesh, outMesh) turns a 2:1-balanced forest into a solver-ready Mesh2D_t. Each leaf becomes an element (leaf-list order); for every leaf face the Stage-4a FaceNeighbor classification drives the emitted connectivity:

  • domain boundary → sideInfo(3)=0, sideInfo(5) = the base element's BC id on that side;
  • same-level leaf → a conforming interior side (sideInfo(3)=neighbour, (4)=10*side+flip, a shared global side id);
  • one-level-finer neighbour → this leaf is the big side of a 2:1 mortar; the two small elements are the finer neighbour node's children on the shared face, with the big edge coordinate [-1,0]/[0,1] mapped to the correct child through the face flip;
  • one-level-coarser neighbour → a small side, filled when its big side is processed.

Mortar sides carry sideInfo(1)=mortar index and sideInfo(3)=sideInfo(5)=0, and the emitted mortarInfo(1:8, :) follows the exact layout of the hand-built SimpleMortarMesh (big elem/side; small elem + 10*side+flip per sub-edge; two sub-edge global side ids), so the existing mortar side-exchange, projection, and flux machinery consume it unchanged. Leaf geometry comes from LeafCoords; the output uses a serial decomposition (Stage 5 will re-partition).

Validated at two levels. A standalone structural test confirms - on an adaptively refined, balanced forest - the side classification is exclusive, conforming sides are reciprocal, every mortar's big and small sides reference the same mortar index, and (geometrically) each big face is exactly tiled by its two small faces. The full-pipeline CI test test/adaptive_mortar_2d.f90 refines a structured mesh, balances, emits the mesh, builds its geometry (strictly positive Jacobians; total area equal to the base mesh), and runs the real SideExchange + MortarExchange on a globally linear field: the external trace matches the interior trace to roundoff on every conforming and 2:1 mortar side - the same criterion the hand-built mortar-mesh tests use, now on a mesh produced entirely by the AMR pipeline.

With this, the serial loop closes: indicator → forest.AdaptFromFlags → (transfer solution) → forest.Balance2to1 → EmitMesh yields a runnable adaptive mesh.

2.10 Adaptation-epoch transfer plan (implemented)#

SELF_TransferPlan_2D is the driver layer of Stage 3: it connects the element-local transfer operators (§2.7) to an actual forest mutation. One adaptation epoch is: snapshot the leaf list (nOld, oldLeaf), mutate the forest (at most one AdaptFromFlags, then any number of refinements — Balance2to1, RefineNode), then BuildTransferPlan(forest, nOld, oldLeaf, plan). The plan records, for every new leaf in leaf-list order (the element ordering EmitMesh produces), where its solution comes from in the old element ordering:

  • copy — the leaf survived unchanged;
  • prolong — the leaf descends from an old leaf; the old polynomial is interpolated down the quadtree path, one step per level, so a fresh child re-refined by balance ripple in the same epoch is handled by depth > 1;
  • restrict — the leaf is (an ancestor of) a coarsened family; the four old children are L2-projected onto their parent, then prolonged down any further steps (depth > 0 occurs when a just-coarsened parent is immediately re-refined by balancing).

Reconstruction after the fact is possible because forest node ids are stable: refinement appends nodes and coarsening only detaches children, whose level/parent/quadrant entries persist. Each new leaf ascends its parent chain until it meets an old leaf or a complete old-leaf family; a snapshot that cannot explain a leaf fails loudly. ApplyTransferPlan executes the plan on nodal data in the MappedScalar2D%interior layout and inherits the §2.7 identities (exact prolongation, conservative restriction, exact refine-coarsen round trips).

test/transfer_plan_2d.f90 validates one epoch that simultaneously coarsens (then re-refines) a family, refines a leaf whose child is refined again (depth-2 prolongation), and lets 2:1 balancing ripple into an untouched root: the classification multiplicities are checked exactly; a bilinear field is reproduced at the emitted new mesh's nodes to roundoff through all three transfer kinds; the Jacobian-weighted global integral of a non-polynomial field is conserved to roundoff; and refine-everything/coarsen-everything across two epochs is the identity.

2.11 Model regrid and the AMR controller (implemented)#

Two pieces close the loop around a live, time-stepping model:

  • DGModel2D%Regrid(mesh, geometry) rebinds a model to a new mesh/geometry pair: the mesh-sized solution storage is reallocated and the boundary-condition registrations and maps are rebuilt (mirroring the mesh-sized portion of Init/Free, including the GPU backend's BC device arrays), while everything that is not mesh-sized is preserved — the time state (t, dt, entropy, IO counter), the time-integrator selection, configuration flags, and model-specific parameters, all of which a fresh Init (intent(out)) would reset. The solution interior is left for the caller to fill via the §2.10 transfer.

  • SELF_AMRController_2D owns the forest, the indicator, and the meshes/geometries it emits (double-buffered), and performs one adaptation epoch per Adapt(model, adapted) call: estimate → cap refine flags at a configurable maxLevel → spread refine flags to face neighbours for nHalo passes (so a feature moving at speed \(c\) stays inside the refined band when the adaptation cadence satisfies \(k\,\Delta t\,c \le n_{halo} h_{fine}\)) → AdaptFromFlags + Balance2to1 (a no-op epoch leaves the model untouched) → BuildTransferPlan + EmitMesh + new SEMQuad → model%Regrid, apply the transferred solution, upload to device. RecommendedTimeStep(dtBase) = dtBase / 2^{MaxLevel} gives the level-based explicit-stability bound to pass to ForwardStep after each epoch — exact for the quadtree, whose children are exact half-scale subdivisions.

test/lineareuler2d_amr_soundwave.f90 runs the full loop on a deliberately under-resolved acoustic pulse (LinearEuler2D, radiation boundaries, RK3): the initial adaptation refines around the pulse up to the level cap; every mid-run adaptation conserves the Jacobian-weighted global integral of each prognostic variable to roundoff; the model's time and parameters survive regridding; the mesh evolves as the wave propagates; and the acoustic energy stays finite and non-increasing across the whole adaptive run.

3. Comparison with Trixi.jl#

Trixi.jl offers two families of indicators that drive both shock capturing and AMR:

  • IndicatorHennemannGassner — the Persson–Peraire modal-energy indicator with the next-shell robustification, mapped through a logistic function to a blending coefficient \(\alpha\in[0,1]\). SELF's indicator uses the same modal-energy quantity (§2.1); it stops at \(\sigma_e = \log_{10} S_e\) and thresholds directly rather than forming \(\alpha\), because AMR needs a discrete refine/keep/coarsen decision, not a continuous blend.
  • IndicatorLöhner — a normalized second-difference (curvature) estimate of a nodal quantity. It is cheaper (no modal transform) but noisier and less tightly coupled to spectral resolution. It is a natural future alternative behind the same flag interface.

Two design points worth noting for anyone extending this:

  • Trixi additionally offers a IndicatorMax-style controller and clip/smoothing passes across element neighbors. SELF's indicator is strictly element-local; neighbor smoothing (to avoid isolated refined elements) belongs in the Stage-4 balancing step, not the trigger.
  • Trixi's ControllerThreeLevel maps indicator values to target refinement levels with hysteresis. SELF's two-threshold refine/keep/coarsen flag is the minimal equivalent; a level-target controller can be layered on top once the tree (Stage 2) exists.

An alternative spectral trigger is the Mavriplis (1994) decay-rate estimator, which least-squares-fits \(\lvert\hat u_p\rvert \sim c\,e^{-\sigma p}\) to the modal coefficients and estimates the truncation error from the fitted decay rate and tail. It is more informative for p-adaptivity (choosing how much to refine) but more fragile than the energy-fraction measure; it is a candidate for a future p/hp-adaptive extension.

4. Staged plan for the remaining machinery#

The current mesh (Mesh2D_t) is statically allocated: elemInfo(1:6,1:nElem), sideInfo(1:5,1:4,1:nElem), nodeCoords(...,1:nElem), and a hand-built mortarInfo table. Dynamic AMR requires making the element set mutable while preserving every invariant the solver relies on. The following stages are each independently reviewable and testable.

Stage 2 — h-refinement mesh mutation#

  • (done) Generate child nodeCoords from the parent geometry map by exact isoparametric subdivision, and build refined connectivity deterministically (SubdivideNodeCoords, RefineConnectivity). Metric terms/Jacobians for children come from the existing geometry routines with the refined nodeCoords as input — no change to the geometry algorithms.
  • (done) Uniform (conforming) refinement end to end (UniformRefineMesh), serial.
  • (2b, done) An explicit quadtree / forest-of-quadtrees parent–child structure (SELF_QuadTreeMesh_2D): each leaf is an active element; adaptive refinement replaces a flagged leaf with four children (each spanning a reference sub-quadrant), coarsening merges four siblings back to their parent, with amortized-capacity node growth and traversal-based leaf enumeration. Driven directly by the Stage-1 indicator flags via AdaptFromFlags.
  • (next) Storage compaction: reclaim nodes orphaned by coarsening (currently the node pool grows monotonically).

Stage 3 — Solution transfer (prolongation / restriction) — done#

Implemented in SELF_SolutionTransfer_2D (see §2.7):

  • Prolongation (parent → 4 children): tensor product of the mortar restriction operator Lagrange%mortarR — exact interpolation of the parent polynomial onto the children.
  • Restriction (4 children → parent): tensor product of the \(L^2\)-projection adjoint Lagrange%mortarP; conservative by construction (\(\sum_k P_k R_k = I\) and discrete conservation), so coarsening preserves cell-integrated quantities.
  • The transfer is a separate operator applied between time steps and does not touch the solver's floating-point reductions.

Stage 4 — Mortar regeneration and 2:1 balancing#

  • (4a, done) Face-neighbour navigation on the forest (FaceNeighbor, ascend/descend quadtree search across the base root connectivity, honouring base side pairings and flips) and 2:1 balance (Balance2to1: neighbours more than one level coarser than a leaf are refined, rippling to a fixed point). See §2.8.
  • (4b, done) EmitMesh (SELF_AdaptiveMesh_2D) rebuilds sideInfo + mortarInfo from the balanced tree and emits a solver-ready Mesh2D_t: every face where FaceNeighbor sees a finer/coarser neighbour becomes a 2:1 mortar (the configuration the solver already handles); same-level faces are conforming; leaf geometry comes from LeafCoords. See §2.9.
  • (next) Optional neighbour smoothing of the trigger flags (avoid isolated refined elements), on top of the balance pass.

Stage 5 — MPI dynamic re-partitioning — implemented (v1)#

Status: implemented. QuadTreeMesh2D%InitGlobal builds the rank-replicated forest from allgathered global base tables; EmitMesh builds the global connectivity/mortar tables on every rank and stores only its contiguous slice of a freshly generated decomposition (sideInfo(3) global ids, global nUniqueSides, fully replicated mortarInfo, exactly the invariants SideExchange/MortarExchange require); the controller allgathers the indicator flags per epoch and migrates the solution through an allgathered global old field applied to the rank-local range (ApplyTransferPlanRange). Validated by test/lineareuler2d_amr_soundwave_mpi.f90 (2 ranks): the global element trajectory and entropy history match the serial run, transfers conserve globally, and a leaf-list checksum confirms forest replication. The point-to-point migration upgrade remains open (v2).

Two observations make a correct first version tractable:

  • The forest's leaf list is already a space-filling curve: root-major depth-first traversal is Morton order within each quadtree, so "SFC partitioning" is just contiguous ranges of the existing leaf list — the same contiguous-ownership model (offsetElem) the domain decomposition already uses.
  • The forest is cheap (a few integers per node), so it can be replicated on every rank. If all ranks apply identical flags, they compute identical adapted/balanced forests, transfer plans, and emitted global connectivity — deterministically, with no communication beyond the flags themselves.

Sub-stages:

  • (5a) Global flags + replicated mutation. The controller allgathers the rank-local indicator flags (by the decomposition's element offsets) into a global per-leaf flag array; every rank then runs the same cap/halo/AdaptFromFlags/Balance2to1 sequence on its forest copy. One small collective per epoch, outside the time-stepping loop.
  • (5b) Multi-rank EmitMesh. Every rank builds the same global sideInfo/mortarInfo (deterministic from the forest), then decomposes it exactly the way the existing mesh constructors do for nRanks > 1, so SideExchange/MortarExchange consume the result unchanged. Repartitioning is implicit: each epoch's emitted mesh is re-decomposed over the new leaf list, so equal-count SFC arcs move with the refinement.
  • (5c) Solution migration. v1: allgatherv the old rank-local solutions into a global old field, then ApplyTransferPlan only for the new rank-local elements. Correct and simple; memory is one global solution copy per rank (fine at single-node scale). The point-to-point upgrade (send exactly the source elements each rank's plan references) is a drop-in replacement behind the same interface.
  • Tests on ≥ 2 ranks: forest determinism across ranks (identical leaf checksums after an epoch), global conservation of the transferred solution (mpi_allreduce), and the AMR soundwave regression run distributed.

Stage 6 — GPU device re-allocation (implemented)#

Both parts are implemented and measured on one MI300X (gfx942, ROCm 6.4.3, exclusive node), ultrasound example, 20 epochs, mean of three runs. The starting point was the "correct but unamortized" form: every adapting epoch freed and re-initialized all model storage, re-uploaded mesh/geometry, and moved the solution through a host round trip.

adaptation total AMR share of epoch loop
before 0.998 s 51.7%
+ 6a device-side transfer 0.830 s 47.0%
+ 6b amortized capacity 0.616 s 39.4%

Time integration is unchanged throughout (0.93-0.95 s), so the whole gain is in adaptation: 38% off the cost of an epoch. On the curved multi-material bone-and-marrow case the AMR share falls from 19.7% to 14.2%.

  • (6a) Device-side transfer — implemented. StageSolutionForTransfer and ApplyTransferPlan are type-bound on the model, using the same backend split as Regrid. The GPU override stages the pre-regrid solution device-to-device and applies the plan in TransferSolution_2D_gpu (src/gpu/SELF_SolutionTransfer.cpp), so a single-GPU adapting run moves no solution data across the host link. The kernel applies only the mortar operator pair of the child on the recorded path rather than prolonging to all four and discarding three, so device and host agree to round-off while conservation stays exact. Multi-rank runs keep the host allgather path, as the sequencing note below anticipated.
  • (6b) Amortized capacity — implemented. The prediction recorded here previously was that "at typical cadences allocation cost is expected to be far below the re-upload and geometry-generation cost". The measurement contradicted it: the free/re-initialize cycle was the single largest component of an adaptation at 51.5%, ahead of both the transfer and geometry regeneration. The cost was not the device allocator (hipMalloc + hipFree were only ~8% of an adaptation) but the host-side work a fresh allocation drags along - allocating, zeroing, rebuilding metadata and equation parsers, and uploading the zeros. Resize on the data classes now reuses storage via rank-1 pools with pointer remapping (see src/SELF_DataPool.f90 for why this avoids changing any kernel stride).

Behaviour change to be aware of. On a GPU build the transferred solution is left on the device, so solution%interior (the host mirror) is stale after Adapt. This matches the rest of the time loop, where the device is authoritative and a caller wanting host data calls UpdateHost() first, as Write_DGModel2D_t does. It used to be incidentally fresh because the transfer ran on the host.

How this scales with refinement depth. Depth is the knob that reaches ultrasound length scales (the base mesh stays 16x16; Lr and the epoch length track maxLevel). Adaptation cost grows faster than time integration, so the AMR share rises with depth - but it rises from a lower base and more slowly than before:

maxLevel f0 elements baseline + 6a/6b + 6c geometry reused
3 0.10 MHz 2,092 51.7% 40.1% 38.1% 72.8%
4 0.20 MHz 4,156 56.8% 43.1% 42.1% 73.0%
5 0.40 MHz 7,528 58.9% 45.9% 45.1% 73.9%
6 0.80 MHz 14,896 61.6% 48.4% 47.6% 74.6%
7 1.60 MHz 30,760 62.9% 50.6% 49.9% 75.4%
8 3.20 MHz 66,616 out of memory 51.8% 51.0% 76.5%

maxLevel 8 could not run at all before this work: it exhausted the device after ~555 adaptations (see the leak below). It now completes, which is the first datapoint at a production-relevant frequency. The reused fraction rises with depth, as expected - the refined annulus is a thinner shell relative to the whole mesh - so geometry reuse pays off more, not less, at production resolution. Per-epoch element counts are identical at every depth before and after 6c, which is the check that the incremental geometry changed no physics.

A leak this work exposed. boundarynormal_gpu was allocated in Init_MappedScalar2D and never freed, leaking 12.8 kB per element per adaptation at N=7, nvar=5 across the model's five MappedScalar2D fields. It aborted a 640-epoch maxLevel-8 run with an out-of-memory error after ~555 adaptations on a 192 GiB device. Fixed, with a device-memory regression test (test/data_2d_device_memory.f90) since nothing in the suite had been watching device memory.

Stage 6c - geometry (implemented)#

Re-profiling after 6a/6b put geometry at the top of an adaptation: 46.8%, ahead of Regrid, split between allocate/zero/free churn (Init_SEMQuad + Free_SEMQuad, 21.6%) and real computation (GenerateFromMesh_SEMQuad, 25.2%). Both are addressed:

  • Persistent geometry. Tensor2D joins the amortized-storage scheme (it was the last class on plain allocate/hipMalloc, and dxds/dsdx are the largest arrays in SEMQuad), SEMQuad gains Resize, and the controller owns two long-lived geometry buffers instead of allocating one per epoch. The nGeo -> N interpolant and node-coordinate staging that GenerateFromMesh used to rebuild on every call are cached, since only the element count varies between epochs.
  • Incremental generation. Most elements do not change in an epoch, and the transfer plan already names them: SELF_TRANSFER_COPY is set only when the new and old element are the same forest node. Because LeafCoords is a pure function of root coordinates, level and quadrant path, and per-element geometry generation is strictly element-local, such an element's whole geometry block is reproducible bit for bit - so it is copied forward rather than recomputed. The changed elements are generated compacted into a scratch geometry and scattered into place, which keeps the existing generation loops untouched. The two buffers alternate so the previous epoch's geometry is readable while the new one is assembled.
  • Multi-rank. sourceElem indexes the global old element list while a rank holds only its slice, so reuse additionally requires the source to be locally owned - a range test against the previous decomposition. On one rank that is always true; on several, migrated elements are regenerated. No communication is added.

Measured: 72.8% of elements reused per epoch at maxLevel 3, rising to 76.5% at maxLevel 8. Adaptation cost 0.616 s -> 0.585 s, and the bone-and-marrow case 14.2% -> 10.7% AMR share.

Correctness is checked numerically rather than argued. test/geometry_2d_reuse.f90 requires every kept element to match its previous-epoch source bit for bit, and adds a positive-Jacobian check plus the discrete metric identity sum_i d(J dsdx(i,j))/dxi_i = 0, which the suite did not previously test. test/geometry_2d_reuse_mpi.f90 covers the multi-rank branch, asserting both that reuse and regeneration actually fire and that the resulting mixed assembly equals a from-scratch generation exactly. Three env-gated switches (SELF_AMR_GEOM_NO_REUSE, SELF_AMR_GEOM_FULL, SELF_AMR_GEOM_VERIFY) allow a suspected geometry problem to be localized in one run without a rebuild.

A correctness bug the amortized scheme introduced#

Worth recording, because the reasoning error is easy to repeat. Stage 6b removed the UpdateDevice from Resize on the grounds that uploading freshly zeroed arrays is waste when the caller overwrites them. That is true of the solution but not of every field, and it silently dropped an invariant Init had provided: after Init a field's device buffer was defined, because the zeros had been uploaded. After Resize it held whatever was previously in that memory, and a fresh hipMalloc is uninitialized.

The low-storage Runge-Kutta update reads its accumulator before writing it - grk = rk_a*grk + dSdt with rk3_a(1) = 0, where grk is workSol%interior_gpu. Multiplying by zero annihilates any finite leftover, which is why this survived initial testing, but 0*NaN and 0*Inf are NaN. A resized workSol buffer that happened to contain a NaN bit pattern poisoned the solution on the next step; the wave died, the indicator found nothing to refine, and the mesh coarsened back to the base mesh. It presented as an intermittent, allocation-history-dependent failure - the same source on the same node passing at one time and failing at another - and no test caught it, because none runs enough GPU adaptation epochs for a stale buffer to land on such a pattern. EnsureDeviceBuffer now zeroes on resize, restoring the invariant unconditionally; a device-side fill costs HBM bandwidth rather than a PCIe transfer, so the point of skipping UpdateDevice is preserved.

The general lesson: "the caller overwrites it before reading" is a per-field claim, not a property of an amortized-storage scheme.

Still open. The geometry copies and the device upload are both still done the simple way: element blocks are copied on the host (roughly eleven array-slice copies per element), and the whole geometry is re-uploaded each epoch regardless of how little changed. That is why 6c yields ~5% rather than the ~35% a naive reading of the 46.8% attribution would predict. The fix is a device-side indexed-copy kernel - structurally the 6a transfer kernel without the operator - plus uploading only the compacted changed elements. There is now a measurement justifying it rather than an estimate. Mesh emission (EmitMesh, ~5% of an adaptation) is also still a full rebuild per epoch.

Driver integration#

Adaptation runs between time steps (e.g. every k steps): estimate → flag → (optional neighbor smoothing) → refine/coarsen tree → transfer solution → regenerate mortars/balance → repartition (MPI) → re-upload (GPU). The time-integration loop itself is unchanged; the solver sees only a (possibly) different, still-valid mesh at the top of the next step.

5. Application roadmap: LinearEuler2D ultrasound point source under AMR#

This section maps the serial AMR machinery above onto a first running application: a single point-source wavelet in the ultrasound frequency range, propagating through a ~1 m × 1 m domain with the 2-D linear Euler model, on a dynamically adapting mesh. Serial CPU and single GPU are the initial targets (Stage 5 MPI repartitioning is not required). Each gap below is an additive, independently testable piece; existing model/mesh interfaces stay untouched.

5.1 What already works (no changes needed)#

  • Flux coupling on adaptive meshes — EmitMesh produces the same mortarInfo layout the LinearEuler2D mortar tests (test/lineareuler2d_mortar_soundwave.f90) already exercise.
  • Per-epoch time step — ForwardStep(tn, dt, ioInterval) takes dt on every call, so a driver that re-computes dt after each adaptation needs no time-integrator changes.
  • Output of a changing mesh — every HDF5 snapshot written by WriteModel carries its own /controlgrid/geometry alongside the solution, so per-snapshot meshes are already representable; pyself reads them file-by-file.
  • Initial condition — SphericalSoundWave (Gaussian pressure pulse) generates an outgoing wavelet whose spectral content is set by the pulse half-width Lr; choosing Lr of a few millimetres puts the dominant wavelength in the ultrasound band with zero model changes.

5.2 Gap 1 — Transfer plan: old-leaf → new-leaf solution mapping — implemented#

Status: implemented in SELF_TransferPlan_2D (see §2.10), validated by test/transfer_plan_2d.f90.

AdaptFromFlags + Balance2to1 mutate the forest but record no correspondence between the pre- and post-adaptation leaf lists, which the Stage-3 transfer operators need. Because node ids are stable in the pool (coarsening only orphans nodes), the plan can be built after the fact from a snapshot of the old leaf array:

  • new leaf is an old leaf → copy;
  • new leaf descends from an old leaf → prolong along the quadtree path (Balance2to1 ripple can refine a fresh child again, so prolongation must handle multiple levels, applying ProlongToChildren per step of the path);
  • new leaf is the parent of four old leaves → restrict (always exactly one level: AdaptFromFlags coarsens one level per call and balancing never coarsens).

Deliverable: a BuildTransferPlan (forest + saved old-leaf list → typed plan) and an ApplyTransfer driver mapping solution%interior(:,:,oldIdx,:) to the new element ordering. Tests: adapt→transfer conservation of ∫u dA (Jacobian-weighted), refine-then-coarsen reversibility through a full plan, and a balanced two-level ripple case.

5.3 Gap 2 — Model regrid: rebinding a live DGModel2D to an emitted mesh — implemented#

Status: implemented as DGModel2D%Regrid + SELF_AMRController_2D (see §2.11), validated by test/lineareuler2d_amr_soundwave.f90. The level-based time step of §5.4 ("now") is RecommendedTimeStep.

DGModel2D storage (7 MappedScalar/Vector objects) is sized by nElem at Init, and Init is intent(out) — it resets t, model parameters (rho0), BC registrations, and the IO counter. Rather than making the model mutable, add an external AMR controller module (e.g. SELF_AMRController_2D) that owns the forest, the indicator, and double-buffered Mesh2D/SEMQuad instances, and performs one adaptation epoch:

  1. indicator%Estimate on the current solution (driving variable: pressure, ivar=3);
  2. flag halo expansion — grow refine flags to face-neighbours of flagged leaves (via FaceNeighbor) so the moving wavefront cannot outrun the refined band between epochs; this is the "neighbour smoothing" already anticipated in §4;
  3. snapshot old leaves → AdaptFromFlags → Balance2to1 → BuildTransferPlan → EmitMesh;
  4. new SEMQuad geometry (Init + GenerateFromMesh), free the old buffer;
  5. save model scalars (t, rho0, integrator choice, flags), Free + Init the model on the new mesh/geometry, restore scalars, apply the transferred solution, UpdateDevice.

Step 5 works on GPU today because a fresh Init allocates correctly sized device arrays; host-side transfer with an UpdateHost/UpdateDevice round-trip per epoch is acceptable at demo cadence (Stage-6 device-side transfer remains the later optimization). Background fields c (var 4) and rho0 (var 5) ride the same prolong/restrict — exact for uniform media.

5.4 Gap 3 — Time-step control#

  • Now (required): level-based global dt. For a quadtree, the fine-level element scale is exactly h_root / 2^maxLevel, so dt_epoch = dt_base / 2^maxLevel with dt_base chosen for the base mesh by the usual explicit-DG bound dt ≈ C·h/(c·N²). Deterministic, free, and no geometry reduction is needed.
  • Later (optional): local time stepping (LTS) — leaves at level ℓ subcycle with dt/2^ℓ. This touches the RK update and requires time-interpolated interface/mortar data between levels, i.e. exactly the time-integration and flux-exchange machinery that is frozen by policy; it needs its own design + review round (call it Stage 7). Cost analysis for this demo says it is not needed to start: with the refined band confined to the wavefront annulus, a global fine dt costs nElem_total × fine-step-count, and most elements are coarse and cheap; LTS buys roughly 2^maxLevel× on the coarse bulk — worth having, not blocking.

5.5 Gap 4 — The example and its CI-scale test — implemented#

Status: implemented as examples/linear_euler2d_amr_ultrasound_pointsource.f90 (water, c₀ = 1500 m/s, f₀ ≈ 100 kHz, 16×16 base at N = 7, level-3 cap, radiation boundaries). The example is registered as a CI test at a 6-epoch (30 µs) default; SELF_AMR_ULTRASOUND_EPOCHS=60 extends it to a full-domain movie run. The generic AMR-loop mechanics are separately covered by test/lineareuler2d_amr_soundwave.f90 (§2.11).

examples/linear_euler2d_amr_pointsource.f90 (plus a reduced test/ variant):

  • Domain [0,1]² m via mesh%StructuredMesh, radiation BCs on all four sides; source at the centre.
  • Medium: water (c = 1500 m/s, rho0 = 1000 kg/m³) with f₀ ≈ 100 kHz → λ = 15 mm (an air / 40 kHz variant, λ ≈ 8.6 mm, also fits but needs one more refinement level).
  • Resolution: base 16×16 (h₀ = 62.5 mm), N = 7, max level 3 (h = 7.8 mm), giving ≈ 15 points per wavelength on the fine level — comfortable for wave propagation; the coarse bulk intentionally under-resolves the front so the indicator must refine to keep σ below threshold.
  • dt ≈ 3×10⁻⁸ s at level 3; an end time of ~0.3 ms (front travels 45 cm) is ~10⁴ steps — seconds-to-minutes serial CPU, trivial on one GPU.
  • Adaptation cadence: regrid every k steps with k·dt·c ≤ one fine element (k ≈ 100 at the numbers above, with the §5.3 halo providing the safety margin).
  • CI assertions: solution NaN-free; entropy finite and non-increasing (upwind flux + radiation BCs are dissipative); refinement actually occurs (forest%MaxLevel() > 0, leaf count grows) and coarsening occurs behind the front (leaf count later shrinks); Jacobian-weighted transfer conservation defect at machine precision per epoch.

5.6 Gap 5 — Visualization: pressure field + mesh skeleton — implemented#

Status: implemented as examples/linear_euler2d_amr_plot.py (h5py + numpy + matplotlib only). Each snapshot's field and geometry are interpolated from the Gauss control points to a uniform per-element grid including the element edges (barycentric Lagrange, exact for the polynomial data), rendered as a filled pressure field with the element-outline wireframe overlaid, one PNG per snapshot plus an MP4 when ffmpeg is available.

A companion examples/linear_euler2d_amr_plot.py (pyself + matplotlib/pyvista) that, per snapshot: renders the pressure field from /controlgrid/solution and overlays the element wireframe traced from each element's four edges in /controlgrid/geometry, then assembles PNG frames into a movie. Because each file carries its own geometry, frames with different element counts need no special handling. This is the artifact that shows the refinement band tracking the expanding wavefront.

5.7 Deferred / follow-on#

  • Time-dependent transducer source. A true point forcing (e.g. Ricker wavelet at f₀) rather than an initial pulse: source2d currently has no access to position or time, so this needs a localized-forcing hook plus per-epoch source relocation (the containing element changes identity on regrid). Physically nicer (continuous-wave and pulse-train experiments); not required for the first demo.
  • LTS (Stage 7) as scoped in §5.4, and device-side transfer (Stage 6).
  • Storage compaction of orphaned forest nodes on long runs (§4, Stage 2 "next").

5.8 Suggested PR sequence#

PR Content Depends on
1 Transfer plan (old→new leaf map, multi-level prolong) + tests —
2 AMR controller (halo flags, regrid orchestration, level-based dt) + adapt-epoch soundwave test 1
3 Ultrasound example, CI-scale test, plotting script, docs 2
4 GPU epoch test; device-side transfer (Stage 6a) and amortized capacity (Stage 6b) - both done, profiling justified them 2
5 (design first) LTS; time-dependent point forcing 3

6. Three-dimensional AMR (octrees + face mortars)#

Everything above describes the 2-D implementation, which came first and remains the more thoroughly exercised of the two. The 3-D stack is a deliberate transcription of it - octrees in place of quadtrees, 2:1 face mortars in place of edge mortars, eight children in place of four - and is implemented and validated end to end. This section records what is the same, what is genuinely different, and what is not done yet.

The design record is AMR (3D) Design; it is a delta document, and where a 3-D module is a mechanical transcription the 2-D design remains the authoritative rationale.

6.1 Status#

Component 2-D 3-D
Modal-decay indicator (CPU + GPU) Implemented Implemented (SELF_RefinementIndicator_3D)
h-refinement primitives, uniform refinement Implemented Implemented (SELF_RefinementPrimitives_3D)
Adaptive forest, level tracking, 2:1 balancing Quadtree Octree (SELF_OctreeMesh_3D)
Nonconforming mesh emission Implemented Implemented (face mortars)
Transfer plan + solution transfer Implemented Implemented (SELF_TransferPlan_3D)
Model regrid + AMR controller Implemented Implemented (SELF_AMRController_3D)
MPI re-partitioning (replicated forest, allgathered migration) Implemented (v1) Implemented (v1)
Device-side solution transfer Implemented (Stage 6a) Implemented (§6.4)
Amortized high-water-mark storage Implemented (Stage 6b) Inherited via the shared data classes
Geometry reuse across an epoch Implemented (Stage 6c) Implemented (SELF_AMR_GEOM_* switches)
Example ultrasound point source linear_euler3d_amr_spherical_soundwave
Coarsen-wake regression Implemented Not yet

The pilot model is LinearEuler3D. Note it carries nvar = 6: u, v, w, P are advanced in time, while c and rho0 are spatially varying but time-constant background fields. They are still solution variables, so the transfer must carry all six - a transfer that moved only the prognostic variables would silently blank the medium on every newly created element.

6.2 What is genuinely different from 2-D#

  • Eight children, not four. transferAxc/Ayc/Azc(1:8) map an octant to its (x,y,z) half, in CGNS corner order: children 1-4 walk the x/y quadrants of the lower-z layer, 5-8 the same quadrants of the upper-z layer.
  • Triple tensor products. Prolongation is U_child = (R_kx ⊗ R_ky ⊗ R_kz) U_parent; restriction is U_parent = Σ_c (P_kx ⊗ P_ky ⊗ P_kz) U_child(c). Each 1-D P_k carries the half-interval Jacobian, so the triple product carries 1/8 rather than 1/4, which is what makes the 3-D restriction conservative. The identities are the same: Restrict(Prolong(u)) = u, and Σ w³ u_parent = (1/8) Σ_c Σ w³ u_child.
  • Cost scales far more steeply. A level of refinement multiplies the element count by 8 rather than 4, and each element carries (N+1)³ nodes rather than (N+1)². Adaptation is therefore a much larger share of a 3-D epoch than of a 2-D one at comparable settings: in the maxLevel-3 case measured in §6.4, adaptation is 89% of the epoch loop before this change and 88% after, against roughly 50% for the 2-D ultrasound case.
  • Balancing is face-based. The forest balances across faces and tolerates 2-level edge and corner jumps, because DG face mortars carry all interface data and hanging edges/corners carry none. EmitMesh guards MaxLevelJump() <= 1.

6.3 The transfer protocol and where the solution lives#

The plan (SELF_TransferPlan_3D) is built on the host after the last forest mutation, and classifies every new leaf as COPY, PROLONG or RESTRICT, with a depth and an octant path for any further descent. It is applied around the regrid as a three-step, type-bound protocol on the model, so the backend split that already selects Regrid selects the transfer too:

call model%StageSolutionForTransfer()   ! preserve the field; Regrid may then free it
call model%Regrid(newMesh,newGeom)
call model%ApplyTransferPlan(plan,interp,eFirst,eLast)

Staging exists because Regrid reallocates (or resizes) the storage the solution lives in. The portable implementation stages into a host array whose lifetime is exactly stage-to-apply; applying without a preceding stage is a hard error, pinned by test/dgmodel3d_guard_transfer_unstaged.f90.

ApplyTransferPlan takes an optional uGlobal, the multi-rank escape hatch. The forest is rank-replicated and migration is gather-then-slice, so on more than one rank the controller allgathers the old field on the host and passes it through uGlobal; each rank then fills exactly its own new element range and elements that changed ranks are migrated by construction. That allgather is inherently a host operation, so the multi-rank branch stays on the portable host transfer on every backend, exactly as in 2-D.

Where the solution lives after Adapt. On a single-rank GPU build the transferred solution is left on the device, so solution%interior (the host mirror) is stale. This matches the rest of the time loop, where the device is authoritative and a caller wanting host data calls UpdateHost() first, as Write_DGModel3D_t does. Before the device transfer existed the mirror happened to be fresh; do not rely on that. On CPU builds the two are the same storage.

6.4 Device-side transfer (implemented)#

Until this landed, a single-rank GPU adaptation epoch moved the entire solution across the host link twice - UpdateHost before the regrid, UpdateDevice after - and ran the triple tensor-product interpolation on the CPU in between. Both steps are now overridden on the GPU backend: staging is a device-to-device copy into a model-owned buffer, and the plan is applied by TransferSolution_3D_gpu (src/gpu/SELF_SolutionTransfer.cpp).

Kernel. One workgroup per new element with (N+1)² threads, where thread (a,b) owns one (a,b) pencil and loops the free index through each of the three directional passes - the same decomposition the 3-D modal indicator uses. An (N+1)³ thread block is not an option: it is 4096 threads at N = 15, past the block limit, and it would push the working buffers into per-thread scratch. The three (N+1)³ working buffers are static __shared__ arrays sized AMR3D_MAXNP³ with AMR3D_MAXNP = 12, i.e. 41 kB reserved per block at any degree; the Fortran caller guards N+1 <= 12 so a larger degree fails loudly rather than overrunning.

A dynamic-shared variant sized to the degree actually in use (12 kB at N = 7, the idiom SELF_MatrixMultiply.cpp uses) was implemented and measured, and was 1-3% slower on a B300 across three runs at two refinement depths - the SM shared budget there is large enough that 41 kB per block does not gate occupancy, so the launch-time sizing was pure overhead. It was therefore not adopted. The argument for it has not gone away on a 64 kB-LDS AMD device, where the static footprint admits one workgroup per CU rather than five; that case has not been measured, and since the alternative on every platform is the host round trip this kernel replaces, the static form cannot regress anything as it stands.

One deliberate divergence from the host. The host descent calls ProlongToChildren, which forms all eight children and discards seven; the kernel applies only the operator triple of the child actually on the recorded path - an eighth of the work per descent step. What is dropped is the seven discarded children, which are independent computations writing disjoint slices, and not any term of the retained value: each contraction sums the same products against the same mortar column in the same ascending index order as the host loop, so the reduction order is preserved. Device and host nonetheless agree to round-off rather than bitwise, because the device compiler contracts these multiply-accumulates into FMAs.

test/solution_transfer_3d_device.f90 pins the result value by value against the portable ApplyTransferPlanRange, on an epoch built to contain all three transfer kinds and prolongation depths 0-2, with a non-polynomial field that differs in every variable. That is the check the AMR regressions cannot make: conservation and entropy non-growth are invariants a kernel with a transposed direction could still satisfy.

Measurements. One B300 (sm_103, CUDA 13, fp64), linear_euler3d_amr_spherical_soundwave, 20 epochs, mean of three runs, before and after differing only in the three files that implement the transfer. tAdapt is reported per adapting epoch, because epochs whose leaf set is unchanged cost only an indicator pass and mixing the two makes the metric depend on how many epochs happened to adapt:

maxLevel elements adapting epochs before after
1 512 1 119.5 ms 101.6 ms -15.0%
2 4,096 10 571.6 ms 516.3 ms -9.7%
3 20,784 20 2182.6 ms 1959.6 ms -10.2%

Time integration is unchanged, as it must be - nothing on this path runs inside ForwardStep. tForwardStep over the same runs was 0.672 -> 0.615 s, 2.152 -> 2.116 s and 5.199 -> 5.176 s, all inside the run-to-run spread (which is under 1% on tAdapt).

Where the saving actually comes from, which is not where you would guess. It is tempting to credit the two eliminated full-field host/device copies, and the 2-D Stage 6a note above leads with them. Do the arithmetic for the maxLevel-3 case: 20,784 elements x 125 nodes x 6 variables x 8 B is a 125 MB field, so the round trip is 250 MB, which at PCIe Gen5 rates is roughly 5 ms - 0.2% of a 2183 ms adaptation. The copies are not the story at this scale.

The 223 ms actually saved is the host-side interpolation: the portable ApplyTransferPlanRange runs the triple tensor-product contractions on one CPU core for every new element, and that work moves onto the GPU. The kernel also does an eighth of the descent work the host does, because it prolongs onto the child on the recorded path instead of forming all eight and discarding seven. Both effects scale with element count, which is why the saving holds at roughly 10% of adaptation across a 40x range in mesh size while the copy time would have faded to nothing.

The remaining ~90% of an adaptation is geometry generation, regrid and the indicator - untouched by this change, and where the next 3-D optimization work belongs.

6.5 What is not done in 3-D#

  • A coarsening regression. 2-D has lineareuler2d_amr_coarsen_wake, which pins the amplitude gate's ability to release the mesh behind a passing front (§2.1.1). There is no 3-D analogue, so 3-D coarsening is exercised only incidentally, by the soundwave regression and by the transfer tests' hand-built epochs.
  • Point-to-point migration (Stage-5 v2). Multi-rank runs still allgather the old field on the host. Until that changes, the device transfer cannot be used on more than one rank, because the migration it would have to feed from is a host operation.
  • A 3-D visualization script. 2-D ships examples/linear_euler2d_amr_plot.py. The 3-D example writes the same self-describing HDF5 snapshots (solution + geometry per file, so the changing mesh needs no special handling downstream) but no renderer is provided.

7. References#

  • P.-O. Persson and J. Peraire, Sub-cell shock capturing for discontinuous Galerkin methods, AIAA 2006-112 (2006).
  • S. Hennemann, A. M. Rueda-Ramírez, F. J. Hindenlang, and G. J. Gassner, A provably entropy stable subcell shock capturing approach for high order split form DG for the compressible Euler equations, J. Comput. Phys. 426, 109935 (2021).
  • C. Mavriplis, Adaptive mesh strategies for the spectral element method, Comput. Methods Appl. Mech. Engrg. 116 (1994) 77–86.
  • M. Schlottke-Lakemper et al., Trixi.jl: Adaptive high-order numerical simulations of hyperbolic PDEs in Julia, and the Trixi.jl IndicatorHennemannGassner / IndicatorLöhner documentation.