On the Entropy Projection and the Robustness of High Order Entropy Stable Discontinuous Galerkin Schemes for Under-Resolved Flows

1 Introduction

Discontinuous Galerkin (DG) schemes have received interest within computational fluid dynamics (CFD) due to their high order accuracy and ability to handle unstructured curved meshes. In particular, there has been interest in DG methods for simulations of under-resolved flows [1–5]. Among such schemes, “entropy stable” DG methods based on a “flux differencing” formulation have received interest due to their robustness with respect to shocks and turbulence [6–9].

Entropy conservative and entropy stable flux differencing schemes were originally formulated for finite difference methods in [10, 11]. They were extended to tensor product grids using discontinuous spectral collocation schemes (also known as discontinuous Galerkin spectral element methods, or DGSEM) [12, 13]. Entropy stable collocation schemes were extended to simplicial meshes in [14, 15] using multi-dimensional summation-by-parts (SBP) operators [16]. Non-collocation entropy stable schemes have also been developed. These schemes began with staggered grid schemes on tensor product grids in [17], which were later extended to simplicial elements in [18]. “Modal” entropy stable DG formulations [19–21] have been utilized to construct a variety of new entropy stable schemes, including Gauss DG methods [22, 23] and reduced order models [24]. We note that under appropriate choices of quadrature, these “modal” formulations reduce to collocation-type entropy stable schemes. Entropy stable schemes have since been extended to an even wider array of discretizations, such as line DG methods, discontinuous Galerkin difference methods, and C⁰ continuous discretizations [25–27].

The main difference between non-collocation and collocation-type entropy stable schemes is the use of transformations between conservative variables and entropy variables together with projection or prolongation operators to facilitate a discrete proof of entropy stability. This is referred to as the “entropy projection” in [19, 25] and as the interpolation or prolongation of entropy variables in [17, 27]. This approach is also equivalent to the mixed formulation of [28]. We will refer to this transformation as the “entropy projection” for the remainder of the paper.

The motivation for introducing the entropy projection has been to enable the use of more accurate quadrature rules or novel basis functions. This has been at the cost of additional complexity and issues related to the sensitivity of the entropy variables for near-vacuum states [19, 27]. To the best of the authors’ knowledge, no inherent advantages in using the entropy projection have been observed in the literature. This paper focuses on the following observation: high order entropy stable schemes based on the entropy projection appear to be more robust than entropy stable collocation schemes for two and three dimensional simulations of under-resolved variable-density fluid flows with small-scale features.

The structure of the paper is as follows: Section 2 reviews mathematical formulations of entropy stable schemes which involve the entropy projection. Section 3 documents the observed difference in robustness for a variety of problems in two and three dimensions, and provides analysis and numerical experiments which support that the primary difference between unstable and stable schemes is the entropy projection. Section 4 conjectures potential explanations for why the entropy projection might improve robustness. We conclude with Section 5, which explores potential applications towards under-resolved flow simulations.

2 Formulation of High Order Entropy Stable Discontinuous Galerkin Schemes

In this section, we provide a brief description of high order entropy stable schemes in 1D. More detailed derivations, multi-dimensional formulations, and extensions to curved grids can be found in [14, 15, 19, 21, 22, 24].

The notation in this paper is motivated by notation in [15, 29]. Unless otherwise specified, vector and matrix quantities are denoted using lower and upper case bold font, respectively. Spatially discrete quantities are denoted using a bold sans serif font. Finally, the output of continuous functions evaluated over discrete vectors is interpreted as a discrete vector.

For example, if $\mathsf{x}$ denotes a vector of point locations, i.e., ${(\mathsf{x})}_i={\mathit{x}}_i$, then $u(\mathsf{x})$ is interpreted as the vector

Similarly, if $\mathsf{u}=u(\mathsf{x})$, then $f(\mathsf{u})$ corresponds to the vector

Vector-valued functions are treated similarly. For example, given a vector-valued function $\mathit{f}:{\mathbb{R}}^n\to {\mathbb{R}}^n$ and a vector of coordinates $\mathsf{x}$, we adopt the convention that ${\left(\mathit{f}(\mathsf{x})\right)}_i=\mathit{f}({\mathit{x}}_i)$.

2.1 Conservation Laws With Entropy

In this section, we review the construction of entropy conservative and entropy stable schemes for a one-dimensional system of nonlinear conservation laws

where s(u) is a source term. We assume the domain is exactly represented by a uniform mesh consisting of non-overlapping intervals D^k, and that the solution u(x) is approximated by degree N polynomials over each element. We also introduce entropy conservative numerical fluxes f_S(u_L, u_R) [30], which are bivariate functions of “left” and “right” states u_L, u_R. In addition to being symmetric and consistent, entropy conservative numerical fluxes satisfy an “entropy conservation” property

here, v_L, v_R are entropy variables evaluated at the left and right states, and ψ(u) denotes the “entropy potential”. Examples of expressions for entropy variables and entropy potentials can be found in [14].

2.2 Collocation Formulations

Degree N entropy stable collocation schemes are typically built from Legendre-Gauss-Lobatto (LGL) quadrature rules with (N + 1) points. Let $\mathsf{x},\mathsf{w}$ denote vectors of quadrature points and weights on the reference interval [−1, 1]. Let ℓ_i(x) denote Lagrange polynomials at LGL nodes, and let $\mathsf{u}$ denote the vector of solution nodal values u(x_i). Define the matrices

here ${\mathsf{V}}_f$ is a face interpolation or extraction matrix which maps from volume nodes to face nodes. Flux derivatives are discretized using a “flux differencing” approach involving summation-by-parts (SBP) operators and entropy conservative fluxes [30]. An entropy stable collocation formulation can now be defined on an element D^k as follows:

where h is the size of the element D^k and ◦ denotes the matrix Hadamard product [10–12].¹ Here, ${\mathsf{f}}^{\ast }$ is a vector which contains numerical fluxes at the left and right endpoints of the interval

where ${\mathsf{u}}_1^{+},{\mathsf{u}}_{N+1}^{+}$ denote exterior nodal values on neighboring elements. If ${\mathit{f}}^{\ast }$ is an entropy conservative flux, then the resulting numerical method is semi-discretely entropy conservative. If ${\mathit{f}}^{\ast }$ is an entropy stable flux (for example, Lax-Friedrichs flux, HLLC, and certain matrix penalizations [14, 31]) then the resulting scheme also dissipates entropy.

2.3 “Modal” Formulations

Degree N entropy stable “modal” DG schemes generalize collocation schemes to arbitrary choices of quadrature. In one dimension, this allow for the use of higher accuracy volume quadratures. In higher dimensions, modal formulations also enable more general choices of surface quadrature. These schemes introduce an additional “entropy projection” step to facilitate the semi-discrete proof of entropy stability or conservation.

We now assume the solution is represented using some arbitrary basis over each element, such that $u(x)={\sum }_j{\mathsf{u}}_j{\varphi }_j(x)$. Let x, w now denote a general quadrature rule with positive quadrature weights. We define quadrature-based interpolation matrices ${\mathsf{V}}_q,{\mathsf{V}}_f$, the mass matrix $\mathsf{M}$, and the modal differentiation matrix $\widehat{\mathsf{Q}}$

We introduce the quadrature-based projection matrix ${\mathsf{P}}_q={\mathsf{M}}^{-1}{\mathsf{V}}_q^T\mathrm{diag}\left(\mathsf{w}\right)$. Using ${\mathsf{P}}_q$ and $\widehat{\mathsf{Q}}$, we can construct quadrature-based differentiation and extrapolation matrices $\mathsf{Q}$, $\mathsf{E}$

To accommodate general quadrature rules which may not include boundary points, we introduce hybridized SBP operators ${\mathsf{Q}}_h$ on the reference interval [−1, 1]

The use of such operators simplifies the implementation for general quadrature rules and nodal sets which do not include boundary nodes [19, 32]. Next, we define ${\mathsf{V}}_h$ as the interpolation matrix to both volume and surface quadrature points

We also introduce the L² projection of the entropy variables and the “entropy projected” conservative variables $\tilde{\mathsf{u}}$

which are defined by evaluating the mapping from conservative to entropy variables $\mathit{u}\left(\mathit{v}\right)$ using the projected entropy variables. Here, $\mathit{v}\left(\mathit{u}\right)$ denotes the mapping from conservative to entropy variables. Note that the projected entropy variables $\mathsf{v}$ is a vector corresponding to modal coefficients, while $\tilde{\mathsf{u}}$ corresponds to point values at volume and face quadrature points.

An entropy stable modal DG discretization over a single element D^k is then

Note that the right hand side formulation is evaluated not using the conservative variables $\mathsf{u}$, but the “entropy projected” conservative variables $\tilde{\mathsf{u}}$.

While we have presented entropy stable DG schemes using a general “modal” DG framework, the formulation reduces to existing methods under appropriate choices of quadrature and basis. For example, specifying LGL quadrature on a tensor product element recovers entropy stable spectral collocation schemes [22]. SBP discretizations without an underlying basis on simplices [14–16] can also be recovered for appropriate quadrature rules by redefining the interpolation and projection matrices ${\mathsf{V}}_q,{\mathsf{P}}_q$ [33].

3 Numerical Comparisons of Collocation and Entropy Projection Schemes

In this section, we will demonstrate numerically that a significant difference in robustness is observed between collocation and entropy projection-based discretizations of the Euler and ideal MHD equations. For the Euler equations, we study the Kelvin-Helmholtz, Rayleigh-Taylor, and Richtmeyer-Meshkov instabilities, and for the MHD equations we study a magnetized Kelvin-Helmholtz instability. All of these examples exhibit small-scale turbulent-like features. Moreover, we observe a difference in robustness between entropy stable collocation and entropy projection-based methods independently of the polynomial degrees, mesh resolutions, and type of mesh (e.g., quadrilateral or triangular). We focus on the following entropy stable DG methods:

• On quadrilateral meshes:

(1) DGSEM: collocation scheme based on the tensor product of one-dimensional (N + 1) point LGL quadrature,

(2) Gauss DG: a “collocation” scheme based on the tensor product of one-dimensional (N + 1) point Gauss quadrature. The entropy projection is used to evaluate interface fluxes [22],

• On triangular meshes:

(1) SBP: a collocation scheme based on multi-dimensional summation-by-parts finite difference operators [14, 16],

(2) Modal: a modal formulation utilizing quadrature rules which exactly integrate entries of the volume and face mass matrices [19].

Remark 1. It is known that the Kelvin-Helmholtz, Rayleigh-Taylor, and Richtmeyer-Meshkov instabilities are notoriously sensitive to initial conditions and discretization parameters, and that numerical schemes may not converge to a unique solution [34, 35]. Instead, this paper focuses on these problems as stress tests of robustness.Unless specified otherwise, all numerical experiments utilize a Lax-Friedrichs interface flux with Davis wavespeed estimate [36]. We also experimented with HLL and HLLC surface fluxes, but did not notice a significant difference. We also note that instead of discontinuous initial conditions, we utilize smoothed approximations for each problem considered here.All experiments are also performed on uniform meshes. For triangular meshes, this mesh is constructed by bisecting each element of a uniform quadrilateral mesh along the diagonal. Unless specified otherwise, all results are produced using the Julia [37] simulation framework Trixi.jl [38, 39]. For most experiments, we utilize an optimized adaptive 4th order 9-stage Runge-Kutta method [40] implemented in OrdinaryDiffEq.jl [41]. The absolute and relative tolerances are set to 10^–7 unless specified otherwise. Scripts generating main results are included in a companion repository for reproducibility [42].We note that the robustness, efficiency, and high order accuracy of both entropy stable DGSEM and entropy stable Gauss DG schemes have been verified in previous works [7–9, 22, 23], and will not be addressed in detail in this paper. However, the difference in robustness between the two methods has not been previously observed in the literature, and will be the focus of this work.

3.1 Euler Equations of Gas Dynamics

We consider first the two and three-dimensional problems for the Euler equations of gas dynamics. The conservative variables for the three-dimensional Euler equations are density, momentum, and total energy, $\mathit{u}=(\rho ,\rho \mathbf{v},E)$, where the vector $\mathbf{v}=(u,v,w)$ contains the velocities in x, y and z, respectively. The flux reads

where $\mathbf{I}$ is the 3 × 3 identity matrix, γ is the heat capacity ratio, and $p=(\gamma -1)\left(E-\rho \Vert \mathbf{v}{\Vert }^2/2\right)$ is the gas pressure. For two-dimensional problems, we neglect the third component of the velocity, w, and $\mathbf{I}$ becomes the 2 × 2 identity matrix.

All the following experiments use the entropy conservative and kinetic energy preserving flux of Ranocha [43, 44]; however, similar results were observed when experimenting with the entropy conservative flux of Chandrashekar [45].

3.1.1 Two Dimensional Kelvin-Helmholtz Instability

We perform additional experiments analyzing the robustness of entropy stable DGSEM and Gauss DG for the Kelvin-Helmholtz instability. The domain is [−1, 1]² with initial condition from [46]:

where B(x, y) is a smoothed approximation to a discontinuous step function

Each solver is run until final time T_final = 15. As can be observed in Figure 1, the solution differs significantly between the N = 3 and N = 7 simulations. This is likely a consequence of the well-known sensitivity of the Kelvin-Helmholtz instability to small perturbations and numerical resolutions [34, 35]. End times for each simulation can be found in Table 1.

FIGURE 1. Snapshots of density for the Kelvin-Helmholtz instability using an entropy stable Gauss DG scheme on uniform quadrilateral meshes.

TABLE 1. End time for simulations of the 3D Kelvin-Helmholtz instability on hexahedral meshes. “Collocation” refers to a nodal DGSEM discretization, while “entropy projection” refers to a method based on Gauss nodes.

3.1.2 Two Dimensional Rayleigh-Taylor Instability

The two-dimensional Rayleigh-Taylor instability generates small-scale flow features through buoyancy or gravity effects [47, 48]. The setup involves a heavy and light fluid suspended above one another separated by a curved interface, and buoyancy or gravity results in displacement of the lighter fluid into the heavier one. This displacement causes velocity shear and the formation of additional Kelvin-Helmholtz instabilities along the interface. The domain is [0, 1/4] × [0, 1].

Let $d_{a,b}(x)=a+\frac{1}{2}(1+\tanh (sx))(b-a)$ denote a smoothed approximation (with slope s) to a discontinuous function with values a for x < 0 and b for x > 0. The initial condition is given by

where $c=\sqrt{\gamma p/\rho }$ is the speed of sound. Here, we borrow from [49] and multiply the y-velocity perturbation by sin(πy)⁶ so that u, v satisfy wall boundary conditions. We also add gravity source terms to the y-momentum and energy equations:

where s = 15. Note that the sign of gravity is such that the light fluid flows up into the heavy fluid. Reflective wall boundary conditions are imposed at all boundaries using mirror states, which results in an entropy stable scheme under the Lax-Friedrichs flux [14, 50]. Figure 2 shows snapshots of the density for a degree N = 3 entropy stable Gauss DG scheme on a mesh of 32 × 128 elements at various times. End times for each simulation can be found in Table 2.

FIGURE 2. Density for a Rayleigh-Taylor instability for a degree N = 3 entropy stable Gauss DG scheme on a mesh of 32 × 128 elements.

TABLE 2. End time for simulations of the Kelvin-Helmholtz instability on quadrilateral and triangular meshes. On quadrilateral meshes, “collocation” refers to a nodal DGSEM discretization, while “entropy projection” refers to a method based on Gauss nodes. On triangular meshes, “collocation” refers to nodal SBP discretization, while “entropy projection” refers to a modal entropy stable DG method.

3.1.3 Two Dimensional Richtmeyer-Meshkov Instability

The Richtmeyer-Meshkov instability generates small-scale flow features by passing a shock over a stratified fluid [47, 51]. The domain for this setup is [0, 40/3] × [0, 40], and the initial density and pressure are given by

where we again set the slope s = 15. The initial velocities are both set to zero, i.e., u, v = 0. We approximate the discontinuous initial condition using smoothed Heaviside functions with a slope of s = 2 due to the size of the domain. Reflective wall boundary conditions are imposed everywhere. Figure 3 shows pseudocolor plots of the density using a degree N = 3 entropy stable Gauss DG on a uniform mesh of 32 × 96 quadrilateral elements. End times for each simulation can be found in Table 3.

FIGURE 3. Density for the Richtmeyer-Meshkov instability using a degree N = 3 entropy stable Gauss DG with 32 × 96 elements. The domain is [0, 40/3] × [0, 40].

TABLE 3. End time for simulations of the Rayleigh-Taylor instability on quadrilateral and triangular meshes. On quadrilateral meshes, “collocation” refers to a nodal DGSEM discretization, while “entropy projection” refers to a method based on Gauss nodes. On triangular meshes, “collocation” refers to nodal SBP discretization, while “entropy projection” refers to a modal entropy stable DG method.

3.1.4 Three-Dimensional Kelvin-Helmholtz Instability

For completeness, we also verify that a difference in robustness is observed for instability-type problems in three dimensions. Due to the high computational cost of entropy stable DG methods on tetrahedral meshes, we restrict ourselves to hexahedral meshes for the following experiments. We adapt the Kelvin-Helmholtz instability to three dimensions using the following initial condition:

where B is defined as in Eq. 5. Table 3 shows the results, which are similar to previous results for the two-dimensional test problems. We note for this example, both the relative and absolute adaptive time-step tolerances were set to 10^–8 instead of 10^–7. This was necessary to avoid crashes for the entropy projection method at degrees N = 6 and N = 7 on the finer N_cells = 32 mesh.

3.2 Ideal GLM-MHD Equations

Next, we consider the ideal GLM-MHD equations. These equations use generalized Lagrange multiplier (GLM) technique to evolve towards a solution that bounds the magnetic field divergence. When the magnetic field divergence is non-zero, the GLM-MHD system requires the use of non-conservative terms to achieve entropy stability and to ensure Galilean invariance in the divergence cleaning technique.

The non-conservative GLM-MHD system without source terms reads

where the state variables are density, momentum, total energy, magnetic field, and the so-called divergence-correcting field, $\mathit{u}=(\rho ,\rho \mathbf{v},E,\mathsf{B},\psi )$, and the vectors $\mathbf{v}=(u,v,w)$ and $\mathsf{B}=(B_1,B_2,B_3)$ contain the velocities and magnetic field in x, y and z, respectively. The flux reads

where $\mathbf{I}$ is again the 3 × 3 identity matrix, γ is the heat capacity ratio, c_h is the hyperbolic divergence-cleaning speed, and $p=(\gamma -1)\left(E-\left(\rho \Vert \mathbf{v}{\Vert }^2-\Vert \mathbf{B}{\Vert }^2-{\psi }^2\right)/2\right)$ is the gas pressure. Finally, the non-conservative term reads

To construct a two-dimensional version of the GLM-MHD system, we replace $\mathsf{I}$ by a rectangular 3 × 2 identity matrix and neglect the flux in z. However, we keep the third component of the velocity and magnetic field because plasma systems admit three-dimensional electromagnetic interactions in two-dimensional problems. For details about the GLM-MHD system, we refer the reader to [52].

The non-conservative GLM-MHD system (Eq. 6) can be discretized using the collocation (Eq. 2) and modal (Eq. 3) formulations by replacing the volume term $\mathsf{F}$ and the surface term ${\mathsf{f}}^{\ast }$ [53]. In the collocation formulation the new terms read

and in the modal formulation they read

In addition to the symmetric two-point flux f_S, we use a non-symmetric two-point term ${\mathbf{\Phi }}^{\diamond }$ to account for the non-conservative term in the equation. The following experiment uses the non-conservative term presented by Rueda-Ramírez et al. [53] and the entropy conservative flux of Hindenlang and Gassner [54], which is a natural extension of the entropy conservative, kinetic energy preserving, and pressure equilibrium preserving Euler flux of Ranocha [43, 44] to the GLM-MHD system.

3.2.1 Two Dimensional Magnetized Kelvin-Helmholtz Instability

To test the robustness of entropy projection schemes for the GLM-MHD system, we propose a modification of the Euler two-dimensional Kelvin-Helmholtz instability of Section 3.1.1. The domain is [−1, 1]² with the initial condition:

where B(x, y) is as defined in Eq. 5. Each solver is run until final time T_final = 15.

For this example, we set c_h as the maximum wave speed in the domain for the initial condition (Eq. 10) and keep it constant throughout the simulation. This standard way of selecting c_h has been shown to control the divergence error efficiently without affecting the time-step size [52, 55]. We observed that smaller values of c_h affect the robustness of the schemes for this problem, and higher values of c_h increase the stiffness of the problem which can also lead to a crash if the tolerance for the adaptive time-stepping method is set too loosely.

Figure 4 shows pseudocolor plots of the density at T = 10 for the magnetized Kelvin-Helmholtz instability problem obtained with the entropy stable Gauss DG using polynomial degrees N = 3 and N = 7 on uniform meshes of 64 × 64 and 32 × 32 quadrilateral elements, respectively. A comparison with Figure 1 shows that the addition of a vertical magnetic field extends the flow features in the y direction and suppresses many of the vortical structures at T = 10. MHD turbulence eventually develops in the domain after T = 10, which leads to later crash times for this example. End times for each simulation can be found in Table 3.

FIGURE 4. Snapshots of density for the magnetized Kelvin-Helmholtz instability using an entropy stable Gauss DG scheme on uniform quadrilateral meshes.

3.3 Overview of Results

Tables 1,2,4,5 show what time the solver ran until for each solver on both quadrilateral and triangular meshes. We observe the pattern that, for degree N > 1, entropy stable methods which utilize the entropy projection appear be more robust than collocation-type schemes. Moreover, this pattern appears to hold independently of the polynomial degree and mesh size.

TABLE 4. End time for simulations of the Richtmeyer-Meshkov instability on quadrilateral and triangular meshes. On quadrilateral meshes, “collocation” refers to a nodal DGSEM discretization, while “entropy projection” refers to a method based on Gauss nodes. On triangular meshes, “collocation” refers to nodal SBP discretization, while “entropy projection” refers to a modal entropy stable DG method.

TABLE 5. End time for simulations of the magnetized Kelvin-Helmholtz instability on quadrilateral and triangular meshes. On quadrilateral meshes, “collocation” refers to a nodal DGSEM discretization, while “entropy projection” refers to a method based on Gauss nodes. On triangular meshes, “collocation” refers to nodal SBP discretization, while “entropy projection” refers to a modal entropy stable DG method.

3.4 Dependence of Robustness on Atwood number

While the numerical results in the previous section indicate a difference between different entropy stable schemes, they do not provide insight into why and when this difference in robustness manifests. The goal of this section is to establish a relationship between robustness, the Atwood number (a measure of the density contrast), and the use of the “entropy projection” in an entropy stable scheme. We restrict our focus to the Kelvin-Helmholtz instability for this section.

The results presented so far are somewhat unexpected, as the robustness of high order entropy stable DG schemes has been documented for a variety of flows where shocks and turbulent features are present [7–9, 13]. In this section, we conjecture that the documented differences in robustness are due to the presence of both small-scale under-resolved features and significant variations in the density. For example, entropy stable DGSEM methods are known to be very robust for the Taylor-Green vortex, where the density is near-constant throughout the duration of the simulation.

We examine the connection between density contrast and robustness by parametrizing the initial condition by the Atwood number. Given a stratified fluid with two densities ρ₁, ρ₂, the Atwood number is defined as

where it is assumed that ρ₂ ≥ ρ₁. For a constant-density flow, A = 0, while A → 1 indicates a flow with very large density contrasts. We investigate the behavior of different entropy stable methods for a version of the Kelvin-Helmholtz instability parametrized by the Atwood number A:

Figure 5 shows the crash times for the Kelvin-Helmholtz instability using various entropy stable solvers at polynomial degrees 3 and 7. For quadrilateral meshes, we utilize entropy stable DGSEM solvers and entropy stable Gauss DG solvers. For triangular meshes, we utilize entropy stable multi-dimensional SBP solvers and entropy stable modal DG solvers. The DGSEM and SBP solvers are collocation-type schemes, while Gauss and modal DG solvers introduce the entropy projection.

FIGURE 5. Final times a solver an until as a function of Atwood number for the Kelvin-Helmholtz instability for DGSEM and various entropy stable solvers. End times less than final time T_final = 10 indicate a crash.

For degree 3 quadrilateral solvers, we utilize a 32 × 32 mesh, while for degree 7 quadrilateral solvers, we utilize a 16 × 16 mesh. The mesh resolution is halved for polynomial degree 7 simulations so that the total number of degrees of freedom is kept constant. For triangular solvers, we again use 32 × 32 and 16 × 16 uniform meshes, but we compare polynomial degrees 3 and 6, as SBP quadrature rules are available only up to degree 6 in Trixi.jl. We run up to time T_final = 10 for A ∈ [0.1, 0.9] and report the times each simulation ran until. For degree N = 3, we observe that schemes which involve the entropy projection runs until the final time T_final = 10. Collocation-type schemes run to completion for low Atwood numbers, but crash earlier and earlier as the Atwood number increases. At degree N = 7, we observe that while both collocation solvers and entropy projection solvers crash at higher Atwood numbers, entropy projection solvers begin crashing at higher Atwood numbers. For example, on quadrilateral meshes, DGSEM crashes around Atwood number 0.3, while Gauss solvers crash around Atwood number 0.7. We note that crash times for entropy projection schemes also tend to depend on the adaptive time-stepping tolerance. For example, for N = 3 and a 32² mesh, Gauss collocation runs stably to T_final = 10 if the absolute and relative tolerances are reduced to 10^–9. The same is not true of entropy stable collocation-type schemes.

To provide another point of comparison, we ran simulations using an entropy stable DGSEM solver with sub-cell finite volume shock capturing [56] with Zhang-Shu positivity-preserving limiting for the density and pressure [57, 58], which we refer to as DGSEM-SC-PP for shock capturing and positivity preservation.² The entropy stable sub-cell finite volume-based shock capturing scheme utilizes a blending coefficient parameter α ≤ α_max [56]. For these experiments, we set α_max = 0.005, which implies that the low order finite volume solution constitutes at most 0.5% of the final blended solution. Despite the fact that this shock capturing is very weak, the resulting solver greatly improves robustness and enables long simulation times: for N = 3 and a 32 × 32 mesh, DGSEM-SC-PP runs stably to time T_final = 10 for Atwood numbers up to 0.99. However, we have also observed that the minimum value of α_max necessary to avoid solver failure depends on the mesh resolution. For example, for N = 3 and a 64 × 64 mesh, we observe that DGSEM-SC-PP with α_max = 0.005 crashes around t = 6.4871.

Remark 2. We note that DGSEM with α_max = 0.005 shock capturing but no positivity preservation is not robust for the Kelvin-Helmholtz instability. For the initial condition (Eq. 4), N = 3, and a 64 × 64 mesh, DGSEM with shock capturing crashes around time t = 4.8891. For N = 7 and a 32 × 32 mesh, DGSEM with shock capturing crashes around time t = 5.0569. In contrast, DGSEM with only positivity preservation results in the simulation stalling due to a very small time-step.

4 The Role of the Entropy Projection

4.1 Is Robustness Due Only to the Entropy Projection?

While the numerical results up to this point indicate that there is a significant difference in robustness for different entropy stable schemes, it is not yet clear that the increased robustness is due to the entropy projection. For example, the numerical experiments in Section 3 compare entropy stable Gauss DG schemes to DGSEM on tensor product meshes and entropy stable “modal” DG methods to SBP schemes on triangular meshes. In both cases, a collocation scheme is compared to a scheme with higher accuracy numerical integration. Thus, it is not immediately clear whether the difference in robustness is due to the entropy projection or other factors such as the quadrature accuracy. We investigate whether the quadrature accuracy has a significant effect on stability by testing two additional variants of entropy stable DGSEM schemes on quadrilateral meshes. These schemes are purposefully constructed to be “bad” methods (in the sense that they introduce additional work without improving the expected accuracy), and are intended only to introduce the entropy projection. Both have quadrature accuracy similar to or lower than entropy stable DGSEM methods.

The first scheme utilizes LGL points for volume quadrature, but utilizes (N + 1) point Clenshaw-Curtis quadrature at the faces. This scheme can be directly derived from a modal formulation and (despite the lower polynomial exactness of Clenshaw-Curtis quadrature) can be shown to be entropy stable on affine quadrilateral meshes using the analysis in [21]. In order to retain entropy stability, the solution must be evaluated using the entropy projection at face nodes. We argue that the use of Clenshaw-Curtis quadrature does not result in a significant increase in quadrature accuracy over LGL quadrature: while Clenshaw-Curtis quadrature has been shown to be similar to Gauss quadrature for integration of analytic functions [61], for lower numbers of points we observe that the accuracy is comparable to LGL quadrature. Moreover, it was argued in [62] that increasing quadrature accuracy only for surface integrals or only for volume integrals does not provide sufficient anti-aliasing. We refer to this method as “DGSEM with face-based entropy projection” in Figure 6.

FIGURE 6. Degree N = 3 and 64 × 64 grid Kelvin-Helmholtz simulations at T = 5. All methods run until T = 25, while DGSEM crashes at T ≈ 3.5.

Remark 3. We note that one can also build an entropy stable scheme from a combination of LGL volume points and Gauss face points. While this method possesses much of the simplicity and advantageous features of entropy stable DGSEM methods while also displaying improved robustness, this method results in a suboptimal rate of convergence by one degree [21].The second scheme we test is similar to the staggered scheme of [17]. However, while the original scheme of Parsani et al. combines degree N Gauss points with degree (N + 1) LGL points, we combine degree N Gauss points with degree N LGL points. This is a “useless” staggering in that it does not increase the accuracy of integration compared with DGSEM, and is intended only to introduce the entropy projection into the formulation.³ We refer to this method as “DGSEM with volume-based entropy projection” in Figure 6.Figure 6 shows snapshots of the density for the Kelvin-Helmholtz instability for a degree N = 3 mesh of 64 × 64 elements for each method. While the plots for the Gauss DG and DGSEM with face-based entropy projection have qualitative similarities, we observe that DGSEM with volume-based entropy projection results in a noisier solution. This may be due to inconsistency in terms of accuracy between the two quadrature rules used (e.g., (N + 1) point LGL and Gauss quadratures). However, all three entropy projection schemes remain stable, and we have verified that they are able to run until T = 25 without crashing.We also compute crash times for each method for the Kelvin-Helmholtz instability with Atwood numbers A ∈ [0.1, 0.9]. These crash times are also compared to crash times of an entropy stable DGSEM method. These computations are performed on both a degree N = 3 mesh of 32 × 32 elements, as well as a degree N = 7 mesh of 16 × 16 elements. Figure 7 plots the crash times for each method. We observe that all schemes which involve the entropy projection run stably for a wider range of Atwood numbers than entropy stable DGSEM, and that this effect becomes even more pronounced for degree N = 7. However, for both the N = 3 and N = 7 experiments, the entropy stable Gauss schemes are stable for the widest ranges of Atwood numbers.These results indicate that incorporating the entropy projection does have a significant effect on the robustness of an entropy stable method, but that the entropy projection is not the only factor which impacts robustness. However, a detailed analysis of factors such as quadrature accuracy is out of the scope of this current work.

FIGURE 7. Final times a solver an until as a function of Atwood number for the Kelvin-Helmholtz instability for DGSEM and different variants of entropy stable solvers based on the entropy projection. End times less than final time T_final = 10 indicate a crash.

4.2 Why Is There a Difference in Robustness for Different Entropy Stable Methods?

While the results from previous sections suggest that the entropy projection plays a role in the robustness of an entropy stable scheme, it is not clear why it plays a role. While we do not have a thorough theoretical understanding of the entropy projection, initial experiments indicate that entropy projection schemes behave most differently from collocation schemes when the solution is either under-resolved or have near-zero density or pressure.

We illustrate the aforementioned behavior of the entropy projection using the one-dimensional compressible Euler equations. The conservative variables for the Euler equations are density, momentum, and total energy (ρ, ρu, E). Let s(u) = log(p/ρ^γ) denote the specific entropy. The entropy variables for the convex entropy S(u) = −ρs(u)/(γ − 1) are given by

Recall that the main steps of the entropy projection are as follows:

(1) Evaluate the entropy variables using degree N polynomial approximations of the conservative variables

(2) Compute the quadrature-based L² projection of the entropy variables to degree N polynomials

(3) Re-evaluate the conservative variables in terms of the projected entropy variables.

These re-evaluated conservative variables are then used to compute contributions from an entropy stable DG formulation.

It was demonstrated numerically in [19] that the entropy projection is high order accurate for sufficiently regular solutions. However, the behavior of the entropy projection was not explored for under-resolved or near-vacuum solution states. We illustrate this behavior using the following solution state:

where p_min > 0 is the minimum pressure, and k is a parameter which controls the frequency of oscillation. As k increases, the solution states in Eq. 11 become more and more difficult to resolve, and as p_min → 0, the solution approaches vacuum and the entropy approaches non-convexity.

Figure 8 illustrates the effect of increasing k and decreasing p_min on the entropy projected conservative variables for a degree N = 2 approximation on a coarse mesh of eight elements. As k increases and the solution becomes under-resolved, the entropy projection develops large jumps at the interface. Similarly, as p_min decreases from 1 to 1/10, the entropy projection develops large jumps at the interface. We note that for both increased k and decreased p_min, spikes do not appear in the interior of the element.

FIGURE 8. Illustration of the effect of larger k (under-resolution) and smaller p_min (near-vacuum state) on the entropy projection. A degree N = 2 approximation and mesh of 8 elements were used.

This indicates that the error in the entropy projection is influenced by both the numerical resolution and how close the entropy is to becoming non-convex. We denote the continuous entropy projection by $\tilde{\mathit{u}}=\mathit{u}\left({\mathrm{\Pi }}_N\mathit{v}({\mathit{u}}_h)\right)$. Then, by the mean value theorem, we can bound the difference between the conservative and entropy-projected variables

where $\frac{\partial \mathit{u}}{\partial \mathit{v}}$ is evaluated at some intermediate state between u_h and $\tilde{\mathit{u}}$. The latter term in the bound $‖\mathit{v}({\mathit{u}}_h)-{\mathrm{\Pi }}_N\mathit{v}({\mathit{u}}_h)‖$ is small when the entropy variables are well-resolved, which we expect to be true when the solution is well-resolved and the mapping between conservative and entropy variables is well-conditioned. Conversely, high frequency components of the solution are often amplified when v(u) is highly nonlinear or the solution is under-resolved (this is the motivation behind filtering for stabilization [63–65]), and we expect $‖\mathit{v}({\mathit{u}}_h)-{\mathrm{\Pi }}_N\mathit{v}({\mathit{u}}_h)‖$ to be large for such settings. The former term $‖\frac{\partial \mathit{u}}{\partial \mathit{v}}‖$ is large when the mapping between conservative and entropy variables is nearly singular, which occurs when the entropy is nearly non-convex (for example, near-vacuum states).

4.2.1 What Role Does Entropy Dissipation Play?

The previous section illustrates that entropy projection schemes are likely to differ from collocation schemes most when the solution is under-resolved or has near-zero density or pressure. Moreover, since the entropy projected variables in Figure 8 display spikes at the interfaces, it seems possible that the entropy projection would change the manner in which entropy dissipative interface dissipation terms are triggered. To test this hypothesis, we compute the evolution of entropy over time for the Kelvin-Helmholtz instability using both entropy stable Gauss DG and DGSEM-SC-PP, which is an entropy stable DGSEM with a shock capturing technique that consists in blending a sub-cell finite volume scheme with the DGSEM in an element-wise manner [56] and Zhang-Shu’s positivity preserving limiter [57, 58]. The blending of the finite volume scheme is capped at 0.5% in order to avoid unnecessary numerical dissipation. We also compare entropy evolution for a scheme that blends a sub-cell finite volume scheme with the DGSEM in a subcell-wise manner [66], which we refer to as DGSEM-subcell. The blending factors are chosen for each node (or subcell) to enforce lower bounds on density and pressure based on the low order solution, ρ ≥ 0.1 ρ^FV, p ≥ 0.1 p^FV. For this choice of lower bound, we observe high order accuracy for a two-dimensional sinusoidal entropy wave [67]. While this scheme is not provably entropy stable, it was demonstrated numerically in [66] that the use of subcell blending factors requires significantly lower levels of limiting compared with an element-wise limiting factor.

Figure 9 shows the evolution of the integrated entropy over the entire domain (which we have shifted to be positive) for the Kelvin-Helmholtz instability. Since periodic boundary conditions are used, the integrated entropy for the semi-discrete formulation can be proven to decrease over time. We observe that all four methods display similar entropy dissipation behavior until time t ≈ 1.2, after which DGSEM shows less entropy dissipation than either Gauss or DGSEM-SC-PP. However, while DGSEM-SC-PP initially dissipates more entropy than Gauss DG, the entropy dissipation for Gauss DG increases and overtakes that of DGSEM-SC-PP around time t ≈ 4. Since entropy dissipation in both Gauss DG and DGSEM-SC-PP schemes is triggered by under-resolved flows (either through a modal indicator or through jump penalization terms) and since the Kelvin-Helmholtz instability generates increasingly small scales at larger times, this suggests that entropy dissipation for Gauss DG may be activated more strongly but at smaller scales than DGSEM-SC-PP. In contrast, Gauss DG dissipates more global entropy than DGSEM-subcell, though DGSEM-subcell eventually catches up to Gauss DG for N = 3.

FIGURE 9. Evolution of entropy over time for the Kelvin-Helmholtz instability.

Our initial hypothesis was that the entropy projection in Gauss DG schemes results in larger interface jumps, which would trigger more entropy dissipation through jump penalization terms. However, this does not appear to be consistent with numerical results for entropy conservative schemes. To test these schemes, we focus on the three-dimensional Taylor-Green vortex. We note that the observed loss of robustness stands in stark contrast to the observed robustness of high order entropy stable and split-form DGSEM for the Taylor-Green vortex [8, 13, 22]. This can be explained by the fact that the density remains near-constant over time for the Taylor-Green vortex; for a Kelvin-Helmholtz initial condition with a constant density, DGSEM runs stably up to final time T = 25 for each of the previous numerical settings. Thus, while the Taylor-Green vortex generates small-scale flow features, it is a more benign test case when evaluating the robustness of high order entropy stable DG schemes.

However, when using a purely entropy conservative scheme (which can be constructed by utilizing entropy conservative interface fluxes), DGSEM methods can display non-robust behavior for the Taylor-Green vortex. We run the Taylor-Green vortex to final time T_final = 20 using a variety of entropy conservative schemes: DGSEM, Gauss DG, as well as an entropy stable C⁰ continuous Galerkin spectral element method (CGSEM) and a periodic finite difference method. We note that, because an entropy conservative scheme can be constructed given any summation-by-parts or skew-symmetric operator [12, 14, 26], we are able to implement an entropy conservative C⁰ continuous spectral element method and periodic finite difference method by constructing global difference operators from the tensor product of one-dimensional operators. These one-dimensional operators are provided by the Julia library SummationByPartsOperators.jl [68].

Table 6 shows the end simulation time for each solver. We observe that again, despite the absence of any entropy dissipation, the Gauss DG solver is more robust than the DGSEM solver. The continuous spectral element solver CGSEM is also significantly more robust than the DGSEM solver, though it does lose robustness at higher orders and finer grid resolutions. We also ran periodic finite difference operators for grids with 4, 6, 8, 10, 12 nodes in each dimension with orders of accuracy 2, 4, 6, 8, 10. We observe that the periodic finite difference operator is as robust as the Gauss DG solver: for every grid resolution and order specified, the finite difference solver ran up to the final time T_final = 20.

TABLE 6. End time for entropy conservative simulations of the Taylor-Green vortex on hexahedral meshes.

These experiments indicate that robustness for schemes involving the entropy projection is not solely due to the entropy dissipative terms. These experiments also show that robustness is improved for CGSEM and periodic finite difference solvers, neither of which contains interface terms. Since these results are on relatively coarse resolutions and utilize an entropy conservative scheme (when most practical schemes are entropy stable), further numerical experiments are necessary to carefully analyze the effect that different discretizations have on robustness.

5 Applications Toward Under-Resolved Simulations

We conclude the paper with a discussion on a comparison between three schemes which include dissipative terms (entropy stable Gauss DG, entropy stable DGSEM-SC-PP, and DGSEM-subcell) for an under-resolved simulation. We run the Kelvin-Helmholtz instability using the initial condition (Eq. 4), but modify the y-velocity perturbation to break symmetry of the resulting flow

We run the simulation up to final time T_final = 25. We use both a degree N = 3 mesh of 64 × 64 elements and a degree N = 7 mesh of 32 × 32 elements, each of which contains the same number of degrees of freedom. Due to the sensitivity of the Kelvin-Helmholtz instability problem and the long time window of the simulation, the results for each scheme are qualitatively very different.

Figures 10, 11 show snapshots of density and pressure for the entropy stable DGSEM-SC-PP and Gauss DG schemes. We observe that in both cases, the flow scales present in the DGSEM-SC-PP scheme are noticeably larger than those observed in the Gauss scheme. This is notable because the DGSEM-SC-PP scheme applies a very small amount of shock capturing: dissipation is added by blending the high order scheme with a low order finite volume scheme, and the amount of the blended low order solution is capped at 0.5%. However, even a small amount of dissipation produces a noticeable change on small-scale features in the resulting flow. We also observe the presence of shocklets or compression waves in the pressure, which mirror observations made in [69].⁴

FIGURE 10. Density and pressure for the Kelvin-Helmholtz instability at T_final = 25 on a N = 3 mesh of 64² elements.

FIGURE 11. Density and pressure for the Kelvin-Helmholtz instability at T_final = 25 on a N = 7 mesh of 32² elements.

For N = 3, the scales observed in DGSEM-subcell scheme are noticeably smaller than those of DGSEM-SC-PP but similar to those of the Gauss DG scheme. For N = 7, the scales observed in the DGSEM-subcell scheme are again smaller than those of DGSEM-SC-PP, but appear to be slightly larger than those of the Gauss DG scheme. To avoid qualitative speculation, we compare these flows by computing the angle-averaged power spectra of the velocity weighted by $\sqrt{\rho }$ at final time T_final = 25 [70, 71]. We follow [3, 7] and generate a grid of uniformly spaced points by evaluating the degree N polynomial solution at (N + 1) equally spaced points along each dimension in the interior of each element of a uniform Cartesian mesh. The power spectra can then be computed from a fast Fourier transform of the resulting data. Figure 12 shows the power spectra, which appear consistent with a k^−7/3 rate of decay from two-dimensional turbulence theory [71]. Moreover, we observe that the entropy stable Gauss DG scheme retains more energetic information than both DGSEM-SC-PP and DGSEM-subcell, though a spurious spike in the energy for Gauss DG schemes is observed near the higher wavenumbers for N = 3.

FIGURE 12. Weighted power spectra for DGSEM-Subcell, entropy stable DGSEM-SC-PP, and entropy stable Gauss DG schemes.

6 Conclusion

This paper shows that for variable density flows which generate small-scale features, there are differences in robustness between entropy stable schemes which incorporate the entropy projection and those which do not. These differences in robustness are observed to depend on the Atwood number (measuring the density contrast) and persist across a range of polynomial degrees, mesh resolutions, and types of discretization. However, the mechanisms behind improved robustness for entropy projection schemes are currently unknown.

We note that any conclusions drawn concerning the robustness of DGSEM and Gauss DG should be restricted to the instability-type problems studied here. These results do not imply that Gauss is uniformly more robust than DGSEM. Moreover, Gauss schemes are more computationally expensive than DGSEM schemes and result in smaller maximum stable timesteps [22, 72–74], so the appropriate scheme will depend on the use case.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/trixi-framework/paper-2022-robustness-entropy-projection.

Author Contributions

All authors contributed to the conception and design of the paper. JC, HR, and AR-R contributed numerical experiments. JC drafted the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

HR was partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2044-390685587, Mathematics Münster: Dynamics-Geometry-Structure. This work has received funding from the European Research Council through the ERC Starting Grant “An Exascale aware and Uncrashable Space-Time-Adaptive Discontinuous Spectral Element Solver for Non-Linear Conservation Laws” (Extreme), ERC grant agreement no. 714487 (GG and AR-R). TW was supported in part by the Exascale Computing Project, a collaborative effort of two U.S. Department of Energy organizations (Office of Science and the National Nuclear Security Administration) responsible for the planning and preparation of a capable exascale ecosystem, including software, applications, hardware, advanced system engineering, and early testbed platforms, in support of the nation’s exascale computing imperative. TW was also supported in part by the JC Faculty Chair in Science at Virginia Tech. JC gratefully acknowledges support from the National Science Foundation under award DMS-CAREER-1943186. GG and ARR acknowledge funding through the Klaus-Tschira Stiftung via the project “HiFiLab.”

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

This work used the Extreme Science and Engineering Discovery Environment (XSEDE) Expanse at the San Diego Supercomputer Center through allocation TG-MTH200014 [75]. This work was performed on the Cologne High Efficiency Operating Platform for Sciences (CHEOPS) at the Regionales Rechenzentrum Köln (RRZK) and on the group cluster ODIN. We thank RRZK for the hosting and maintenance of the clusters. The authors also thank Fabian Föll for providing the initial condition for the Richtmeyer-Meshkov instability.

Footnotes

¹Since the entries of $\mathsf{F}$ are vector-valued, the Hadamard product $\left(\mathsf{Q}-{\mathsf{Q}}^T\right)◦\mathsf{F}$ should be understood as each scalar entry of $\left(\mathsf{Q}-{\mathsf{Q}}^T\right)$ multiplying each vector-valued entry of $\mathsf{F}$.

²For DGSEM-SC-PP, we utilize a 4-stage 3rd order adaptive strong stability preserving (SSP) Runge-Kutta time-stepping method [59, 60] with stepsize controller and efficient implementation of [40], which is necessary to ensure fully discrete positivity.

³This scheme can also be derived by beginning with an entropy stable DGSEM scheme and replacing the diagonal LGL mass matrix with the fully integrated dense mass matrix computed using Gauss quadrature. The resulting scheme can be made entropy stable by evaluating the spatial formulation using the entropy projection. More specifically, the appropriate entropy projection for this setting interpolates the entropy variables at Gauss nodes, then interpolates to LGL nodes.

⁴We note that these “shocklets” are not strictly shock waves, as the flow is not supersonic.

References