Procedures

Resolve the communicator a new decomposition will run on and register it with the process-wide lifecycle bookkeeping.

Perform one adaptation epoch on the model (see the module documentation). On return, adapted reports whether the mesh changed; when it did, the model is already rebound to the new mesh with the solution transferred (conservatively), and the caller should re-evaluate its time step (RecommendedTimeStep) before the next ForwardStep. When the leaf set is unchanged the model is untouched.

Mutate the forest from a per-leaf flag array (indexed over the current leaves in this%leaf order, e.g. the flag array produced by SELF_RefinementIndicator_2D):

Register boundary conditions for the advection-diffusion model. Hyperbolic: mirror (no-normal-flow) wall condition on the solution. Parabolic: no-stress condition that zeros the normal gradient at the boundary, giving zero diffusive flux through walls. This is stable with the Bassi-Rebay (BR1) method.

Register boundary conditions for the advection-diffusion model. Hyperbolic: mirror (no-normal-flow) wall condition on the solution. Parabolic: no-stress condition that zeros the normal gradient at the boundary, giving zero diffusive flux through walls. This is stable with the Bassi-Rebay (BR1) method.

Register the (CPU) no-normal-flow and radiation boundary conditions. GPU builds call this parent and then overwrite both registrations with the device kernels in AdditionalInit_LinearEuler2D.

Register the (CPU) radiation boundary condition. GPU builds overwrite this registration with the device kernel in AdditionalInit_LinearEuler3D.

Adds a pressure-balanced warm bubble perturbation.

Adds a pressure-balanced warm bubble perturbation.

Allgather an integer array with perElem entries per element from the decomposition's contiguous rank-local element ranges into the global element ordering.

Allgather a real(prec) array with perElem entries per element from the decomposition's contiguous rank-local element ranges into the global element ordering.

\addtogroup SELF_SupportRoutines @{ \fn AlmostEqual Compares two floating point numbers and determines if they are equal (to machine precision).

Execute a transfer plan on nodal element data: uNew(:,:,li,:) receives old element data copied, prolonged (exact interpolation), or restricted (conservative L2 projection) according to plan entry li. interp must be the solution interpolant the data lives on (its mortar operators drive the transfer). Runs once per adaptation epoch - not a per-time-step hot path - so clarity is preferred over fused loops here.

Apply the transfer plan on the device (Stage 6a), writing solution%interior_gpu directly and moving no solution data across the PCIe/xGMI link.

Transfer the staged pre-regrid solution onto the regridded mesh through plan, filling the rank-local element range [eFirst,eLast] of the new solution.

Execute the contiguous sub-range eFirst..eLast of a transfer plan: uNew(:,:,k,:) receives the data of new element eFirst+k-1. uOld is the full (global) old field; the output is only the requested slice. This is the decomposed-mesh (AMR Stage 5) entry point: each rank passes its own contiguous range of the new element ordering and a gathered global old solution, and fills exactly its rank-local storage.

Enforce the 2:1 balance condition: no leaf face may separate elements differing by more than one refinement level. Iterates to a fixed point - in each sweep, any leaf whose equal-or-larger neighbour is a leaf two or more levels coarser triggers refinement of that coarser neighbour; refinement can ripple, so sweeps repeat until nothing changes. The leaf set is rebuilt on return.

LLF Riemann flux on GPU — fully device-resident.

LLF Riemann flux on GPU — fully device-resident.

LMARS interface flux on GPU.

LMARS (Low-Mach Approximate Riemann Solver, Chen et al. 2013) interface flux. No hydrostatic pressure split: gravity is folded into SourceMethod via the Souza non-conservative form using the geopotential carried in the state vector (variable index 5).

LMARS interface flux on GPU. No hydrostatic pressure split: gravity is handled by the Souza non-conservative source term using the geopotential carried as state variable index 6.

LMARS (Low-Mach Approximate Riemann Solver, Chen et al. 2013) interface flux. No hydrostatic pressure split: gravity is folded into SourceMethod via the Souza non-conservative form using the geopotential carried in the state vector (variable index 6).

Fill mesh%nodeCoords(:,:,:,e) for one element. Place the four corner nodes from the file's node table, then perform linear (nGeo=1) placement or transfinite interpolation (nGeo>1) using the (possibly curved) edge curves. For straight sides the edge curve is filled in here by linear interpolation between corners at Chebyshev-Gauss-Lobatto parametric coordinates.

Fill mesh%nodeCoords(:,:,:,:,e) for one hexahedral element by transfinite (Coons) interpolation of its six face grids. Faces flagged in the mesh file use the file's face-point grids; the remaining faces are bilinear patches of their four corner nodes evaluated at Chebyshev-Gauss-Lobatto parametric coordinates. Edge curves are extracted from the face grids (preferring a flagged face when an edge borders one flagged and one bilinear face) and the standard Boolean-sum formula x = Pxi + Peta + Pzeta - PxiPeta - PxiPzeta - PetaPzeta + PxiPeta*Pzeta combines face, edge, and corner contributions. For an element with all-straight faces this reduces to trilinear interpolation of the eight corners.

Axis-aligned bounding box of each element. The box is taken over the element's boundary nodes (geometry%x%boundary), which are evaluated at the element edges (reference coordinate = +/-1) and therefore reach the true element boundary and corners. Using the interior nodes alone is not sufficient: for Gauss quadrature the interior nodes are strictly inside [-1,1], so a box built from them under-covers the element by a fixed fraction of its size and rejects points lying on element edges/corners. The boundary nodes give the correct extent for both Gauss and Gauss-Lobatto node sets.

Axis-aligned bounding box of each element, taken over the boundary nodes (geometry%x%boundary) which are evaluated on the element faces (reference coordinate = +/-1) and so reach the true element boundary, edges and corners. See BuildElementBBoxes_2D for why the interior nodes alone are insufficient (Gauss interior nodes under-cover the element).

Fill newGeom for the emitted mesh, reusing the previous epoch's geometry for every element that did not change and generating only the rest (AMR Stage 6c). nReused reports how many elements were copied rather than computed.

Build the aggregated halo exchange tables for the mesh this domain decomposition belongs to.

Build the nodal->modal transform Pmodal for a 1-D degree-N interpolant whose nodes are controlPoints(1:N+1). The transform is the exact inverse of the L2-normalized Legendre Vandermonde V(i,p) = Ltilde_{p-1}(x_i); Pmodal(ii,p) stores V^{-1}(p,ii) so that uhat(p) = sum_ii Pmodal(ii,p) * u(ii). The Vandermonde inverse is formed in double precision for conditioning and cast back to the working precision.

Build the old->new transfer plan for the adaptation epoch that took the forest from the leaf configuration (nOld, oldLeaf) - a snapshot of (forest%nLeaves, forest%leaf) taken before mutating - to its current leaf configuration. See the module documentation for the allowed mutations within one epoch. The plan is valid until the forest is mutated again.

Base method for calculating entropy of a model When this method is not overridden, the entropy is simply set to 0.0. When you develop a model built on top of this abstract class or one of its children, it is recommended that you define a convex mathematical entropy function that is used as a measure of the model stability.

GPU implementation of the EC-DG 2-D tendency.

Computes du/dt = source - EC-DG flux divergence.

GPU-resident tendency for ESAtmo2D.

ESAtmo2D tendency = EC inviscid pipeline (parent) + optional constant-coefficient Laplacian diffusion (BR1 weak-form DG).

GPU-resident tendency for ESAtmo3D. The inviscid pipeline is identical to ECDGModel3D's GPU CalculateTendency; if either nu or kappa is positive, the constant-coefficient Laplacian divergence (BR1 weak-form) is then accumulated into fluxDivergence before forming dSdt.

ESAtmo3D tendency = EC inviscid pipeline (parent) + optional constant-coefficient Laplacian diffusion (BR1 weak-form DG).

Convert a (signed, possibly out-of-range) real cell index to an integer cell index clamped to [0, nCells-1].

\addtogroup SELF_SupportRoutines @{ \fn CompareArray Compares to INTEGER arrays and determines if they are identical.

Copy whole-element geometry blocks from src into myGeom: element srcIdx(k) of src becomes element dstIdx(k) of myGeom, for k = 1..n (AMR Stage 6c).

Returns the two corner-node IDs delimiting a given local side using SELF's 2D side convention: Side 1 South = [CN1, CN2] Side 2 East = [CN2, CN3] Side 3 North = [CN4, CN3] Side 4 West = [CN1, CN4]

Count the rank-local elements currently carrying the requested flag value (SELF_AMR_REFINE / SELF_AMR_KEEP / SELF_AMR_COARSEN). Provided as a convenience for drivers and tests; a global count across MPI ranks is the caller's responsibility.

GPU-resident fill of diffFlux%boundaryNormal with the SIPG-stabilised BR1 diffusive flux: f = -coeff(avg_grad . n)nmag + tau(uL - uR)nmag

Fill diffFlux%boundaryNormal with the SIPG-stabilised BR1 flux:

GPU-resident fill of diffFlux%boundaryNormal with the SIPG-stabilised BR1 diffusive flux: f = -coeff(avg_grad . n)nmag + tau(uL - uR)nmag tau = eta_penaltycoeff(N+1)^2/length_scale, computed on the host.

Fill diffFlux%boundaryNormal with the SIPG-stabilised BR1 flux:

GPU-resident fill of diffFlux%interior with the constant-coefficient Laplacian flux F_d(iVar) = -coeff(iVar) * d(s_iVar)/dx_d.

Fill diffFlux%interior with the constant-coefficient Laplacian flux at every interior node:

GPU-resident fill of diffFlux%interior with the constant-coefficient Laplacian flux F_d(iVar) = -coeff(iVar) * d(s_iVar)/dx_d, where coeff is 0 for rho, nu for momentum, kappa for rho*theta.

Fill diffFlux%interior with the constant-coefficient Laplacian flux at every interior node:

GPU override of DiracDelta_2D_Points_t. When the scalar is a device-resident MappedScalar2D, the kernel writes scalar%interior_gpu directly (host array left stale; caller may UpdateHost). For non-GPU scalar subtypes, defers to the inherited host implementation.

Scatter a discrete Dirac delta of unit strength (S = 1) onto a 2D MappedScalar, one variable per stored point. Variable p of scalar receives the delta associated with point p.

GPU override of DiracDelta_3D_Points_t. See DiracDelta_2D_Points.

Scatter a discrete Dirac delta of unit strength (S = 1) onto a 3D MappedScalar, one variable per stored point. See DiracDelta_2D_Points_t for the full specification; this is the direct 3D analogue.

GPU implementation of the reference-element split-form divergence. df must be a device pointer to a scalar field of size (N+1)(N+1)nElem*nVar.

Computes the split-form (two-point) divergence of a 2-D vector field in the reference element (computational coordinates).

GPU implementation of the reference-element split-form divergence. df must be a device pointer to a scalar field of size (N+1)^3 * nElem * nVar.

Computes the split-form (two-point) divergence of a 3-D vector field in the reference element (computational coordinates).

Create a six-element mesh with two 2:1 mortar interfaces, used to validate the mortar machinery in configurations the SimpleMortarMesh cannot reach:

Build outMesh (a conforming-or-mortar Mesh2D_t) from a 2:1-balanced forest. baseMesh is the mesh the forest was initialised from (supplies BC metadata and the communicator; on nRanks > 1 the forest must be rank-replicated so every rank emits identical global tables). The forest must already be balanced (MaxLevelJump <= 1); EmitMesh does not mutate it.

Grow the node arrays (amortized doubling) so at least need nodes fit.

High-water-mark device allocation (AMR Stage 6b). Grows ptr to hold neededBytes, reusing the existing allocation when it already does, and records the capacity in allocBytes. Contents are not preserved across a growth - every caller rewrites the buffer before reading it.

Grow the device buffers to hold the current logical arrays, reusing existing allocations when the bytes already fit (AMR Stage 6b).

Grow each device buffer to hold the current logical array, reusing the existing allocation when the bytes already fit (AMR Stage 6b). A device pointer carries no shape, so unlike the host side no remapping is needed - only a byte-capacity check.

Grow each device buffer to hold the current logical array, reusing the existing allocation when the bytes already fit (AMR Stage 6b/6c). A device pointer carries no shape, so only byte capacity has to be tracked.

Grow each device buffer to hold the current logical array, reusing the existing allocation when the bytes already fit (AMR Stage 6b). A device pointer carries no shape, so only byte capacity has to be tracked.

Guarantee that pool has at least needed elements, growing it by poolGrowth when it must reallocate. Existing contents are NOT preserved: every caller re-establishes the array contents after resizing (the solution through the AMR transfer, the geometry through GenerateFromMesh, everything else by being written before it is read), so copying the old data would be wasted bandwidth.

Prepare the cached GenerateFromMesh scratch for a mesh with nGeo/quadrature and nElem elements (AMR Stage 6c). The nGeo -> N interpolant is built once and reused; the node coordinate staging is resized, so an adapting run stops rebuilding either one per epoch.

Quadratic entropy: eta(u) = u^2 / 2

Quadratic entropy: eta(u) = u^2 / 2

Mathematical entropy: total energy density (kinetic + internal).

Mathematical entropy: total energy density (kinetic + internal).

Acoustic energy for the PML model. The background density is the scalar this%rho0 (the PML variant does not carry a per-node background density; its variable 5 is the auxiliary phi_u, not rho0), and the sound speed is taken from s(4). This overrides the parent entropy_func, which reads rho0 from s(5).

The entropy function is the sum of kinetic and internal energy For the linear model, this is

The entropy function is the sum of kinetic and internal energy For the linear model, this is

GPU path: one device thread per element computes the modal transform and energy ratios, writes the indicator and flag on device, then copies both back to the host so the (host-side) mesh-adaptation logic can consume them directly.

Compute the per-element modal-energy indicator sigma_e and refine/keep/coarsen flag from the nodal solution field. ivar selects the driving variable in [1,solution%nVar]; passing SELF_AMR_ALLVARS (=0) reduces the indicator over all variables by taking, per element, the largest (least smooth) smoothness ratio S_e before the log10.

Evaluate a 2D MappedScalar at the cached points on the device. The caller is responsible for hipMalloc-ing values_dev with capacity nPointsnVarprec bytes. Layout is column-major (nPoints, nVar).

Evaluate a 2D MappedScalar at all located points by tensor-product Lagrange interpolation at the stored reference coordinates. Points with elements(p) == 0 receive a value of zero. When LocatePoints has cached the per-point basis at the matching polynomial degree, this routine reuses it and skips Lagrange-polynomial evaluation altogether.

Evaluate a 3D MappedScalar at all located points by tensor-product Lagrange interpolation. Points with elements(p) == 0 receive zero. When LocatePoints has cached the basis at the matching degree, this routine reuses it.

Find the equal-or-larger face neighbour of node across its local side s using the classic quadtree ascend/descend search. Returns: nbr - neighbour node id (0 = physical domain boundary). It is either a LEAF at any level <= level(node), or an INTERNAL node at exactly level(node) (meaning the shared face is subdivided on the neighbour side, i.e. finer neighbours exist). ns - the neighbour's local side that faces node. nf - the flip between the two shared edges (0 same direction, 1 reversed), inherited from the base-mesh face where the search crosses a root boundary. With this, a 2:1 hanging face is exactly "nbr is a leaf with level(nbr) = level(node)-1", and finer neighbours are exactly "nbr is internal".

Interior flux. Variables 1-4 use the parent linear Euler 2D flux; auxiliary variables 5-7 carry zero flux in both directions and are evolved purely by the PML source term.

\addtogroup SELF_SupportRoutines @{ \fn ForwardShift Shift an array integers by one index forward, moving the last index to the first.

Forward steps the model using the associated tendency procedure and time integrator

Release the forest, the indicator, and any controller-emitted mesh/geometry. The caller-owned base mesh/geometry are untouched. A model still pointing at a controller-emitted mesh must not be used after this call.

Free the 1D DG model, including GPU BC arrays.

Free the 2D DG model, including GPU BC arrays.

Free the 3D DG model, including GPU BC arrays.

Free EC Advection 2D, including GPU BC arrays.

Free EC Advection 3D, including GPU BC arrays.

Frees all memory (host and device) associated with an instance of the Lagrange class

Frees all memory (host and device) associated with an instance of the Lagrange class

Frees all memory (host and device) associated with an instance of the Lagrange_t class

Generate geometry for nElem elements directly from their mesh node coordinates (AMR Stage 6c). GenerateFromMesh is a thin wrapper over this.

Returns the node associated with the given bcid. If the bcid is not found, a null pointer is returned.

Returns the current simulation time stored in the model % t attribute

Map a global element position (gx,gy,gz), with gx in [1,nX] etc., to the tile-ordered element id used throughout the structured mesh. This is the inverse of the (i,ti) decomposition: gx = i + nxpertile*(ti-1).

Heap sort of the halo side list by ascending 64-bit key, permuting the companion arrays alongside the keys. O(n log n), in place; the relative order of equal keys is irrelevant because global side ids are unique.

GPU-accelerated mirror BC for 2D EC Advection.

Mirror boundary condition on the solution: sets the exterior state equal to the interior state. With LLF this zeroes the upwind dissipation at the wall.

Prescribed-zero boundary state for the tracer.

No-normal-flow (wall) boundary condition.

GPU-accelerated no-normal-flow BC for 2D Entropy-Stable Atmosphere. Reflects normal momentum, mirrors density, rho*theta, and Phi.

No-normal-flow boundary condition for 2D linear Euler equations. Reflects the velocity vector about the boundary normal while preserving pressure, sound speed, and background density.

GPU-accelerated no-normal-flow BC for 2D Linear Euler.

No-normal-flow BC for the PML-augmented linear Euler model. Variables 1-4 are treated identically to the parent LinearEuler2D no-normal-flow BC. Auxiliary variables 5-7 are given a zero exterior state; since they carry zero Riemann flux the exterior value is mathematically irrelevant, but zeroing it keeps the boundary state clean for diagnostics.

No-normal-flow boundary condition for 2D linear shallow water equations. Reflects the velocity vector about the boundary normal while preserving the free surface elevation.

Radiation BC: zero acoustic perturbation in the exterior state; the sound speed (variable 4) and background density (variable 5) are copied from the interior side so the Riemann solver sees a consistent c and rho0 (impedance-matched, non-reflecting outflow).

GPU-accelerated radiation BC for 2D Linear Euler.

Radiation BC for the PML-augmented linear Euler model: zero acoustic perturbation in the exterior state, sound speed copied from interior so the Riemann solver sees a consistent c, and auxiliary variables set to zero.

GPU-accelerated mirror BC for 3D EC Advection.

Mirror boundary condition on the solution: sets the exterior state equal to the interior state. With LLF this zeroes the upwind dissipation at the wall.

Prescribed-zero boundary state for the tracer.

No-normal-flow (wall) boundary condition.

GPU-accelerated no-normal-flow BC for 3D Entropy-Stable Atmosphere. Reflects normal momentum, mirrors density and rho*theta.

Radiation BC: zero acoustic perturbation in the exterior state; the sound speed (variable 5) and background density (variable 6) are copied from the interior side so the Riemann solver sees a consistent c and rho0.

Attach the controller to an initialized model. The model's current mesh becomes the forest's base mesh (level 0); its geometry interpolant drives the indicator and the solution transfer. Thresholds are the sigma = log10 modal-energy-ratio cut-offs of the refinement indicator (refineThreshold > coarsenThreshold; see SELF_RefinementIndicator_2D). ivar is the driving solution variable (or SELF_AMR_ALLVARS). maxLevel >= 0 caps the refinement depth; nHalo >= 0 sets the refine-flag halo width in elements.

Initialize the 1D DG model, then upload BC element/side arrays to GPU.

Initialize the 2D DG model, then upload BC element/side arrays to GPU.

Initialize the 3D DG model, then upload BC element/side arrays to GPU.

Initialize EC Advection 2D, then upload BC element/side arrays to GPU.

Initialize EC Advection 3D, then upload BC element/side arrays to GPU.

Initialize an instance of the Lagrange class On output, all of the attributes for the Lagrange class are allocated and values are initialized according to the number of control points, number of target points, and the types for the control and target nodes. If a GPU is available, device pointers for the Lagrange attributes are allocated and initialized.

Initialize an instance of the Lagrange class On output, all of the attributes for the Lagrange class are allocated and values are initialized according to the number of control points, number of target points, and the types for the control and target nodes. If a GPU is available, device pointers for the Lagrange attributes are allocated and initialized.

Initialize an instance of the Lagrange_t class On output, all of the attributes for the Lagrange_t class are allocated and values are initialized according to the number of control points, number of target points, and the types for the control and target nodes. If a GPU is available, device pointers for the Lagrange_t attributes are allocated and initialized.

Allocate host + device storage for nPoints in nDim dimensions. The basis-cache device buffers are allocated lazily in UpdateDevice once the polynomial degree N is known (it is set by LocatePoints).

Allocate storage for a point cloud of size nPoints in nDim dimensions. nDim must be 2 or 3.

Initialize the forest with one root per base-mesh element (all leaves at level 0). The base geometry (node coordinates) is copied so leaf geometry can be regenerated after any amount of refinement without holding a reference to the mesh. Requires a single-rank mesh (a decomposed mesh only stores its local elements); a rank-replicated forest over a decomposed base is built by gathering the global tables and calling InitGlobal.

Allocate the indicator for an interpolant of degree interp%N and nElem elements and precompute the nodal->modal transform matrix from the interpolant control points.

Allocate the interior array for a 2-D two-point vector field. The interior array has rank 6 with layout (n,i,j,nEl,nVar,idir).

Allocate the interior array for a 3-D two-point vector field. The interior array has rank 7 with layout (n,i,j,k,nEl,nVar,idir).

Build a rank-replicated forest over a decomposed base mesh: allgather the global node coordinates, side table, and material ids (by the decomposition's contiguous element ownership) and initialize the forest from the global tables. Collective over the mesh communicator; runs once, at controller initialization.

Initialize the forest directly from GLOBAL base-mesh tables (one root per global base element, all leaves at level 0). This is the initialization path for a rank-replicated forest over a decomposed base mesh (AMR Stage 5): every rank passes the same gathered tables and holds an identical forest. rootNbr carries global element ids (0 = physical boundary), rootNbrSide/rootFlip decode the base sideInfo(4) pairing, rootBC the base boundary-condition id per side, and rootMaterial the base material id per element.

Invert the n-by-n matrix A by Gauss-Jordan elimination with partial pivoting, in double precision. Used once at initialization for the small (N+1) Legendre Vandermonde; A is a well-conditioned Vandermonde in an orthonormal basis, so a direct solve is appropriate.

Physical geometry-node coordinates of leaf leafIndex, produced by repeated exact isoparametric subdivision of its root element along the quadtree path. geomInterp must be a degree-nGeo Lagrange interpolant on the mesh's geometry (quadrature) nodes - the same interpolant SELF_MeshRefinement_2D builds. Level 0 leaves return the root geometry directly.

Run the host spatial-hash + Newton search (inherited from Points_t), then mirror per-point state to the device.

Locate each stored physical point inside the 2D SEMQuad geometry. On exit, elements(p) holds the (rank-local) element id and coordinates(p,:) holds the (s,t) reference coordinates, both for points that resolve. Points that resolve to no element are marked with elements(p) = 0.

Locate each stored physical point inside the 3D SEMHex geometry. On exit, elements(p) and coordinates(p,1:3) hold the located element id and reference (s,t,u) for points that resolve; elements(p) = 0 otherwise.

Numerically stable logarithmic mean (Ismail-Roe 2009 / Ranocha 2018).

Numerically stable logarithmic mean (Ismail-Roe 2009 / Ranocha 2018).

Size the backing pools for (Np,nVar,nElem) and remap the public arrays onto them at the exact logical shape. Bounds-remapping a rank-1 contiguous target to a rank-4 pointer is what keeps every extent and stride identical to a plain allocate while allowing the pool underneath to be larger and reused. See SELF_DataPool.

Size the backing pools for (Np,nVar,nElem) and remap the public arrays onto them at the exact logical shape. See SELF_DataPool for why this is done with pointer remapping rather than by over-allocating the element dimension.

Size the backing pools for (Np,nVar,nElem) and remap the public arrays onto them at the exact logical shape. See SELF_DataPool for why this is done with pointer remapping rather than by over-allocating the element dimension.

Scan the mesh boundary condition IDs and populate the elements/sides arrays for each registered boundary condition.

Scan the mesh sideInfo and populate the elements/sides arrays for each registered boundary condition.

Scan the mesh sideInfo and populate the elements/sides arrays for each registered boundary condition.

Computes the divergence of a 2-D vector using the weak form On input, the attribute of the vector is assigned and the attribute is set to the physical directions of the vector. This method will project the vector onto the contravariant basis vectors.

Computes the divergence of a 3-D vector using the weak form On input, the attribute of the vector is assigned and the attribute is set to the physical directions of the vector. This method will project the vector onto the contravariant basis vectors.

Computes the divergence of a 3-D vector using the weak form On input, the attribute of the vector is assigned and the attribute is set to the physical directions of the vector. This method will project the vector onto the contravariant basis vectors.

Calculates the gradient of a function using the weak form of the gradient and the average boundary state. This method will compute the average boundary state from the and attributes of

Calculates the gradient of a function using the weak form of the gradient and the average boundary state. This method will compute the average boundary state from the and attributes of

GPU implementation of the physical-space split-form divergence for a 2-D two-point vector on a curvilinear mesh.

Computes the physical-space divergence of a 2-D split-form vector field on a curvilinear mesh.

GPU implementation of the physical-space split-form divergence for a 3-D two-point vector on a curvilinear mesh.

Computes the physical-space divergence of a 3-D split-form vector field on a curvilinear mesh.

Calculates the gradient of a function using the strong form of the gradient in mapped coordinates.

Calculates the gradient of a function using the strong form of the gradient in mapped coordinates.

Calculates the gradient of a function using the strong form of the gradient in mapped coordinates.

Calculates the gradient of a function using the strong form of the gradient in mapped coordinates.

Highest refinement level among the active leaves.

Largest refinement-level difference across any leaf face (0 on a conforming or uniformly refined forest, 1 on a 2:1-balanced adaptive forest). Because FaceNeighbor returns the equal-or-larger neighbour, every level difference is observed from the finer leaf as a coarser leaf neighbour; internal (finer) neighbours contribute nothing from this side.

GPU implementation of the mortar exchange (see the base class for the algorithm) : traces are staged, reoriented, restricted, and projected in device memory with the SELF_Mortar kernels.

Fills the extBoundary attribute on all sides that participate in a 2:1 nonconforming (mortar) interface.

GPU implementation of the vector mortar exchange; the kernels treat the (variable, direction) pairs as 2*nvar independent trace lines.

Fills the extBoundary attribute on all sides participating in a 2:1 nonconforming (mortar) interface; vector analogue of the scalar MortarExchange (see MappedScalar2D_t for the algorithm description).

GPU implementation of MortarFluxCollect (see the base class for the algorithm and conservation statement). Stages the small sides' boundaryNormal traces in the mortar buffer, then overwrites the big side's integrand on device.

Replaces the big-side boundaryNormal trace on each mortar interface with the L2 projection of the two small sides' boundaryNormal traces.

Post the aggregated halo exchange: one MPI_Irecv/MPI_Isend pair per neighboring rank, carrying every (side,variable,component) boundary trace shared with that rank in a single packed device buffer. The boundary array is laid out with the component index outermost, so the pack/unpack kernels treat the vector as 2*nvar scalar variables. Packed buffers are allocated on first use; the shared side tables are built by SideExchange before this is called.

Post the aggregated halo exchange: one MPI_Irecv/MPI_Isend pair per neighboring rank, carrying every (side,variable,component) boundary trace shared with that rank in a single packed device buffer. The boundary array is laid out with the component index outermost, so the pack/unpack kernels treat the vector as 3*nvar scalar variables. Packed buffers are allocated on first use; the shared side tables are built by SideExchange before this is called.

GPU-resident analogue of the base-class mortar message posting; messages are posted on device memory (GPU-aware MPI), following MPIExchangeAsync.

Posts the point-to-point messages required for mortar interfaces whose big and small elements reside on different ranks. The big-side rank sends its edge trace to each remote small-side rank and receives each remote small-side trace; tags are built from the sub-edge global side ids, following the conforming SideExchange tag convention.

GPU-resident mortar message posting for vector data; messages are posted on device memory (GPU-aware MPI), one per variable and physical direction.

Vector analogue of the scalar mortar exchange message posting; each physical direction of each variable is exchanged as its own message.

One-directional messages for MortarFluxCollect on device memory : each remote small side sends its boundaryNormal trace to the big side's rank.

Posts the one-directional messages for MortarFluxCollect : each remote small side sends its boundaryNormal trace to the big side's rank.

Newton inverse-map for the 2D SEM element iEl: solve X(xi) = xTarget for the reference coordinate xi in R^2, where X is the high-order interpolant. The Jacobian dX/dxi is taken from geometry%dxds (the covariant basis tensor), interpolated to xi via Lagrange basis.

Newton inverse-map for the 3D SEM element iEl. See NewtonInverse_2D.

Select the geometry buffer to fill this epoch: whichever of the two is not currently active, so the previous epoch's geometry stays intact and readable (AMR Stage 6c).

L2-normalized Legendre polynomial Ltilde_p(x) = L_p(x)*sqrt((2p+1)/2) on [-1,1], evaluated with the standard three-term recurrence. The normalization gives int_{-1}^{1} Ltilde_p Ltilde_q dx = delta_pq.

No-stress boundary condition for the BR1 diffusive flux. Reflects the interior gradient so that the normal component of the averaged gradient vanishes at the boundary: sigma_ext = sigma_int - 2 (sigma_int . nhat) nhat This gives zero diffusive flux through the wall and is unconditionally stable for Bassi-Rebay.

Parabolic boundary condition: zero diffusive flux normal to the wall (no-stress for momentum, no-heat-flux for rho*theta).

GPU-accelerated parabolic no-stress / no-heat-flux BC for 2D Entropy-Stable Atmosphere.

No-stress boundary condition for the BR1 diffusive flux. Reflects the interior gradient so that the normal component of the averaged gradient vanishes at the boundary: sigma_ext = sigma_int - 2 (sigma_int . nhat) nhat This gives zero diffusive flux through the wall and is unconditionally stable for Bassi-Rebay.

Parabolic boundary condition: zero diffusive flux normal to the wall (no-stress for momentum, no-heat-flux for rho*theta).

GPU-accelerated parabolic no-stress / no-heat-flux BC for 3D Entropy-Stable Atmosphere. Reflects the normal component of the solution gradient at every wall node so that BR1 averaging gives avgGrad . n = 0 (zero diffusive flux through the wall) for every variable.

Populate the elements and sides arrays for a registered boundary condition. Called after scanning the mesh to determine which faces belong to each bcid.

PostStepHook is a template routine invoked by the time integrators once immediately after each completed time step (after this%t has advanced by this%dt). The default implementation is a no-op.

PreStepHook is a template routine invoked by the time integrators once immediately before each time step is taken (before any Runge-Kutta stages of that step). The default implementation is a no-op.

PreTendencyHook is a template routine that is used to house any additional calculations that you want to execute at the beginning of the tendency calculation routine. This default PreTendencyHook simply returns back to the caller without executing any instructions

Prolong (interpolate) a parent element's nodal solution onto its four children. uChildren(:,:,:,c) is the field on child c (ordering 1=SW,2=SE,3=NE,4=NW). Exact for the parent's degree-N polynomial representation.

.true. if child c's side s is interior to its parent (its neighbour across s is a sibling).

The local side directly across an element from side s.

Sibling child index obtained by reflecting c across side s (swap the y-half for the horizontal faces S/N, the x-half for the vertical faces E/W).

Sub-position (1 or 2, in the positive direction of side s) of child c along side s; 0 if c does not touch side s.

Reader for HOHQMesh text mesh files in the ISM and ISM-MM formats. The format is auto-detected from the first line: * Line equal to "ISM-MM" (trimmed) => ISM-MM with per-element material name strings and a 4-int count line that includes an unused nEdges field (the ISM-MM writer in HOHQMesh does NOT emit an edge block). * Anything else is treated as plain ISM: the first line is itself the count line "nNodes nElems polyOrder" and there are no per-element material names.

Reader for HOHQMesh 3-D (hexahedral) text mesh files in the ISM and ISM-MM formats, as written by HOHQMesh's WriteISMHexMeshFile (Source/3DSource/Mesh3DOutputMethods.f90). Unlike the 2-D writer, the 3-D writer emits NO format header line: the first line is always the count line "nNodes nElems polyOrder". The two variants differ only in the per-element corner-node line: * ISM : 8 corner-node ids * ISM-MM : 8 corner-node ids followed by a material-name string The variant is auto-detected from the presence of the 9th token.

Recompute the active leaf set by depth-first traversal from the roots. Traversal (rather than a scan of all nodes) is what makes orphaned nodes left behind by coarsening invisible.

Level-based explicit-stability time step: refinement halves the element scale per level, so a time step dtBase that is stable on the base (level-0) mesh scales to dtBase / 2**MaxLevel on the current forest. Deterministic and exact for the quadtree (child elements are exact half-scale subdivisions); no geometry reduction is needed.

Build the connectivity of a uniformly (2:1 in each direction) refined 2-D mesh from its base connectivity. Every base element p is split into four children (global id 4*(p-1)+c). The scheme is fully deterministic integer bookkeeping - no coordinate hashing, no flip recomputation, no MPI:

Subdivide one leaf node into four children. Does nothing (with a warning) if the node is already refined. The caller is responsible for rebuilding the leaf list afterwards (or calling AdaptFromFlags, which does so).

Register a boundary condition function with the given bcid and bcname. If the bcid is already registered, the function pointer is updated. The elements and sides arrays are not allocated here; call PopulateBoundaries after scanning the mesh.

GPU regrid: release the device copies of the old boundary-condition element/side arrays, rebuild the model storage and BC maps on the new mesh (base Regrid), then upload the new BC arrays to the device - the same device bookkeeping Init/Free do around the base Init/Free.

Rebind a live model to a new mesh/geometry pair (AMR regrid). The mesh-sized solution storage is reallocated and the boundary-condition registrations and maps are rebuilt for the new mesh, while everything that is not mesh-sized is preserved: the time state (t, dt, entropy, IO counter), the time-integrator selection, configuration flags, and any model-specific parameters (Init is intent(out) and would reset all of these). nvar/nstepped are unchanged - the model solves the same equations on a new mesh.

Retire one decomposition from the process-wide lifecycle bookkeeping and, if it was the last one and SELF initialized MPI, finalize.

Base method for reporting the entropy of a model to stdout. Only override this procedure if additional reporting is needed. Alternatively, if you think additional reporting would be valuable for all models, open a pull request with modifications to this base method.

Base method for reporting the entropy of a model to stdout. Only override this procedure if additional reporting is needed. Alternatively, if you think additional reporting would be valuable for all models, open a pull request with modifications to this base method.

Base method for reporting the entropy of a model to stdout. Only override this procedure if additional reporting is needed. Alternatively, if you think additional reporting would be valuable for all models, open a pull request with modifications to this base method.

Method that can be overridden by users to report their own custom metrics after file io

Method that can be overridden by users to report their own custom metrics after file io

This method can be used to reset all of the boundary elements boundary condition type to the desired value.

This method can be used to reset all of the boundary elements boundary condition type to the desired value.

This method can be used to reset all of the boundary elements boundary condition type to the desired value.

Rebind to a new element count, reusing host pools and device buffers where they fit (AMR Stage 6b), and without uploading zeros the way Init does.

Rebind to a new element count, reusing storage where it fits (AMR Stage 6b). The mortar staging buffer is sized by mesh%nMortars rather than nElem, and nMortars changes with the mesh independently of the element count, so it is invalidated here and lazily re-created at the correct size by the next mortar exchange.

Rebind to a new element count, reusing host pools and device buffers where they fit (AMR Stage 6b). Inherits the Vector2D resize (host pools plus the five device buffers) and adds this class's mortar buffers, which are sized by mesh%nMortars rather than nElem: nMortars changes with the mesh independently of the element count, so a stale buffer would silently under-size the next mortar exchange. They are invalidated here and lazily re-created at the correct size on next use.

Rebind to a new element count, reusing storage where it fits (AMR Stage 6b). See Resize_MappedScalar2D_t for why the mortar staging buffer is invalidated rather than resized.

Rebind to a new element count, reusing host pools and device buffers where they fit. Deliberately does NOT call UpdateDevice: Init uploads freshly zeroed arrays, which in the adaptive loop is pure waste because the solution transfer overwrites them immediately.

Rebind a live object to a new element count, reusing the existing storage when it fits (AMR Stage 6b). Unlike Init - which is intent(out) and therefore resets the whole object, reallocates, zeroes and, on GPU builds, uploads those zeros - this preserves metadata and equation parsers (neither depends on nElem) and touches storage only when it must grow.

Rebind a live geometry to a new element count, reusing storage where it fits (AMR Stage 6c). This replaces the Free + Init cycle the adaptive loop performed on a freshly allocated SEMQuad every epoch, which threw away exactly the amortization Stage 6b introduced: each member Free released its pools and device buffers, and each Init reallocated, zeroed, rebuilt metadata and equation parsers, and uploaded the zeros.

Rebind to a new element count, reusing host pools and device buffers where they fit. Deliberately does NOT call UpdateDevice: Init uploads freshly zeroed arrays, which is pure waste in the adaptive loop because the arrays are rewritten immediately after.

Rebind a live object to a new element count, reusing existing storage when it fits (AMR Stage 6b/6c). Unlike Init - which is intent(out), so it resets the object, reallocates, zeroes and reconstructs the 4*nVar equation parsers - this preserves metadata and parsers (neither depends on nElem) and touches storage only when it must grow.

Rebind to a new element count, reusing host pools and device buffers where they fit. Deliberately does NOT call UpdateDevice: Init uploads freshly zeroed arrays, which is pure waste in the adaptive loop because the fields are rewritten immediately after.

Rebind a live object to a new element count, reusing existing storage when it fits (AMR Stage 6b). Unlike Init - which is intent(out), so it resets the object, reallocates, zeroes and reconstructs the equation parsers - this preserves metadata and parsers (neither depends on nElem) and touches storage only when it must grow. Preserving the parsers also avoids repeating the EquationParser construction that Init performs 2*nVar times as an amdflang workaround.

Read the Stage 6c geometry debug switches once.

Restrict (L2-project) the solution on four children back onto their parent element. Conservative, and the exact left inverse of ProlongToChildren.

Local Lax-Friedrichs (Rusanov) Riemann flux for linear advection. Entropy-stable: provides symmetric dissipation at element interfaces. Uses the per-face spectral radius lambda = |u.n|, matching upwind (Godunov) for linear advection: tangential faces (u.n=0) contribute zero flux, avoiding spurious dissipation across element interfaces whose face normal is perpendicular to the velocity.

Local Lax-Friedrichs (Rusanov) Riemann flux. Provided as a fallback; the model overrides BoundaryFlux directly with the LMARS solver.

Impedance-matched Riemann flux for acoustic variables 1-4, with zero flux returned for the auxiliary variables 5-7. The acoustic formula is identical to the parent LinearEuler2D model; see riemannflux2d_LinearEuler2D_t for the derivation.

Characteristic-decomposition (impedance-matched) interface flux for linear acoustics with possibly discontinuous sound speed.

Local Lax-Friedrichs (Rusanov) Riemann flux for linear advection. Uses the per-face spectral radius lambda = |u.n|, matching upwind (Godunov) for linear advection: tangential faces (u.n=0) contribute zero flux, avoiding spurious dissipation across element interfaces whose face normal is perpendicular to the velocity.

Local Lax-Friedrichs (Rusanov) Riemann flux.

Impedance-matched (characteristic/Godunov) Riemann flux, identical in form to the 2-D model's. The interface states are resolved with the per-side acoustic impedances Z = rho0*c (each side using its own background density rho0 and sound speed c), so material interfaces in a heterogeneous field are handled with the physically correct transmission/reflection:

Apply boundary conditions for the solution on GPU. Syncs boundary data from device, applies host-side BC dispatch (periodic defaults + registered BCs), then syncs back to device.

Apply boundary conditions for the solution. Periodic boundaries are set as the default; registered boundary conditions overwrite specific endpoints.

Apply registered boundary conditions for the solution. Each boundary condition method loops over its own boundary faces.

Apply registered boundary conditions for the solution.

Set the description of the ivar-th variable

Set the constant-coefficient Laplacian diffusion coefficients (kinematic momentum diffusivity and thermal diffusivity, both in m^2/s) and the dimensionless SIPG jump penalty. Setting nu or kappa > 0 enables the gradient pipeline so that the diffusive flux methods receive solutionGradient.

Set the constant-coefficient Laplacian diffusion coefficients (kinematic momentum diffusivity and thermal diffusivity, both in m^2/s) and the dimensionless SIPG jump penalty. Setting nu or kappa > 0 enables the gradient pipeline so that the diffusive flux methods receive solutionGradient.

Sets the equation parser for the ivar-th variable

Sets the equation parser for the idir direction and ivar-th variable

Sets the equation parser for the idir direction and ivar-th variable

Apply gradient boundary conditions on GPU. Syncs gradient boundary data from device, applies host-side BC dispatch (periodic defaults + registered BCs), then syncs back to device.

Apply boundary conditions for the solution gradient. Periodic boundaries are set as the default; registered boundary conditions overwrite specific endpoints.

Apply registered boundary conditions for the solution gradient. Each boundary condition method loops over its own boundary faces.

Apply registered boundary conditions for the solution gradient.

Initialise a hydrostatically balanced atmosphere with uniform potential temperature theta0, zero velocity, and the geopotential Phi = g*y carried as state variable index 5.

Initialise a hydrostatically balanced atmosphere with uniform potential temperature theta0, zero velocity, and the geopotential Phi = g*z carried as state variable index 6.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Sets the this % interior attribute using the eqn attribute, geometry (for physical positions), and provided simulation time.

Set the name of the ivar-th variable

Five conserved variables: (rho, rhou, rhov, rhotheta, Phi), where Phi = gy is the geopotential. Phi has zero flux (volume and surface) so its tendency is identically zero; it is carried in the state vector solely so that the Souza et al. (2023) non-conservative gravity flux differencing in SourceMethod can read it node-locally.

Six conserved variables: (rho, rhou, rhov, rhow, rhotheta, Phi), where Phi = g*z is the geopotential. Phi has zero flux (volume and surface) so its tendency is identically zero; it is carried in the state vector solely so that the Souza et al. (2023) non-conservative gravity flux differencing in SourceMethod can read it node-locally.

Populate the per-node sigma_x, sigma_y fields. Only nodes inside elements whose material name starts with this%pml_material_prefix receive non-zero damping; all other nodes are forced to zero so the interior solution remains unmodified.

Copy user-supplied physical coordinates into the cloud. xIn must be shape (nPoints, nDim) matching the init dimensions.

Sets the model % t attribute with the provided simulation time

Update the refine/coarsen thresholds without rebuilding the transform.

Sets the time integrator method, using a character input

Set the units of the ivar-th variable

Complete the aggregated halo exchange started by SideExchangeStart: wait on the persistent requests and unpack the received traces into extBoundary, applying side flips. The host blocks in MPI_Waitall (which also drives MPI progress) while previously launched device kernels continue to execute. No-op if no exchange is in flight.

Complete the aggregated halo exchange started by SideExchangeStart: wait on the persistent requests and unpack the received traces into extBoundary, applying side flips. The host blocks in MPI_Waitall (which also drives MPI progress) while previously launched device kernels continue to execute. No-op if no exchange is in flight.

Begin the aggregated halo exchange and launch the local (same-rank) side exchange kernel.

Begin the aggregated halo exchange and launch the local (same-rank) side exchange kernel.

Create the smallest 2:1 nonconforming (mortar) mesh: one "big" element of size 2dx x 2dx whose east edge is shared with the west edges of two "small" dx x dx elements. The mesh is conforming everywhere except at the single mortar interface.

No source term — upload the zero-initialised host array to device.

No source term — upload the zero-initialised host array to device.

Souza et al. (2023) non-conservative gravity flux differencing on GPU. The geopotential lives at solution(:,:,:,5); the source for rho*v is computed via the SBP-EC two-point form using log-mean density and the contravariant metric.

Souza et al. (2023) non-conservative gravity flux differencing.

Souza et al. (2023) non-conservative gravity flux differencing on GPU. The geopotential lives at solution(:,:,:,:,6); the source for rho*w is computed via the SBP-EC two-point form using log-mean density and the contravariant metric. Fully device-resident.

Souza et al. (2023) non-conservative gravity flux differencing.

Hu (2001) unsplit PML source term, evaluated per node. In the interior (sigma_x = sigma_y = 0) this leaves the acoustic variables untouched and the auxiliaries integrate q in time but never re-enter the dynamics (since the coupling coefficient sigma_xsigma_y is zero). Inside the PML, the (sigma_x + sigma_y) damping plus the sigma_xsigma_y*phi term produce the correct perfectly-matched behaviour for the linear Euler 2D system.

This subroutine sets the initial condition for a weak blast wave problem. The initial condition is given by

This subroutine sets the initial condition for a weak blast wave problem. The initial condition is given by

Stage the pre-regrid solution on the DEVICE (Stage 6a), replacing the base implementation's device-to-host copy of the whole field. Regrid subsequently releases solution%interior_gpu, so the field is copied device-to-device into a buffer this model owns; ApplyTransferPlan then reads it from there.

Preserve the current solution ahead of a regrid, so that Regrid may release the storage it lives in. Pair with ApplyTransferPlan, which consumes the staged copy:

Isoparametric subdivision of one element's geometry node coordinates into its four children. The parent geometry is the degree-nGeo Lagrange interpolant through parentCoords (sampled at geomInterp's control points, i.e. the mesh geometry nodes); each child node coordinate is that interpolant evaluated at the corresponding point of the parent reference square. Exact for any polynomial geometry of degree <= nGeo and for any control-node type, so straight-sided and curved (isoparametric) elements are both handled without approximation.

Arithmetic-mean two-point flux for linear advection. Entropy-conserving with respect to eta(u) = u^2/2.

Entropy-conserving two-point flux function.

Souza et al. (2023, JAMES) entropy-conservative two-point flux for 2-D compressible Euler in (rho, rhov, rhotheta) variables with p = p0(rhoRd*theta/p0)^gamma.

Arithmetic-mean two-point flux for linear advection. Entropy-conserving with respect to eta(u) = u^2/2.

Entropy-conserving two-point flux function (3-D). flux(ivar, d) is the d-th physical component (d=1:x, d=2:y, d=3:z). Override in concrete models. Stub returns zero.

Souza et al. (2023, JAMES) entropy-conservative two-point flux for compressible Euler in (rho, rhov, rhotheta) variables with p = p0(rhoRd*theta/p0)^gamma.

Contravariant EC two-point flux on GPU — fully device-resident.

Contravariant EC two-point flux on GPU — fully device-resident.

Computes pre-projected SCALAR contravariant two-point fluxes for all node pairs and stores them in twoPointFlux%interior(n,i,j,iel,ivar,r).

Computes pre-projected SCALAR contravariant two-point fluxes for all node pairs, following Trixi.jl for curved meshes. Each direction r uses the correct partner AND the correct averaged metric Ja^r.

Souza et al. (2023) entropy-conservative two-point flux on GPU.

Pre-projected scalar contravariant two-point Souza et al. (2023) EC flux. For each node pair (a, b) along reference direction r:

Souza et al. (2023) entropy-conservative two-point flux on GPU. Fully device-resident.

Pre-projected scalar contravariant two-point Souza et al. (2023) EC flux. For each node pair (a, b) along reference direction r:

Create a fully triply-periodic structured hexahedral mesh and store it in SELF's unstructured mesh format. Element geometry and ordering are identical to UniformStructuredMesh; the only difference is connectivity. The faces on the six domain boundaries are wired as interior faces whose neighbor is the element on the opposite side of the domain. This realises the triply periodic box T^3 = [0,Lx] x [0,Ly] x [0,Lz] required by, e.g., the Arnold-Beltrami-Childress (ABC) flow benchmark.

\addtogroup SELF_SupportRoutines @{ \fn UniformPoints Generates a REAL(prec) array of N points evenly spaced between two points.

Produce meshOut, a uniformly refined copy of meshIn: every element is split into four children. Child geometry is generated by exact isoparametric subdivision of the parent geometry (so the refined mesh represents the identical curved domain), and connectivity / flips are inherited deterministically from the base mesh. meshOut is a fully-formed, conforming Mesh2D_t with 4x the elements, ready for geometry generation and time stepping.

Create a structured mesh and store it in SELF's unstructured mesh format. The mesh is created in tiles of size (tnx,tny). Tiling is used to determine the element ordering.

Create a structured mesh and store it in SELF's unstructured mesh format. The mesh is created in tiles of size (tnx,tny,tnz). Tiling is used to determine the element ordering.

Copy elements, coordinates, and the per-point Lagrange basis cache from host to device. Lazily (re)allocates the cache buffers if their degree changed since the previous call.

Computes a solution update as , where dt is either provided through the interface or taken as the Model's stored time step size (model % dt)

Computes a solution update as , where dt is either provided through the interface or taken as the Model's stored time step size (model % dt)

Computes a solution update as , where dt is either provided through the interface or taken as the Model's stored time step size (model % dt)

Computes a solution update as , where dt is either provided through the interface or taken as the Model's stored time step size (model % dt)

Computes a solution update as , where dt is either provided through the interface or taken as the Model's stored time step size (model % dt)

Computes a solution update as , where dt is either provided through the interface or taken as the Model's stored time step size (model % dt)

Push the geometry the solver kernels read to the device. Mirrors the uploads that GenerateFromMesh's own path performs, for use when geometry was assembled by element copy rather than generated (AMR Stage 6c).

Generate the geometry the original way and report, per quantity, the largest discrepancy against the incrementally built one. Diagnostic only; gated by SELF_AMR_GEOM_VERIFY.

Writes the metadata to a HDF5 file using the fields : * /metadata/{group}/name/{varid} * /metadata/{group}/description/{varid} * /metadata/{group}/units/{varid}