TrixiAtmo.jl API

TrixiAtmo

TrixiAtmo.jl is a simulation package for atmospheric models based on Trixi.jl

See also: trixi-framework/TrixiAtmo.jl

FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerEnergyEquationsWithGravity3D)

Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using an estimate c of the speed of sound.

References:

• Xi Chen et al. (2013) A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian Coordinate DOI: 10.1175/MWR-D-12-00129.1

AbstractCovariantEquations{NDIMS, 
                           NDIMS_AMBIENT, 
                           GlobalCoordinateSystem,
                           NVARS} <: AbstractEquations{NDIMS, NVARS}

Abstract type used to dispatch on systems of equations in covariant form, in which fluxes and prognostic variables are stored and computed in terms of their contravariant components defining their expansions in terms of the local covariant tangent basis. The type parameter NDIMS denotes the dimension of the manifold on which the equations are solved, while NDIMS_AMBIENT is the dimension of the ambient space in which such a manifold is embedded. Some references on discontinuous Galerkin methods in covariant flux form are listed below:

• M. Baldauf (2020). Discontinuous Galerkin solver for the shallow-water equations in covariant form on the sphere and the ellipsoid. Journal of Computational Physics 410:109384. DOI: 10.1016/j.jcp.2020.109384
• M. Baldauf (2021). A horizontally explicit, vertically implicit (HEVI) discontinuous Galerkin scheme for the 2-dimensional Euler and Navier-Stokes equations using terrain-following coordinates. Journal of Computational Physics 446:110635. DOI: 10.1016/ j.jcp.2021.110635
• L. Bao, R. D. Nair, and H. M. Tufo (2014). A mass and momentum flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere. A mass and momentum flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere. Journal of Computational Physics 271:224-243. DOI: 10.1016/j.jcp.2013.11.033

When using this equation type, functions which are evaluated pointwise, such as fluxes, source terms, and initial conditions take in the extra argument aux_vars, which contains the geometric information needed for the covariant form. The type parameter GlobalCoordinateSystem specifies the global coordinate system used to define the covariant tangent basis, and may be either GlobalCartesianCoordinates or GlobalSphericalCoordinates. The GlobalCoordinateSystem type parameter also specifies the coordinate system with respect to which the initial condition should be prescribed.

ChristoffelSymbolsAutodiff()

Struct used for multiple dispatch on functions that compute the Christoffel symbols. This option uses ForwardDiff.jl to compute

\[\Gamma_{jk}^i = \frac{1}{2}G^{il}\big(\partial_j G_{kl} + \partial_k G_{jl} - \partial_l G_{jk}\big)\]

using forward-mode automatic differentiation.

ChristoffelSymbolsCollocationDerivative()

Struct used for multiple dispatch on functions that compute the Christoffel symbols. Letting $I^N$ denote the degree $N$ polynomial interpolation operator collocated with the scheme's quadrature nodes, this option computes the Christoffel symbols at each quadrature node using the approximation

\[\Gamma_{jk}^i \approx \frac{1}{2}G^{il}\big(\partial_j I^N G_{kl} + \partial_k I^N G_{jl} - \partial_l I^N G_{jk}\big).\]

CompressibleEulerEnergyEquationsWithGravity2D(gamma)

The compressible Euler equations with gravity

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho e \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e +p) v_1 \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ (\rho e +p) v_2 \end{pmatrix} = \begin{pmatrix} 0 \\ - \rho \frac{\partial}{\partial x} \phi \\ - \rho \frac{\partial}{\partial y} \phi \\ - \rho v_1 \frac{\partial}{\partial x} \phi - rho v_2 \frac{\partial}{\partial y} \phi \end{pmatrix}\]

for an ideal gas with ratio of specific heats gamma in two space dimensions. Here, $\rho$ is the density, $v_1$, $v_2$ are the velocities, $e$ is the specific total energy rather than specific internal energy, $phi$ is the gravitational potential, and

\[p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right)\]

the pressure.

CompressibleEulerEnergyEquationsWithGravity3D(gamma)

The compressible Euler equations with gravity

\[\frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \\ \rho e \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ \rho v_1 v_3 \\ ( \rho e +p) v_1 \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ \rho v_1 v_3 \\ ( \rho e +p) v_2 \end{pmatrix} + \frac{\partial}{\partial z} \begin{pmatrix} \rho v_3 \\ \rho v_1 v_3 \\ \rho v_2 v_3 \\ \rho v_3^2 + p \\ ( \rho e +p) v_3 \end{pmatrix} = \begin{pmatrix} 0 \\ -\rho \frac{\partial}{\partial x} \phi \\ -\rho \frac{\partial}{\partial y} \phi \\ -\rho \frac{\partial}{\partial z} \\ -\rho v_1 \frac{\partial}{\partial x} \phi -\rho v_2 \frac{\partial}{\partial y} \phi -\rho v_3 \frac{\partial}{\partial z} \phi \end{pmatrix}\]

for an ideal gas with ratio of specific heats gamma in three space dimensions. Here, $\rho$ is the density, $v_1$, $v_2$, $v_3$ the velocities, $e$ the specific energy rather than specific internal energy, $\phi$ is the gravitational potential, and

\[p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2+v_3^2) \right)\]

the pressure.

CompressibleRainyEulerEquations2D{RealT <: Real} <:
    AbstractCompressibleRainyEulerEquations{2, 9}

The compressible Euler equations in two dimensions with gravity and separate densities for dry air ($\rho_d$), moist air ($\rho_m$), and rain (\rho_r).

\[\begin{alignat}{4} \partial_t \, &&(\rho_d) &+ \nabla \cdot (\rho_d\, v) &&= & 0\,,\\ \partial_t \, &&(\rho_m) &+ \nabla \cdot (\rho_m\, v) &&= & -S_r \,,\\ \partial_t \, &&(\rho_r) &+ \nabla \cdot (\rho_r\, (v-v_r\, e_z)) &&= & S_r\,,\\ \partial_t\,&&(\rho v) &+ \nabla \cdot (\rho v \otimes v - \rho_r\, v_r \,v\otimes e_z + p\cdot\textrm{Id}) &&= & -\rho\, g\, e_z\,,\\ \partial_t\, &&(E) &+ \nabla \cdot ((E + p)\,v - (c_l\, T + 1/2\, \langle v, v\rangle)\, \rho_r\, v_r\, e_z) &&= &\,\,-\rho\, g\, \langle e_z, v\rangle\,. \end{alignat}\]

Here, moist air is the sum of vapor ($\rho_v$) and cloud water ($\rho_c$), which are given implicitly by the nonlinear system

\[\begin{align} (c_{vd}\, \rho_d + c_{vv}\, \rho_v + c_l\, \rho_l)\,T + \rho_v\, L_{\textrm{ref}} + 1/2\,\rho\,\langle v, v\rangle - E &= 0\,,\\ \min \left( \frac{e_s(T)}{R_v\,T}, \rho_m\right) - \rho_v &=0\,,\\ \rho_m - \rho_v - \rho_c &=0\,, \end{align}\]

where $e_s(T)$ is the modeled saturation vapor pressure, $v_r$ the modeled terminal rain, fall velocity, and $S_r$ a modeled source term for rain conversion.

Reference

S. Doppler et al. (2024). A discontinuous Galerkin approach for atmospheric flows with implicit condensation. Journal of Computational Physics. 499:112713. DOI: 10.1016/j.jcp.2023.112713

CovariantLinearAdvectionEquation2D{GlobalCoordinateSystem} <:  
    AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3}

Denoting the covariant derivative by $\nabla_j$ and summing over repeated indices, a variable-coefficient linear advection equation can be defined on a two-dimensional manifold in three-dimensional ambient space as

\[\partial_t h + \nabla_j (hv^j) = 0.\]

We treat this problem as a system of equations in which the first variable is the scalar conserved quantity $h$, and the second two are the contravariant components $v^1$ and $v^2$ used in the expansion

\[\vec{v} = v^i \vec{a}_i = v^1 \vec{a}_1 + v^2 \vec{a}_2,\]

where $\vec{a}_1 = \partial \vec{x} / \partial \xi^1$ and $\vec{a}_2 = \partial \vec{x} / \partial \xi^2$ are the so-called covariant basis vectors, and $\xi^1$ and $\xi^2$ are the local reference space coordinates. The velocity components are spatially varying but assumed to be constant in time, so we do not apply any flux or dissipation to such variables. The resulting system is then given on the reference element as

\[J \frac{\partial}{\partial t} \left[\begin{array}{c} h \\ v^1 \\ v^2 \end{array}\right] + \frac{\partial}{\partial \xi^1} \left[\begin{array}{c} J h v^1 \\ 0 \\ 0 \end{array}\right] + \frac{\partial}{\partial \xi^2} \left[\begin{array}{c} J h v^2 \\ 0 \\ 0 \end{array}\right] = \left[\begin{array}{c} 0 \\ 0 \\ 0 \end{array}\right],\]

where $J = \lVert\vec{a}^1 \times \vec{a}^2 \rVert$ is the area element. Note that the variable advection velocity components could alternatively be stored as auxiliary variables, similarly to the geometric information.

CovariantShallowWaterEquations2D{GlobalCoordinateSystem} <:  
    AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3}

Denoting the covariant derivative by $\nabla_j$ and summing over repeated indices, the shallow water equations can be expressed on a two-dimensional surface in three-dimensional ambient space as

\[\begin{aligned} \partial_t h + \nabla_j (hv^j) &= 0,\\ \partial_t (hv^i) + \nabla_j \tau^{ij} + gh G^{ij}\partial_j b &= -fJ G^{ij}\varepsilon_{jk} hv^k, \end{aligned}\]

where $h$ is the geopotential height (equal to the total geopotential height $H$ for zero bottom topography), $v^i$ and $G^{ij}$ are the contravariant velocity and metric tensor components, $g$ is the gravitational acceleration, $b$ is the bottom topography, $f$ is the Coriolis parameter, $J$ is the area element, $\varepsilon$ is the Levi-Civita symbol, and $\partial_j$ is used as a shorthand for $\partial / \partial \xi^j$. The contravariant momentum flux tensor components are given by

\[\tau^{ij} = hv^i v^j + \frac{1}{2}G^{ij}gh^2.\]

The covariant shallow water equations with constant bottom topography can be formulated on the reference element as a system of conservation laws with a source term (implemented in the exported function source_terms_geometric_coriolis), as given by

\[J \frac{\partial}{\partial t} \left[\begin{array}{c} h \\ hv^1 \\ hv^2 \end{array}\right] + \frac{\partial}{\partial \xi^1} \left[\begin{array}{c} J h v^1 \\ J \tau^{11} \\ J \tau^{12} \end{array}\right] + \frac{\partial}{\partial \xi^2} \left[\begin{array}{c} J h v^2 \\ J \tau^{21} \\ J \tau^{22} \end{array}\right] = J \left[\begin{array}{c} 0 \\ -\Gamma^1_{jk}\tau^{jk} - f J \big(G^{12}hv^1 - G^{11}hv^2\big) \\ -\Gamma^2_{jk}\tau^{jk} - f J \big(G^{22}hv^1 - G^{21}hv^2\big) \end{array}\right].\]

Note that the geometric contribution to the source term involves the Christoffel symbols of the second kind, which can been expressed in terms of the covariant metric tensor components $G_{ij}$ as

\[\Gamma_{jk}^i = \frac{1}{2}G^{il}\big(\partial_j G_{kl} + \partial_k G_{jl} - \partial_l G_{jk}\big).\]

References

• M. Baldauf (2020). Discontinuous Galerkin solver for the shallow-water equations in covariant form on the sphere and the ellipsoid. Journal of Computational Physics 410:109384. DOI: 10.1016/j.jcp.2020.109384
• L. Bao, R. D. Nair, and H. M. Tufo (2014). A mass and momentum flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere. A mass and momentum flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere. Journal of Computational Physics 271:224-243. DOI: 10.1016/j.jcp.2013.11.033

GlobalCartesianCoordinates()

Struct used for dispatch, specifying that the covariant tangent basis vectors should be defined with respect to a global Cartesian coordinate system.

GlobalSphericalCoordinates()

Struct used for dispatch, specifying that the covariant tangent basis vectors should be defined with respect to a global spherical coordinate system.

MetricTermsCovariantSphere(christoffel_symbols = ChristoffelSymbolsAutodiff())

Struct specifying options for computing geometric information for discretizations in covariant form based on an exact representation of the spherical geometry. Currently, the only field is christoffel_symbols, specifying the approach used to compute the Christoffel symbols, for which the options are ChristoffelSymbolsAutodiff or ChristoffelSymbolsCollocationDerivative.

MetricTermsCrossProduct()

Struct used for multiple dispatch on functions that compute the metric terms. When the argument metric_terms is of type MetricTermsCrossProduct, the contravariant vectors are computed using the cross-product form.

MetricTermsInvariantCurl()

Struct used for multiple dispatch on functions that compute the metric terms. When the argument metric_terms is of type MetricTermsInvariantCurl, the contravariant vectors are computed using the invariant curl form.

References

• Kopriva, D. A. (2006). Metric identities and the discontinuous spectral element method on curvilinear meshes. Journal of Scientific Computing 26, 301-327. DOI: 10.1007/s10915-005-9070-8
• Vinokur, M. and Yee, H. C. (2001). Extension of efficient low dissipation high order schemes for 3-D curvilinear moving grids. In Caughey, D. A., and Hafez, M. M. (eds.), Frontiers of Computational Fluid Dynamics 2002, World Scientific, Singapore, pp. 129–164. DOI: 10.1142/9789812810793_0008

NonlinearSolveDG

Newton method, which can be called in every stage via callbacks.

Parameters

• residual::Function: function evaluating the residual
• jacobian::Function: function evaluating the Jacobian of the residual
• variables_index_vector::Vector{Int64}: vector of indices of entries of the solution vector the Newton method operates on
• tolerance::Real: tolerance for termination of the Newton method
• max_iterations::Int64: maximal number of iterations of the Newton method

ShallowWaterEquations3D(; gravity, rotation_rate = 0, H0 = 0)

Rotating shallow water equations (SWE) in three space dimensions to be solved on curved manifolds (e.g., on a spherical shell). The equations are given by

\[\begin{aligned} \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h v_1) + \frac{\partial}{\partial y}(h v_2) + \frac{\partial}{\partial z}(h v_3) &= s_h \\ \frac{\partial}{\partial t}(h v_1) + \frac{\partial}{\partial x}\left(h v_1^2 + \frac{g}{2}h^2\right) + \frac{\partial}{\partial y}(h v_1 v_2) + \frac{\partial}{\partial z}(h v_1 v_3) g h \frac{\partial b}{\partial x} &= s_{hv_1} \\ \frac{\partial}{\partial t}(h v_2) + \frac{\partial}{\partial x}(h v_1 v_2) + \frac{\partial}{\partial y}\left(h v_2^2 + \frac{g}{2}h^2\right) + \frac{\partial}{\partial z}(h v_2 v_3) + g h \frac{\partial b}{\partial y} &= s_{hv_2} \\ \frac{\partial}{\partial t}(h v_3) + \frac{\partial}{\partial x}(h v_1 v_3) + \frac{\partial}{\partial y}(h v_2 v_3) + \frac{\partial}{\partial z}\left(h v_3^2 + \frac{g}{2}h^2\right) + g h \frac{\partial b}{\partial z}&= s_{hv_3}. \end{aligned}\]

The unknown quantities of the SWE are the water height $h$ and the velocities $\mathbf{v} = (v_1, v_2, v_3)^T$. The gravitational acceleration is denoted by g.

The 3D Shallow Water Equations (SWE) extend the 2D SWE to model shallow water flows on 2D manifolds embedded within 3D space. To confine the flow to the 2D manifold, a source term incorporating a Lagrange multiplier is applied to the momentum equations using the function source_terms_lagrange_multiplier. This term effectively removes momentum components that are normal to the manifold, ensuring the flow remains constrained within the 2D surface.

To incorporate the effect of the rotation of the manifold, use the function source_terms_coriolis, which adds the necessary Coriolis source terms to the momentum equations assuming a rotation around the $z$ axis with a rotation rate in radians per time unit given by rotation_rate. To incorporate both Coriolis forces and the Lagrange multiplier terms, use source_terms_coriolis_lagrange_multiplier.

The additional quantity $H_0$ is also available to store a reference value for the total water height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness.

In addition to the unknowns, TrixiAtmo.jl currently stores the bottom topography values at the approximation points despite being fixed in time. This is done for convenience of computing the bottom topography gradients on the fly during the approximation as well as computing auxiliary quantities like the total water height $H$ or the entropy variables. This affects the implementation and use of these equations in various ways:

• The flux values corresponding to the bottom topography must be zero.
• The bottom topography values must be included when defining initial conditions, boundary conditions or source terms.
• AnalysisCallback analyzes this variable.
• Trixi.jl's visualization tools will visualize the bottom topography by default.

References:

• J. Coté (1988). "A Lagrange multiplier approach for the metric terms of semi-Lagrangian models on the sphere". Quarterly Journal of the Royal Meteorological Society 114, 1347-1352. DOI: 10.1002/qj.49711448310
• F. X. Giraldo (2001). "A spectral element shallow water model on spherical geodesic grids". DOI: 10.1002/1097-0363(20010430)35:8<869::AID-FLD116>3.0.CO;2-S

SplitCovariantShallowWaterEquations2D{GlobalCoordinateSystem} <:  
    AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3}

Alternative flux formulation of CovariantShallowWaterEquations2D given on the reference element by

\[J \frac{\partial}{\partial t} \left[\begin{array}{c} h \\ hv^1 \\ hv^2 \end{array}\right] + \frac{\partial}{\partial \xi^1} \left[\begin{array}{c} J h v^1 \\ J h v^1 v^1 \\ J h v^2 v^1 \end{array}\right] + \frac{\partial}{\partial \xi^2} \left[\begin{array}{c} J h v^2 \\ J h v^1 v^2 \\ J h v^2 v^2 \end{array}\right] + J \left[\begin{array}{c} 0 \\ \Upsilon^1 \\ \Upsilon^2 \end{array}\right] = J\left[\begin{array}{c}0 \\ s^1 \\ s^2 \end{array}\right].\]

In the above, the non-conservative differential terms in the momentum equations are given by

\[\Upsilon^i = \frac{1}{2}hv^j\big(G^{ik}\partial_j v_k - \partial_j v^i\big) + ghG^{ij}\partial_j (h + b),\]

where we allow for a variable bottom topography defined by $h_s$, and the algebraic momentum source terms implemented in source_terms_geometric_coriolis are given by

\[s^i = -\frac{1}{2}\big(\Gamma_{jk}^i hv^j v^k - G^{ik}\Gamma_{jk}^lh v^j v_l \big) - f JG^{ij}\varepsilon_{jk} hv^k.\]

In the above, we employ the same notation as in CovariantShallowWaterEquations2D (including summation over repeated indices) and note that the covariant velocity components are given by $v_i = G_{ij} v^j$. To obtain an entropy-conservative scheme with respect to the total energy

\[\eta = \frac{1}{2}h(v_1 v^1 + v_2v^2) + \frac{1}{2}gh^2 + ghb,\]

this equation type should be used with volume_flux = (flux_ec, flux_nonconservative_ec).

flux_chandrashekar(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEnergyEquationsWithGravity3D)

Entropy conserving two-point flux by

• Chandrashekar (2013) Kinetic Energy Preserving and Entropy Stable Finite Volume Schemes for Compressible Euler and Navier-Stokes Equations DOI: 10.4208/cicp.170712.010313a

flux_ec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature2D)

Entropy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_ec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature3D)

Entropy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_ec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperatureWithGravity2D)

Entropy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_ec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperatureWithGravity3D)

Entropy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_ec(u_ll, u_rr, aux_vars_ll, aux_vars_rr,
                           orientation::Integer,
                           equations::SplitCovariantShallowWaterEquations2D)

Symmetric part of an entropy-conservative flux for the shallow water equations in covariant form. Note that this does not include the pressure term or the non-symmetric curvature correction term. When used with flux_nonconservative_ec for the nonconservative volume and surface terms, this flux recovers the formulation described in the following paper for the special case of the Euclidean metric $G_{ab} = \delta_{ab}$:

• N. Wintermeyer, A. R. Winters, G. J. Gassner, and D. A. Kopriva (2017). An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry. Journal of Computational Physics 300:240-242. DOI: 10.1016/j.jcp.2017.03.036

flux_ec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature1D)

Entropy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_ec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperatureWithGravity1D)

Entropy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_fjordholm_etal(u_ll, u_rr,
                    normal_direction::AbstractVector,
                    equations::ShallowWaterEquations3D)

Total energy conservative (mathematical entropy for shallow water equations). When the bottom topography is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. For well-balancedness in the volume flux use flux_wintermeyer_etal.

Details are available in Eq. (4.1) in the paper:

• Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011) Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography DOI: 10.1016/j.jcp.2011.03.042

flux_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction,
                    equations::CompressibleEulerEnergyEquationsWithGravity2D)

Kinetic energy preserving two-point flux by

• Kennedy and Gruber (2008) Reduced aliasing formulations of the convective terms within the Navier-Stokes equations for a compressible fluid DOI: 10.1016/j.jcp.2007.09.020

flux_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction,
                    equations::CompressibleEulerEnergyEquationsWithGravity3D)

Kinetic energy preserving two-point flux by

• Kennedy and Gruber (2008) Reduced aliasing formulations of the convective terms within the Navier-Stokes equations for a compressible fluid DOI: 10.1016/j.jcp.2007.09.020

flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer,
                                    equations::ShallowWaterEquations3D)
flux_nonconservative_fjordholm_etal(u_ll, u_rr,
                                    normal_direction::AbstractVector,
                                    equations::ShallowWaterEquations3D)

Non-symmetric two-point surface flux discretizing the nonconservative (source) term of that contains the gradient of the bottom topography ShallowWaterEquations3D.

This flux can be used together with flux_fjordholm_etal at interfaces to ensure entropy conservation and well-balancedness.

Further details for the original finite volume formulation are available in

• Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011) Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography DOI: 10.1016/j.jcp.2011.03.042

and for curvilinear 2D case in the paper:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036

flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                      equations::ShallowWaterEquations3D)
flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                      normal_direction::AbstractVector,
                                      equations::ShallowWaterEquations3D)

Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography ShallowWaterEquations3D.

For the surface_flux either flux_wintermeyer_etal or flux_fjordholm_etal can be used to ensure well-balancedness and entropy conservation.

Further details are available in the papers:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036
• Patrick Ersing, Andrew R. Winters (2023) An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on curvilinear meshes DOI: 10.48550/arXiv.2306.12699

flux_ranocha(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerEnergyEquationsWithGravity2D)

Entropy conserving and kinetic energy preserving two-point flux by

• Hendrik Ranocha (2018) Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods for Hyperbolic Balance Laws PhD thesis, TU Braunschweig

See also

• Hendrik Ranocha (2020) Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators Proceedings of ICOSAHOM 2018

flux_ranocha(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleEulerEnergyEquationsWithGravity3D)

Entropy conserving and kinetic energy preserving two-point flux by

• Hendrik Ranocha (2018) Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods for Hyperbolic Balance Laws PhD thesis, TU Braunschweig

See also

• Hendrik Ranocha (2020) Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators Proceedings of ICOSAHOM 2018

flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction,
                equations::CompressibleEulerEnergyEquationsWithGravity2D)

This flux is is a modification of the original kinetic energy preserving two-point flux by

• Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) Kinetic energy and entropy preserving schemes for compressible flows by split convective forms DOI: 10.1016/j.jcp.2018.08.058

The modification is in the energy flux to guarantee pressure equilibrium and was developed by

• Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) Preventing spurious pressure oscillations in split convective form discretizations for compressible flows DOI: 10.1016/j.jcp.2020.110060

flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction,
                equations::CompressibleEulerEnergyEquationsWithGravity3D)

This flux is is a modification of the original kinetic energy preserving two-point flux by

• Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018) Kinetic energy and entropy preserving schemes for compressible flows by split convective forms DOI: 10.1016/j.jcp.2018.08.058

The modification is in the energy flux to guarantee pressure equilibrium and was developed by

• Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020) Preventing spurious pressure oscillations in split convective form discretizations for compressible flows DOI: 10.1016/j.jcp.2020.110060

flux_wintermeyer_etal(u_ll, u_rr,
                      normal_direction::AbstractVector,
                      equations::ShallowWaterEquations3D)

Total energy conservative (mathematical entropy for shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a volume_flux. For the surface_flux either flux_wintermeyer_etal or flux_fjordholm_etal can be used to ensure well-balancedness and entropy conservation.

Further details are available in Theorem 1 of the paper:

• Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry DOI: 10.1016/j.jcp.2017.03.036

P4estMeshCubedSphere2D(trees_per_face_dimension, radius;
                        polydeg, RealT=Float64,
                        initial_refinement_level=0, unsaved_changes=true,
                        p4est_partition_allow_for_coarsening=true,
                        element_local_mapping=false)

Build a "Cubed Sphere" mesh as a 2D P4estMesh with 6 * trees_per_face_dimension^2 trees.

The mesh will have no boundaries.

Arguments

• trees_per_face_dimension::Integer: the number of trees in the two local dimensions of each face.
• radius::Integer: the radius of the sphere.
• polydeg::Integer: polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree.
• RealT::Type: the type that should be used for coordinates.
• initial_refinement_level::Integer: refine the mesh uniformly to this level before the simulation starts.
• unsaved_changes::Bool: if set to true, the mesh will be saved to a mesh file.
• p4est_partition_allow_for_coarsening::Bool: Must be true when using AMR to make mesh adaptivity independent of domain partitioning. Should be false for static meshes to permit more fine-grained partitioning.
• element_local_mapping::Bool: option to use the alternative element-local mapping from Appendix A of Guba et al. (2014). If set to true, the four corner vertex positions for each element will be obtained through an equiangular gnomonic projection (Ronchi et al. 1996), and the tree node coordinates within the element (i.e. the field tree_node_coordinates) will be obtained by first using a bilinear mapping based on the four corner vertices, and then projecting the bilinearly mapped nodes onto the spherical surface by normalizing the resulting Cartesian coordinates and scaling by radius. If set to false, the equiangular gnomonic projection will be used for all tree node coordinates.

P4estMeshQuadIcosahedron2D(trees_per_face_dimension, radius;
                           polydeg, RealT=Float64,
                           initial_refinement_level=0, unsaved_changes=true,
                           p4est_partition_allow_for_coarsening=true)

Build a quad-based icosahedral mesh as a 2D P4estMesh with 60 * trees_per_face_dimension^2 trees (20 triangular faces of the icosahedron, each subdivided into 3 parent quads, each of which subdivided into trees_per_face_dimension^2 trees).

The node coordinates of the trees will be obtained using the element-local mapping from Appendix A of Guba et al. (2014). See P4estMeshCubedSphere2D for more information about the element-local mapping.

The mesh will have no boundaries.

Arguments

• trees_per_face_dimension::Integer: the number of trees in the two local dimensions of each parent quad.
• radius::Integer: the radius of the sphere.
• polydeg::Integer: polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree.
• RealT::Type: the type that should be used for coordinates.
• initial_refinement_level::Integer: refine the mesh uniformly to this level before the simulation starts.
• unsaved_changes::Bool: if set to true, the mesh will be saved to a mesh file.
• p4est_partition_allow_for_coarsening::Bool: Must be true when using AMR to make mesh adaptivity independent of domain partitioning. Should be false for static meshes to permit more fine-grained partitioning.

     clean_solution_lagrange_multiplier!(u, equations::ShallowWaterEquations3D, normal_direction)

Function to apply Lagrange multiplier discretely to the solution in order to constrain the momentum to a 2D manifold.

The vector normal_direction is perpendicular to the 2D manifold. By default, it is the normal contravariant basis vector.

contravariant2global(u, aux_vars, equations)

Transform the vector u of solution variables with the momentum or velocity given in terms of local contravariant components into the global coordinate system specified by the GlobalCoordinateSystem type parameter in AbstractCovariantEquations. u is a vector type of the correct length nvariables(equations). Notice the function doesn't include any error checks for the purpose of efficiency, so please make sure your input is correct. The inverse conversion is performed by global2contravariant.

examples_dir()

Return the directory where the example files provided with TrixiAtmo.jl are located. If TrixiAtmo.jl is installed as a regular package (with ]add Trixi), these files are read-only and should not be modified. To find out which files are available, use, e.g., readdir:

Examples

readdir(examples_dir())

flux_ec_rain(u_ll, u_rr, orientation_or_normal_direction,
             equations::CompressibleRainyEulerEquations2D)

Entropy-conserving flux including rain, derived in: Fabian Höck A Discontinuous Galerkin Method for Moist Atmospheric Dynamics with Rain Master's thesis, University of Cologne, 2025

flux_etec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature2D)

Entropy and total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_etec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature3D)

Entropy and total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_etec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperatureWithGravity2D)

Entropy and total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_etec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperatureWithGravity3D)

Entropy and total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_etec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature1D)

Entropy and total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_etec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperatureWithGravity1D)

Entropy and total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_nonconservative_artiano_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEnergyEquationsWithGravity2D)

Well-balanced gravity term for a constant potential temperature background state for the CompressibleEulerEnergyEquationsWithGravity2D developed by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_nonconservative_artiano_etal(u_ll, u_rr,
								  normal_direction::AbstractVector,
								  equations::CompressibleEulerEnergyEquationsWithGravity3D)

Well-balanced gravity term for constant potential temperature background state by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_nonconservative_artiano_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D)

Well-balanced gravity term for constant potential temperature background state by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_nonconservative_artiano_etal(u_ll, u_rr,
								  normal_direction::AbstractVector,
								  equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D)

Well-balanced gravity term for constant potential temperature background state by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

fluxnonconservativeartianoetal(ull, u_rr, orientation::Integer, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D)

Well-balanced gravity term for constant potential temperature background state by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_nonconservative_ec(u_ll, u_rr, aux_vars_ll, aux_vars_rr,
                        orientation::Integer,
                        equations::SplitCovariantShallowWaterEquations2D)

Non-symmetric part of an entropy-conservative flux for the shallow water equations in covariant form, consisting of pressure and bottom topography terms, as well as a curvature correction term. This can be used for both the volume and surface terms.

flux_nonconservative_souza_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEnergyEquationsWithGravity2D)

Kinetic and potential energy preserving (KPEP) gravity term for the CompressibleEulerEnergyEquationsWithGravity2D developed by

• Souza et al. The Flux-Differencing Discontinuous {G}alerkin Method Applied to an Idealized Fully Compressible Nonhydrostatic Dry Atmosphere [DOI: 10.1029/2022MS003527] (https://doi.org/10.1029/2022MS003527)

flux_nonconservative_souza_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEnergyEquationsWithGravity3D)

• Souza et al. The Flux-Differencing Discontinuous {G}alerkin Method Applied to an Idealized Fully Compressible Nonhydrostatic Dry Atmosphere [DOI: 10.1029/2022MS003527] (https://doi.org/10.1029/2022MS003527)

flux_nonconservative_souza_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D)

• Souza et al. The Flux-Differencing Discontinuous {G}alerkin Method Applied to an Idealized Fully Compressible Nonhydrostatic Dry Atmosphere [DOI: 10.1029/2022MS003527] (https://doi.org/10.1029/2022MS003527)

flux_nonconservative_souza_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D)

• Souza et al. The Flux-Differencing Discontinuous {G}alerkin Method Applied to an Idealized Fully Compressible Nonhydrostatic Dry Atmosphere [DOI: 10.1029/2022MS003527] (https://doi.org/10.1029/2022MS003527)

fluxnonconservativesouzaetal(ull, u_rr, orientation::Integer, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D)

• Souza et al. The Flux-Differencing Discontinuous {G}alerkin Method Applied to an Idealized Fully Compressible Nonhydrostatic Dry Atmosphere [DOI: 10.1029/2022MS003527] (https://doi.org/10.1029/2022MS003527)

flux_nonconservative_surface_simplified(u_ll, u_rr, aux_vars_ll, aux_vars_rr,
                                        orientation::Integer,
                                        equations::SplitCovariantShallowWaterEquations2D)

For bottom topography which is continuous across element interfaces, we can significantly simplify the nonconservative surface terms, such that only the pressure contribution remains. In such cases, this flux is equivalent to flux_nonconservative_ec when used as a surface flux, but should not be used as a volume flux.

flux_nonconservative_waruszewski_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEnergyEquationsWithGravity2D)

Well-balanced gravity term for an isothermal background state for the CompressibleEulerEnergyEquationsWithGravity2D developed by

• Maciej Waruszewski and Jeremy E. Kozdon and Lucas C. Wilcox and Thomas H. Gibson and Francis X. Giraldo (2022) Entropy stable discontinuous {G}alerkin methods for balance laws in non-conservative form: Applications to the {E}uler equations with gravity DOI: 10.1016/j.jcp.2022.111507

The well-balancedness on curvilinear coordinates was proven by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_nonconservative_waruszewski_etal(u_ll, u_rr,
									 normal_direction::AbstractVector,
									 	equations::CompressibleEulerEulerEquationsWithGravity3D)

Well-balanced gravity term for isothermal background state

• Maciej Waruszewski and Jeremy E. Kozdon and Lucas C. Wilcox and Thomas H. Gibson and Francis X. Giraldo (2022), Entropy stable discontinuous {G}alerkin methods for balance laws in non-conservative form: Applications to the {E}uler equations with gravity DOI: 10.1016/j.jcp.2022.111507

The well balanced on curvilinear coordinates was proven by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_nonconservative_waruszewski_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D)

Well-balanced gravity term for isothermal background state

• Maciej Waruszewski and Jeremy E. Kozdon and Lucas C. Wilcox and Thomas H. Gibson and Francis X. Giraldo (2022) Entropy stable discontinuous {G}alerkin methods for balance laws in non-conservative form: Applications to the {E}uler equations with gravity DOI: 10.1016/j.jcp.2022.111507

The well balanced on curvilinear coordinates was proven by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_nonconservative_waruszewski_etal(u_ll, u_rr,
									 normal_direction::AbstractVector,
									 	equations::CompressibleEulerPotentialTemperatureEquationsWithGravity3D)

Well-balanced gravity term for isothermal background state

• Maciej Waruszewski and Jeremy E. Kozdon and Lucas C. Wilcox and Thomas H. Gibson and Francis X. Giraldo (2022), Entropy stable discontinuous {G}alerkin methods for balance laws in non-conservative form: Applications to the {E}uler equations with gravity DOI: 10.1016/j.jcp.2022.111507

The well balanced on curvilinear coordinates was proven by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

fluxnonconservativewaruszewskietal(ull, u_rr, orientation::Integer, equations::CompressibleEulerPotentialTemperatureEquationsWithGravity1D)

Well-balanced gravity term for isothermal background state

• Maciej Waruszewski and Jeremy E. Kozdon and Lucas C. Wilcox and Thomas H. Gibson and Francis X. Giraldo (2022), Entropy stable discontinuous {G}alerkin methods for balance laws in non-conservative form: Applications to the {E}uler equations with gravity DOI: 10.1016/j.jcp.2022.111507

The well balanced on curvilinear coordinates was proven by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_tec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature2D)

Total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_tec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature3D)

Total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_tec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperatureWithGravity2D)

Total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_tec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperatureWithGravity3D)

Total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_tec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature1D)

Total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

flux_tec(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquationsPotentialTemperature1D)

Total energy conservative two-point flux by

• Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025) Structure-Preserving High-Order Methods for the Compressible Euler Equations in Potential Temperature Formulation for Atmospheric Flows (https://arxiv.org/abs/2509.10311)

global2contravariant(u, aux_vars, equations)

Transform the vector u of solution variables with momentum or velocity components given with respect to the global coordinate system into local contravariant components. The global coordinate system is specified by the GlobalCoordinateSystem type parameter in AbstractCovariantEquations. u is a vector type of the correct length nvariables(equations). Notice the function doesn't include any error checks for the purpose of efficiency, so please make sure your input is correct. The inverse conversion is performed by contravariant2global.

have_aux_node_vars(equations)

Trait function determining whether equations requires the use of auxiliary variables. Classical conservation laws such as the CompressibleEulerEquations2D do not require auxiliary variables. The return value will be True() or False() to allow dispatching on the return type.

initial_condition_barotropic_instability(x, t, equations)

Barotrotropic instability initiated by a perturbation applied to a mid-latitude jet. The velocity field is a purely zonal flow, given as a function of the latitude $\theta$ as

\[v_\lambda(\theta) = \begin{cases} u_0 \exp(-4 / (\theta_1 - \theta_0)^{-2})\exp((\theta - \theta_0)^{-1} (\theta - \theta_1)^{-1}), & \quad \theta_0 < \theta < \theta_1, \\ 0 & \quad \text{otherwise}, \end{cases}\]

where $u_0 = 80 \ \mathrm{m}/\mathrm{s}$, $\theta_0 = \pi/7$, and $\theta_1 = \pi/2 - \theta_0$. The background geopotential height field is given by

\[h_0(\theta) = 10158 \ \mathrm{m} - \frac{a}{g} \int_{-\pi/2}^\theta v_\lambda(\theta')\big(2\Omega\sin\theta' + v_\lambda(\theta')\tan\theta' / a \big)\, \mathrm{d}\theta',\]

where $a = 6.37122 \times 10^3\ \mathrm{m}$ is the Earth's radius, $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$ is the Earth's gravitational acceleration, and $\Omega = 7.292 \times 10^{-5}\ \mathrm{s}^{-1}$ is the Earth's rotation rate. The perturbation is then added to obtain the following geopotential height field:

\[h(\lambda, \theta) = \begin{cases} h_0(\theta) + \delta h \cos\theta \exp(-(\lambda/\alpha)^2) \exp(-((\theta_2 -\theta)/\beta)^2), & \quad -\pi < \lambda < \pi,\\ h_0(\theta), & \quad \text{otherwise}, \end{cases}\]

where $\lambda$ is the longitude coordinate, and we take $\alpha = 1/3$, $\beta = 1/15$, and $\delta h = 120 \ \mathrm{m}$. This problem was proposed in the following paper:

• J. Galewsky, R. K. Scott, and L. M. Polvani (2004). An initial-value problem for testing numerical models of the global shallow-water equations. Tellus A 56.5:429–440. DOI: 10.3402/tellusa.v56i5.14436

initial_condition_gaussian(x, t, equations)

This Gaussian bell case is a smooth initial condition suitable for testing the convergence of discretizations of the linear advection equation on a spherical domain of radius $a = 6. 37122 \times 10^3\ \mathrm{m}$, representing the surface of the Earth. Denoting the Euclidean norm as $\lVert \cdot \rVert$, the initial height field is prescribed as a function of the position $\vec{x}$ relative to the centre of the Earth by

\[h(\vec{x}) = h_0 \exp \Big(-b_0 \big(\lVert \vec{x} - \vec{x}_0 \rVert / \lVert \vec{x} \rVert\big)^2 \Big),\]

where $h_0 = 1 \times 10^3\ \mathrm{m}$ is the height of the bell, $b_0 = 5$ is the width parameter, and $\vec{x}_0$ is the position of the centre of the bell, which is initialized at a longitude of $3\pi/2$ and a latitude of zero. The velocity field corresponds to a solid body rotation with a period of 12 days at an angle of $\alpha = \pi/4$ from the polar axis. Denoting $\vec{\omega}$ as the corresponding angular velocity vector, the velocity is therefore initialized as

\[\vec{v}(\vec{x}) = \vec{\omega} \times \vec{x}.\]

This problem is adapted from Case 1 of the test suite described in the following paper:

• D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224. DOI: 10.1016/S0021-9991(05)80016-6

initial_condition_geostrophic_balance(x, t, equations)

Steady geostrophic balance for the spherical shallow water equations, corresponding to a purely zonal velocity field given as a function of the latitude $\theta$ by $v_\lambda(\theta) = v_0 \cos\theta$, where we define $v_0 = 2\pi a / (12 \ \mathrm{days})$ in terms of the Earth's radius $a = 6.37122 \times 10^3\ \mathrm{m}$. The height field then varies with the latitude as

\[h(\theta) = \frac{1}{g} \Big(gh_0 - \Big(a \Omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta\Big),\]

where $gh_0 = 2.94 \times 10^4 \ \mathrm{m}^2/\mathrm{s}^2$, $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and $\Omega = 7.292 \times 10^{-5}\ \mathrm{s}^{-1}$. This problem corresponds to Case 2 of the test suite described in the following paper:

• D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224. DOI: 10.1016/S0021-9991(05)80016-6

initial_condition_isolated_mountain(x, t, equations)

Zonal flow over an isolated mountain with a profile given in terms of the latitude $\lambda$ and longitude $\theta$ as

\[h_s(\lambda,\theta) = h_{s0} (1 - \sqrt{\min(R^2, (\lambda-\lambda_0)^2 + (\theta-\theta_0)^2)}/R),\]

where $h_{s0} = 2000 \ \text{m}$, $\lambda_0 = -\pi/2$, $\theta_0 = \pi/6$, and $R =\pi/9$. The initial velocity field is given by $v_\lambda(\theta) = v_0 \cos\theta$, where $v_0 = 20 \ \mathrm{m/s}$, and the total geopotential height $H = h+h_s$ is given by

\[H(\theta) = H_0 - \frac{1}{g}\Big(a \Omega v_0 + \frac{1}{2} v_0^2\Big)\sin^2\theta,\]

where $H_0 = 5960 \ \mathrm{m}$, $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, and $\Omega = 7.292 \times 10^{-5}\ \mathrm{s}^{-1}$. To use this test case with SplitCovariantShallowWaterEquations2D, the keyword argument auxiliary_field = bottom_topography_isolated_mountain should be passed into the SemidiscretizationHyperbolic constructor. This problem corresponds to Case 5 of the test suite described in the following paper:

• D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224. DOI: 10.1016/S0021-9991(05)80016-6

initial_condition_rossby_haurwitz(x, t, equations)

Rossby-Haurwitz wave case for the spherical shallow water equations, where the zonal and meridional velocity components are given, respectively, as functions of the longitude $\lambda$ and latitude $\theta$ by

\[\begin{aligned} v_\lambda(\lambda,\theta) &= a \omega \cos \theta+a K \cos ^{R-1} \theta \left(R \sin ^2 \theta-\cos ^2 \theta\right) \cos (R \lambda),\\ v_\theta(\lambda,\theta) &= -a K R \cos ^{R-1} \theta \sin \theta \sin (R \lambda), \end{aligned}\]

where $\omega = K = 7.848 \times 10^{-6} \ \mathrm{s}^{-1}$ and $R = 4$ are given constants, and $a = 6.37122 \times 10^3\ \mathrm{m}$ is the Earth's radius. Taking $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\Omega = 7.292 \times 10^{-5} \ \mathrm{s}^{-1}$, and $h_0 = 8000 \ \mathrm{m}$ and defining the functions

\[\begin{aligned} A(\theta) &= \frac{\omega}{2}(2 \Omega+\omega) \cos^2 \theta + \frac{1}{4} K^2 \cos^{2 R} \theta\Big((R+1) \cos^2\theta +\left(2 R^2-R-2\right) - \big(2 R^2 / \cos^2 \theta\big) \Big), \\ B(\theta) &= \frac{2(\Omega+\omega) K}{(R+1)(R+2)} \cos ^R \theta\big((R^2+2 R+2) - (R+1)^2 \cos^2 \theta\big), \\ C(\theta) &= \frac{1}{4} K^2 \cos^{2 R} \theta\big((R+1) \cos^2 \theta-(R+2)\big), \end{aligned}\]

the initial height field is given by

\[h(\lambda,\theta) = h_0 + \frac{a^2}{g}\Big(A(\theta) + B(\theta)\cos(R\lambda) + C(\theta)\cos(2R\lambda) \Big).\]

This problem corresponds to Case 6 of the test suite described in the following paper:

• D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber (1992). A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224. DOI: 10.1016/S0021-9991(05)80016-6

initial_condition_unsteady_solid_body_rotation(x, t, equations)

Unsteady analytical solution to the spherical shallow water equations, corresponding to a solid body rotation with a prescribed bottom topography. Assuming the domain to be a sphere of radius $a = 6.37122 \times 10^3\ \mathrm{m}$, letting $\vec{x}$ denote the position relative to the centre of the Earth, and letting $\vec{e}_x$, $\vec{e}_y$, and $\vec{e}_z$ denote the Cartesian basis vectors, we define the rotating frame

\[\vec{b}_x(t) = \cos(\Omega t)\vec{e}_x + \sin(\Omega t)\vec{e}_y, \quad \vec{b}_y(t) = -\sin(\Omega t)\vec{e}_x + \cos(\Omega t)\vec{e}_y, \quad \vec{b}_z(t) = \vec{e}_z,\]

as a function of time $t$. We also define the associated coordinate transformation

\[\vec{\varphi}(\vec{x},t) = (\vec{x} \cdot \vec{b}_x(t)) \vec{e}_x + (\vec{x} \cdot \vec{b}_y(t)) \vec{e}_y + (\vec{x} \cdot \vec{b}_z(t)) \vec{e}_z\]

as well as a fixed axis $\vec{c} = -\sin(\alpha)\vec{e}_x + \cos(\alpha)\vec{e}_y$, where $\Omega = 7.292 \times 10^{-5}\ \mathrm{s}^{-1}$ is the Earth's rotation rate, and we take $\alpha = \pi/4$. For a bottom topography prescribed as

\[h_s(\vec{x}) = \frac{1}{2g}(\vec{\Omega} \cdot \vec{x})^2\]

where $\vec{\Omega} = \Omega\vec{e}_z$ and $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$ are the Earth's axis of rotation and gravitational acceleration, respectively, the time-dependent velocity field is given as

\[\vec{v}(\vec{x},t) = v_0\, \vec{\varphi}(\vec{c},t) \times \vec{x}/\lVert \vec{x} \rVert,\]

and the total geopotential height $H = h+b$ is given by

\[H(\vec{x},t) = \frac{1}{2g}\left(\big(v_0 \,\vec{\Omega}\cdot \vec{x} - \vec{\varphi}(\vec{c},t) \cdot \vec{x}/\lVert \vec{x} \rVert \big)^2 + (\vec{\Omega} \cdot \vec{x})^2 + 2k_1\right),\]

where $v_0 = 2\pi a / (12 \ \mathrm{days})$, $k_1 = 133681 \ \mathrm{m}^2/\mathrm{s}^2$, and $\lVert \cdot \rVert$ denotes the Euclidean norm. To use this test case with SplitCovariantShallowWaterEquations2D, the keyword argument auxiliary_field = bottom_topography_unsteady_solid_body_rotation should be passed into the SemidiscretizationHyperbolic constructor. This analytical solution was derived in the following paper:

• M. Läuter, D. Handorf, and K. Dethloff (2005). Unsteady analytical solutions of the spherical shallow water equations. Journal of Computational Physics 210:535–553. DOI: 10.1016/j.jcp.2005.04.022

source_terms_coriolis(u, du, x, t, 
                      equations::ShallowWaterEquations3D,
                      normal_direction)

Source term function to apply the Coriolis force with an angular velocity of equations.rotation_rate around the $z$ axis.

The vector normal_direction is perpendicular to the 2D manifold. By default, it is the normal contravariant basis vector.

source_terms_coriolis_lagrange_multiplier(u, du, x, t,
                                          equations::ShallowWaterEquations3D,
                                          normal_direction)

Computes the Coriolis source term (source_terms_coriolis) and the Lagrange multiplier source term (source_terms_lagrange_multiplier).

source_terms_lagrange_multiplier(u, du, x, t,
                                 equations::ShallowWaterEquations3D,
                                 normal_direction)

Source term function to apply a Lagrange multiplier to the semi-discretization in order to constrain the momentum to a 2D manifold.

The vector normal_direction is perpendicular to the 2D manifold. By default, it is the normal contravariant basis vector.

transform_initial_condition(initial_condition, equations)

Takes in a function with the signature initial_condition(x, t, equations) which returns an initial condition given in terms of global Cartesian or zonal/meridional velocity components, and returns another function initial_condition_transformed(x, t, equations) or initial_condition_transformed(x, t, aux_vars, equations) which returns the same initial data, but transformed to the appropriate prognostic variables used internally by the solver. For the covariant form, this involves a transformation of the global velocity components to contravariant components using global2contravariant as well as a conversion from primitive to conservative variables. For standard Cartesian formulations, this simply involves a conversion from primitive to conservative variables. The intention here is to have a set of test cases (for example, initial_condition_gaussian) for which the initial condition is prescribed using a standardized set of primitive variables in a global coordinate system, and transformed to the specific prognostic variables required for a given model.