Skip to content

Adaptive Mesh Refinement (AMR)#

This page describes SELF's approach to adaptive mesh refinement for 2-D discontinuous Galerkin spectral element models. It documents the refinement trigger that is implemented today and lays out a staged design for the remaining mesh-adaptation machinery so that the invasive parts (dynamic mesh mutation, solution transfer, MPI re-partitioning, and GPU re-allocation) can be reviewed and landed incrementally.

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.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)

! solution : the model's MappedScalar2D (or any Scalar2D); ivar selects the driving
! variable, or SELF_AMR_ALLVARS (=0) reduces over all variables (most conservative).
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)
n = amr%CountFlagged(SELF_AMR_REFINE)

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. 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.