Sec. Computational Physics
Volume 8 - 2020 | https://doi.org/10.3389/fphy.2020.500690
A Novel Robust Strategy for Discontinuous Galerkin Methods in Computational Fluid Mechanics: Why? When? What? Where?
- 1Division of Mathematics, Department of Mathematics and Computer Science, Center for Data and Simulation Science, University of Cologne, Cologne, Germany
- 2Division of Computational Mathematics, Department of Mathematics, Linköping University, Linköping, Sweden
In this paper we will review a recent emerging paradigm shift in the construction and analysis of high order Discontinuous Galerkin methods applied to approximate solutions of hyperbolic or mixed hyperbolic-parabolic partial differential equations (PDEs) in computational physics. There is a long history using DG methods to approximate the solution of partial differential equations in computational physics with successful applications in linear wave propagation, like those governed by Maxwell’s equations, incompressible and compressible fluid and plasma dynamics governed by the Navier-Stokes and the Magnetohydrodynamics equations, or as a solver for ordinary differential equations (ODEs), e.g., in structural mechanics. The DG method amalgamates ideas from several existing methods such as the Finite Element Galerkin method (FEM) and the Finite Volume method (FVM) and is specifically applied to problems with advection dominated properties, such as fast moving fluids or wave propagation. In the numerics community, DG methods are infamous for being computationally complex and, due to their high order nature, as having issues with robustness, i.e., these methods are sometimes prone to crashing easily. In this article we will focus on efficient nodal versions of the DG scheme and present recent ideas to restore its robustness, its connections to and influence by other sectors of the numerical community, such as the finite difference community, and further discuss this young, but rapidly developing research topic by highlighting the main contributions and a closing discussion about possible next lines of research.
1 A Brief Introduction to DG
The first discontinuous Galerkin (DG) type discretisation is either attributed to Reed and Hill in 1973  for an application to steady state scalar hyperbolic linear advection to model neutron transport, or to Nitsche in 1971  who introduced a discontinuous finite element method (FEM) to solve elliptic problems with non-conforming approximation spaces. It was however a series of papers by Cockburn and Shu et al. starting 20 years later [3–6] that introduced the modern form of the so-called Runge-Kutta DG scheme. They applied the method especially to nonlinear hyperbolic problems such as the compressible Euler equations on unstructured simplex grids with slope limiting to capture shocks. Bassi and Rebay were the first that extended the DG method to the compressible Navier-Stokes equations . They used a fully discontinuous ansatz based on a mixed variational formulation, where they rewrote the second order partial differential equation (PDE) into a first order system. The resulting DG formulation requires numerical fluxes for the advective as well as for the diffusive part. Although the methods gave reasonable results for the compressible Navier-Stokes equations, an analysis of the method in Arnold and Brezzi et al. [8, 9] applied to pure elliptic problems revealed how to improve the method in terms of convergence rate, adjoint consistency, and stability. Since the introduction of its modern form, the DG method has been applied and advanced by many researchers across different scientific disciplines around the world. The DG method is used in a wide range of applications such as compressible flows [10–12], electromagnetics and optics [13–16], acoustics [17–21], meteorology [22–25], and geophysics [26, 27]. The first book available on DG was basically a collection of papers . Since then, many different text books on DG are available focusing on theoretical developments as well as specific implementation details, e.g., [29–31].
One of the main applications of DG methods is the discretisation of nonlinear advection-diffusion problems of the form
where u is the vector of conserved quantities, e.g., the mass, momentum, or energy. The vector
with the fluxes fi in each spatial direction xi, i = 1,2,3. The viscous fluxes are denoted by
thus modeling parabolic effects, e.g., heat conduction. The problem is typically defined on a given spatial domain
The DG scheme is based on a Galerkin type weak formulation. For the sake of simplicity, we drop the viscous second order terms in what follows and refer to, e.g. , for a complete description of the advection-diffusion case. To construct the approximation space of the DG method, the domain is split into non-overlapping elements
where P(N) is the number of polynomial basis functions depending on the polynomial degree N. The time dependent polynomial coefficients are Uj(t), and
Similarly, for block vectors,
We obtain the weak formulation by multiplying each equation by a polynomial test function
If we choose the test function
that depends on the two values at the element interface, i.e. on the value inside the element U and outside from the neighbor element U+. Typically, the numerical surface fluxes are constructed from (approximate) Riemann solvers, e.g., . With the surface numerical flux, we arrive at the semi-discrete DG formulation of the advection problem in weak form
or if we apply integration-by-parts once more to the volume terms it becomes the so-called strong DG formulation
where the surface contribution is a penalty between the numerical surface flux and the normal flux evaluated from the interior of an element. Note, the weak and strong form DG formulations are equivalent .
There are still many critical decisions necessary before the DG formulation produces an algorithm that can be implemented. The type of element shape needs to be decided as well as which polynomial basis. For example, modal polynomial basis functions for tetrahedral elements or nodal tensor product polynomials for hexahedral elements. In addition, the surface integral and the volume integral needs to be discretized. In most cases, the integrals are approximated with numerical quadrature rules, e.g., high order Gauss-Legendre quadrature and cubature. Many variants and detailed descriptions can be found in text books on DG methods and their implementation, e.g., [28–31]. It is important to note that these choices involving the element type, basis functions, and approximation of inner products all have a major impact on the performance of the resulting DG scheme in terms of computational complexity and robustness due, e.g., to the presence of spurious oscillations near discontinuities that result in unphysical solution states (like negative density or pressure) or aliasing instabilities. Many mechanisms exist in the DG community to combat spurious oscillations (i.e., shock capturing) such as slope [3, 5, 35] or WENO [36, 37] limiters, filtering [29, 38, 39], finite volume sub-cells [40–42], MOOD-type limiting [43–47], or artificial viscosity [48, 49]. The issue of shock capturing will not be discussed further. However, aliasing errors, how they arise within the DG method and strategies to remove said errors and increase robustness will be discussed at length in this article.
With the above decisions, we arrive at the generic semi-discrete ordinary differential equation (ODE) form of the DG scheme, which can be integrated in time with an appropriate high order explicit or implicit ODE solver, e.g., [50–54].
The resulting DG method is high order accurate and has excellent dispersion and dissipation behavior, e.g., [55, 56]. Furthermore, due to its compact stencil (only interface neighbor data is needed) the DG scheme is well known for its excellent parallel scaling, e.g., [57, 58], and its ability to handle unstructured and non-conforming grids, e.g., [16, 54, 59–63]. These nice properties of the DG methodology are one reason why more and more researchers apply and extend the DG methodology to many different problem setups in computational physics. However, DG is not the perfect discretisation and there are unfortunately some issues that necessitate detailed analysis and discussion.
The remainder of this review article gives the answers to why we need novel developments, Section 2, when the novel developments started, Section 3, what the key ideas of these novel strategies are, Section 4, and where there are still open questions toward future research directions, Section 5.
2 Why Do We Need a Novel Robust Strategy?
Throughout the analysis and discussions in this manuscript we describe different types of stability for a numerical approximation. Principally, we concentrate on the stability and boundedness of the spatial DG discretization.
2.1 On the L2-Stability of the DG Method
It is easy to show that the DG scheme is L2-stable for linear advection problems with constant coefficients due to its Galerkin nature, e.g., as a special case in Ref. 64. As a brief illustrative example let’s consider the one-dimensional scalar linear advection model
where f(u) = au with constant velocity a. The respective strong form DG scheme reads
We get the discrete evolution of the L2-norm
Assuming time continuity, the first term reduces to
The volume term contribution is
which shows that the volume contribution can be lifted to the boundary. In total we have
The discrete evolution of the L2-norm only depends on the choice of the numerical flux function f*. A simple choice would be the central flux
Summing over all elements E in the domain to get the total L2-norm and assuming periodic boundary conditions, we get
For linear problems basically all DG variants are stable, but what about nonlinear problems? To address nonlinear problems, there was a crucial contribution by Jiang and Shu in 1994 , who demonstrated that for scalar nonlinear hyperbolic problems, the DG method is L2-stable provided: 1) Exact evaluation of all integrals are used; 2) Entropy stable numerical surfaces fluxes are used at the element interfaces. The L2-stability result with conditions 1) and 2) also extends to symmetric variable coefficient hyperbolic systems . Again, as a brief example to illustrate the important steps of the analysis, we consider a scalar one-dimensional problem Eq. 11 with a simple quadratic flux function
The resulting discrete evolution of the L2-norm is
which agin only depends on the choice of the numerical surface flux function f*. Note, that for the central flux choice
Following the same arguments as above for the linear advection problem, we sum over all of the elements in the domain and obtain L2-stability for the nonlinear scalar hyperbolic problem for quadrature rules with at least 3N−1 integration precision.
Unfortunately, even ignoring the practical issues with the assumption of exact integration for a moment, the results of Jiang and Shu cannot be directly extended to general nonlinear hyperbolic systems, e.g., the compressible Euler equations. A key step in the analysis of Jiang and Shu is that the test functions ϕ in the DG formulation Eq. 10 are replaced with the discrete DG solution
which itself is a linear combination of the test functions
2.2 On the Entropy Stability of the DG Method
In the previous section it was shown for a scalar nonlinear conservation law that the DG solution U is bounded in the L2-norm if a special choice of the numerical surface flux is chosen, such as the one in Eq. 20. We conjecture that an analogous statement of nonlinear stability should be true for systems of nonlinear conservation laws. However, in general, stability in the L2-norm is insufficient to exclude unphysical phenomena such as expansion shocks . To remove the possibility of such phenomena we generalize the notion of a stability estimate for nonlinear problems.
Before we discuss the mathematics of a general nonlinear hyperbolic system, a detour is taken to examine an important underlying physical principle. In particular, we introduce concepts from thermodynamics, which is a branch of physics relating the heat, temperature, or entropy of a given physical system to energy and work. The laws of thermodynamics are some of the most fundamental laws in all of physics. This is because they play an important role in describing how, and predicting why, physical systems behave and evolve the way that we observe them. Moreover, thermodynamics provides fundamental rules to decide how a physical system cannot behave. That is, what type of solution behavior is physically meaningful and what is not. From a mathematical point of view, we note that satisfying the second law of thermodynamics is not enough to guarantee uniqueness of the PDE solution and that conditions for uniqueness are an active topic of research in the analysis of said PDEs, see for example [69–72].
The first law of thermodynamics concerns the conservation of the total energy in a closed system. The second law of thermodynamics states that the entropy of a closed physical system tends to increase overtime and, importantly, that it cannot shrink. The laws of thermodynamics must be satisfied simultaneously at all times, otherwise a mathematical solution can exhibit strange and obviously incorrect behavior. For example, a fluid that only conserves its total energy but does not take care on the entropy, i.e. satisfying the first law of thermodynamics but not the second law, could transfer all of its internal energy into kinetic energy. The result would be a very fast, but very cold jet of air. Such a flow configuration has never been observed in nature. This discrepancy is removed when incorporating the second law of thermodynamics where the transfer of energies are regulated. For reversible processes the entropy remains constant over time (isentropic) and the time derivative of the total system entropy is zero. For irreversible processes the entropy increases and the time derivative is positive. Solution dynamics where the total system entropy shrinks in time are never observed and deemed unphysical.
To discuss these ideas in a mathematical context consider a system of nonlinear hyperbolic conservations laws
where we take the viscous flux components in Eq. 1 to be zero. Typically, the diffusion terms are dissipative in nature and are mostly uncritical. A prototypical example of a purely hyperbolic system of nonlinear conservation laws is the compressible Euler equations modeling inviscid gas dynamics. A smooth solution that satisfies the system of PDEs Eq. 23 corresponds to a reversible process. One of the difficulties, either analytically or numerically, of nonlinear hyperbolic PDEs is that the solution may develop discontinuities (e.g., shocks) regardless of the continuity of the initial conditions . A discontinuous solution of Eq. 23 corresponds to an irreversible process and dissipates entropy.
To mathematically account for possible discontinuous solutions, system Eq. 23 is considered in its weak form. Just as in the Galerkin discretisation in Section 1, the weak form of the PDE is found by multiplying the governing equations by a smooth test function
Another form of the conservation law is its integral form that, under the assumption of differentiable fluxes, arises from Gauss’ theorem
and holds for arbitrary control volumes
It turns out that there are additional conserved quantities, e.g., the entropy, which are not explicitly built into the nonlinear hyperbolic system Eq. 23 but are still a consequence of the PDE. In order to reveal this auxiliary conservation law we define a convex (mathematical) entropy function s = s(u) that is a scalar function and depends nonlinearly on the conserved variables u. This allows the definition of a new set of entropy variables
that provides a one-to-one mapping between the conservative variable space and entropy space . If we contract the nonlinear hyperbolic system Eq. 23 from the left with the entropy variables w we have
From the definition of the entropy variables and assuming continuity in time, we know that
Further, each of the flux vectors in the coordinate directions xi must satisfy a compatibility condition
where fis, i = 1,2,3 is a corresponding entropy flux . We point out that a chain rule is the linchpin of the manipulations Eqs. 28 and 29 to move from the space of conservative state variables into the space of entropy variables. In the continuous setting this is not an issue under certain continuity assumptions. However, in a numerical setting it is extraordinarily difficult (or even impossible) to recover the chain rule with discrete differentiation, e.g., . We postpone the discussion on this issue and what it means for a high order DG numerical approximation to Section 4.
By definition of a mathematical entropy, contracting the nonlinear hyperbolic system into entropy space, as in Eq. 27, and assuming the solution is smooth (i.e. a reversible process) results in the auxiliary conservation law for the entropy
The corresponding integral form of the entropy conservation is given by
for an arbitrary volume
for discontinuous solutions.
As an illustrative example for the contraction of a nonlinear hyperbolic system of PDEs into entropy space, we consider the compressible Euler equations of gas dynamics in one spatial dimension
with the ideal gas assumption
where γ is the adiabatic constant. The convex entropy function s(u) for the compressible Euler equations is not unique . However, a common choice for the mathematical entropy function is the scaled negative thermodynamic entropy [80–84]
with the corresponding entropy flux fs = v s(u). From this definition of the mathematical entropy function we get the entropy variables
and see that each component of the entropy variables is a highly nonlinear function of the state vector components. Regardless, the variables Eq. 36 contract the one dimensional compressible Euler equations into entropy space and, when integrated over the domain
The discrete equivalent of the entropy inequality Eq. 32 is referred to as entropy stability. It is a generalization of the L2-stability statement to systems of nonlinear hyperbolic PDEs, e.g., [79, 86]. An additional requirement built into the entropy stability condition is that the fluxes
At present, we restrict ourselves to one of the key ingredients in the analysis to derive a nonlinear entropy stability estimate for general nonlinear hyperbolic systems, see e.g., [78, 86]. It is natural to develop an entropy stable DG approximation because the continuous and discrete analysis both rely on a weak form of the governing equations. However, for entropy stability the nonlinear system is not multiplied by the solution u as was the case for the L2-stability analysis. Instead the equation is multiplied with the entropy variables wEq. 26, which are nonlinear functions of the state u, see e.g., the compressible Euler equations in (36). Thus, a direct combination of this approach with the analysis of Jiang and Shu from Section 2.1 is not possible. For a polynomial DG ansatz U, the discrete entropy variables W = W(U) are no longer polynomials of degree N and do not belong to the space of test functions
To overcome the issue with the test function space and to enable an entropy stability estimate for the DG method, Hughes et al.  as well as Hildebrand and Mishra et al. [92–94] used a space-time DG approach with an ansatz directly written in terms of the entropy variables. The idea is to make the DG ansatz in entropy space, i.e. to approximate the entropy variables
with a polynomial of degree N. Ignoring the time discretisation for brevity, the DG formulation changes to
which shows that the scheme is still formulated in conservative form, however all the conserved variables u now depend implicitly on the polynomial approximation of the entropy variables W. As this approach is naturally implicit, a straightforward and elegant extension of the scheme in time is to use a temporal DG scheme on top of the spatial DG scheme, resulting in a fully implicit space-time DG formulation. Hildebrand and Mishra proved that the resulting discretisation is entropy stable provided: i) Exact evaluation of all integrals; ii) Entropy stable numerical fluxes at the spatial surface integrals and upwind fluxes (due to causality) in the temporal surface integrals are used. These conditions on the space-time DG approximation are very similar to those imposed by Jiang and Shu for the scalar nonlinear hyperbolic case discussed for the nonlinear Burgers’ equation.
Unfortunately, the assumption 1) on exact integration is extremely difficult to guarantee and implement in practical simulations, which we describe in more detail in the next subsection. Without proper exact quadratures of the integral terms, the chain rule does not hold discretely. Such inexact approximations of functions are referred to as aliasing errors. These can occur in realistic simulations and might cause instabilities. Thus, robustness is still an issue and ad hoc dealiasing or stabilization mechanisms, e.g., artificial diffusion  are necessary.
2.3 The Unpleasant Role of Numerical Integration, Nonlinearities, Variational Crimes, and Aliasing
As described above it is necessary to discretize the variational formulation and make certain choices in the DG approximation such as the polynomial basis functions and, especially, the discrete integration to approximate the surface and volume integrals. Unfortunately, only in the rarest and simplest cases is it possible to avoid these discretisation steps and use an exact evaluation of the integrals. Hence, the notion of “variational crimes” is introduced to express the steps necessary to turn the formulation into an actual algorithm that can be implemented.
One of the biggest problems when discrediting nonlinear advection-diffusion problems is that in many interesting cases, the nonlinearity is non-polynomial. Our exemplary problem, the compressible Euler equations depend on the mass density u1, the momentum density in the x direction u2, and energy density u3. Often, these are also denoted as
This is important in the context of the DG discretisation because if the variables u1 and u2 are polynomials of degree N, the velocity v is not a polynomial, but a rational function. This occurs not only for the velocity but also for other quantities that are needed to evaluate the advective fluxes
If we focus for instance on the strong form DG volume integral from Eq. 38, we see that the core part to evolve the DG solution in time is an L2 projection of the flux divergence onto the polynomial basis
Unfortunately, these aliasing issues are not simply an abstract and “ugly” theoretical oddity without practical consequences. On the contrary, aliasing plays an important role when using DG methods for realistic complex applications to model nonlinear phenomena. It is worth pointing out that one of the advantages of DG, it’s very low dissipation errors, are in this particular point of view also it’s biggest problem. Due to the inherent low numerical dissipation in a high order DG method, there is no in-built self-defence against the aliasing issues and any instability that they may create. A repercussion of this fact is that it has become naturalized in the numerics community that especially the high order variants of the DG method, with very low dissipation errors, have robustness issues in practical applications. For instance DG approximations of the compressible Euler and Navier-Stokes equations are known to sometimes fail due to aliasing instabilities, e.g., . This instability can manifest itself through the observation that the kinetic energy artificially grows in the simulation, while the inner energy decreases. Note, that the total energy is conserved by construction with the DG method; however, this exchange of kinetic and internal energy is unphysical and violates the second law of thermodynamics and is purely a result of the variational crimes (inexact integration).
An obvious solution to these problematic variational crimes, nonlinearities, and aliasing is to decrease their deleterious effects as much as possible. While technically unavoidable in the strictest mathematical sense, it is possible to increase the amount of Legendre-Gauss quadrature nodes to evaluate all integrals “consistently” such that the inexactness errors are on the order of machine precision, see e.g., Kirby et al. . This approach is quite effective and immediately has a positive stabilizing effect on many applications with nonlinear PDEs, see e.g., [39, 96]. However, it is clear that the computational complexity drastically increases when arbitrarily increasing the number of quadrature nodes. Hence, one often tunes the increased number of Legendre-Gauss quadrature nodes and takes as many as needed to make a simulation stable—which is, of course, unsatisfying and ad hoc, as it highly depends on the particular problem setup. In comparison, another approach is to not directly fight the variational crimes themselves, but the consequences they induce. This is achieved by applying well designed filters to the solution, with the purpose to clip out higher order aliasing content in an effort to decrease its effect, e.g., . It is needless to say that this filtering approach is also very ad hoc and depends on many parameters that need tuning depending on the particular problem one wishes to solve with the high order DG scheme.
While these ad hoc stabilization techniques are reinforced by little to no mathematical analysis or rigor, the prevailing consensus in the DG community has been that, in practice, they work reasonably well. In that, consistent integration (sometimes referred to as over-integration) and filtering increase the robustness of high order DG methods to a level such that they could be applied to model challenging physical problems allowing them to shine with their high order accuracy and low dispersion and dissipation errors. Especially in the turbulence community, many research groups started to apply high order DG methods with stabilization in the context of implicit large eddy simulation with excellent results competitive with others from the broader numerics community, e.g., [97, 98]. However, in Moura et al.  it was reported that certain configurations of the DG method for the inviscid Taylor Green vortex problem kept crashing, even when drastically increasing the number of quadrature nodes in the surface and volume integrals. In fact, the amount of quadrature points was increased up to the point where the DG scheme was no longer computationally feasible, but the simulations still crashed. The inviscid Taylor Green vortex setup in this case was used to investigate the case of a very high Reynolds number flow with severe under resolution common in realistic turbulence setups. These findings were also verified by the authors of this review article and have a strong consequence for the DG community. While the ad hoc stabilization techniques were “good enough” in the sense that they helped to make the DG scheme run for a broad range of interesting problems, this approach is apparently not bullet proof. Further, it is impossible to tell a priori for which cases the stabilization will work and for which cases it will not.
This one example, where the high order DG scheme was not stable and could not finish the simulation illustrates, fundamentally, that the removal of aliasing and variational crimes cannot be reliably done in an ad hoc fashion. Instead we need a better understanding of these aliasing errors and how they can be removed from inexact and/or under resolved discretisations. Furthermore, we require a mathematically sound approach to address these aliasing errors in the DG approximation. That is, we need a novel strategy to design robust high order DG methods to approximate the solution of nonlinear advection-diffusion systems.
3 When Did the Novel Development Start?
In 2013, two landmark results completely reshaped the development of the DG method moving forward. First, in his PhD thesis, Fisher extended the work on entropy stable schemes LeFloch and Rhode  as well as the high order entropy stable schemes of LeFloch et al.  to the summation-by-parts (SBP) finite difference framework with high order boundary closures in Refs. 88 and 101 respectively. LeFloch et al. and Fisher et al. found that the entropy stability estimates for low-order FV methods, developed by Tadmor , can be extended to high order accuracy. Whereas the high order reconstruction of LeFloch et al. was for periodic domains (i.e. without considering finite domain boundaries), the SBP finite difference framework includes special boundary closures and are applicable for finite domains. Kreiss et al. [102–105] introduced the SBP finite difference framework to specifically mimic integration-by-parts. Integration-by-parts is a valuable tool for the construction of stability estimates. Further discussion on SBP is given by, e.g., Olsson [106, 107], Strand , Nordström  and Svärd and Nordström .
To briefly introduce the main ideas of the classic SBP finite difference framework, we consider a discretisation in one spatial dimension on a finite interval E = [−1,1]. Within this interval, we consider a set of N+1 regular grid nodes xj that include the boundaries x0 = 1 and xN = 1. On this grid, a continuous function u(x,t) is represented as the grid node values
The discrete integration and differentiation need to be compatible for a SBP operator to satisfy the property
Grouping terms and using that
which is a discrete approximation of the integration-by-parts formula for the corresponding functions u and v
hence the name summation-by-parts.
With the SBP property, it is directly possible to show L2-stability of the finite difference scheme for constant coefficient linear advection problems. Starting again as an example with the scalar problem Eq. 11 and the linear flux f = au, we get the following SBP finite difference semi-discretisation
As stated above, one motivation of the SBP framework is to mimic the energy analysis of finite element discretisations like that found in the DG analysis presented above for linear advection. We proceed and first get a corresponding Galerkin type variational form by multiplying with the discrete integration matrix
This form is valid for all grid functions VT multiplied from the left and hence is a direct approximation of the variational Galerkin form, e.g., Eq. 12. We point out that this finite difference approximation ignores the surface terms that are specific to the discontinuous Galerkin scheme. Here, the arbitrary grid function VT takes the role of the test function ϕ. Hence, we can mimic the next step in the analysis, i.e. replacing the test function with the DG solution
The volume term can be reformulated with the SBP property to move the discrete derivative onto the test function and generate boundary data
From this step we see that the contribution of the volume terms in the SBP finite difference scheme can be again lifted to the interval boundaries
but this time without any assumption on a necessary quadrature rule precision. In fact, this is general and holds for all diagonal norm SBP finite difference operators. Again, with the assumption of time continuity and periodic boundary conditions (UN=U0), it follows that the discrete L2-norm of the SBP finite difference solution
We emphasize again that neither the reconstruction techniques of LeFloch et al., nor the SBP finite difference framework as a whole, depend on integration or exact evaluation of integrals. Thus, in contrast to the DG stability results discussed above, the stability results obtained for SBP finite differences by Fisher et al. do not assume exact evaluation of any integrals. Thus, such methods yield efficient algorithms with feasible implementations that have provable stability estimates.
The second important result was separately discovered in 2013 by Gassner . He realized that the base operators of the nodal discontinuous Galerkin spectral element method (DGSEM) have the diagonal norm SBP property as long as the collocation nodes
In the nodal DGSEM-LGL framework, similar to the finite difference framework, the solution coefficients of the DG polynomial are nodal values
which have the Kronecker delta property that
where we used the short hand notation for the spatial derivative of the Lagrange basis
As was shown in Ref. 111, these two discrete operators are again compatible and provide the SBP property
which means that the DGSEM-LGL operators belong to the class of diagonal norm SBP operators. This simple property of one-dimensional discrete integration-by-parts is the basis for a whole polynomial spectral calculus  that includes, for instance, discrete version of Gauss’ law on curvilinear grids in three spatial dimensions.
Returning to the discussion on stability, the LGL quadrature rule with N+1 points has an integration precision of 2N−1. Thus, the DGSEM-LGL is stable for scalar linear advection as shown above. However, the DGSEM-LGL is not stable for nonlinear problems, e.g., for the quadratic flux function discussed in Section 2.1 where an integration precision of 3N−1 is necessary i.e., exact integration of the volume terms. But using the SBP property of the DGSEM-LGL operators, it is possible to apply ideas similar to Fisher et al. and construct a novel DGSEM with LGL quadrature, that is discretely L2-stable for the nonlinear Burgers’ equation, without the assumption on exact evaluation of the integrals . These first results have been extended and compounded upon for the compressible Euler equations [114–118], the shallow water equations [63, 119–121], the compressible Navier-Stokes equations [32, 122, 124], non-conservative multi-phase problems , magnetohydrodynamics [125, 126], relativistic Euler , relativistic magnetohydrodynamics , the Cahn-Hilliard equations , incompressible Navier-Stokes (INS) , and coupled Cahn-Hilliard and INS  among many other complex PDE models and DG discretization types e.g., .
4 What Is the Key Idea?
To recap the discussion, there are several obstacles that make it difficult to obtain entropy stability estimates for high order DG methods in the case of a general nonlinear system of hyperbolic PDEs: 1) The assumption of exact evaluation of integrals is unfeasible in practice; 2) We need to contract with entropy variables w that nonlinearly depend on the conservative quantities u, which means that we need to replace the test functions by a projection (or interpolation) of w(u) ; 3) We need to satisfy a discrete version of the chain rule to contract the flux divergence into entropy space, i.e. a discrete version of
In what follows, the key ideas to resolve all three issues are presented. We focus on issue 1) and consider a scalar nonlinear problem with quadratic fluxes discussed above first. Then we ramp-up the complexity in the second subsection and discuss how to extend the novel approach to general systems and how to resolve all issues (i)-(iii).
4.1 On the Conservative Form, Split Forms and Skew-Symmetry
To illustrate the general idea of a split form and how to incorporate it into a high order DG approximation we examine our simple scalar nonlinear hyperbolic conservation law, the Burgers’ equation. We start with the conservative form
that can be rewritten into its advective form
which is equivalent in the continuous case for smooth solutions. We can also consider an equivalent combination of the two forms
While in the continuous case with smooth solutions all of these forms are equivalent, it is important to note that in the discrete case this is not true. Considering the DGSEM operators with N+1 LGL nodes, the volume terms for the conservative form are given by
where U is the vector of values of u at the LGL nodes,
Only for polynomial functions u with degree
In all other cases, i.e. in the general case for arbitrary nodal vectors U, we get
The discrete forms are different because of different aliasing errors. Whereas the conservative form computes a discrete derivative of U2, the second form computes a “clean” derivative of U, but on the other hand needs to compute the product of two functions U and
We note that the idea of split formulations was already introduced in the spectral community to develop stable numerical methods for the incompressible Navier-Stokes equations e.g., , but is especially prominent in the finite difference fluid dynamics community e.g., [134–139]. Split formulations are used as a built-in dealiasing mechanism to stabilize numerical methods e.g., . Combinations of different forms of the advective terms of the compressible Euler equations yield finite difference approximations that are more robust than the standard conservative ones. In a perfect world, it would be desirable if we could choose the split form parameter α such that the different aliasing errors cancel exactly. Unfortunately, in general, it is not possible to cancel the aliasing errors for each grid node; however, what we will show next is that it is possible to cancel the aliasing errors in a way to get a global L2-stability estimate similar to the estimates by Jiang and Shu , but without the assumption of exact integration.
First, we derive the strong DG formulation of the split form of Burgers’ Eq. 58 by multiplying with a test function ϕ, integrating over a grid cell E, inserting a numerical flux function f* at the grid cell interface to account for the discontinuous nature of our ansatz, and use integration-by-parts again to arrive at
where the flux is f(U) = U2/2. We consider specifically the DGSEM-LGL variant to get discrete operators that satisfy the SBP property with diagonal norm (mass matrix)
Next, we follow the standard procedure to derive an L2-stability estimate by multiplying the scheme with the DG solution and use the quadrature rule to numerically integrate over the element i.e., we multiply by
where we used the fact that
the evolution of the discrete L2-norm in the grid cell E. Next, we focus on the volume terms. Note, that in the analysis of Jiang and Shu exact integration was assumed to contract the volume contribution to the surface. This is very important, as it allows direct control over the stability of the scheme with the choice of the numerical interface flux F*. Without the assumption of exact integration however, we look at the influence of the choice of the split form parameter α instead. We realize that the second term in the volume integral can be transposed
The term with the boundary evaluation matrix
This estimate is now analogous to the one with exact integration Eq. 19 and hence, with the same arguments, the choice of the numerical flux functions as
gives again a discrete stability estimate
for the split form DGSEM-LGL with
We note that this particular choice of numerical flux function exactly preserves the discrete L2-norm. If one considers non-smooth solutions this choice would be inappropriate as for e.g., shocks, because the L2-norm needs to decrease as u2 is a mathematical entropy for Burgers’ equation. For scalar equations, s(u) = u2/2 is the square entropy (which leads to an L2-stability estimate e.g., [99, 141]) that gives the simple entropy variables
As an example, a simple choice of an entropy dissipative numerical flux function is that of Rusanov
An important question is how we can extend this approach to general nonlinear systems. Following the same ideas, it was possible to derive split forms for the shallow water equations [63, 119, 142], a simplified version of the compressible Euler equations. There are many split forms for compressible Euler that, for instance, give kinetic energy preserving properties e.g., [143, 144]. However, up to now, no split form for the compressible Euler equations is known that gives the desired discrete entropy stability estimate. The problems are issues ii) and iii) mentioned in the beginning of the section, where we need the discrete chain rule property to contract the volume terms to the surface. Concluding this subsection we revisit the derivations of Burgers’ equation and make two important observations.
First, with the proper choice of
Skew-symmetry is strongly connected to entropy, see e.g., Tadmor . Multiplying the spatial derivative term by u as in the L2-stability analysis gives
which shows that the skew-symmetric form gives a product-rule type form in the stability analysis that can directly be contracted to the divergence form i.e., contracts to the surface when integrating. In fact, for this simple problem, the chain rule needed for contraction reduces to the simpler product rule. Analogously, we get for the discrete skew-symmetric volume terms of the DGSEM-LGL
Thus, in our derivation, we already used a specific discrete version of the chain rule (product rule) to get our estimate. The question is, how to extend this idea to the general case?
The second important observation pioneered for SBP schemes by Fisher in his PhD thesis  (in the spirit of earlier work by LeFloch et al. ) is that the particular skew-symmetric volume terms
In combination with the first observation we get the property that this new discrete derivative operator (or divergence operator in the multi-dimensional case) satisfies
We will see in the next subsection, how these observations guide the path to discrete stability estimates for general hyperbolic PDE systems.
4.2 On the Discrete Entropy Stability of the DGSEM-LGL
We have demonstrated how to build a high order skew-symmetric DG approximation of the scalar nonlinear Burgers’ equation. To do so required a very particular discrete derivative operator Eq. 78 that was the key to restore discrete entropy stability. We now discuss how to extend the split form approach to general systems of nonlinear hyperbolic conservation laws. For general nonlinear systems, it is unclear how to explicitly construct the split form to obtain a discrete chain rule property. In particular, the compatibility condition on the physical fluxes obtained when one contracts into entropy space Eq. 29 that we reproduce here, assuming one spatial dimension, due to their pertinence in the present discussion
As previously indicated, the chain rule is either unfeasible or even impossible to directly recover with discrete differentiation. With this in mind we apply the product rule to this compatibility condition between the physical flux f and the entropy flux fs to find
A principle motivation for this manipulation is because it is far easier to recover the product rule discretely than it is the chain rule. That is, we already have a particular discrete equivalent for the product rule if the discrete derivative matrix
Next, we aim to find a discrete version of the new compatibility condition Eq. 81 following the ideas of Tadmor . Tadmor analyzed low order FV schemes and developed conditions on the numerical surface flux to derive a discretely entropy conserving scheme. In the context of the low order FV methodology, our unknowns in the elements are mean values that are naturally discontinuous across grid cell interfaces. As mentioned above, the idea to resolve these discontinuities with numerical flux functions was also used in the construction of the DG approximation. Consider the contribution to the compatibility condition Eq. 81 at an arbitrary interface. It depends on the discrete values in the current cell and the direct neighbor of that cell, denoted again with a “+”. We approximate all derivatives with first order differences and define Tadmor’s entropy conservation condition on the numerical surface flux function
which matches the particular entropy conservative flux derived in Section 4.1. We note again that the entropy conservative flux is symmetric in its arguments U+ and U, and is consistent to the physical flux in the sense that for the same arguments we recover the PDE flux
However, for systems of nonlinear hyperbolic conservation laws Eq. 83 is a single algebraic condition for a system vector of unknown flux quantities. Therefore, care must be taken to define an entropy conservative numerical flux function that remains physically consistent. That being said, the entropy conservation condition on the numerical surface flux Eq. 83 is an incredibly powerful statement. Provided we know an explicit form of the entropy variables, the physical flux, and the entropy flux we can define an appropriate numerical flux that ensures entropy consistency for a low order FV numerical approximation. A general form for such a numerical flux was developed by Tadmor  defined as a phase integral
To evaluate the phase integral form of the numerical flux function requires a certain quadrature rule that defines a path through phase space. Though theoretically useful, this phase integral form is computationally prohibitive for practical simulations, even for low order numerical approximations. However, over the past 20 years “affordable” versions of the entropy conservative flux function
The key to these numerically tractable versions of the entropy conservative numerical flux function is to evaluate the components of the physical flux at various mean states between U+ and U. Note that these mean states can take on incredibly complex forms that depend on the arithmetic mean, the product of arithmetic means, or more uncommon quantities like the logarithmic mean. Complete details on the derivation of such numerical flux functions can be found in e.g., [81, 83, 146, 147].
For illustrative purposes we summarize the specific form of an entropy conservative numerical flux for the one dimensional compressible Euler equations with the ideal gas assumption due to Chandrashekar . First, we introduce notation for the arithmetic mean and the logarithmic mean for two quantities a and a+:
Also, we introduce a variable proportional to the inverse of the temperature
which simplifies the form of the entropy variables Eq. 36 to be
Then, applying the entropy conservation condition Eq. 83 and many algebraic manipulations later, we arrive at an analytical expression of an entropy conservative numerical flux for the compressible Euler equations
So far, the discussion on entropy conservative numerical approximations has all been in the context of low order finite volume methods. It is possible to create a high order entropy aware scheme with ENO  or WENO type reconstructions . However, as mentioned above, a critical and remarkable result of Fisher’s work is that a low order finite volume entropy conservative scheme can be extended to an arbitrarily high order accurate spatial scheme, when it is based on diagonal norm SBP operators . As was described for Burgers’ equation in the last two observations in Section 4.1, the crucial part is to move the contribution to the entropy production from the volume terms to the surface via a discrete version of the chain rule. Once stability is only governed by surface contributions, it is possible to control the stability with a proper choice of numerical interface flux function. In this sense, the analysis of the high order scheme reduces to a similar problem as for the low order finite volume method Eq. 83 and the well-established theoretical analysis tools and results can be reused. It is worth mentioning that
We return to a strong form DG approximation of a nonlinear hyperbolic system of conservation laws that takes a form identical to Eq. 38 but now only considered in one spatial dimension
Here the standard tools of a nodal DG approximation have been applied:
1. The solution and physical fluxes are approximated with polynomials.
2. Any integrals in the variational formulation are approximated with high order LGL quadrature.
3. The interpolation and quadrature nodes are collocated.
Importantly, these steps mean that the discrete DG differentiation matrix is a SBP operator. Furthermore, we use the entropy conservative numerical surface flux at the interface and the particular discrete derivative projection
Now, if we take the test function
We note that this remarkable property holds for general nonlinear systems with available entropy estimate and a corresponding low order entropy conserving flux
which is the discrete analogue of the integral form of the entropy conservation law discussed in Section 2.2. The resulting DGSEM-LGL is entropy conservative by construction and it is important to note we have assumed no exactness on the integration. From this baseline entropy conservative numerical scheme, that does not dissipate entropy by construction, we can create a high order DGSEM-LGL that enforces the entropy inequality Eq. 32. We do so by introducing dissipation at the element interfaces via the choice of the numerical surface flux function, e.g., the Rusanov flux Eq. 72. More complex dissipation techniques are also available that dissipate solution information according to the different wave strengths with complete details found in, e.g., [152–156]. We finally note that this discussion was restricted to one spatial dimension for the sake of convenience and simplicity. Extensions to general three dimensional curvilinear coordinate systems are available, see e.g., [32, 101, 116, 125, 157, 158] for details.
4.3 Validation of Robustness and Application to Space Physics of the Entropy Stable DGSEM
In this subsection, we demonstrate two exemplary simulation results of a DG scheme based on the key ideas outlined above. The general split form DGSEM with LGL nodes on three dimensional curvilinear hexahedral unstructured meshes is implemented in the open source software FLUXO (project-fluxo/fluxo at github), written in modern Fortran with a special emphasis on massively parallel CPU based hardware. The main focus of the software is on compressible Navier-Stokes and visco-resistive MHD equations. Time integration of the semi-discrete form is done with a fourth order accurate low storage Runge-Kutta method of Carpenter and Kennedy .
For the validation of the robustness we revisit an important numerical contribution of Moura et al. . They were the first to report of a test case that the DG scheme with (numerical) exact integration was not able to run, demonstrating, that further improvement on the robustness of the DG methodology was necessary. For this validation test, we consider the compressible Navier-Stokes equations (viscous case) or the compressible Euler equations (inviscid case, basically setting the viscosity parameter to zero). The considered problem is the Taylor-Green vortex in a fully periodic domain, which serves as a test case for a fully periodic turbulent box [0,1]3, that starts with a smooth initial velocity field
and transitions to turbulence during its temporal evolution until it reaches a state similar to homogeneous turbulence. The initial density is uniform
For the robustness investigation, we consider an inviscid flow (with the viscosity parameter is zero) and focus on three particular setups, where N = 1,3,7 with number of elements NQ = 56,28,14 respectively. This ensures that for all three computations the overall number of DOF per conserved quantity is equal, about 1.4 million. There are many more investigations of different configurations presented in Ref. 160, but they all demonstrate the same behavior: while the low order variants N = 1,3 seem to be relatively robust with full integration, the higher the polynomial degree, the less stable the DG method becomes. And for the case N = 7 with NQ = 14 the simulation crashed at about a simulation time of tcrash = 8.4, even when increasing the quadrature nodes from 83 up to 323 = 32768 per element. In contrast, the novel entropy stable DGSEM with standard LGL nodes runs all configurations without crashing.
Furthermore, it is possible to run this challenging test case even without any artificial dissipation, i.e. with the
FIGURE 1. The estimate of the numerical Reynolds number of an under resolved simulation with polynomial degree N = 7 and eight elements (643 DOF). The physical Reynolds number of the Taylor-Green vortex setup is Re = 1600. The entropy conservative (EC) scheme retains the physical Reynolds number remarkably well and is virtually dissipation free. The entropy stable (Rusanov) variant clearly introduces stabilizing dissipation as soon as the scales can no longer be resolved starting at about t = 3.
We note, that this test case would crash for the standard DG scheme, however for the presented novel DGSEM-LGL with the proper discrete chain rule, it runs with the dissipative Rusanov numerical flux
For an exemplary application we consider a complex test case from space physics. We focus on the electrodynamic and plasma interaction of the moon Io with the strong magnetic field of Jupiter. Io is embedded in a dense plasma torus, induced by the magnetosphere of Jupiter, and it exhibits interesting plasma flow characteristics containing steep gradients and discontinuities . The general problem setup is illustrated in Figure 2. Neglecting neutral density, relativistic, viscous, resistive, and Hall effects, this MHD flow within Io's plasma torus can be modeled with the ideal MHD equations. For such magnetohydrodynamic flows, the entropy stable DGSEM solver in FLUXO uses a hyperbolic divergence cleaning mechanism to enforce the divergence-free constraint on the magnetic field variables
FIGURE 2. Io’s electrodynamic interaction is unique due to its fast rotation and the influence of the strong magnetic field of Jupiter. Plasma interactions with Io’s atmosphere lead to mass loss in the form of ions and neutrons. These neutrons then ionize through radiative effects and accumulate around Io forming a plasma torus. Consequently, this flow of magnetized plasma past the obstacle Io, combined with atmospheric interactions, are the engine behind Io’s plasma interactions with Jupiter’s magnetosphere.
As Io orbits Jupiter, plasma from the torus streams in and forms a local ionosphere that induces polarisation charges and modifies the electric field, thus changing the local Lorentz force and damping electron and ion flow close to Io. The plasma flow is strongly reduced inside Io’s (very weak) atmosphere, which can be modeled by incorporating a neutral collision source term to the ideal MHD system, Saur et al. [162, 164],
with the collision frequency
FIGURE 3. Standing Alfvén current tube magneto-spheric disturbances, termed Alfvén wings, that spread and propagate away along field lines in both directions. The development of the Alfvén wings are observable in the regions with decreasing plasma bulk velocity v1 and perturbed magnetic field B1 north and south of Io in the xz plane. It is known that the Alfvén current tubes are bent back by a constant angle with respect to the unperturbed background magnetic field. Further, the perturbation of the magnetic field is positively correlated with the perturbation of the velocity field in the northern Alfvén wing and negatively correlated in the southern wing .
In the dimensionless computational domain we scale the radius of the atmosphere to take a spherical shape with radius one and locate the center of the sphere at the origin
The ionospheric processes uses the exponential blending of the collision frequency above with a dilatation factor d = 150/1820. The initial conditions for the flow are taken as
that will evolve to a final time T = 5. The gas constant is taken to be
In anticipation to capture the relevant physical interactions at the sphere as well as the development of the Alfvén wings best, we exploit the geometric flexibility of the entropy stable DGSEM solver and divide the computational domain into an unstructured, curvilinear mesh presented in Figure 4. Within each element we use a polynomial order of N = 3.
FIGURE 4. Computational mesh for the Alvén wing test problem built from 14,336 curved elements. At the origin the elements are curved to capture the spherical shape of the neutral gas cloud
In Figure 5 we show a 2D-slice of the B1 and v1 components at y = 0 of the entropy stable approximation at the final time which presents the numerically generated Alfvén wings from the entropy stable DGSEM. It also demonstrates the expected positive correlation of the B1 variable with the velocity variable v1 in the northern Alfvén wing and the negative correlation in the southern wing. Moreover, we consider profile slices in these B1 and v1 along the line z = 5 and compare the results to a solution computed by the open source software ZEUS (www.astro.princeton.edu/jstone/zeus.html) in Figure 6. ZEUS is a FV solver written in spherical coordinates that uses explicit time integration. For the presented comparison, ZEUS used 10 million grid cells (DOF) in total whereas the entropy stable DGSEM solver used approximately 14,336 elements with N=3 polynomials and a total of about 1 million DOF. The reduction of DOFs also translated in a nearly 10 fold reduction of overall CPU time when both codes were run in parallel using MPI on 100 cores. This increased efficiency of the high order DGSEM-LGL, the increased robustness due to discrete entropy stability, and the geometrical flexibility are several advantages of this novel DG framework.
FIGURE 5. Alfvén wings numerically computed with the entropy stable DGSEM for the plasma interaction of a spherical gas cloud. The snapshot is a slice in the xz plane at y = 0 at the final time T = 5. The polynomial order in each spatial direction was N = 3 in each spatial direction. As expected, the Alfvén wings evolve from the northern and southern poles of the neutral gas cloud and are bent back by a constant angle with respect to the background magnetic field. This bending was taken into account in the construction of the curvilinear mesh.
FIGURE 6. One dimensional visualization of the Alfvén wing solution along the line z = 5 from the xz plane slice at y = 0. A comparison is performed between the entropy DG solver with N = 3 in each spatial direction and the first order finite volume solver ZEUS. The entropy stable DG approximation uses 90% fewer DOF compared to the 10 million DOF used for the ZEUS computation. Qualitatively, the solutions are very similar.
5 Where to Go Next?
The response of the DG community to the split form DGSEM with LGL quadrature on tensor-product hexahedra has been astounding. However, naturally, there are still many limitations of this method. Some that have been recently addressed and many others still open. So far, we discussed semi-discrete DGSEM-LGL variants with tensor product expansions on possible curvilinear unstructured hexahedral meshes. Direct extensions of this variant include non-conforming meshes [166, 167], moving meshes [168–170], different related versions such as e.g., the line DG method , and a fully discrete space-time approach without the assumption on time continuity [172–175]. An exciting recent development are explicit modified Runge-Kutta methods that retain the semi-discrete entropy stability estimates .
A downside of LGL is that in comparison to the Legendre-Gauss (LG) points, the accuracy in dispersion and dissipation is lower, see e.g., . Unfortunately, LG points do not include the boundary nodes, hence they do not directly satisfy the classic SBP property and the presented developments cannot be directly applied to this case. However, there are several developments where the framework was extended to construct entropy stable variants with LG nodes. One is based on a staggered grid approach by Parsani et al. . The authors used the standard LG nodes to span the solution, however to compute the discrete derivative operator, they interpolate onto a higher order staggered LGL grid where they can use the SBP property. Then they ensure that the back-projection is still entropy stable to retain the stability estimate on top of the accuracy of the LG nodes. Another approach is presented by e.g., Ortleb  where the author directly constructed a scheme that preserves the kinetic energy with LG nodes. Although SBP could not directly be applied, the difference lies only in the boundary operator
A generalization onto multi-dimensional domains is given in e.g., , termed multidimensional SBP operators. This is an interesting development, as it resolves another limitation of the classic DGSEM-LGL. In some applications, it is favourable to have more flexibility when generating meshes for geometries with complex shapes. In this case, meshes with simplex element types such as triangles and tetrahedra or even hybrid element types such as prisms and pyramids are desired. We will mention some recent developments and extensions here, but stress that this list is not complete and there are many more. Returning to the multidimensional SBP framework, this approach’s strength is that it can be used to construct stable methods on simplex meshes e.g., by Chan [115, 180], Chen and Shu , Hicken et al. , and Crean et al. . Another interesting approach to generate DG scheme on general meshes is presented by Chan . He shows that a special projection directly with the entropy variables collocated at the grid nodes can give a SBP type property that can be used to construct stable discretisations.
As stated the response of the DG community has been astonishing with an explosion of developments, extensions, and new insights. However, there are still many unresolved issues that need to be researched in the future to further evolve high order DG methods into a viable tool for computational physics. The most important problem is still robustness. Although entropy stability significantly improves the stability of the scheme in many applications, there are still situations where no significant gain in robustness of the DG scheme can be observed [183, 184]. A possible reason could be, for instance, the incorrect choice of mathematical entropy function. While there is typically only one physical (thermodynamic) entropy, there are many mathematical ones that often lead to a stability estimate in corresponding norms of the solution. So what are the important entropy quantities to consider? What about e.g., kinetic energy, cross helicity and related quantities? For example, in turbulence, the kinetic energy and proper prediction of its behavior seems to play an important role [144, 161, 185, 186]. Besides getting discrete entropy stability estimates, it is possible to use the split form DGSEM-LGL approach from Section 4.1 to construct DG schemes that discretely preserve the kinetic energy and are discretely compatible for the inner and the total energy, e.g., . It could be shown in Refs. 160 and 161 that such kinetic energy preserving DG schemes behave favourably for the simulation of compressible turbulence in comparison to the standard DG, especially in combination with subgrid turbulence models. However, just as in the development of entropy conservative (and in turn entropy stable) numerical flux functions there is nonuniqeuness. There exist many possible solutions to create a numerical flux that is e.g., kinetic energy preserving or entropy conservative or both e.g., [81, 83, 144, 185, 187]. Moreover, the discrete behavior of the kinetic energy as it evolves in time for under-resolved turbulent flows is quite different even between fluxes that are all provably kinetic energy preserving on paper [116, 187, 188]. So, an important question for the future is thus: What are the important quantities not only from a mathematical, but also from a physical point of view?
A fundamental issue in this context are problem setups that involve physical discontinuities e.g., shock waves in the compressible Euler or ideal MHD equations. Discontinuities trigger another instability, inherent in high order methods: the Gibbs phenomenon i.e., numerical oscillations. These oscillations can be devastating as under- or overshoots can cause non-physical state solutions e.g., negative density or pressure. Hence, positivity is a necessary criterion for all numerical methods when simulating such problems. However, up to this date there is still not enough research into the topic of entropy stability and positivity e.g., . It is worth pointing out, that mathematically, the entropy function is only well-defined for positive solutions and, hence, is strongly connected to positivity. Generalizing this discussion, it is evident that entropy stability is “not enough” as a property for the numerical method. We need more properties, such as e.g., positivity. However, this is also where we reach uncharted research territory as even for many continuous problems e.g., the compressible Navier-Stokes equations it is up to this point unclear to show positivity even for the model itself.
The overview in this work focused on the volume contributions and the underlying tools (physical and mathematical) which led to the entropy stable DG method. However, the contributions at the physical boundaries have been ignored. Properly posing the boundary conditions to be entropy stable for a given model, like the compressible Euler or Navier-Stokes equations, has been considered [189–195], but this is remains an active area of research particularly because the treatment and behavior of the solution (whether on the continuous or discrete level) is directly related to the validity of a mathematical PDE model and directly tied into issues of well-posedness.
Concluding, we are currently at an exciting development stage with high order DG methods, where we can mimic important continuous stability estimates by careful construction of discrete operators. However, practical simulations show that these are not enough for the most complex problems that we desire to simulate. Plus, our numerical schemes and their properties are very close to the current analytical knowledge we have about the physical models. It is very hard to progress with the numerical developments further than what is analytically known: Which properties are important? How do you show positivity? Or a more general question: How does one prove physicality of the solutions for a given PDE model? It seems that the answers can only be given in close collaboration of researchers from physics and mathematics.
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.
GG was supported by the European Research Council (ERC) under the European Union’s Eights Framework Program Horizon 2020 with the research project Extreme, ERC Grant Agreement No. 714487. Support for open access publication was provided by Linköping University.
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.
This work was partially performed on the Cologne High Efficiency Operating Platform for Sciences (CHEOPS) at the Regionales Rechenzentrum Köln (RRZK) at the University of Cologne. The authors thank Marvin Bohm for his contributions to the space physics discussion.
2. Nitsche JA. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. Abh Math Sem Univ Hamburg (1971) 36:9–15. doi:10.1007/BF02995904
3. Cockburn B, Hou S, Shu C-W. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV: the multidimensional case. Math Comput (1990) 54:545–81. doi:10.2307/2008501
4. Cockburn B, Shu CW. The Runge-Kutta local projection -discontinuous Galerkin method for scalar conservation laws. Rairo-Matehmatical Model Num Anal-Modelisation Mathematique et Anale Num (1991) 25:337–61. doi:10.1051/m2an/1991250303371
7. Bassi F, Rebay S. A high order accurate discontinuous finite element method for the numerical solution of the compressible Navier-stokes equations. J Comput Phys (1997) 131:267–79. doi:10.1006/jcph.1996.5572Get
14. Deng S, Cai W, Astratov V. Numerical study of light propagation via whispering gallery modes in microcylinder coupled resonator optical waveguides. Optic Express (2004) 12:6468–80. doi:10.1364/opex.12.006468
15. Kopriva DA, Woodruff SL, Hussaini MY. Discontinuous spectral element approximation of Maxwell’s Equations In: B Cockburn, G Karniadakis, and C-W Shu, editors. Proceedings of the international symposium on discontinuous Galerkin methods. New York: Springer-Verlag (2000) p. 355–61.
16. Kopriva DA, Woodruff SL, Hussaini MY. Computation of electromagnetic scattering with a non-conforming discontinuous spectral element method. Int J Numer Methods Eng (2002) 53:105–22. doi:10.1002/nme.394
21. Wilcox LC, Stadler G, Burstedde C, Ghattas O. A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media. J Comput Phys (2010) 229:9373–96. doi:10.1016/j.jcp.2010.09.008
22. Bonev B, Hesthaven JS, Giraldo FX, Kopera MA. Discontinuous Galerkin scheme for the spherical shallow water equations with applications to tsunami modeling and prediction. J Comput Phys (2018) 362:425–48. doi:10.1016/j.jcp.2018.02.008
24. Giraldo F, Restelli M. A study of spectral element and discontinuous Galerkin methods for the Navier–Stokes equations in nonhydrostatic mesoscale atmospheric modeling: equation sets and test cases. J Comput Phys (2008) 227:3849–77. doi:10.1016/j.jcp.2007.12.009
25. Restelli M, Giraldo F. A conservative discontinuous Galerkin semi-implicit formulation for the Navier-Stokes equations in nonhydrostatic mesoscale modeling. SIAM J Sci Comput (2009) 31:2231–57. doi:10.1137/070708470
26. Fagherazzi S, Furbish D, Rasetarinera P, Hussaini MY. Application of the discontinuous spectral Galerkin method to groundwater flow. Adv Water Res (2004a) 27:129–40. doi:10.1016/j.advwatres.2003.11.001
28. Cockburn B, Karniadakis G, Shu C-W. The development of discontinuous Galerkin methods. In: B Cockburn, G Karniadakis, and C-W Shu, editors. Proceedings of the inter)national symposium on discontinuous Galerkin methods. New York: Springer-Verlag (2000). p. 3–50.
38. Bohm M, Schermeng S, Winters AR, Gassner GJ, Jacobs GB. Multi-element SIAC filter for shock capturing applied to high-order discontinuous Galerkin spectral element methods. J Sci Comput (2019) 81:820–44.
40. Hennemann S, Rueda-Ramírez AM, Hindenlang FJ, Gassner GJ (2020). A provably entropy stable subcell shock capturing approach for high order split form DG for the compressible Euler equations. J Comput Phys 109935. doi:10.1016/j.jcp.2020.109935
43. Dumbser M, Loubère R. A simple robust and accurate a posteriori sub-cell finite volume limiter for the discontinuous Galerkin method on unstructured meshes. J Comput Phys (2016) 319:163–99. doi:10.1016/j.jcp.2016.05.002
44. Dumbser M, Zanotti O, Loubère R, Diot S. A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws. J Comput Phys (2014) 278:47–75. doi:10.1016/j.jcp.2014.08.009
45. Fambri F. Discontinuous Galerkin methods for compressible and incompressible flows on space–time adaptive meshes: toward a novel family of efficient numerical methods for fluid dynamics. Arch Comput Methods Eng (2020) 27:199–283. doi:10.1007/S11831-018-09308-6
46. Giri P, Qiu J. A high-order Runge-Kutta discontinuous Galerkin method with a subcell limiter on adaptive unstructured grids for two-dimensional compressible inviscid flows. Int J Numer Methods Fluid (2019) 91:367–94. doi:10.1002/fld.4757
47. Zanotti O, Fambri F, Dumbser M, Hidalgo A. Space–time adaptive ADER discontinuous Galerkin finite element schemes with a posteriori sub-cell finite volume limiting. Comput. Fluids (2015) 118:204–24. doi:10.1016/j.compfluid.2015.06.020
48. Ching EJ, Lv Y, Gnoffo P, Barnhardt M, Ihme M. Shock capturing for discontinuous Galerkin methods with application to predicting heat transfer in hypersonic flows. J Comput Phys (2019) 376:54–75. doi:10.1016/j.jcp.2018.09.016
53. Gassner G, Haas M. An explicit high order accurate predictor-corrector time integration method with consistent local time-stepping for discontinuous Galerkin schemes. AIP Conf Proceed (2009) 1168:1188–91. doi:10.1063/1.3241276
54. Kopera MA, Giraldo FX. Analysis of adaptive mesh refinement for IMEX discontinuous Galerkin solutions of the compressible Euler equations with application to atmospheric simulations. J Comput Phys (2014) 275:92–117. doi:10.1016/j.jcp.2014.06.026
55. Gassner G, Kopriva DA. A comparison of the dispersion and dissipation errors of Gauss and Gauss-Lobatto discontinuous Galerkin spectral element methods. SIAM J Sci Comput (2010) 33:2560–79. doi:10.1137/100807211
57. Altmann C, Beck AD, Hindenlang F, Staudenmaier M, Gassner GJ, Munz C-D. An efficient high performance parallelization of a discontinuous Galerkin spectral element method. In: R Keller, D Kramer, and J-P Weiss, editors Facing the multicore-challenge IIIof lecture notes in computer science. Springer Berlin Heidelberg (2013), Vol. 7686, p. 37–47.
58. Hindenlang F, Gassner GJ, Altmann C, Beck A, Staudenmaier M, Munz C-D. Explicit discontinuous Galerkin methods for unsteady problems. Comput Fluids (2012) 61:86–93. doi:10.1016/j.compfluid.2012.03.006
59. Dumbser M, Käser M, Toro EF. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes–V. local time stepping and p-adaptivity. Geophys J Int (2007) 171:695–717. doi:10.1111/j.1365-246X.2007.03427.x
62. Warburton T. Application of the discontinuous Galerkin method to Maxwell’s equations using unstructured polymorphic -finite elements. In: B Cockburn, G Karniadakis, and C-W Shu, editors. Proceedings of the international symposium on discontinuous Galerkin methods. New York: Springer-Verlag (2000). p. 451–8.
63. Wintermeyer N, Winters AR, Gassner GJ, Kopriva DA. An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry. J Comput Phys (2017) 340:200–42. doi:10.1016/j.jcp.2017.03.036
64. Kopriva DA, Gassner GJ. An energy stable discontinuous Galerkin spectral element discretization for variable coefficient advection problems. SIAM J Sci Comput (2014) 34:A2076–A2099. doi:10.1137/130928650
70. Klingenberg C, Markfelder S. Non-uniqueness of entropy-conserving solutions to the ideal compressible MHD equations. In: A Bressan, M Lewicka, D Wang, and Y Zheng, editors. Hyperbolic problems: theory, numerics, applications. AIMS on Applied Mathematics (2020), Vol. 10, p. 491–8.
75. Lax PD. Weak solutions of nonlinear hyperbolic conservation equations and their numerical computation. Commun Pure Appl Math (1954) 7:159–93. 10.1002/cpa.3160070112
79. Ranocha H. Mimetic properties of difference operators: product and chain rules as for functions of bounded variation and entropy stability of second derivatives. BIT Num Math (2019) 59:547–63. doi:10.1007/s10543-018-0736-7
81. Chandrashekar P. Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier-Stokes equations. Commun Comput Phys (2013) 14:1252–86. doi:10.4208/cicp.170712.010313a
91. Hughes TJ, Franca L, Mallet M. A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics. Comput Methods Appl Mech Eng (1986) 54:223–34. doi:10.1016/0045-7825(86)90127-1
93. Hiltebrand A, Mishra S. Entropy stability and well-balancedness of space-time DG for the shallow water equations with bottom topography. Netw Heterogeneous Media (2016) 11:145–62. doi:10.3934/nhm.2016.11.145
96. Mengaldo G, Grazia DD, Moxey D, Vincent PE, Sherwin SJ. Dealiasing techniques for high-order spectral element methods on regular and irregular grids. J Comput Phys (2015) 299:56–81. doi:10.1016/j.jcp.2015.06.032
97. Beck AD, Bolemann T, Flad D, Frank H, Gassner GJ, Hindenlang F, et al. High-order discontinuous Galerkin spectral element methods for transitional and turbulent flow simulations. Int J Numer Methods Fluid (2014) 76:522–48. doi:10.1002/fld.3943
98. Persson P-O, Bonet J, Peraire J. Discontinuous Galerkin solution of the Navier-Stokes equations on deformable domains. Comput Methods Appl Mech Eng (2009) 198:1585–95. doi:10.1016/j.cma.2009.01.012
104. Kreiss H-O, Scherer G. Finite element and finite difference methods for hyperbolic partial differential equations. Math Aspect Finite Elem Part Diff Equat (Elsevier) (1974) 195–212. doi:10.1016/B978-0-12-208350-1.50012-1
105. Kreiss H-O, Scherer G. On the existence of energy estimates for difference approximations for hyperbolic systems (1977). Tech. Rep. Uppsala:Department of Scientific Computing. Uppsala University.
111. Gassner G. A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods. SIAM J Sci Comput (2013) 35:A1233–A1253. doi:10.1137/120890144
113. Kopriva DA. A polynomial spectral calculus for analysis of DG spectral element methods. In Bittencourt ML, Dumont NA, and Hesthaven, JS, editors. Spectral and high order methods for partial differential equations ICOSAHOM 2016. Cham:Springer International Publishing (2017) p. 21–40.
114. Carpenter M, Fisher T, Nielsen E, Frankel S. Entropy stable spectral collocation schemes for the Navier-Stokes equations: discontinuous interfaces. SIAM J Sci Comput (2014) 36:B835–B867. doi:10.1137/130932193
116. Gassner G, Winters A, Kopriva DA. Split form nodal discontinuous Galerkin schemes with summation-by-parts property for the compressible Euler equations. J Comput Phys (2016a) 327:39–66. doi:10.1016/j.jcp.2016.09.013
118. Gouasmi A, Duraisamy K, Murman SM. Formulation of entropy-stable schemes for the multicomponent compressible Euler equations. Comput Methods Appl Mech Eng (2020) 363:112912. doi:10.1016/j.cma.2020.112912
119. Gassner GJ, Winters AR, Kopriva DA. A well balanced and entropy conservative discontinuous Galerkin spectral element method for the shallow water equations. Appl Math Comput (2016b) 272(2):291–308. doi:10.1016/j.amc.2015.07.014
122. Parsani M, Carpenter M, Nielsen E. Entropy stable discontinuous interfaces coupling for the three-dimensional compressible Navier-Stokes equations. J Comput Phys (2015a) 290:132–8. doi:10.1016/j.jcp.2015.02.042
123. Parsani M, Carpenter M, Nielsen E. Entropy stable wall boundary conditions for the three-dimensional compressible Navier-Stokes equations. J Comput Phys (2015b) 292:88–113. doi:10.1016/j.jcp.2015.03.026
125. Bohm M, Winters AR, Gassner GJ, Derigs D, Hindenlang F, Saur J. An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. Part I: theory and numerical verification. J Comput Phys (2018) 422:108076. doi:10.1016/j.jcp.2018.06.027
129. Manzanero J, Rubio G, Kopriva DA, Ferrer E, Valero E. A free–energy stable nodal discontinuous Galerkin approximation with summation–by–parts property for the Cahn–Hilliard equation. J Comput Phys (2020c) 403:109072. doi:10.1016/j.jcp.2019.109072
130. Manzanero J, Rubio G, Kopriva DA, Ferrer E, Valero E. An entropy–stable discontinuous Galerkin approximation for the incompressible Navier–Stokes equations with variable density and artificial compressibility. J Comput Phys (2020a) 408:109241. doi:10.1016/j.jcp.2020.109241
131. Manzanero J, Rubio G, Kopriva DA, Ferrer E, Valero E. Entropy–stable discontinuous Galerkin approximation with summation–by–parts property for the incompressible Navier–Stokes/Cahn–Hilliard system. J Comput Phys (2020b) 408:109363. doi:10.1016/j.jcp.2020.109363
132. Chen T, Shu C-W. Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws. J Comput Phys (2017) 345:427–61. doi:10.1016/j.jcp.2017.05.025
134. Ducros F, Laporte F, Soulères T, Guinot V, Moinat P, Caruelle B. High-order fluxes for conservative skew-symmetric-like schemes in structured meshes: application to compressible flows. J Comput Phys (2000) 161:114–39. doi:10.1006/jcph.2000.6492
136. Kennedy CA, Gruber A. Reduced aliasing formulations of the convective terms within the Navier–Stokes equations for a compressible fluid. J Comput Phys (2008) 227:1676–700. doi:10.1016/j.jcp.2007.09.020
141. LeFloch PG, Mohammadian M. Why many theories of shock waves are necessary: kinetic functions, equivalent equations, and fourth-order models. J Comput Phys (2008) 227:4162–89. doi:10.1016/j.jcp.2007.12.026
142. Winters AR, Gassner GJ. A comparison of two entropy stable discontinuous Galerkin spectral element approximations for the shallow water equations with non-constant topography. J Comput Phys (2015) 301:357–76. doi:10.1016/j.jcp.2015.08.034
145. Tadmor E. Perfect derivatives, conservative differences and entropy stable computation of hyperbolic conservation laws. Disc Contin Dyn Syst-A (2016) 36:4579–98. 10.3934/dcds.2016.36.4579
146. Fjordholm US, Mishra S, Tadmor E. Well-balanced and energy stable schemes for the shallow water equations with discontiuous topography. J Comput Phys (2011) 230:5587–609. doi:10.1016/j.jcp.2011.03.042
148. Fjordholm US, Mishra S, Tadmor E. Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws. SIAM J Numer Anal (2012) 50:544–73. doi:10.1137/110836961
150. Ranocha H. Shallow water equations: split-form, entropy stable, well-balanced, and positivity preserving numerical methods. GEM-International Journal on Geomathematics (2017) 8:85–133. doi:10.1007/s13137-016-0089-9
152. Barth TJ. Numerical methods for gasdynamic systems on unstructured meshes. In: D Kröner, M Ohlberger, and C Rohde, editors An introduction to recent developments in theory and numerics for conservation laws of lecture notes in computational science and engineering. Springer Berlin Heidelberg (1999), Vol. 5, p. 195–285.
153. Derigs D, Winters AR, Gassner GJ, Walch S. A novel averaging technique for discrete entropy-stable dissipation operators for ideal MHD. J Comput Phys (2017) 330:624–32. doi:10.1016/j.jcp.2016.10.055
155. Ray D, Chandrashekar P. An entropy stable finite volume scheme for the two dimensional Navier–Stokes equations on triangular grids. Appl Math Comput (2017) 314:257–86. doi:10.1016/j.amc.2017.07.020
156. Winters AR, Derigs D, Gassner GJ, Walch S. A uniquely defined entropy stable matrix dissipation operator for high Mach number ideal MHD and compressible Euler simulations. J Comput Phys (2017) 332:274–89. doi:10.1016/j.jcp.2016.12.006
158. Rojas D, Boukharfane R, Dalcin L, Fernández DCDR, Ranocha H, Keyes DE, et al. On the robustness and performance of entropy stable collocated discontinuous Galerkin methods. J Comput Phys (2020) 109891. doi:10.1016/j.jcp.2020.109891
160. Winters AR, Moura RC, Mengaldo G, Gassner GJ, Walch S, Peiro J, et al. A comparative study on polynomial dealiasing and split form discontinuous Galerkin schemes for under-resolved turbulence computations. J Comput Phys (2018) 372:1–21. doi:10.1016/j.jcp.2018.06.016
162. Saur J, Neubauer FM, Strobel DF, Summers ME. Three-dimensional plasma simulation of Io’s interaction with the Io plasma torus: asymmetric plasma flow. J Geophys Res: Space Physics (1999) 104:25105–26. doi:10.1029/1999JA900304
164. Saur J, Neubauer FM, Connerney J, Zarka P, Kivelson MG. Plasma interaction of Io with its plasma torus. In: Bagenal F, Dowling, TE, and McKinnon, WB, editors. Jupiter. The Planet, Satellites and Magnetosphere. vol. 1. Cambridge: Cambridge University Press (2004). p. 537–60.
166. Friedrich L, Winters AR, Fernández DCDR, Gassner GJ, Parsani M, Carpenter MH. An entropy stable non-conforming discontinuous Galerkin method with the summation-by-parts property. J Sci Comput (2018) 77:689–725. doi:10.1007/s10915-018-0733-7
167. Kopriva DA, Winters AR, Bohm M, Gassner GJ. A provably stable discontinuous Galerkin spectral element approximation for moving hexahedral meshes. Comput Fluids (2016) 139:148–60. doi:10.1016/j.compfluid.2016.05.023
168. Schnücke G, Krais N, Bolemann T, Gassner GJ. Entropy stable discontinuous Galerkin schemes on moving meshes for hyperbolic conservation laws. J Sci Comput (2020) 82:1–42. doi:10.1007/s10915-020-01171-7
169. Yamaleev NK, Fernandez DCDR, Lou J, Carpenter MH. Entropy stable spectral collocation schemes for the 3-D Navier-Stokes equations on dynamic unstructured grids. J Comput Phys (2019) 399:108897. doi:10.1016/j.jcp.2019.108897
171. Friedrich L, Schnücke G, Winters AR, Fernández DCDR, Gassner GJ, Carpenter MH. Entropy stable space–time discontinuous Galerkin schemes with summation-by-parts property for hyperbolic conservation laws. J Sci Comput (2019) 80:175–222. doi:10.1007/s10915-019-00933-2
174. Parsani M, Boukharfane R, Nolasco IR, Fernández DCDR, Zampini S, Hadri B, et al. High-order accurate entropy-stable discontinuous collocated Galerkin methods with the summation-by-parts property for compressible CFD frameworks: scalable SSDC algorithms and flow solver. J Comput Phys (2020) 424:109844. doi:10.1016/j.jcp.2020.109844
175. Ranocha H, Sayyari M, Dalcin L, Parsani M, Ketcheson DI. Relaxation Runge–Kutta methods: fully discrete explicit entropy-stable schemes for the compressible Euler and Navier–Stokes equations. SIAM J Sci Comput (2020) 42:A612–A638. doi:10.1137/19M1263480
176. Parsani M, Carpenter MH, Fisher TC, Nielsen EJ. Entropy stable staggered grid discontinuous spectral collocation methods of any order for the compressible Navier–Stokes equations. SIAM J Sci Comput (2016) 38:A3129–A3162. doi:10.1137/15m1043510
181. Hicken JE, Fernández DCDR, Zingg DW. Multidimensional summation-by-parts operators: general theory and application to simplex elements. SIAM J Sci Comput (2016) 38:A1935–A1958. doi:10.1137/15m1038360
182. Crean J, Hicken JE, Fernández DCDR, Zingg DW, Carpenter MH. Entropy-stable summation-by-parts discretization of the Euler equations on general curved elements. J Comput Phys (2018) 356:410–38. doi:10.1016/j.jcp.2017.12.015
185. Jameson A. Formulation of kinetic energy preserving conservative schemes for gas dynamics and direct numerical simulation of one-dimensional viscous compressible flow in a shock tube using entropy and kinetic energy preserving schemes. J Sci Comput (2008) 34:188–208. doi:10.1007/s10915-007-9172-6
186. Kuya Y, Kawai S. A stable and non-dissipative kinetic energy and entropy preserving (KEEP) scheme for non-conforming block boundaries on Cartesian grids. Comput Fluids (2020) 200:104427. doi:10.1016/j.compfluid.2020.104427
187. Ranocha H. Entropy conserving and kinetic energy preserving numerical methods for the Euler equations using summation-by-parts operators. In: SJ Sherwin, D Moxey, J Peiró, PE Vincent, and C Schwab, editors Spectral and high order methods for partial differential equations ICOSAHOM 2018. Cham: Springer International Publishing (2020). p. 525–35.
189. Dalcin L, Rojas D, Zampini S, Fernández DCDR, Carpenter MH, Parsani M. Conservative and entropy stable solid wall boundary conditions for the compressible Navier–Stokes equations: adiabatic wall and heat entropy transfer. J Comput Phys (2019) 397:108775. doi:10.1016/j.jcp.2019.06.051
191. Hindenlang FJ, Gassner GJ, Kopriva DA. Stability of wall boundary condition procedures for discontinuous Galerkin spectral element approximations of the compressible Euler equations. In: SJ Sherwin, D Moxey, J Peiró, PE Vincent, and C Schwab, editors. Spectral and high order methods for partial differential equations ICOSAHOM 2018. Springer International Publishing (2020). p. 3–19.
192. Parsani M, Carpenter MH, Nielsen EJ. Entropy stable wall boundary conditions for the three-dimensional compressible Navier–Stokes equations. J Comput Phys (2015c) 292:88–113. doi:10.1016/j.jcp.2015.03.026
Keywords: discontinuous Galerkin method, robustness, split form, dealiasing, summation-by-parts, second law of thermodynamics, entropy stability
Citation: Gassner GJ and Winters AR (2021) A Novel Robust Strategy for Discontinuous Galerkin Methods in Computational Fluid Mechanics: Why? When? What? Where?. Front. Phys. 8:500690. doi: 10.3389/fphy.2020.500690
Received: 01 October 2019; Accepted: 26 November 2020;
Published: 29 January 2021.
Edited by:Christian F. Klingenberg, Julius Maximilian University of Würzburg, Germany
Reviewed by:Hendrik Ranocha, King Abdullah University of Science and Technology, Saudi Arabia
Francesco Fambri, Max Planck Institute for Plasma Physics (IPP), Germany
Copyright © 2021 Gassner and Winters. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Andrew R. Winters, firstname.lastname@example.org