Investigating Inherent Numerical Stabilization for the Moist, Compressible, Non-Hydrostatic Euler Equations on Collocated Grids

Abstract

This study investigates inherent numerical dissipation due to upwind fluxes and reconstruction strategies for collocated Finite-Volume integration of the Euler equations. Idealized supercell simulations are used without any explicit dissipation. Flux terms are split into: mass flux, pressure, and advected quantities. They are computed with the following upwind strategies: central, advectively upwind, and acoustically upwind. This is performed for third and ninth-order-accurate reconstructions with and without Weighted Essentially Non-Oscillatory limiting. Acoustic-only upwinding for pressure and mass flux terms and advective-only upwinding for advected quantities is the most flexible simplification found. It reduces data movement and computations. Assuming a constant speed of sound in acoustic upwinding gives similar results to using the true speed of sound. Dissipation from upwind adapts automatically to grid spacing, time step, reconstruction accuracy, and flow smoothness. While stability is maintained even at 21st-order spatial accuracy, there is a limit to the spatial order of accuracy for which upwinding alone can create a realizable solution in the conditions of this study. Convex combinations of upwind and central solutions for flux terms also reduced dissipation, but as the central proportion grows, solutions become physically unrealizable. The range of length scales of the kinetic energy spectra can be extended along k^−5/3 to smaller spatial scales by reducing dissipation either with higher-order reconstructions or using convex combinations of upwind and central fluxes. However, not all extensions of the length scale range along k^−5/3 exhibit physically realizable solutions, even though the spectra appear to be physical.

Introduction

Viscosity in the atmosphere is small. Therefore, viscous dissipation does not dominate small, resolved scales for the grid spacings used for atmospheric models performing numerical weather prediction and climate simulations. For this reason, the inviscid Euler equations have a physical motivation when seeking a numerical approximation to the fluid dynamics of the gas-only portion of these models—often labeled as the “dynamical core” or “dycore.” It is important to note that though inherent numerical dissipation can approximate some portions sub-grid-scale mixing, it cannot account for dynamics that involve interactions between phase change, sub-grid-scale eddy dynamics, and/or radiation (Jablonowski & Williamson, 2011; Pressel et al., 2017).

All numerical approximations to the Euler equations need some form of dissipation or stabilization in order to avoid accumulation of energy at the smallest scales, especially in 3-D where the energy cascade involves transfer of energy downscale. This is because viscosity is not an inherent part of the Euler equation set, and this stabilization can be thought of as being related to accounting for a portion of sub-grid-scale eddy dynamics that eventually dissipate Kinetic Energy (KE) at Kolmogorov microscales. This provides a link between inherent numerical dissipation/stabilization of the Euler equations and Implicit Large Eddy Simulation approaches (Margolin & Wyszogrodzki, 2007).

Stabilization can come from several parts of the numerical discretization, including (but certainly not limited to):

• Grid staggering choices (Fox-Rabinovitz, 1994; Staniforth & Thuburn, 2012);

• Geometric structure preserving approaches arising from the use of exterior calculus and quasi-Hamiltonian formulations (Cotter & Shipton, 2012; Eldred & Randall, 2017; Eldred et al., 2019; Gassner et al., 2016; Waruszewski et al., 2022);

• Explicit post-hoc dissipation in the forms of viscosity, hyper-viscosity, divergence damping, sub-gird-scale eddy models and/or averaged turbulence models (Dennis et al., 2012; Jablonowski & Williamson, 2011; Santiago & Martilli, 2010; Stoll et al., 2020);

• Dissipative interpolations such as monotone (including first-order interpolation), low-order, odd-order, and non-oscillatory approaches (Colella & Woodward, 1984; Liu et al., 1994); and

• Preferring upwind information either when interpolating or when reconciling disparate interpolated state or flux vectors at the same spatial location (e.g., with a Riemann solver) (Godunov & Bohachevsky, 1959; Roe, 1981).

Many grid staggering as well as geometric structure preserving approaches (these are related but, to a degree, separable strategies) are engineered to avoid certain unstable linear modes, providing inherent stabilization at least with respect to those modes (Eldred & Le Roux, 2018, 2019; Eldred & Randall, 2017; Eldred et al., 2019; Staniforth & Thuburn, 2012). Additionally, they can provide stability through the preservation of global solution invariants such as total energy and total potential enstrophy (Eldred & Randall, 2017) in an (implied) flux-form. Collocated grids/“A-grids” are grids for which all model variables are expressed in the same physical location for a given grid point or cell, and these can be expressed as point values, cell averages, or element-based spectral representations.

There are pros and cons for staggered and collocated grids regarding accuracy, stability, and computational efficiency in various contexts; but the goal of this paper is not to investigate that comparison. Rather, the goal is to better understand the stabilization effects on collocated grids that come from two sources in particular: (a) preferring upwind information for cell-edge fluxes; and (b) dissipative interpolation.

Naive discretizations on collocated grids typically suffer from unstable linear modes that require some means of dissipation, whether implicit to the numerical discretization or explicitly added after the fact. This is not to say that modes do not also exist in some staggered grid applications, for example, (Adcroft et al., 1999). For many collocated grid approaches, some means of preferring upwind information in the discretization is used to provide stabilization, and this is most commonly done when computing cell-edge fluxes in a Finite-Volume (FV) approach (Godunov & Bohachevsky, 1959; Roe, 1981).

Most element-based methods, particularly those with weak-form spatial integrals like Discontinuous Galerkin, can use upwind information at element edges as well (Cheruvu et al., 2007; Giraldo & Restelli, 2008; M. R. Norman, 2015; M. R. Norman & Finkel, 2012). However, element-based methods (especially at high spatial orders of accuracy) derive relatively little of their inherent numerical dissipation from element edge state estimates, meaning they can use relatively diffusive Riemann solvers without consequence (Nair et al., 2005). This is because the discrete update of element-based methods, especially at high-order, is typically dominated by a spatial integral over the element's internal Degrees Of Freedom (e.g., a “stiffness matrix”) that usually has no notion of upwind information. Therefore, FV methods are of particular interest here because in the absence of source terms, the entire FV update is formed by cell-edge fluxes, each of which feels the influence of upwind information present in fluxes.

There are some attractive properties of preferring upwind information in FV fluxes that motivate further investigation. While this is not a unique observation (Jablonowski & Williamson, 2011), it will be shown that preferring upwind information for cell-edge fluxes leads to a dissipation mechanism that automatically responds to spatial order of accuracy, grid spacing, time step size, and local flow smoothness. Also, there is nuance to the idea of “upwind.” For instance, do all variables in the flux vector use upwind information, or do only some variables use upwind information? For a given variable, is the upwind preference relative to advective dynamics, acoustic dynamics, or both?

It is also important to note that this study is confined to subsonic flows such as those observed in typical weather and climate modeling. With supersonic flows, subtle choices regarding how to introduce upwind preferences in the discretization can have profound effects such as the issue with entropy violation in rarefactions in the original Roe scheme and other solution artifacts at hypersonic speeds (Qu et al., 2021). However, with subsonic flows, the impacts of subtle upwind choices do not impact the solution nearly as much. For instance, how a flux Jacobian is averaged at a cell interface or which variables are reconstructed (state, primitive, or characteristic variables) have relatively small impacts on the physical realism of the solution at subsonic speeds.

There are numerous spatially collocated discretizations available for the compressible, non-hydrostatic Euler equations (Ahmad & Lindeman, 2007; Dennis et al., 2012; Giraldo & Restelli, 2008; M. R. Norman, 2021; M. R. Norman et al., 2011). This study uses a FV discretization cast in conservation form (or flux form) on a collocated/A-grid.

While ideas of stabilization and dissipation can be hard to separate in practice, in this study emphasis is placed on the notion of stabilization and the concept of producing a physically realizable flow, which is easier to define than the concept of a realistic flow, particularly in idealized settings and with chaotic flows in which initially small differences diverge quickly to disparate solutions. This study investigates how inherent numerical dissipation mechanisms stabilize the flow to create a realizable solution. A realizable solution is defined by the following constraints:

1. The solution is free from floating point values of NaN (Not a Number), infinity, and negative infinity;

2. The solution does not produce obviously unrealistic wind speeds. In this case, we will determine 150 m per second or greater to be physically unrealizable (in the context of idealized supercell flow);

3. The solution does not contain permanent and global oscillations in any coordinate direction for any variable (meaning transient, local oscillations responding to discontinuities are permitted). Specifically, this entails 2Δx oscillations that span an entire dimension and do not abate over time. These show up in KE spectra plots as well as large peaks at small wavelengths;

4. Tracer, density, and thermodynamic quantities remain non-negative.

5. The solution does not a contain significant degree of gridpoint scale noise (often evident in vertical velocity) located far away from convection cells.

   • (a)

     Constraint number five, above, is the most subjective of the constraints, but there will be solutions in this study that clearly show a level of gridscale noise in vertical velocity that cannot be realized by the physics of fluid flow. As a clear example of this last constraint being violated, consider comparing Figure 7a against the ground truth high-resolution simulation of Figure 7c. Figure 7a exhibits gridpoint-scale convection regions far away from where convection should be physically happening.

There are five hypotheses being evaluated in this study:

1. H1: There are simplifications to the fully upwind Riemann solver of (M. R. Norman, 2021) (and fully upwind solvers in general) that result in physically realizable solutions;

2. H2: Upwind dissipation has different effects for different types of terms in the flux vector;

3. H3: There are noticeable differences in preferring upwind information according to advective velocities versus acoustic dynamics;

4. H4: There is an upper limit to spatial order of accuracy that can produce physically realizable solutions with preferring upwind information in fluxes being the sole means of stabilization; and

5. H5: Physically realizable discretizations can be obtained by convex combinations of upwind estimates and central (non-upwind) estimates of flux terms.

The Moist, Compressible, Non-Hydrostatic Euler Equations

In this work, we will study a commonly-used equation set for atmospheric models: the moist compressible (non-hydrostatic) Euler equations using the constant κ approximation (Eldred et al., 2022), simplified here by assuming a flat, doubly periodic domain with no rotation and a simple geopotential ϕ = gz. These represent a moist extension of the equations used in (M. R. Norman, 2021): [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]where ρ is total density; u, v, and w are the wind velocities in the x, y, and z directions, respectively; θ is the virtual potential temperature, ρ_v is water vapor density, ρ_ψ is the density of a tracer substance that adds mass but due to the constant κ approximation does not contribute to the thermodynamics (such as precipitated species or cloud liquid); ρ_d is the density of dry air; p is total pressure assumed to be the sum of dry air and vapor pressures; T is the air temperature, R_d and R_v are dry air and water vapor ideal gas constants, γ = c_p,d/c_v,d ≈ 1.4, is the ratio of specific heats for dry air, and $${C}_{0}={\left({R}_{d}{p}_{0}^{-{R}_{d}/{c}_{p,d}}\right)}^{\gamma }$$; p₀ = 10⁵ Pa; g = 9.8 m s⁻² is acceleration due to gravity; c_p,d is the specific heat of dry air at constant pressure; and c_v,d is the specific heat of dry air at constant volume. It is important to note that the virtual potential temperature used here is not same as the (moist) potential temperature. In this equation set, microphysical transitions between moist phases as well as fallout and aggregation of precipitated species is assumed to be operator split from the dynamics so that there are no latent heating or precipitation fallout/aggregation source terms while integrating the moist Euler equations.

The key feature of this equation set is the constant κ approximation, which decouples trace species such as ρ_v and ρ_ψ from the thermodynamics. This allows pressure to be expressed solely in terms of total density ρ and virtual potential temperature θ using the dry ideal gas constants. Further, full mass loading is used in this equations set, meaning momentum feels the effects of tracer constituents that do not contribute to total pressure in this model but do have mass contribution. Immersed and precipitated liquid and ice species are assumed to have no contributions to total pressure.

The equation of state for pressure can be recast in terms of its differential as follows: [Image Omitted. See PDF]where $${c}_{s}=\sqrt{\gamma p/\rho }$$ is the speed of acoustic propagation. This differential form of the equation of state will be used in several places later in this study.

Defining Upwind Notions via the Diagonalization

When only considering advective dynamics, the notion of “upwind” direction is trivially defined as the opposite direction of the wind velocity, $$\vec{u}$$. However, acoustics propagate on the boundary of an expanding distorted sphere in 3-D, meaning upwind dynamics technically take in all spatial directions at the same time. In this case, it is less clear how to define what the upwind direction means and how to compute an upwind state. Since this study uses dimensionally split dynamics, in 1-D, acoustic dynamics propagate in both directions at the same time (at the boundaries of an expanding line). Therefore, before computing an upwind state, the hyperbolic equation set must first be cast into its “characteristic” form, which describes a set of waves traveling with finite and (generally) non-zero velocities. In characteristic form, the notion of upwind once again becomes trivial to define.

Full Diagonalization

One can recast the Equation set (1) in a slightly different form: [Image Omitted. See PDF]where $$\mathbf{f}={\left[\rho u,\rho {u}^{2}+p,\rho uv,\rho uw,\rho u\theta \right]}^{\top }$$, $$\mathbf{g}={\left[\rho v,\rho vu,\rho {v}^{2}+p,\rho vw,\rho v\theta \right]}^{\top }$$, and $$\mathbf{h}={\left[\rho w,\rho wu,\rho wv,\rho {w}^{2}+p,\rho w\theta \right]}^{\top }$$.

In this form, one can use the derivative chain rule to explore the hyperbolic behavior of the system in each spatial direction. Considering only the x-direction for now and neglecting the source terms gives: [Image Omitted. See PDF]where A = ∂f/∂q is the “flux Jacobian” matrix. One can diagonalize the matrix: A = RΛL, where R is the matrix whose columns are right eigenvectors, Λ is a diagonal matrix whose diagonal components are eigenvalues, and L = R⁻¹ is the matrix whose rows are left eigenvectors. Left-multiplying the equation set by L and assuming the flux Jacobian is locally constant in time and uniform in space gives (also known as “locally frozen”): [Image Omitted. See PDF]where $$\mathbf{w}=L\mathbf{q}={\left[{w}_{1},{w}_{2},{w}_{3},{w}_{4},{w}_{5}\right]}^{\top }$$ are scalar “characteristic variables,” and λ_p is the pth diagonal component of Λ. These describe a decoupled set of transport equations, transporting characteristic variables, w_p, at finite velocities, λ_p, making the equation set “hyperbolic.” In the x-direction, the Jacobian and its diagonalization for the main dynamics variables ρ, ρu, ρv, ρw, and ρθ are: [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

From the eigenvalues, there are three waves traveling at velocity u, one traveling at velocity u − c_s, and one traveling at velocity u + c_s, where $${c}_{s}=\sqrt{\gamma p/\rho }$$ is the speed of sound. The two waves containing the speed of sound have dynamics that mix advective and acoustic dynamics, and it can be seen that these incorporate acoustics traveling in both directions by decomposing acoustics into separate modes, each with their own velocity. The upwind direction for each wave is the negative of the eigenvalue or wave velocity, λ_p.

Defining the Upwind State

Based on the characteristic form of the equation set, one can then define an upwind state at a cell interface based on the left- and right-hand estimates of the state from interpolation centered on the two adjacent cells. These states are denoted q⁻ and q⁺, respectively. For each wave, we compute the upwind characteristic variable using the characteristic velocity, λ_p: [Image Omitted. See PDF]where ℓ_p is the pth row of L. The vector of upwind characteristic variables, w_upw, is transformed back into an upwind state vector using the right Eigen-matrix: q_upw = Rw_upw. In this study, A = RΛL is computed using $$\left({\mathbf{q}}^{-}+{\mathbf{q}}^{+}\right)/2$$.

Separating Out Acoustics

It is possible to separate the flux Jacobian into two matrices that, when summed together, give the original Jacobian. There is not a unique solution for this, but one approach that yields physically intuitive results is to construct a fast, acoustic Jacobian, A_acoust, by setting u → 0, v → 0, and w → 0 in A.

Fast, Acoustic Jacobian

Setting u → 0, v → 0, and w → 0 in A yields the following Jacobian matrix and associated Eigenvalues/characteristic wave velocities: [Image Omitted. See PDF]

From this, it is clear that the resulting dynamics propagate at purely acoustic speeds. It is also clear that when including all dynamical variables, this system is only weakly hyperbolic, meaning some of the waves have zero speed for all fluid states. However, by writing out the associated equations from the acoustic-only flux Jacobian and using the implications of that equation set, one can reduce the purely acoustic equations to only two variables. [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

The simplification in Equation 15 is due to the equation of state's differential form in Equation 4. Equation 16 tells us that potential temperature does not change in time during purely acoustic dynamics. Using the fact that during purely acoustic motions, θ does not change in time, one can take the equation of state's differential with respect to time to get: [Image Omitted. See PDF]

Combining this relationship with the conservation of mass, the acoustic dynamics can be simplified to a system of only two variables: momentum and pressure: [Image Omitted. See PDF]

This gives a strongly hyperbolic and purely acoustic diagonalization: [Image Omitted. See PDF]

What this simple system provides is a convenient way of computing a purely acoustically upwind state for pressure and momentum, which makes it easy to test the effects of different types of upwinding for different terms in the flux vector. This system matches physical intuition that acoustics are a fast adjustment between pressure and density mediated by momentum. Also, equation set (18) can be further simplified (assuming a constant speed of sound with respect to time for the moment) to the second-order wave equation for pressure: $${\partial }_{tt}p={c}_{s}^{2}{\partial }_{xx}p$$, demonstrating that this represents simple linear acoustic propagation.

It is also worth noting that Equation 18 is also identical to applying a “relaxation technique” to acoustic-only equation for normal momentum, $${\partial }_{t}\left(\rho u\right)+{\partial }_{x}p=0$$, assuming a maximum eigenvalue that is the speed of sound, c_s (Banda & Seaid, 2005). The relaxation approach could be thought of as a generalized hyperbolic equation set that gives rise to what are essentially Lax-Friedrichs or sometimes Local Lax-Friedrichs (LLF) fluxes. Solving for upwind fluxes using equation set (18) verifies that this system uses LLF fluxes of the form: [Image Omitted. See PDF]where $$\mathbf{q}={\left[p,\rho u\right]}^{\top }$$ and $$\mathbf{f}={\left[{c}_{s}^{2}\rho u,p\right]}^{\top }$$, a “+” superscript denotes the interpolation at the cell interface from the cell on the right, and a “−” superscript denotes the interpolation at the cell interface from the cell on the left.

The first term of Equation 20 in involves a straightforward central estimate of the flux with no upwind preference. The second term of Equation 20 has the form of a second-order viscosity term cast as an interface flux, but with notable differences from a standard second-order viscosity term. First, the coefficient of viscosity automatically scales proportionally to $${\left({\Delta }x\right)}^{2}{\Delta }{t}^{-1}$$ such that it automatically responds to changes in grid spacing and time step size. Second, it is responding to interpolations of the state at the cell interface rather than the cell-averaged state in adjacent cells (the latter being equivalent to first-order interpolation). This, in turn, has two effects. The magnitude of viscosity decreases as the accuracy of the interpolation increases because, in general, higher-order, and more accurate interpolation leads to smaller differences in state estimates at cell edges. Also, the magnitude of viscosity decreases as the smoothness of the flow increases from the same argument. Therefore, the preference for upwind information at cell interfaces provides an adaptive dissipation mechanism that responds to the accuracy of interpolation, the smoothness of the flow, the grid spacing, and the time step.

Yet a further possible simplification to equation set (18) is that one could assume a speed of sound that is constant in time and either uniform in space or a time-constant function of height only. The only simplification tested here is a constant speed of sound in some of the experiments to determine its effects on stabilization in the overall scheme.

Slow, Advective Dynamics

For slow, advective dynamics, all information propagates in the direction of the fluid velocity u, v, and w in the x, y, and z directions, respectively. Therefore, determining the upwind state in terms of advective dynamics is trivial. In fact, it does not actually need to be performed in a Riemann solver separate from reconstruction/interpolation. It is trivial to show that shifting a stencil with an odd number of cell averages for constraints (and therefore presumably a reconstruction with an odd order of accuracy) to be centered about the upwind cell adjacent to the interface yields the same results as computing two estimates at a cell interface and choosing the upwind value.

Connections to Existing Riemann Solvers

The idea of splitting a Jacobian into fast and slow portions in order to handle them separately has been explored in various contexts in literature before. It is commonly done for semi-implicit or implicit-explicit (IMEX) integration of low-Mach flows where advective dynamics are performed time-explicitly and acoustic dynamics are performed time-implicitly (Gardner et al., 2018; Giraldo et al., 2013; Weller et al., 2013). The means of splitting fast and slow equation terms is varied in the literature. It is also common in low-Mach and all-Mach Riemann solvers (Chen et al., 2013; Li & Chen, 2019; Liou, 2006). Typically this involves controlling dissipation differently for acoustic and advective terms, and that has a relationship with the convex combinations of upwind and central fluxes, particularly when performed only for advective terms, in Section 3.4.

Finite-Volume Discretization With Cell-Centered Reconstructions

The equations are discretized using dimensional splitting in time, meaning the equation set can be reduced to a generalized 1-D set of conservation laws of the form: [Image Omitted. See PDF]

This is integrated spatially over a local cell domain $${{\Omega }}_{i}\in \left[{x}_{i-1/2},{x}_{i+1/2}\right]$$ where x_i±1/2 = x_i ± Δx_i/2: [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]where $${\widetilde{\mathbf{f}}}_{i}\left(x,t\right)$$ and $${\widetilde{\mathbf{s}}}_{i}\left(x,t\right)$$ are spatial reconstructions of the vectors f and s, respectively, over the ith cell domain, Ω_i, (the time dimension is left continuous for now), and $$\mathcal{R}\left({\mathbf{f}}^{-},{\mathbf{f}}^{+}\right)$$ is a Riemann solver to reconcile the multi-valued flux at a cell interface due to samples from disparate reconstructions at the same location in space. The cells fully span the spatial domain without overlap.

The variables explicitly reconstructed include ρ′, u, v, w, and $${\left(\rho \theta \right)}^{\prime }$$, where a primed quantity denotes perturbation from a hydrostatic background profile. Then, hydrostasis is added at interpolation locations after reconstructions to give ρ, u, v, w, and ρθ at cell edges. MPI exchanges of cell-edge reconstructions are performed, and then the Riemann solver is performed. The Jacobians (and thus the upwind directions) are determined as the linear average of the two states at the cell edge. The interface averaged variables ρ, u, v, w, and θ are computed first. Then, pressure and speed of sound are computed from those averaged variables to provide the Jacobian diagonalizations needed for upwinding.

In this formulation, the inherent dissipative properties of the scheme are largely confined to how the fluxes are computed, which mainly involves the choice of spatial reconstructions and the choice of Riemann solver. Recall that time is left continuous for now. Since upwinding is a primary object of investigation in this study, focus is placed there. The flux vector will be broken into three terms: (a) the pressure; (b) the direction-normal mass flux; and (c) the advected quantities. Each of these can use: (a) a centrally-averaged value with no upwinding; (b) upwinding according to the direction-normal fluid velocity (i.e., advectively upwinded); (c) upwinding according to purely acoustic dynamics in equation set (18); and (d) upwinding according to the full, native flux Jacobian's diagonalization accounting for both acoustic and advective waves. While this, at face value, would lead to many possible combinations, these can be reduced in practice because most of the configurations lead to either unstable or clearly unphysical solutions.

This decomposition of the flux terms is motivated by the flux Jacobian's diagonalization, which decomposes well into acoustically dominated waves, which primarily involve the mass fluxes and pressure terms; and advectively dominated waves, which primarily involve the advected quantities. A further constraint on how one can decompose the flux vector into distinct types of terms is that the same mass flux should ideally be shared among all equations in the model to maintain mass consistency with how density is evolved. That way, when dividing by density to obtain the desired quantity from its mass-weighted value, a consistent value is obtained.

Reconstruction Strategy

Variation within cells is reconstructed using Vandermonde interpolation of the surrounding cell averages. Odd-order-accurate interpolation is natural for cell-centered reconstructions because directionally unbiased, cell-centered stencil sizes are naturally odd. Odd-ordered interpolations have a hyper-viscous term leading the error term, which results in a small amount of dissipation in the resulting interpolation—the strength of which depends upon the order of accuracy. After reconstructing variation, point values of the state and tracers are interpolated at the cell edges for input into the Riemann solver. Interpolations are performed with and without Weighted Essentially Non-Oscillatory (WENO) limiting, and the WENO limiting used is from (M. R. Norman, 2021) for only the scalar terms (density, potential temperature, and tracers).

Momenta interpolations are left unlimited, which is often done in weather and climate simulations. This is because many if not most physics parameterizations do not involve wind velocities, and those that do tend to be dissipative in nature such as sub-grid-scale turbulence parameterizations and surface layer mixing. Therefore, it is of greater importance to limit scalars because oscillations and positivity violations in scalars tend to create poor results from physical parameterizations. Further, since this is solely subsonic flow, it is not necessary to handle phenomena such as shocks in the velocities with limiting.

For accuracy, only perturbations from a hydrostatic background state are interpolated for density, ρ, and for mass-weighted potential temperature, ρθ. An analytically derived background hydrostatic profile is subtracted from cell averages and added to state estimates at the cell edges. Also, density is divided out of momenta and tracer quantities prior to interpolation and multiplied back afterward to avoid the influence of the non-linear decay of density with increasing height while interpolating with polynomials.

Time Stepping

An optimal time-explicit third-order, three stage Strong Stability Preserving Runge-Kutta (SSPRK3) time integrator (Gottlieb et al., 2009) is used with a maximum CFL value of 0.6 in all tests with all spatial dimensionalities. In practice, the SSPRK3 time integrator retains its advertised SSP CFL value of one only when the spatial operator is first-order accurate. When the spatial operator is higher-order accurate, numerical experiments performed by the authors in a 1-D advection context show that the maximum SSP CFL value of the combined space-time operator typically drops to two-thirds of the time stepper value when used with a first-order spatial scheme (i.e., the values given in literature), at least in the context of the spatial discretizations used here. This time stepping choice introduces a measure of dissipation in itself. Further, the time step is calculated assuming a maximum wave speed of 350 + 80 m per second for estimated acoustic and wind speeds, respectively.

To achieve dimensional splitting, the x, y, and z-direction solves are all computed independently inside each Runge-Kutta stage. This has been shown, even in non-orthogonal coordinates, to introduce relatively little error, especially in spatially less smooth regimes that are common in weather and climate forecasting (Katta et al., 2015).

The Flux Terms Being Tested

The test code computes the three types of flux terms (pressure, normal mass fluxes, and advected quantities shown in green, blue, and red below, respectively) using each of the four aforementioned upwinding techniques: [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

To recall, the four methods of flux calculation from two estimates at the same cell edges are:

1. Average of left and right states (central)

2. Upwind according to advective wind velocity (advec)

3. Upwind according to the full, native Jacobian (full)

4. Upwind according to acoustics only (acoust)

   • The fourth treatment above is applicable to only the mass flux and pressure terms, while the other three treatments are applicable to all terms

Numerical Experiments

All numerical experiments in this study use an idealized supercell initialization to create challenging flows with a combination of convectively unstable (in the troposphere) and convectively stable (at and above the tropopause) vertical layers and strong wind velocities and shears. Note that no explicit forms of dissipation are used in any of these numerical experiments, meaning there is no explicit viscosity, hyperviscosity, divergence damping, sponge layer, radiating boundary conditions, or any other treatment that could be interpreted as explicit numerical damping. All dynamics generated by the model (including acoustic and gravity waves generated by unbalanced initial conditions) are trapped inside the domain with only inherent numerical dissipation used to maintain stability. That is not to say that this treatment is necessarily a physical choice but that emphasis is desired on the stabilization properties of upwind preference and dissipative interpolation. Removing all explicit dissipation strategies is a helpful way to investigate this in a more isolated context.

This initialization is described in Section 3.1. Simple “Kessler” microphysics (specifically, that used in Zarzycki et al. (2016)) are used in this test case, which evolve water vapor, cloud liquid, and precipitation liquid. Kessler microphysics involves a saturation adjustment between water vapor and cloud liquid along with appropriate latent heating and cooling as well as precipitation auto-conversion from cloud liquid and precipitation fallout. The Kessler code provided (https://github.com/ClimateGlobalChange/DCMIP2016/blob/master/interface/kessler.f90) by the Dynamical Core Model Inter-comparison Project (DCMIP) was refactored into portable C++ for use in this study. The entire model is GPU-capable, written using a portable C++ call Yet Another Kernel Launcher (M. Norman et al., 2022).

There are four sets of numerical experiments performed in this study:

1. Investigating the different upwind strategies for the three different flux terms (Hypotheses H1, H2, and H3)

2. Investigating the effects of assuming a time-constant speed of sound for acoustically upwind pressure and mass flux (Hypothesis H1);

3. Investigating the effects of using convex combinations of upwind and central fluxes (Hypothesis H5); and

4. Investigating the upper bounds of spatial order that produce physically realizable simulations with upwind dissipation (Hypothesis H4).

For many of these investigations, KE power spectra will be plotted as a function of spatial wavenumber. In this study, the KE is not mass weighted, meaning it is computed as: [Image Omitted. See PDF]

The spectral power is computed as the square of the absolute value of each complex coefficient of the forward real-to-complex 2-D FFT computed on the x-y plane and averaged over the vertical levels. The spectra are also averaged over the last 8 hr of simulation at half hour increments. The wavenumber is computed as $$k=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}$$, and spectra with identical wavenumber are averaged for that wavenumber. Finally, spectra and corresponding wavenumbers are binned into equally spaced bins, and the mean of those bins are used for plotting.

Initializing an Idealized Supercell Test Case

The idealized supercell data used for initialization in this study is very similar to that which is defined in Klemp et al. (2015), Zarzycki et al. (2016) but for planar geometry instead of spherical geometry, with no Coriolis force, and with a slightly simplified implementation procedure. A model domain of 100 × 100 × 20 km is used, a tropopause height of z_trop = 12 km is assumed, and a temperature profile is assumed to be a linear interpolation between $$T\left(z=0\right)=300\,\mathrm{K}$$, $$T\left(z={z}_{\mathit{trop}}\right)=213\,\mathrm{K}$$, and $$T\left(z=20\,\text{km}\right)=213\,\mathrm{K}$$. The temperature is initialized as constant above the tropopause. A surface pressure of 10⁵ Pa is assumed. The relative humidity is assumed to be: [Image Omitted. See PDF]with the additional constraint that the dry mixing ratio is capped at a maximum value of $${q}_{v,\max }=\max \left({\rho }_{v}/{\rho }_{d}\right)=0.014$$ (the effects of which are applicable to the surface boundary layer). Dry mixing ratio is calculated from saturation mixing ratio and relative humidity using the following definition for saturation mixing ratio: [Image Omitted. See PDF]

Using this information, one can combine the equation of state and hydrostasis to determine a hydrostatically balanced pressure profile: [Image Omitted. See PDF] [Image Omitted. See PDF]

This implies: [Image Omitted. See PDF]

This equation is integrated using quadrature within small intervals designed to create quadrature point values of pressure useful for further integration needed to compute initial cell averages.

To be more clear about the initialization procedure, please refer to the following algorithm:

• For each cell

  • –

    For each interval between GLL quadrature points within the cell

    • ∗

      For each GLL quadrature point within the interval:

    • - Compute specified temperature, T

    • - Compute dry pressure, p_d based on an exact integration of the combined equation of state and hydrostatic relationship for dry air:

      [Image Omitted. See PDF]

    • - Compute saturation mixing ratio based on temperature, T, and dry pressure, p_d. Note that this involves a simplification that saturation mixing ratio is relative to dry pressure rather than full pressure, but it produces a very similar vertical profile as the original test case. Moist addition to pressure is only of the order of 1%.

    • - Compute specified relative humidity, RH

    • Limit the RH if necessary to keep the mixing ratio capped at a maximum value of 0.014

    • - Compute and store the right-hand-side of $$\left(33\right)$$ in preparation for integrating full pressure via quadrature

• –

  Integrate moist pressure to the next GLL quadrature point within the cell using sub-quadrature points in this interval

For each cell

• –

  For each quadrature point within the cell

  • ∗

    Use specified temperature, relative humidity, and integrated hydrostatic moist pressure to compute mixing ratio, dry density, water vapor density, and potential temperature

• –

  Integrate quadrature points within the cell to compute an integrated initial average value for the cell for each of the state variables

Beyond this hydrostatically balanced, convectively unstable background profile initialization, a linear vertical profile for x-direction wind velocity, u, from −15 m s⁻¹ at the surface to 15 m s⁻¹ at z = 5 km is initialized. Wind in the y and z directions is initialized to zero. Also, an initial thermal perturbation is added to realize the convective instability of the initial profile via an initial vertical lift. This perturbation is defined as: [Image Omitted. See PDF] [Image Omitted. See PDF]and $$\left({x}_{0},{y}_{0},{z}_{0}\right)=\left(50\,\text{km},50\,\text{km},1.5\,\text{km}\right)$$ and $$\left({x}_{R},{y}_{R},{z}_{R}\right)=\left(10\,\text{km},10\,\text{km},1.5\,\text{km}\right)$$.

This will initiate convection at the horizontal domain center that then splits into two cells that propagate north and south in the y dimension and propagate slowly in the positive x-direction. Traditionally, this is simulated for two model hours, but in many simulations in this study, it will be simulated for one model day. Figure 1 shows volume renderings of clouds, precipitation, and 5 km vertical velocity after 30, 60, 90, and 120 min of simulation using 1,000 × 1,000 × 200 cells.

Maintaining Convective Instability

One further addition is made during model simulation. At each model time step, the model is nudged equally for all model columns such that the model average column maintains the same profile that was initialized. Since this nudges the overall model mean column, this does not affect variability between columns or perturbations from the background state toward which it is being nudged. The nudging is performed with a time scale of 900 s. What this does, in effect, is maintains the convective instability that the supercell convection would typically stabilize and forces the model into a convective equilibrium wherein the instability never dissipates. Note that this does not produce a damping effect in any sense. Rather, it continually excites instability in the model domain to maintain strong convection during the entire simulation. However, even with the maintenance of instability, the results look similar to existing supercell test cases that do not use this nudging approach after two model hours.

Investigating Upwind Flux Strategies

As mentioned before, the three types of flux terms (pressure, mass flux, and advected quantities) from Section 2.6 are each being tested with the four different notions of upwinding: central, advectively upwind, acoustic-only upwind, and fully upwind. Acoustic-only upwind fluxes only apply to mass flux and pressure, not to advected quantities. Also, fully upwind fluxes gives extremely similar results to using acoustic-only fluxes for pressure and mass flux and using advectively upwind fluxes for advected quantities (called an “acoust-acoust-advect” configuration from here out to denote the treatment for pressure, mass flux, and advected quantities, in that order). As evidence for this, please see Figure 2 the KE spectra below comparing the two approaches at ninth-order accuracy with no limiting. This similarity holds for third-order accuracy and for WENO limiting as well. In this similarity, the “acoust-acoust-advect” approach is slightly more dissipative.

Because of how similar the full Riemann solver is to using acoustic-only upwinding for pressure and mass flux terms and advective-only upwinding for advected terms, only three upwinding strategies are tested in this experiment for mass flux and pressure (central, advectively upwind, and acoustically upwind); and only two strategies are tested for advected quantities (central and advectively upwind). This means there are 3 × 3 × 2 = 18 Riemann solver combinations to test. A given combination will be described by the upwind treatment in the order of pressure, mass flux, and advected quantity—meaning that “acoust-advect-central” would mean acoustic-only upwinding is used for the pressure terms, advective-only upwinding is used for the mass flux terms, and central values are used for the advected quantities.

Identifying Unstable and Highly Unphysical Configurations

The model's stability using various flux upwinding options is tested for four different reconstruction strategies in Table 1: third-order and ninth-order accuracies both with and without WENO limiting for scalar quantities. Supercell simulations were run for one model day using 50 × 50 × 40 cells. “Unstable” means infinities appeared in the solution, “unphysical” means persistent, visible, and global 2Δx oscillations appeared in the solution for at least one direction and variable. “Stable” means the simulation completed, and oscillations (if any) are either non-persistent and local to a discontinuity or hard enough to discern visibly that further investigation is worthwhile. Not every simulation labeled as “stable” will result in a physically realizable solution.

Table 1 Table of the Stability of Various Upwind Strategies for Various Flux Terms for Different Spatial Reconstruction Strategies by Running a 3-D Idealized Supercell Simulation on a Grid of 50 × 50 × 40 Cells for One Model Day (86,400 s)

Pressure	Mass flux	Advected	nolim9	weno9	nolim3	weno3
acoust	acoust	advect	Stable	Stable	Stable	Stable
acoust	acoust	central	Unphysical	Unstable	Stable	Stable
acoust	advect	advect	Unphysical	Unstable	Unphysical	Unphysical
acoust	advect	central	Unstable	Unstable	Unphysical	Unphysical
acoust	central	advect	Stable	Stable	Stable	Stable
acoust	central	central	Unphysical	Unstable	Unphysical	Unphysical
advect	acoust	advect	Unphysical	Unstable	Unphysical	Unphysical
advect	acoust	central	Unstable	Unstable	Unstable	Unstable
advect	advect	advect	Unstable	Unstable	Unstable	Unstable
advect	advect	central	Unstable	Unstable	Unstable	Unstable
advect	central	advect	Unstable	Unstable	Unstable	Unstable
advect	central	central	Unstable	Unstable	Unstable	Unstable
central	acoust	advect	Unphysical	Unphysical	Stable	Unphysical
central	acoust	central	Unstable	Unstable	Unstable	Unstable
central	advect	advect	Unstable	Unstable	Unstable	Unstable
central	advect	central	Unstable	Unstable	Unstable	Unstable
central	central	advect	Unphysical	Unstable	Unphysical	Unphysical
central	central	central	Unstable	Unstable	Unstable	Unstable

Note. “Unstable” means the solution reached infinity. “Unphysical” means the solution contains strong, visible, and global 2 Δx oscillations in at least one direction. “Stable” means that oscillations in the solution (if any) are local to a discontinuity and not persistent. “Advected” refers to advected quantities in the flux terms. “acoust,” “advect,” and “central” indicate acoustic-only, advective-only, and central strategies for determining upwind interface fluxes.

From Table 1, it is clear that the stability of different upwinding options generally improves when using a more dissipative third-order reconstruction, which matches intuition. When going from non-limited interpolation to WENO interpolation, the stability actually decreases, which may go against intuition for some. Keep in mind, however, that a significant weakness of the WENO approach is that in the presence of non-centered discontinuities, it relaxes to polynomials that are biased toward a spatial direction and are no longer centered on the cell being reconstructed. This increases variation in the Vandermonde basis functions over that cell domain. In the one-sided case, where the cell being reconstructed is at the very edge of the stencil of constraint, this can even lead to instability, since the sampling location is not well-bounded by the domain of constraint in a FV context (see Figure 3 of M. R. Norman and Nair (2018)).

From Table 1, the largest number of stable configurations occur when the pressure term is acoustically upwind, indicating that the pressure term is likely the most important for stability. Notably, there are no stable configurations that use advectively upwind values for either pressure or mass flux terms. Arguably, this is because whenever pressure or mass flux use an upwind value with regard to advection, for subsonic flow they inevitably are using a downwind value for one of the acoustic modes. Since these terms are acoustically active, this likely causes a destabilizing effect. There are stable configurations for all reconstruction strategies when using a central approach to the mass fluxes along with appropriate upwind definitions for the pressure terms and advected quantities.

Further Determining Physically Realizable Configurations

Each stable simulation from the previous experiment is now simulated using 250 × 250 × 50 cells for one model day. Figure 3 shows KE spectra for all reconstruction strategies using an “acoust-acoust-advect” upwind strategy. Since WENO limiting is only used for scalars and not for winds, it is unsurprising that the damping effects on KE spectra are relatively small. It is also clear from the plot that using ninth-order accuracy maintains KE in the smaller wavelengths better than using third-order accuracy, which is an expected result. Each of these configurations show no accumulation of energy at the smallest wavelengths, which confirms that for all reconstruction strategies, this upwinding strategy gives physically realizable results. Visual inspection of each of these reveals that the solutions are also free of significant gridscale noise.

In Figures 4a and 4b, the KE spectra for the “acoust-central-advect” upwinding strategy is plotted at ninth and third orders of accuracy, respectively, using the “acoust-acoust-advect” upwinding strategy previously deemed as stable as a baseline. When WENO limiting is used, clearly there is a very rapid and unphysical increase in energy at the very smallest 2Δx wavelengths. When WENO limiting is not used, however, the increase in energy at the 2Δx scale is much smaller and less concerning at ninth-order accuracy and is hardly noticeable at all at third-order accuracy. Visual inspection shows that both the non-limited third-order and ninth-order solutions with “acoust-central-advect” upwinding are also free of gridpoint noise and can therefore be classified as physically realizable. The WENO-limited “acoust-central-advect” third and ninth-order simulations do not show visible gridscale noise, but the KE spectra are unphysical enough that they fail one of the criteria and should be considered physically unrealizable.

In Figure 4c, the KE spectra of the “acoust-acoust-central” upwinding strategy is plotted at third-order accuracy with and without WENO limiting; and the “central-acoust-advect” upwinding strategy is plotted at third-order accuracy without limiting—together, encapsulating the remaining configurations deemed stable enough to warrant further investigation in Section 3.2.1. A baseline of ninth-order accuracy with an “acoust-acoust-advect” upwinding strategy is plotted as well. Both third-order “acoust-acoust-central” simulations not only demonstrate a highly unphysical increase in energy at smaller scales (along with a decrease in energy at larger scales), but they also reveal significant grid point noise that is clearly not physically realizable. It isn't immediately clear what causes the decrease in energy at larger scales for the acoust-acoust-central scheme, but regardless, it is not physically realizable. The third-order unlimited “central-acoust-advect” simulation, however, has a much more reasonable looking KE spectra, and visual inspection reveals the solution is free of significant gridscale noise, meaning that configuration should be deemed physically realizable.

Summary of Physically Realizable Configurations

In summary, the following configurations are found to be physically realizable from the previous experiments:

• “Acoust-acoust-advect” upwinding is realizable for all reconstruction strategies

• “Acoust-central-advect” upwinding is realizable only for non-limited reconstruction strategies

• “Central-acoust-advect” upwinding is realizable only for third-order accuracy with no limiting

Therefore, if the developer restricts themselves to low-order accuracy, then the pressure terms can be calculated in a central fashion. If the developer restricts themselves to not use WENO limiting, then the mass flux terms can be calculated in a central fashion. Otherwise, the developer should prioritize acoustic upwinding for pressure and mass flux terms and advective upwinding for advected quantities. Note that using the full upwind definition from the diagonalization of the full Jacobian can be more or less substituted with acoustic-only upwinding for pressure and mass flux and advective-only upwinding for advected quantities with the minor caveat that the latter is slightly more diffusive (albeit simpler and cheaper to implement).

Assuming a Constant Speed of Sound

A potentially helpful simplification to determining the acoustically upwind states of pressure and mass flux terms would be to assume that the speed of sound is constant. The upwind values, expanded from the diagonalization of the equation set (18), for the acoustic fluxes are given below with an assumption of a “locally frozen” true acoustic speed, c_s: [Image Omitted. See PDF] [Image Omitted. See PDF]where $${c}_{s}^{\star }$$ is an assumed constant speed of sound in the acoustic-only flux Jacobian in Equation 18 and its subsequent diagonalization. If this were confined to the acoustic-only equation system, the momentum would be multiplied by the true squared speed of sound, $${c}_{s}^{2}$$, to constitute the flux of the pressure time tendency. In that form, it gives the LLF flux, and in that form, a larger magnitude of c_s would lead to larger dissipation.

However, in the form above, it is not entirely clear what a model developer should expect to encounter from different sound speed assumptions in the Riemann solver. The pressure term appears to exhibit greater dissipation with increasing assumed speed of sound, but the momentum term appears to exhibit less dissipation with increasing assumed speed of sound. Therefore, numerical experiments are performed to gather information about the effects of smaller and larger magnitudes for the assumed speed of sound in the acoustic-only flux Jacobian.

In Figure 5, the ninth-order accurate solution with no limiter is computed using the “acoust-acoust-advect” upwinding strategy with a variety of constant sound speeds along with the actual sound speed as a reference. Using a constant speed of sound (in meters per second) of 150 was unstable, but stable simulations were found in at least the domain of $$\left[200,500\right]$$ meters per second. The primary finding is that assuming a constant speed of sound has relatively little effect on the KE spectra overall. Using 200 m per second decreased dissipation at small scales and increased dissipation at large scales but only by relatively small amounts. Using 300–400 m per second generally gave results close to using the actual speed of sound. Using 500 m per second gave the most KE power at larger scales while hovering in the mean of the constant sound speed experiments in the small scales.

In summary, the user can assume a constant speed of sound within at least the bounds of $$\left[200,500\right]$$ meters per second in the acoustic-only upwind calculations for pressure and mass flux while retaining numerical stability and with relatively small effects on the resulting solution. Choosing a realistic value between 300 and 400 m per second gives a result most similar to using the actual sound speed.

Using Convex Combinations of Upwind and Central Fluxes

Using acoustic-only upwinding for pressure and mass flux terms and advective-only upwinding for advected quantities is a stable and physically realizable approach to the compressible, non-hydrostatic Euler equations. It stands to reason that using a convex combination of central and upwind fluxes could be stable and physically realizable with a reduced level of dissipation. The ultimate effect of this approach will be to reduce the magnitude of dissipation created by the preference for upwind information. This idea is somewhat related to the idea of using anti-dissipation to restore the contact wave of a Riemann solver for shock tube simulations such as for the “HLLC” Riemann solver (Toro et al., 1994). The difference is that this technique is done for specific flux terms rather than certain characteristic variables.

The parameter, $${\alpha }_{up}\in \left[0,1\right]$$, will be used to denote the proportion of upwind information used in this convex sum for pressure and mass flux terms such that 1 − α_up is the proportion of central information used. Therefore, α_up = 1 means only upwind information is used, and α_up = 0 means only central estimates are used. The parameter, $${\beta }_{up}\in \left[0,1\right]$$, will be used to denote the proportion of upwind information used for advected quantities and has the same behavior as α_up.

The effects of convex sums of upwind and central interface fluxes are determined through numerical experiments using 250 × 250 × 50 cells with ninth-order accuracy and WENO limiting used for scalars but not momenta. Using α_up = 0 (with any value for β_up) and β_up = 0 (with any value for α_up) both provide unstable solutions reaching infinity values. When α_up = 0, this is equivalent to a “central” approach for pressure and mass flux. When β_up = 0, this is equivalent to a “central” approach for advected quantities. Numerical experiments were performed with a tensored grid of $${\alpha }_{up},{\beta }_{up}\in \left\{0.2,0.4,0.6,0.8,1.0\right\}$$, giving a total of 25 numerical experiments, all of which completed without infinity values or visible persistent global oscillations. Note that α_up = β_up = 1 gives the same solution as the “acoust-acoust-advect” upwind strategy at ninth-order with WENO limiting for scalar quantities.

Figure 6a shows the effects of convolving in central fluxes for advected quantities, and Figure 6b shows the effects of convolving in central fluxes for pressure and mass flux terms (these are treated together in this section). Figure 6c shows the effects of convolving in central fluxes for all terms equally. In each of these figures, the expected effect occurs that KE is preserved at larger magnitudes at smaller scales with increasing portions of the fluxes being central estimates. What was not necessarily expected, however, is that in all experiments that avoid infinity values, the KE spectra asymptote toward extending the k^−5/3 model commonly observed for 3-D turbulence. This should not necessarily be extended to declaring that the increased resolution at small scales is physical, however. The KE spectra, at least, appear to follow physical intuition, but that does not guarantee that it follows the k^−5/3 model for physical reasons. In fact, further inspection will reveal that for at least some experiments, it is not for physically realizable reasons.

On close inspection, many of the spectra exhibit a slight uptick in energy at the very smallest ≈2Δx wavelengths admitted by the model's grid, which is not a desirable attribute, though the magnitude is small. Previous results have shown that small upticks at 2Δx can still provide what appear to be physically realizable results. Another interesting feature is that incorporating central fluxes for advected quantities has a significantly greater effect on the KE spectra than incorporating central fluxes for mass flux and pressure terms.

Figure 7a shows the velocity at 5 km with α_up = β_up = 0.2 with 250 × 250 × 50 cells at 2 hr model time. Figure 7c shows the same with fully upwind fluxes with 1,000 × 1,000 × 200 cells at 2 hr model time as a high-resolution comparison. Clearly, with mostly central fluxes, even though the KE spectra follow what appears to be a physically realistic model, there are strong, spurious vertical velocities far away from where the supercell is located, indicating a lack of physical realism. This is an example of a physically unrealizable amount of grid scale noise. A plot of α_up = β_up = 0.8 in Figure 7b shows a much more realistic looking flow without those spurious vertical velocities. This means there is a degree to which one can include central fluxes for anti-dissipative effects and increased resolution at smaller scales. However, there is a limit to how much this should be done in practice.

Investigating the Limits of Spatial Order of Accuracy

A last set of experiments seeks to test the hypothesis that there is a limit to how high-order-accurate a simulation can be before upwind fluxes alone lead to a physically realizable and stable solution. Therefore, non-limited simulations at 11th, 15th, and 21st orders of accuracy are performed with the “acoust-acoust-advect” upwinding strategy on a grid of 250 × 250 × 50 cells.

The KE spectra of these simulations are plotted in Figure 8. The KE spectra show an alarming increase in energy at small wavelengths, showing that higher-order accuracy has a clear impact on the apparent resolution of the simulation, even though the data being simulated is known to have a continuity only between $${\mathcal{C}}^{1}-{\mathcal{C}}^{2}$$. This behavior is, at first glance, encouraging because in an ideal world, the k^−5/3 cascade has no physical reason for dropping off at small wavelengths to a cascade that is much steeper (as seen in most if not all stable numerical simulations). However, seeking a fully extended energy cascade (i.e., delaying the scale at which numerical dissipation takes over) cannot be a goal in itself. One must ensure that this cascade is extended due to physically realizable reasons and is free of gridscale noise.

Vertical velocity at 5 km is plotted in Figures 9a–9d for 9th, 11th, 15th, and 21st orders of accuracy, respectively. Again, the idea of having gridscale noise that is too large and frequent to be physically realizable is a somewhat subjective notion. However, the solution at 21st-order accuracy has clearly reached that point where it is not physically realizable. The solutions at 11th and 15th orders of accuracy are less clear in this regard. So this confirms Hypothesis H5, though it is not clear at exactly what order of accuracy this occurs. As was the case with taking convex sums of upwind and central flux estimates, it is certainly interesting that even with physically unrealizable results, the KE energy cascade continues to extend along the k^−5/3 model quite closely.

Discussion and Implications Hypothesis H1

Hypothesis H1 was that there are simplifications to the fully upwind Riemann solver of (M. R. Norman, 2021) (and fully upwind solvers in general) that result in physically realizable results; and this hypothesis was confirmed. While others were found, the most reliable of the simplifications is the strategy of discretizing mass flux and pressure using acoustic-only upwind estimates according to the simplified acoustic system in equation set (18) and advected quantities using advective-only upwind estimates.

There are other configurations that only work for non-limited reconstructions or only work at lower orders of accuracy (e.g., Third-order). However, from a practical point of view, they are not as valuable because only pressure or mass flux can be computed using a central estimate. With that constraint, still at least one of those types of terms must be computed upwind according to acoustic dynamics. Given that cell edge estimates of both pressure and mass flux are required in order to get an acoustically upwind flux of either, then one might as well discretize both in a manner that is acoustically upwind. Computationally, very little is saved by only discretizing one of them in an acoustically upwind manner.

From an implementation perspective, there is quite a bit to gain from this simplification to the fully upwind Riemann solver. As mentioned earlier, one can get an advectively upwind interpolation of a variable while it is being reconstructed merely by biasing the stencil of reconstruction to be centered on the upwind cell adjacent to a cell edge. This has three computational benefits for advected quantities:

1. Only one dot product with respect to the stencil is needed instead of two, which reduces the required computations;

2. Only one kernel of computation is needed instead of two since one no longer needs to reconcile two values at the same cell edge location; and

3. Because no reconciliation is needed at cell edges, advected quantities do not need to be communicated over MPI with other domain in preparation for calculating upwind fluxes.

Since all but two of the evolved quantities are classified as “advected quantities”, this means that there is a significant reduction in data volume for the cell-edge MPI data exchange in preparation for computing upwind fluxes (particularly when a large number of advected tracers are used). This also reduces the amount of data accessed in the upwind calculation kernel: only two variables, pressure and mass flux, are accessed in that kernel now. On Graphics Processing Units, this reduces the cost of the upwind flux calculation kernel because data accesses are the most expensive operations (relative to computations). In parallel decompositions, this reduces the amount of data that needs to be transferred over MPI.

In fact, one could choose to redundantly reconstruct pressure and mass flux cell edge estimates for one extra cell per horizontal direction at domain decomposition edges (at the expense of communicating one extra cell per direction for just those two variables in the MPI halo exchange). If this is done, then the upwind calculation can be done in a single kernel and without an extra stage of MPI communication. Certainly for simulations with a large number of tracers, this would provide a clear benefit in most computational contexts.

One further simplification to a fully upwind Riemann solver and to the above-mentioned simplified strategy is that one could assume a constant speed of sound in the acoustic-only upwind calculation for pressure and mass flux. This, in effect, mostly just avoids a square root calculation, but for some, it may be a worthwhile effort. Choosing a realistic constant speed of sound in the upwind calculation of 300–400 m per second produced results quite close to using the actual speed of sound. For time-implicit discretizations that set up a linear system, this strategy can also further simplify the solver as well.

Hypotheses H2 and H3

Hypothesis H2 was that upwind dissipation has different effects for different types of terms in the flux vector; and this hypothesis was confirmed. Hypothesis H3 was that there are noticeable differences between preferring upwind information according to advective velocity and according to acoustic dynamics. H3 was also confirmed in experiments.

There were no physically realizable configurations that used central estimates of advected quantities, which seems to indicate that upwinding for advected quantities is the most important for the moist Euler equations with microphysics. This likely is influenced by the inclusion of microphysics, which create strong interactions with the quality of advection. There were also no physically realizable simulations that compute pressure or mass flux terms in an advectively upwind manner. This is likely because an advectively upwind direction for subsonic flow is, by necessity, downwind for one of the acoustic modes. Finally, there were some physically realizable configurations using central differencing for pressure or mass flux (but not for both at the same time). These, however, came with restrictions on the reconstruction operator—that it must not include WENO limiting or that the order of accuracy is limited.

Hypothesis H4

Hypothesis H4 was that there is an upper limit to the spatial order of accuracy that can produce physically realizable results with only upwind preference for dissipation; and this was confirmed as well—though the order of accuracy at which the solution becomes unrealizable is not entirely clear. Certainly at 21st-order accuracy, the solution is not realizable as evidenced by strong gridscale noise that is not related to supercell convection and does not appear in a high-resolution baseline simulation. The gridscale noise gradually increases with the order of accuracy, but defining exactly which order the noise is unacceptable is a subjective exercise. This is certainly not to say that there aren't other Riemann solvers or dissipation/stabilization choices that could stabilize a solution at extremely high orders of accuracy. But in the context of these numerics and Riemann solver choices, the Riemann solver alone is not enough to create physically realizable solutions at 21st-order accuracy.

Hypothesis H5

Hypothesis H5 was that physically realizable discretizations can be obtained by using convex combinations of upwind and central fluxes as a means of introducing anti-dissipation to varying extents. This hypothesis was confirmed—though as with Hypothesis H4, the point at which solutions become unrealizable is somewhat subjective because the dominant artifact in the solution is a gradual increase of gridscale noise as the proportion of central fluxes increases.

This hypothesis being confirmed has some important implications along with H4 as well in that it provides a fairly easily tunable set of parameters to control the dissipation of the numerical solution. And wherever there are tunable parameters, one can conceive of a potential Machine Learning approach that can automatically adjust these parameters to optimize a chosen cost function. Having two (and possibly more) floating point numbers in the domain of $$\left[0,1\right]$$ that control dissipation monotonically is a powerful opportunity to use ML in the development of numerical algorithms for fluid dynamics to achieve desired properties.

Using ML-based optimization of free parameters in fluid simulations is not a new idea (Bar-Sinai et al., 2019). Also, optimization of free parameters for non-oscillatory reconstruction techniques (Kossaczká et al., 2021) likely has as much or more promise as the proportion of central fluxes to convolve into the final flux. These can also be used in tandem to control dissipation in both the reconstruction and upwinding phases of a numerical algorithm.

A Final Note on Spectra

If one considers the Euler equation (or the Navier-Stokes equations far away from dissipation-dominated scales), physical intuition says that in the absence of other dynamics, there is no reason to expect a deviation from a KE cascade that follows 3-D turbulence when the scales are small enough. This can lead a numerical developer to seek out a KE cascade that maintains a k^−5/3 slope all the way to 2Δx scales while avoiding the much steeper numerical dissipation regime experienced by most if not all stable numerical simulations. Coarsening a very high-resolution simulation can show this effect. However, this study has shown clearly that one can indeed extend the KE spectra along k^−5/3 toward the smallest scales and do so with a very physically unrealizable simulation. Therefore, care should always be taken to ensure that the solution, though following a physical looking KE energy cascade, is behaving in a way that matches basic physical intuition such as the constraints used in this study.

Conclusions and Future Work

The goal of this paper is to gain a better understanding of inherent numerical dissipation in collocated discretizations of the moist, non-hydrostatic, compressible Euler equations with reacting microphysics—particularly from using upwind information for cell-edge fluxes and from using varying levels of dissipation in spatial reconstruction operators. This is explored through a set of five main hypotheses that are tested using numerical experiments, all of which were confirmed.

Upwind fluxes have been used in many numerical implementations of the Euler and Navier-Stokes equations as a dissipation mechanism that inherently adapts to changing grid spacing, time step size, reconstruction accuracy, and flow smoothness. During these numerical experiments, a simplification to a fully upwind solution was detailed that results in tangible reductions in computations and data transfers when compared to a fully upwind approach. The bounds of spatial order of accuracy in the absence of explicit dissipation were also explored, and it was found that though extremely high-order accuracy (21st-order in this case) does not produce infinity values or global, persistent oscillations, it does produce unacceptable levels of gridscale noise inconsistent with physically realizable dynamics.

For future work, it would be interesting to understand these ideas in more realistic contexts. These include extending to spherical and non-orthogonal geometries, the inclusion of highly variable terrain, the use of alternative time stepping operators such as non SSP semi-discrete operators and ADER time stepping procedures. Another area of future work would be to better understand how a model might produce a physical looking energy cascade with physically unrealizable solutions and how to detect when this happens in situ as a warning to the developer that the numerics likely need to change.

The experiments in this study have also given rise to continuous parameters that monotonically affect dissipation and could be tuned by Machine Learning training approaches to optimize developer-determined cost functions.

Acknowledgments

This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This research was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration. Notice: This manuscript has been authored by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan). This article has been co-authored by an employee of National Technology and Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns right, title and interest in and to the article and is responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-plan.

Data Availability Statement

The code used in this study, which is sufficient for reproduction of all results, can be obtained here: https://zenodo.org/record/7781644. The citation is M. Norman (2023).