Procedures

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.

Adds a pressure-balanced warm bubble perturbation.

Adds a pressure-balanced warm bubble perturbation.

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

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.

Axis-aligned bounding box of geometry%x%interior nodes for each element.

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.

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]

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

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

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

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.

\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

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

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

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 density, pressure, and sound speed.

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

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.

GPU-accelerated radiation BC for 2D Linear Euler.

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.

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.

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

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

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.

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.

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.

PreTendency 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 PreTendency simply returns back to the caller without executing any instructions

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.

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.

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.

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.

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.

Uses a local lax-friedrich's upwind flux The max eigenvalue is taken as the sound speed

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.

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

Sets the time integrator method, using a character input

Set the units of the ivar-th variable

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.

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

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:

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

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)

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