Entropy-Bounded Discontinuous Galerkin Scheme for Euler Equations
EEntropy-Bounded Discontinuous Galerkin Scheme for Euler Equations
Yu Lv ∗ , Matthias Ihme Department of Mechanical Engineering, Stanford University, Stanford, CA 94305, USA
Abstract
An entropy-bounded Discontinuous Galerkin (EBDG) scheme is proposed in which the solution is regular-ized by constraining the entropy. The resulting scheme is able to stabilize the solution in the vicinity ofdiscontinuities and retains the optimal accuracy for smooth solutions. The properties of the limiting op-erator according to the entropy-minimum principle are proofed analytically, and an optimal CFL-criterionis derived. We provide a rigorous description for locally imposing entropy constraints to capture multiplediscontinuities. Significant advantages of the EBDG-scheme are the general applicability to arbitrary high-order elements and its simple implementation for two- and three-dimensional configurations. Numericaltests confirm the properties of the scheme, and particular focus is attributed to the robustness in treatingdiscontinuities on arbitrary meshes.
Contents1 Introduction 22 Governing equations 33 Discontinuous Galerkin discretization 44 Entropy principle and entropy-bounded discontinuous Galerkin method 6 L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134.5 Numerical analysis of the limiting operator L . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 ∗ Corresponding author
Email addresses: [email protected] (Yu Lv), [email protected] (Matthias Ihme) a r X i v : . [ m a t h . NA ] N ov Generalization to multi-dimension and arbitrary elements 15 J ∂k The stabilization of solutions near flow-field discontinuities remains an open problem to the discontinuousGalerkin (DG) community. Considerable progress has been made on the development of limiters for two-dimensional quadrilateral and triangular elements. These limiters can be categorized into three classes.Methods that limit the solution using information about the slope along certain spatial directions [1, 2] fallin the first class. The second class of limiters extends this idea by limiting based on the moments of thesolution [3, 4], and schemes in which the DG-solution is projected onto a WENO [5, 6, 7] or Hermit WENO(HWENO) [8] representation fall in the last category. Although these limiters show promising results forcanonical test cases on regular elements and structured mesh partitions, the following two issues related topractical applications have not been clearly answered: • How can discontinuous solutions be regularized on multi-dimensional curved high-order elements? • How can non-physical solutions that are triggered by strong discontinuities and geometric singularitiesbe avoided?The present work attempts to simultaneously address both of these questions.2ecently, positivity-preserving DG-schemes have been developed for the treatment of flow-field discon-tinuities, and relevant contributions are by Zhang and Shu [9, 10, 11]. The positivity preserving methodprovides a robust framework with provable L -stability, preventing the appearance of negative pressure anddensity. Resulting algorithmic modifications are minimal, and these schemes have been used in simulationsof detonation systems with complex reaction chemistry [12, 13].Motivated by these attractive properties, the present work aims at developing an algorithm that avoidsnon-physical solutions on arbitrary elements and multi-dimensional spatial representations. The resultingscheme that will be developed in this work has the following properties: First, by invoking the entropyprinciple, solutions are constrained by a local entropy bound. Second, a general implementation on arbitraryelements is proposed without restriction to a specific quadrature rule. Third, the entropy constraint isimposed on the solutions through few algebraic operations, thereby avoiding the computationally expensiveinversion of a nonlinear system. Fourth, a method for the evaluation of an optimal CFL-criterion is derived,which is applicable to general polynomial orders and arbitrary element types.The remainder of this paper has the following structure. The governing equations and the discretizationare summarized in the next two sections. The entropy-bounded DG (EBDG) formulation is presented inSec. 4, and the derivation of the CFL-constraint and the limiting operator are presented. This analysisis performed by considering a one-dimensional setting, and the generalization to multi-dimensional andarbitrary elements is presented in Sec. 5. Section 6 is concerned with the evaluation of the entropy-boundedDG-scheme, and a detailed description of the algorithmic implementation is given in Sec. 7. The EBDG-method is demonstrated by considering several test cases, and the accuracy and stability are examined inSec. 8. The paper finishes with conclusions.
2. Governing equations
We consider a system of conservation equations, ∂ U ∂t + ∇ · F = 0 in Ω , (1)where the solution variable U : R × R N d → R N v and the flux term F : R N v → R N v × N d . Here, N d denotesthe spatial dimension and N v is the dimension of the solution vector. For the Euler equations, U and F takethe form: U ( x, t ) = ( ρ, ρu, ρe ) T , (2 a ) F ( U ) = ( ρu, ρu ⊗ u + p I , u ( ρe + p )) T , (2 b )where t is the time, x ∈ R N d is the spatial coordinate vector, ρ is the density, u ∈ R N d is the velocity vector, e is the specific total energy, and p is the pressure. Equation (1) is closed with the ideal gas law: p = ( γ − (cid:18) ρe − ρ | u | (cid:19) , (3)3n which γ is the ratio of specific heats, which, for the present work, is set to a constant value of γ = 1 . | · | to represent the Euclidean norm. With this, we define the localmaximum characteristic speed as, ν = | u | + c with c = (cid:114) γpρ , (4)where c is the speed of sound.Because of the presence of discontinuities in the solution of Eq. (1), we seek a weak solution that satisfiesphysical principles. This is the so-called entropy solution. By introducing U as a convex function of U with U : R N v → R , Lax [14] showed that the entropy solution of Eq. (1) satisfies the following inequality: ∂ U ∂t + ∇ · F ≤ , (5)where F : R N v → R N d is the corresponding flux of U . The consistency condition between Eqs. (1) and (5)requires [14]: (cid:18) ∂ U ∂ U (cid:19) T ∂ F ∂ U = ∂ F ∂ U . (6)The weak solution of Eq. (1) that satisfies this additional condition for the pair ( U , F ) is called an entropysolution. With this definition, Eq. (5) is commonly called entropy inequality or entropy condition, and U iscalled the entropy variable. A familiar example for gas-dynamic applications is to relate U to the physicalentropy s with: s = ln( p ) − γ ln( ρ ) + s , (7)where s is the reference entropy. The corresponding definition of the entropy variable and its flux in thecontext of the Euler system is [15]: ( U , F ) = ( − ρs, − ρsu ) . (8)Note that Eq. (7) directly provides a constraint on the positivity of pressure p and density ρ .
3. Discontinuous Galerkin discretization
We consider the problem to be posed on the domain Ω with boundary ∂ Ω. A mesh partition is definedas Ω = ∪ N e e =1 Ω e , where Ω e corresponds to a discrete element of this partition. The edge of element Ω e isdefined as ∂ Ω e . In order to distinguish different sides of the edge, the superscripts “+” and “ − ” are usedto denote the interior and exterior, respectively. We define a global space of test functions as V = ⊕ N e e =1 V e , V e = span { ϕ n (Ω e ) } N p n =1 , (9)where ϕ n is the n th polynomial basis, and N p is the number of bases. On the space V e we seek an approximatesolution to Eq. (1) of the form: U (cid:39) U = ⊕ N e e =1 U e , U e ∈ V e , (10)4here the solution vector U e on each individual element takes the general form U e ( x, t ) = N p (cid:88) m =1 (cid:101) U e,m ( t ) ϕ m ( x ) , (11)and the unknown vector of basic coefficients (cid:101) U e,m ∈ R N v × N p is obtained from the discretized weak solutionof Eq. (1): d (cid:101) U e,m dt (cid:90) Ω e ϕ n ϕ m d Ω − (cid:90) Ω e ∇ ϕ n · F ( U e ) d Ω + (cid:90) ∂ Ω e ϕ + n (cid:98) F ( U + e , U − e , (cid:98) n) d Γ = 0 , (12) ∀ ϕ n with n = 1 , . . . , N p . The numerical Riemann flux (cid:98) F is evaluated based on the states at both sides of theinterface ∂ Ω e and the outward-pointing normal vector (cid:98) n. It is of interest to note that for the particular caseof N p = 1 and ϕ = 1, the weak form reduces to the classical first-order finite-volume (FV) discretization. Itcan also be seen that the DG-scheme does not rely on a specific type of basis functions. Since the followingderivation is based on this mathematical property, we introduce the following lemma. Lemma 1
A polynomial P , P ( x ) = N p (cid:88) m =1 (cid:101) P m ϕ m ( x ) for x ∈ Ω e , (13) with a set of polynomial bases { ϕ m ( x ) , m = 1 , . . . , N p } , can be exactly interpolated by a Lagrangian polyno-mial of N p points { y n ∈ Ω e , n = 1 , . . . , N p } under the condition that (cid:2) ϕ m ( y n ) (cid:3) is non-singular: P ( x ) = N p (cid:88) n =1 P ( y n ) φ n ( x ) for x ∈ Ω e . (14) Proof: By equating Eqs. (13) and (14), and comparing terms, it follows that ϕ m ( x ) = N p (cid:88) n =1 ϕ m ( y n ) φ n ( x ) or [ ϕ m ( x )] = [ ϕ m ( y n )] [ φ n ( x )] , (15) where we use [ · ] to denote a tensor or a vector. Since [ ϕ m ( y n )] is non-singular, Eq. (15) can be inverted: [ φ n ( x )] = [ ϕ m ( y n )] − [ ϕ m ( x )] . Remark 1
The significance of this lemma is that it provides a description to convert any basis set toa Lagrangian basis set with N p interpolation points, as long as they are located at general positions. Tofacilitate the following derivation, we choose the points y n with n = 1 , . . . , N p from the N q quadraturepoints[16]. According to the accuracy requirement of the quadrature scheme for Eq. (12), N q ≥ N p is alwaytrue. . Entropy principle and entropy-bounded discontinuous Galerkin method In this section, we review the entropy principle by considering a three-point FV-setting. Then, we willexplore how to extend this principle to a DG-scheme, which leads to the concept of entropy boundedness.In order to enable the implementation of this concept, two important ingredients will be discussed, namelya time-step constraint and a limiting operator. After conducting numerical analyses by considering a one-dimensional configuration, we will extend the entropy boundedness to multi-dimensional and arbitraryelement types. The dimensional generality, geometric adaptability and simple implementation are majoradvantages of the resulting entropy-bounded DG-method.
To illustrate the entropy principle, we consider a local Lax-Friedrichs flux, which can be written as: (cid:98) F ( U L , U R , (cid:98) n) = 12 ( F ( U L ) + F ( U R )) · (cid:98) n − λ ( U R − U L ) , (16)and λ ≥ max k ∈{ L,R } ν ( U k )is the dissipation coefficient. Note that this flux function satisfies consistency: (cid:98) F ( U, U, (cid:98) n) = F ( U ) · (cid:98) n,conservation: (cid:98) F ( U L , U R , (cid:98) n) = − (cid:98) F ( U R , U L , − (cid:98) n), and Lipschitz-continuity. In the following, we consider thesimplest case of DGP0 scheme, with N p = 1, in a one-dimensional setting. This formulation is consistentwith the classical three-point FV-discretization. For x ∈ Ω e = [ x e − / , x e +1 / ], the discretized solution toEq. (1) can be written as: (cid:101) U e ( t + ∆ t ) = (cid:101) U e − ∆ th (cid:16) (cid:98) F ( (cid:101) U e , (cid:101) U e − , −
1) + (cid:98) F ( (cid:101) U e , (cid:101) U e +1 , (cid:17) , where (cid:101) U e is the basis coefficient, which is identical to the piecewise constant approximation to the exactsolution in Ω e . In the following, we introduce (cid:101) U ∆ te to denote the solution vector (cid:101) U e ( t + ∆ t ), and use thesuperscript ∆ t to denote a temporally updated quantity at t +∆ t . With the numerical flux given in Eq. (16),this discretization preserves the positivity of pressure and density under the CFL-condition [9, 17]:∆ tλh ≤ . (17)In addition, it was discussed in [17] that Eq. (17) satisfies the discrete minimum entropy principle proposedby Tadmor [15], s ( (cid:101) U ∆ te ) ≥ s e ( t ) = min j ∈{ e − ,e,e +1 } s ( (cid:101) U j ) . (18)6o show this property, we can rewrite Eq. (17) and split (cid:101) U e ( t + ∆ t ) into two parts. For x ∈ Ω e , this iswritten as (cid:101) U e ( t + ∆ t ) = 12 (cid:16) (cid:101) U ∆ te,p + (cid:101) U ∆ te,p (cid:17) , (19 a ) (cid:101) U ∆ te,p = (cid:101) U e − ∆ th (cid:16) F ( (cid:101) U e +1 ) − λ e +1 / (cid:101) U e +1 − F ( (cid:101) U e ) + λ e +1 / (cid:101) U e (cid:17) , (19 b ) (cid:101) U ∆ te,p = (cid:101) U e + ∆ th (cid:16) F ( (cid:101) U e − ) + λ e − / (cid:101) U e − − F ( (cid:101) U e ) − λ e − / (cid:101) U e (cid:17) , (19 c )where (cid:101) U ∆ te,p and (cid:101) U ∆ te,p can be viewed as the P0-approximations to the solutions of the hyperbolic systems(under the CFL constraint of Eq. (17)) ∂ U ∂t + (cid:0) F (cid:48) ( U ) − λ e +1 / I (cid:1) ∂ U ∂x = 0 , (20 a ) ∂ U ∂t + (cid:0) F (cid:48) ( U ) + λ e − / I (cid:1) ∂ U ∂x = 0 , (20 b )with the exact (Godunov) flux. If we denote the exact solutions to Eq. (20 a ) and (20 b ) as U p ( x, t + ∆ t )and U p ( x, t + ∆ t ), respectively, then their P0-approximations in Ω e yield (cid:101) U ∆ te,p = h (cid:82) x e +1 / x e − / U p ( x, t + ∆ t ) dx and (cid:101) U ∆ te,p = h (cid:82) x e +1 / x e − / U p ( x, t + ∆ t ) dx . Both equation systems are obtained by imposing a constant shifton the characteristic speeds without modifying the characteristic variables. With these modifications, allcharacteristics in Eq. (20 a ) are right-running while those in Eq. (20 b ) are left-running. The correspondingentropy inequalities take the form ∂ U ∂t + ∂∂x (cid:0) F − λ e +1 / U (cid:1) ≤ , (21 a ) ∂ U ∂t + ∂∂x (cid:0) F + λ e − / U (cid:1) ≤ . (21 b )Without loss of generality, we now consider Eq. (21 a ) and integrate over [ t, t + ∆ t ] × [ x e − / , x e +1 / ],resulting in the following expression: x e +1 / (cid:90) x e − / U (cid:0) U p ( x, t + ∆ t ) (cid:1) dx − x e +1 / (cid:90) x e − / U ( (cid:101) U e ) dx + t +∆ t (cid:90) t (cid:0) F ( U ( x e +1 / , t )) − λ e +1 / U ( U ( x e +1 / , t )) (cid:1) dt − t +∆ t (cid:90) t (cid:0) F ( U ( x e − / , t )) − λ e − / U ( U ( x e − / , t )) (cid:1) dt ≤ . (22)Recognizing that all characteristics are right-running, the temporal integral can be evaluated exact since U ( x e − / , t ) = (cid:101) U e and U ( x e +1 / , t ) = (cid:101) U e +1 under the condition of Eq. (17). Then by utilizing the convexityof U with respect to U , the following estimate is obtained: U ( (cid:101) U ∆ te,p ) = U (cid:32) h (cid:90) x e +1 / x e − / U p ( x, t + ∆ t ) dx (cid:33) ≤ h (cid:90) x e +1 / x e − / U ( U p ( x, t + ∆ t )) dx ≤ U ( (cid:101) U e ) + ∆ th (cid:16) F ( (cid:101) U e ) − λ e +1 / U ( (cid:101) U e ) (cid:17) − ∆ th (cid:16) F ( (cid:101) U e +1 ) − λ e +1 / U ( (cid:101) U e +1 ) (cid:17) . U , F ), given in Eq. (8), it follows s ( (cid:101) U ∆ te,p ) ≥ ρ e ρ ∆ te,p (cid:20) − ∆ th (cid:0) λ e +1 / − u e (cid:1)(cid:21) s ( (cid:101) U e ) + ρ e +1 ρ ∆ te,p ∆ th (cid:0) λ e +1 / − u e +1 (cid:1) s ( (cid:101) U e +1 ) . (23)The constraint (17) ensures that the coefficients in front of s ( (cid:101) U e ) and s ( (cid:101) U e +1 ) are positive and sum to unityaccording to Eq. (19 b ). From these arguments directly follows: s ( (cid:101) U ∆ te,p ) ≥ min { s ( (cid:101) U e ) , s ( (cid:101) U e +1 ) } , (24)and s ( (cid:101) U ∆ te,p ) ≥ min { s ( (cid:101) U e − ) , s ( (cid:101) U e ) } . (25)Combining these two relations with the quasi-concavity of the entropy s (Lemma 2.1 of [11]), the discreteminimum entropy principle of Eq. (18) is obtained.The result above is obtained for a one-dimensional setting. However, in the following we derive arotational version of this entropy principle for multi-dimensional cases by following the idea of [17]. Byconsidering an arbitrary face with a normal (cid:98) n, we can define a tangential vector (cid:98) t. By projecting thevelocity in Cartesian coordinates onto this local coordinate ( (cid:98) n , (cid:98) t) with the following mapping, u → ( u n , u t ) T , with u n = u · (cid:98) n , u t = u · (cid:98) t , with which, the conservative variables and flux terms can be written as follows U = ( ρ, ρu n , ρu t , ρe n , ρe t ) T , (26 a ) F = ( ρu n , ρu n + p, ρu t u n , u n ( ρe n + p ) , ρu n e t ) T , (26 b )and e n = pρ ( γ −
1) + 12 u n , e t = 12 u t . For this augmented system, it can be seen that the variations of density, pressure and entropy are all governedby a 1D reduced system that is parallel to (cid:98) n: U = ( ρ, ρu n , ρe n ) T , (27 a ) F = ( ρu n , ρu n + p, u n ( ρe n + p )) T , (27 b )and its discretized version is identical to Eq. (17). Thus, the solution preserves the positivity of density andpressure and satisfies the entropy principle under the CFL-constraint, Eq. (17), with λ ≥ max j ∈{ e − ,e,e +1 } ν ( U j ).We now conclude the above analysis with the following lemma as a critical element for the subsequentderivation. 8 emma 2 For a three-point system defined on R N d , the solution along an arbitrary direction (cid:98) n, (cid:101) U ∆ te = (cid:101) U e + ∆ th (cid:16) (cid:98) F ( (cid:101) U e , (cid:101) U e − , − (cid:98) n ) + (cid:98) F ( (cid:101) U e , (cid:101) U e +1 , (cid:98) n ) (cid:17) , (28) with the flux function (cid:98) F specified in Eq. (16), preserves the positivity of density and pressure, and satisfiesthe entropy principle: s ( (cid:101) U ∆ te ) ≥ min j ∈{ e − ,e,e +1 } s ( (cid:101) U j ) , (29) under the CFL condition: ∆ tλh ≤ , λ ≥ max j ∈{ e − ,e,e +1 } ν ( (cid:101) U j ) , This three-point system is consistent with that used by Zhang and Shu [9, 11]. The difference is that weintroduce a local entropy bound s e at time t instead of a global entropy bound that is derived from theinitial conditions min x ∈ Ω s ( x, To robustly capture shocks while retaining the high-order benefit of the DG-scheme, sub-cell shockresolution is required [20]. We now extend the discussion by considering a high-order DG-solution withsub-cell representation. In each DG-cell, the whole solution is approximated by a function space. However,there is no guarantee that the high-order ( N p >
1) solution obeys the physical entropy principle. This isthe reason that DG suffers from numerical instability in the vicinity of discontinuities. To suppress theseinstabilities, one approach is to consider imposing constraints based on the behavior of the entropy solution.The positivity-preserving DG-method [9, 10] is a successful example for this approach. Based on the entropyprinciple, Eq. (18), Zhang and Shu [11] extended their implementation to an entropy-based constraint. Here,we propose a general framework that is based on the entropy principle, and major differences and advantageshave been highlighted in Sec. 1.We define the constraint for the high-order DG-scheme as follows: ∀ x ∈ Ω e , s ( U ∆ te ( x )) ≥ min { s ( U ( y )) | y ∈ Ω e ∪ ∂ Ω − e } ≡ s e ( t ) . (30)In this equation, the right-hand-side sets an entropy bound for an element-local solution in Ω e ; with this, werefer to a DG-solution as entropy-bounded if it satisfies this principle. s e ( t ) is a local estimate for the trueentropy minimum in Ω e , | s e ( t ) − min x ∈ Ω e s ( U ( x )) | ∼ O ( h k ), where k is the local order of accuracy. Besides that, s e ( t ) is bounded if the entropy is bounded at the domain boundaries, s e ( t ) ≥ min x ∈ Ω s ( U ( x, t = 0)) = s , where s is the minimum entropy at the initial condition. By imposing this constraint, we expect that the sub-cell9 = t n t = t n + Δ t s e0 (t n ) s e0 (t n ) Ω e Ω e -1 Ω e +1 U s(U) L U Figure 1: (Color online) Schematic of entropy-bounding of the EBDG scheme.
DG-solution is regularized, avoiding the appearance of non-physical solutions. This idea is illustrated in Fig.1. At time level t , s e ( t ) is calculated and used to set a reference bound for the solution at the next step, U ∆ te . If U ∆ te yields entropy undershoot with respect to s e ( t ), it will be modified to satisfy the constraintof Eq. (30). In order to implement this regularization for a high-order DG scheme, the following aspectsrequire addressing:i To impose Eq. (30) on the DG-solution, we introduce a limiting operator L . The regularized solution,denoted by L U ∆ te , requires that s ( L U ∆ te ( x )) ≥ s e ( t ) ∀ x ∈ Ω e . In the following, we relax this condition,and impose Eq. (30) only on the set of quadrature points, D , that are involved in solving the weakform in Eq. (12).ii Guaranteeing that the constraint (30) is always imposed requires the existence of the operator L . Asufficient condition for this is that the element-averaged solution is entropy-bounded, s ( U ∆ te ) ≥ s e ( t ).Enforcing this condition relies on the selection of a proper CFL-condition, and this analysis will bedeveloped in Sec. 4.3 for a one-dimensional system. Subsequently, this analysis is then extended inSec. 5.2 to general multi-dimensional elements.iii Algorithmic details on the implementation of the operator L constraining the element-local DG-solution are discussed in Sec. 4.4.iv The evaluation of the lower bound s e ( t ) that is necessary to constrain the entropy solution is given inSec. 6. The objective now is to extend the analysis for DGP0 to a DG-scheme with high-order polynomial repre-sentations. Consider a one-dimensional domain in which the element Ω e is centered at x e , and a quadrature10ule with weights w q and (cid:80) N q q =1 w q = 1. These quadrature weights are evaluated at the quadrature points x q ∈ [ x e +1 / , x e − / ]. The discretized cell-averaged solution U e is defined as: U e = 1 h (cid:90) Ω e U e dx, (31)(for the P0-case discussed above, U e = (cid:101) U e ), which can be further expanded by a quadrature rule withsufficient accuracy: U e = N q (cid:88) q =1 w q U e ( x q ) , = N q (cid:88) q =1 (cid:0) w q − θ l φ q ( x e − / ) − θ r φ q ( x e +1 / ) (cid:1) U e ( x q ) + θ l U e ( x e − / ) + θ r U e ( x e +1 / ) , = N q (cid:88) q =1 θ q U e ( x q ) + θ l U e ( x e − / ) + θ r U e ( x e +1 / ) , (32)where the first line utilizes the exactness of the quadrature rule, the second line utilizes Lemma 1, and thethird line defines θ q = w q − θ l φ q ( x e − / ) − θ r φ q ( x e +1 / ). Under the condition that θ r,l > θ q ≥
0, thelast line of Eq. (32) is a convex combination. Since the quadrature weights w q are positive, the existenceof θ r,l is guaranteed through the condition w q ≥ θ l φ q ( x e − / ) + θ r φ q ( x e +1 / ). If φ q ( x e ± / ) > θ r,l isconstrained as (0 , min q w q / max { φ q ( x e ± / ) } ]. If some of φ q ( x e ± / ) are negative, they are not essential insetting the upper bound for θ r,l . Remark 2
In the following, θ r,l will be related to a CFL-constraint. To obtain an optimal CFL-number,the largest value of θ r,l needs to be found. This can be formulated as a maximization problem subject to theconstraints, θ r,l > and θ q ≥ . For illustration, we fully discretize Eq. (12) using a forward Euler time integration scheme and insert theresults from Eq. (32). The element-averaged solution in Ω e is then updated as: U ∆ te = U e − ∆ th (cid:16) (cid:98) F ( U e ( x e − / ) , U e − ( x e − / ) , −
1) + (cid:98) F ( U e ( x e +1 / ) , U e +1 ( x e +1 / ) , (cid:17) , = N q (cid:88) q =1 θ q U e ( x q ) + θ l U e ( x e − / ) − ∆ th (cid:16) (cid:98) F ( U e ( x e − / ) , U e − ( x e − / ) , −
1) + (cid:98) F ( U e ( x e − / ) , U ∗ e , (cid:17) + θ r U e ( x e +1 / ) − ∆ th (cid:16) (cid:98) F ( U e ( x e +1 / ) , U ∗ e , −
1) + (cid:98) F ( U e ( x e +1 / ) , U e +1 ( x e +1 / ) , (cid:17) , (33)where U ∗ e = 12 (cid:0) U e ( x e − / ) + U e ( x e +1 / ) (cid:1) − λ (cid:0) F ( U e ( x e +1 / )) − F ( U e ( x e − / )) (cid:1) (34)11s introduced to simplify subsequent analyses. Note that U ∗ e ensures the validity of the second equality inEq. (33) with a λ that is defined in the following lemma. We can see that Eq. (33) contains two three-pointsystems discussed in Sec. 4.1. To guarantee that U ∆ te is entropy-bounded, it is necessary that these systemsconform to the entropy principle of Eq. (18). This leads to the following lemma. Lemma 3
For a one-dimensional DG-system, the element-averaged solution satisfies the entropy principle s ( U ∆ te ) ≥ s e ( t ) = min { s ( U ( y )) | y ∈ Ω e ∪ ∂ Ω − e } , (35) under the CFL-constraint ∆ tλh ≤
12 min { θ l , θ r } , (36) where λ is the maximum wave speed that is evaluated over the set of point-wise solutions λ ≥ max U ν ( U ) with U ∈ { U e − ( x e − / ) , U e ( x e ± / ) , U e +1 ( x e +1 / ) } , (37) given the conditions θ r,l > and θ q ≥ .Proof: First, we need to show that U ∗ e satisfies the discretized entropy principle of Eq. (18). This is indeed thecase since U ∗ e is essentially the left-hand-side of Eq. ( b ) with time step taking ∆ t = h/ (2 λ ) , correspondingto the upper bound of the CFL-constraint. Since U ∗ e is an entropy solution, ν ( U ∗ e ) is bounded by λ , whichmakes the Lax-Friedrichs flux (cid:98) F involving U ∗ e in Eq. (33) valid according to the definition in Eq. (16).Therefore, we have s ( U ∗ e ) ≥ min (cid:8) s ( U e ( x e − / )) , s ( U e ( x e +1 / )) (cid:9) , Second, we reformulate Eq. (33) as U ∆ te = N q (cid:88) q =1 θ q U e ( x q ) + θ l U ∆ te,p + θ r U ∆ te,p , (38) in which U ∆ te,p and U ∆ te,p are the two updated solutions of the three-point system. Their definitions are readilyobtained by comparing Eqs. (38) and (33). The given constraints, θ r,l > and θ q ≥ , guarantee that theform of the convex combination in Eq. (38) always holds. According to Lemma 2, it follows s ( U ∆ te,p ) ≥ min (cid:8) s ( U ∗ e ) , s ( U e ( x e − / )) , s ( U e − ( x e − / )) (cid:9) ,s ( U ∆ te,p ) ≥ min (cid:8) s ( U ∗ e ) , s ( U e ( x e +1 / )) , s ( U e +1 ( x e +1 / )) (cid:9) , under the given CFL-constraint, Eq. (36). Combining this with the quasi-concavity of entropy [11], it follows s ( U ∆ te ) ≥ min (cid:110) s ( U e ( x q )) , s ( U ∆ te,p ) , s ( U ∆ te,p ) (cid:111) , ≥ min { s ( U ( y )) | y ∈ Ω e ∪ ∂ Ω − e } . emark 3 Equation (36) ensures the positivity of p ( U ∆ te ) and ρ ( U ∆ te ) . In this context, we emphasize that the CFL-constraint (36) provides a description for the entropy bounded-ness and does not conflict with the general CFL-constraint for linear stability, CFL L . To distinguish bothconstraints, here we use CFL EB to denote the CFL-number for guaranteeing the entropy boundedness. Ingeneral, the time step has to be selected to satisfy both criteria. Equation (36) shows that CFL EB dependson the value min { θ l , θ r } , and a rigorous evaluation for this will be given below.Although we considers the specific case of a forward Euler time discretization scheme, all the derivationand conclusions are directly applicable to any explicit Runge-Kutta (RK) methods with positive coefficients,since the RK-solution is a convex combination of solutions obtained from several forward Euler sub-steps.In practice, RK-methods are preferred as DG time-integration schemes due to their compatible stabilityproperties [21]. L Following Lemma 3, the entropy constraint is imposed on the set of quadrature points, x ∈ D ⊂ Ω e . Forthe one-dimensional case, D is: D = { x e ± / , x q , q = 1 , . . . , N q ( N q ≥ N p ) } . (39)In the following, we are concerned with the construction of a limiting operator L , such that ∀ x ∈ D , s ( L U ∆ te ( x )) ≥ s e ( t ) , (40)Since the operator L is applied at the end of each sub-iteration, we will omit the superscript ∆ t in thesubsequent analysis. According to the entropy definition (7), Eq. (40) can be written as: p ( L U e ( x )) ≥ exp( s e ) ρ γ ( L U e ( x )) ∀ x ∈ D . (41)To define the operator L , we follow the work of Zhang and Shu [9, 11], and introduce a linear scaling: L U e = U e + ε ( U e − U e ) . (42)The parameter ε is then determined by substituting Eq. (42) into Eq. (41) and by applying Jensen’s in-equality: p (cid:0) (1 − ε ) U e + εU e (cid:1) ≥ (1 − ε ) p ( U e ) + εp ( U e ) , ≥ exp( s e ) (cid:2) (1 − ε ) ρ γ ( U e ) + ερ γ ( U e ) (cid:3) , (43) ≥ exp( s e ) ρ γ (cid:0) (1 − ε ) U e + εU e (cid:1) . Solving for ε gives ε = ττ − [ p ( U e ) − exp( s e ) ρ γ ( U e )] with τ = min (cid:26) , min x ∈D { p ( U e ( x )) − exp( s e ) ρ γ ( U e ( x )) } (cid:27) , (44)13hich is subject to the conditions ρ ( U e ) > , p ( U e ) > exp( s e ) ρ γ ( U e ) . (45)These conditions are automatically guaranteed through the CFL EB -constraint of Lemma 3. While thepositivity condition for pressure is embedded in Eq. (43), the positivity of density must be imposed for all x ∈ D before L is applied, and the methodology for this is presented in [10].Compared to the limiting operator, presented in [11], the herein proposed method is substantially sim-plified. Specifically, the step for imposing the positivity of pressure is avoided; in addition, ε is obtainedfrom an algebraic relation, and does not require a computationally expensive Newton iteration. It is alsonoted that the operator L contains the positivity-preserving limiter as a special case, which is obtained bysetting s e → −∞ . L In this section, numerical properties of the limiting operator are examined.
Conservation.
Integrating Eq. (42) over Ω e , (cid:90) Ω e L U e dx = (1 − ε ) (cid:90) Ω e U e dx + εU e (cid:90) Ω i dx = (cid:90) Ω e U e dx , confirms that the limiting operator preserves the conservation properties of the solution vector. Stability.
Since the positivity of density and pressure is preserved at the quadrature points, L is L -stable,which was shown in [9, 11]. Here, we extend this stability analysis and evaluate the L -stability. Byconsidering a periodic domain and taking the L -norm of Eq. (42) we obtain: || L U e || = (cid:90) Ω e [ U e + ε ( U e − U e )] dx , (46 a )= (1 − ε ) (cid:90) Ω e U e dx − ε ( ε − (cid:90) Ω e U e dx , (46 b ) ≤ (1 − ε ) (cid:90) Ω e U e dx − ε ( ε − (cid:90) Ω e U e dx , (46 c ) ≤ || U e || . (46 d )After integrating over the entire domain, we obtain || L U || ≤ || U || , which shows that L does not affect the stability of the DG-discretization. Further, since L constrains pressureand density, λ in Eq. (36) provides a robust CFL-criterion, without the need for arbitrarily reducing ∆ t toincrease the stability region. 14 ccuracy. In regions where the solution is smooth, we assume that the weak solution before limiting hasoptimal accuracy: || U − U || Ω ≤ C h p +1 , and that undershoots in entropy remain small, so that ε ∼ O ( h p ). Thus, the error is estimated as follows: || L U − U || = (cid:88) e || L U e − U e || , (47 a )= (cid:88) e || ε ( U e − U e ) + (1 − ε )( U e − U e ) || , (47 b ) ≤ (cid:88) e (cid:0) ε || ( U e − U e ) || + 2(1 − ε ) || U e − U e || (cid:1) , (47 c ) ≤ C h p +2 , (47 d )where for simplicity, we introduce U e to denote the element-wise representations to U . Here we use the factthat U e is locally a first-order approximation to U e , U e = U e + C e O ( h ).In the vicinity of a discontinuity, the DG-solution looses its regularity so that the convergence ratereduces to first-order: || U − U || Ω ≤ C h. Triggered by spurious sub-cell solutions, the entropy undershoot can be very large, so that ε ∼ O (1). Byrepeating the above argument, we obtain an estimate for the accuracy of the discontinuous solution: || L U − U || Ω ≤ C h. The accuracy arguments given here are substantiated through numerical tests in Sec. 8.
5. Generalization to multi-dimension and arbitrary elements
The entropy-bounded DG scheme that was presented for one-dimensional systems in the previous sectioncan be generalized to arbitrary elements in multi-dimensions. This extension is the subject of the followinganalysis.Since EBDG does not rely on a specific quadrature rule, any quadrature method can be used as longas it accurately integrates the problem and ensures the positivity of the quadrature weights. The limitingprocedure requires the definition of a new set of quadrature points D for the general multi-dimensionalsetting. The selection of these points is given in the next section. The extension to arbitrary elementsrequires special consideration of the CFL EB number. To present a general formulation for multi-dimensional configurations, we first introduce necessary no-tations to describe general elements with curvatures. For this, we define a geometric mapping function15 : R N d → R N d on a reference element Ω r e , such that x = Φ( r ) maps any point r ∈ Ω r e onto x ∈ Ω e , and J = [ ∂x/∂r ] is the geometric Jacobian. With these specifications, we can write the discretized state vectoras: U e ( x, t ) = N p (cid:88) m =1 (cid:101) U e,m ( t ) ϕ m ( r ) , x = x ( r ) ∈ Ω e , ∀ r ∈ Ω r e (48)The mapping function is commonly parameterized by a polynomial function x ( r ) = (cid:80) N g m =1 (cid:101) x m ϕ gm ( r ), where ϕ gm ( r ) is a Lagrangian interpolation and N g is the number of geometric bases used to represent Ω e . Sincethe reference element is regular, we can use a subspace of r to parameterize the element edges. Therefore,to parameterize the k th edge of Ω e we define g k = P k ( r ) ∈ R N d − such that ∀ r ∈ ∂ Ω r e,k , r = P − k ( g k ), inwhich P − k is the pseudo-inverse of P k . For the physical element, the edge can be represented as: ∂ Ω e,k = (cid:8) x ∈ Ω e | x = Φ( r ) , r = P − k ( g k ) ∈ ∂ Ω r e,k (cid:9) . (49)The integral in Eq. (12) is evaluated using multi-dimensional quadrature rules. Considering the com-plexity of the dimensionality, here we follow the quadrature convention that is (cid:80) N q v =1 w v = V r e (the volumeof Ω r e ) and (cid:80) N kq q =1 w k,q = S r e,k (the area of ∂ Ω e,k ). With these preliminaries, we can evaluated any volumeintegral in Eq. (12) as: (cid:90) Ω e f ( x ) dx = (cid:90) Ω r e |J ( r ) | f ( x ( r )) dr = N q (cid:88) v =1 |J ( r v ) | f ( x ( r v )) w v . (50)The surface integral of a scalar function on ∂ Ω e,k can be written as: (cid:90) ∂ Ω e,k f ( x ) d Γ = (cid:90) ∂ Ω r e,k f ( x ( g k )) |J ∂k | dg k = N ∂q,k (cid:88) q =1 |J ∂k ( g k,q ) | f ( x ( g k,q )) w k,q , (51)and the surface integral of a vector-valued function is evaluated as: (cid:90) ∂ Ω e,k f ( x ) · (cid:98) n d Γ = (cid:90) ∂ Ω r e,k f ( x ( g k )) · J ∂k dg k = N ∂q,k (cid:88) q =1 |J ∂k ( g k,q ) | f ( x ( g k,q )) · (cid:98) n( g k,q ) w k,q , (52)where J ∂k is the surface Jacobian, and (cid:98) n refers to the unit vector J ∂k / |J ∂k | . Note that for a general element,the quadrature expression might be subject to a tiny numerical error bounded by O ( h p +1 ), given theaccuracy requirement for integrating Eq. (12).With the above notation, we are now able to define the set of quadrature points D for general curvedelements: D = N ∂ (cid:91) k =1 { g k,q , q = 1 , . . . , N ∂q,k } (cid:91) { r v , v = 1 , . . . , N q } , (53)where N ∂ is the number of element edges (which is equal to the number of neighbor elements). In thiscontext, it is noted that D includes all quadrature points that are involved in the integration. With thisspecification of D , the limiting operator L , developed in Sec. 4.4, can be directly extended to arbitrary16lements on multi-dimensional configurations. In the following, a CFL EB -constraint is derived that extendsthe results of Lemma 3, thereby ensuring the existence of the general limiter L . Following the same approach as for the one-dimensional derivation in Sec. 4.3, the element-averagedsolution of U ∆ te is evaluated as: U ∆ te = U e − ∆ tV e N ∂ (cid:88) k =1 (cid:90) ∂ Ω e,k (cid:98) F (cid:0) U + e , U − e , (cid:98) n (cid:1) d Γ , (54 a )= N q (cid:88) v =1 |J ( r v ) | w v V e U e ( r v ) − N ∂ (cid:88) k =1 N ∂q,k (cid:88) q =1 ∆ t |J ∂k ( g k,q ) | w k,q V e (cid:98) F (cid:0) U + e ( r ( g k,q )) , U − e ( r ( g k,q )) , (cid:98) n( g k,q ) (cid:1) , (54 b )= N q (cid:88) v =1 θ v U e ( r v ) + N ∂ (cid:88) k =1 N ∂q,k (cid:88) q =1 (cid:2) θ k,q U + e ( r ( g k,q )) − ∆ tζ k,q (cid:16) (cid:98) F (cid:0) U + e ( r ( g k,q )) , U − e ( r ( g k,q )) , (cid:98) n( g k,q ) (cid:1) + (cid:98) F (cid:0) U + e ( r ( g k,q )) , U ∗ e , − (cid:98) n( g k,q ) (cid:1)(cid:17)(cid:105) , (54 c )where θ v = | J ( r v ) | w v V e − N ∂ (cid:88) k =1 N ∂q,k (cid:88) q =1 θ k,q φ v ( r ( g k,q )) (55)is introduced to decompose the volumetric quadrature to obtain U + e ( r ( g k,q )). For notational simplification,we defined ζ k,q = |J ∂k ( g k,q ) | w k,q V e , (56)so such that S e = (cid:80) N ∂ k =1 (cid:80) N ∂q,k q =1 ζ k,q is equal to the surface area of Ω e , and (cid:80) N ∂ k =1 (cid:80) N ∂q,k q =1 ζ k,q (cid:98) n( g k,q ) = 0since Ω e has a closed surface. To apply the results from the three-point system to the multi-dimensionalconfiguration, we introduce the auxiliary variable U ∗ e : U ∗ e = N ∂ (cid:88) k =1 N ∂q,k (cid:88) q =1 ζ k,q S e (cid:20) U + e ( r ( g k,q )) − λ ∗ F ( U + e ( r ( g k,q ))) · (cid:98) n( g k,q ) (cid:21) . (57)It can be shown that U ∗ e is essentially the solution to the following equation: N ∂ (cid:88) k =1 N ∂q,k (cid:88) q =1 ζ k,q (cid:98) F (cid:0) U + e ( r ( g k,q )) , U ∗ e , − (cid:98) n( g k,q ) (cid:1) = 0 , (58)subject to a preselected dissipation coefficient λ ∗ , so that the equality in Eq. (54 c ) holds true. Here, weevaluate λ ∗ from the following relation: λ ∗ = τ max (cid:8) ν ( U ) | U ∈ { U + e ( r ( g k,q )) } (cid:9) , τ = max (cid:110)(cid:112) N d , (cid:112) γ ( γ − (cid:111) (59)and the rationale for this selection is provided later in Remark 4. To prove that U ∆ te is entropy bounded,we present the following lemma. 17 emma 4 U ∗ e in Eq. (57) satisfies s ( U ∗ e ) ≥ min { s ( U ) | U ∈ { U + e ( r ( g k,q )) }} .Proof: For notational simplification, we combine the indices k and q into a single index j , and we denote thetotal number of surface quadrature points on ∂ Ω e by N tot , N tot = (cid:80) N ∂ k =1 N ∂q,k . Considering (cid:80) N tot j =1 ζ j (cid:98) n ( d ) j = 0 and ζ j > , the d th components of the surface-normal vectors have different signs. To denote each component,we introduce the superscript ( d ) . By sorting (cid:98) n ( d ) j so that the first N > tot vector components are positive. Thefollowing statement is true for any d : N > tot (cid:88) j =1 ζ j (cid:98) n ( d ) j = − N tot (cid:88) j = N > tot +1 ζ j (cid:98) n ( d ) j = N par (cid:88) n =1 l n , (60) where l n introduces a partition as illustrated in Fig. 2 and N par is the dimension of this partition. Withthis, we are able introduce a variable mapping, U s + n = U + e ( r j ) , if j − (cid:88) i =1 ζ i (cid:98) n ( d ) i < n (cid:88) i =1 l i ≤ j (cid:88) i =1 ζ i (cid:98) n ( d ) i ,U s − n = U + e ( r j ) , if − j − (cid:88) i = N > tot +1 ζ i (cid:98) n ( d ) i < n (cid:88) i =1 l i ≤ − j (cid:88) i = N > tot +1 ζ i (cid:98) n ( d ) i . With this, Eq. (57) is equivalent to: U ∗ e = 1 S e N tot (cid:88) j =1 ζ j U + e ( r j ) − N tot (cid:88) j =1 ζ j λ ∗ F ( U + e ( r j )) · (cid:98) n j , = N tot (cid:88) j =1 ζ j S e (cid:32) − N d (cid:88) d =1 | (cid:98) n ( d ) j |√ N d (cid:33) U + e ( r j ) + N d (cid:88) d =1 N tot (cid:88) j =1 S e (cid:32) ζ j | (cid:98) n ( d ) j |√ N d U + e ( r j ) − ζ j (cid:98) n ( d ) j λ ∗ F ( d ) ( U + e ( r j )) (cid:33) , = N tot (cid:88) j =1 ζ j S e (cid:32) − N d (cid:88) d =1 | (cid:98) n ( d ) j |√ N d (cid:33) U + e ( r j ) + N d (cid:88) d =1 S e √ N d N par (cid:88) n =1 l n U ∗∗ d,n , where we introduce U ∗∗ d,n = 12 ( U s + n + U s − n ) − √ N d λ ∗ (cid:16) F ( d ) ( U s + n ) − F ( d ) ( U s − n ) (cid:17) , which takes the same form as the left-hand-side of Eq. (34). Note that U ∗∗ d,n is essentially expressed in aone-dimensional setting along x ( d ) . Therefore, one can follow the same argument used for Eq. (34) inLemma 3 to verify that s ( U ∗∗ d,n ) ≥ min (cid:8) s ( U s ± n ) (cid:9) ≥ min (cid:8) s ( U ) | U ∈ { U + e ( r j ) } , j = 1 , . . . , N tot (cid:9) , with the given form of λ ∗ in Eq. (59). As given above, U ∗ e is a convex combination; by using the quasi-concavity of entropy [11], we conclude that s ( U ∗ e ) > min (cid:8) s ( U ) | U ∈ { U + e ( r j ) } , j = 1 , . . . , N tot (cid:9) . nj n j j n j +sign signDefined partition: (d)(d) Figure 2: Illustration of the partition introduced in Lemma 4.
Remark 4
Note that the maximum characteristic speed of U ∗∗ d,n is bounded, ν ( U ∗∗ d,n ) ≤ ν ( U + e ( r )) . Accordingto the combination law given in Appendix A, we have ν ( U ∗ e ) ≤ (cid:112) γ ( γ −
1) max { ν ( U + e ( r )) } ≤ λ ∗ . Withthis, we are now able to show that the maximum characteristic speed of U ∗ e is bounded by the chosen valueof λ ∗ , so that the Lax-Friedrichs flux (cid:98) F ( U + e ( r ) , U ∗ e , − (cid:98) n ) in Eq. (54) is valid according to the definition ofEq. (16). To enforce the entropy boundedness, the decomposition of U ∆ te in Eq. (54) is required to be convex. Thiscan be satisfied under the following condition: θ v ≥ , ∀ v = 1 , . . . , N q ,θ k,q > , ∀ ( k, q ) , q = 1 , . . . , N ∂k , k = 1 , . . . , N ∂ , (61)With this, the entropy boundedness of U ∆ te is shown by the following lemma. Lemma 5
For a general DG element, the element-averaged solution is entropy bounded, s ( U ∆ te ) ≥ s e ( t ) = min { s ( U ( y )) | y ∈ Ω e ∪ ∂ Ω − e } , (62) under the condition that Eq. (61) holds and that the following constraint is fulfilled: ∆ tλ ≤
12 min (cid:26) θ k,q ζ k,q (cid:27) , ∀ ( k, q ) , q = 1 , · · · , N ∂k , k = 1 , · · · , N ∂ , (63) where λ ≥ max { ν ( U ) | U ∈ { U ± e ( r ( g k,q )) } and λ ≥ λ ∗ .Proof: The proof follows Lemma 3, utilizing Lemma 2 and the quasi-concavity of entropy. Note that Lemma 5 does not rely on any assumption regarding the dimensionality or shape of the finiteelement, and is therefore general. Another observation is that Eq. (63) essentially provides an estimate forCFL EB that is only a function of the geometry of the element. For practical applications, we require the19ight-hand-side of Eq. (63) to be as large as possible so that larger time steps can be taken. This can beachieved by solving a convex optimization problem:maximize (cid:18) min (cid:26) θ k,q ζ k,q (cid:27)(cid:19) , (64)subject to Eq. (61)where θ k,q , ζ k,q are properties of the geometry alone. This problem can be solved for each individual elementas a pre-processing step prior to the simulation. Another way to interpret the expression is to identity alength scale from the right-hand-side of Eq. (63), for which the CFL EB number can be explicitly defined.For this, L e = min V e / |J ∂k ( g k,q ) | is used as a characteristic length for Ω e . Hence,min (cid:26) θ k,q ζ k,q (cid:27) ≥ L e min (cid:26) θ k,q w k,q (cid:27) and an alternative expression to Eq. (64) ismaximize min (cid:26) θ k,q w k,q (cid:27) , (65)subject to Eq. (61) , where the optimal solution is the value of CFL EB . With this, the CFL-constraint can be written as∆ tλL e ≤
12 CFL EB , (66)which is used in the following numerical tests. The factor of 1/2 is a consequence of the Riemann fluxformulation. For some of the most relevant element types with regular shapes, the value of CFL EB has beencalculated and listed in Table 1 for different polynomial orders. In practice, we found that the bound inEq. (66) leads to a conservative estimate for the time step. Considering the computation of efficiency andthe constraint for the linear stability by [21], this condition is relaxed and we consider 0 . EB for thefollowing numerical experiments.
6. Evaluation of entropy bound
In this section, we propose an approach for evaluating the entropy bound s e ( t ) to answer the fourthimplementation problem listed in Sec. 4.2. Obviously, the most accurate way for evaluating a lower boundof entropy is to use Newton’s method. However, this approach can significantly impair the efficiency, sincesearching the minimum on a multi-dimensional high-order element is intractable in terms of computationalcost. To overcome this issue, we propose the following two approaches: • User-defined global bound . The first strategy is to let the user specify a global entropy bound, whichis then kept constant and used everywhere in the computational domain. Although this approach is20 able 1: Summary of quadrature orders and optimal CFL numbers for different types of elements. Quadrature rule (QR) ap-plied: Line, Quadrilateral and Brick: tensor-product Gauss-Legendre; Triangle: Dunavant [22]; Tetrahedron: Zhang, et al. [23].(note that Dunavant’s triangle rule includes negative weights for 3rd-and 7th-order quadrature, therefore, only quadrature ruleswith positive weights are used with one extra order.)
Element Order QR on ∂ Ω e QR on Ω e CFL EB (1, 0)(0, 0) p = 1 / 3 0.5 p = 2 / 5 0.167 p = 3 / 7 0.123 p = 4 / 9 0.073 (0, 1)(0, 0) (1, 1)(1, 0) p = 1 3 3 0.25 p = 2 5 5 0.083 p = 3 7 7 0.062 p = 4 9 9 0.036 (0, 1)(0, 0) (1, 0) p = 1 3 4 0.135 p = 2 5 5 0.067 p = 3 7 8 0.058 p = 4 9 9 0.033 (1,1,1)(0,1,1)(1,0,1)(0,0,1) (1,1,0)(0,1,0)(1,0,0)(0,0,0) p = 1 3 3 0.167 p = 2 5 5 0.056 p = 3 7 7 0.041 p = 4 9 9 0.024 (1,0,0)(0,1,0)(0,0,0)(0,0,1) p = 1 4 3 0.066 p = 2 5 5 0.035 p = 3 8 7 0.015 p = 4 9 9 0.013simple and robust, it is not optimal. It is suitable for certain problems with a well-defined entropybound. As example, for a supersonic flow over an airfoil, the free stream entropy can be used toimpose this bound. However, for more complex cases with multiple discontinuities that include severalentropy jumps, such a constant bound is not able to enforce the constraint for all elements. Note thatthis approach recovers the positivity constraint of Zhang and Shu [9, 10] in the limit of s e ( t ) → −∞ . • Estimate of local entropy bound . This strategy imposes an entropy bound for each element anddynamically updates s e ( t ) during the simulation. Instead of relying on a sophisticated search algorithm, s e ( t ) can be approximately evaluated by reusing available information on quadrature points. According21o the definition of s e ( t ), Eq. (30), we also need to consider the set of quadrature points on ∂ Ω − e ,denoted as D − . Therefore, the estimation for s e ( t ) is obtained according to the following formulation, E s e ( t ) = min min x ∈D − s ( U ( x )) , s m − min x ∈D ,x (cid:54) = x m {| x m − x |}| x m − x n | (cid:0) s n − s m (cid:1) (67)where we introduce x m and x n to denote the locations of the minimum and maximum entropy values,respectively, s m = s ( U e ( x m )) = min x ∈D s ( U e ( x )) ,s n = s ( U e ( x n )) = max x ∈D s ( U e ( x )) . Although this estimate is simple and inexpensive, one has to realize that any extrapolation in thevicinity of discontinuities becomes dangerous due to the spurious behavior of the sub-cell solution.However, it can be resolved by referring to the entropy bounds around Ω e at the last time step, s e ( t ) = max (cid:26) E s e ( t ) , min k ∈N e ∪{ e } s k ( t − ∆ t ) (cid:27) , (68)where N e refers to the set of the indices of all neighbor elements of Ω e that share a common edge.For practical tests, we found that the above strategy can be applied in a combined way. Specifically, Eq.(67) is used for initializing the simulation, and Eq. (68) is then applied during the subsequent simulation.22 . Algorithmic implementation Algorithm 1 provides a description of the implementation details of the EBDG scheme.
Algorithm 1:
Implementation of EBDG scheme.
Pre-computation of CFL condition : For each element, solve Eq. (64); alternatively, take CFL EB from Table 1 and compute L e (recommended for simplicity) Initialization : Initialize solution vector U ( x,
0) = U while t ≤ t end dofor each element do Estimate entropy bound s e ( t ) according to Eqs. (67) and (68)Find λ and estimate time step size ∆ t according to Eq. (63) end Find the minimum permissible time step ∆ t min over all elements; for each stage k of a Runge-Kutta integration scheme dofor each element do step 1: Update solution vector U k +1 = U k + ∆ t min R k ( R refers to the residual)step 2: Apply L on U k +1 with s e ( t ) according to Eqs. (42) and (44) endend Advance time t = t + ∆ t min end8. Results and numerical test cases In the following, EBDG is applied to a series of test cases to demonstrate the performance of this method.We begin by considering one-dimensional configurations to confirm the high-order accuracy and essentialconvergence properties. This is followed by two- and three-dimensional cases with specific emphasis onapplications to unstructured meshes and general curved elements.
The first case considers a one-dimensional periodic domain x ∈ [0 ,
1] with smooth initial conditions: ρ ( x,
0) = 1 + 0 . πx ) ,u ( x,
0) = 1 ,p ( x,
0) = 1 . The accuracy is examined by considering different spatial resolutions and polynomial orders. For eachpolynomial order, the CFL number is assigned to 0.8CFL EB , in which CFL EB is taken from Table 1.23nitially, s is set to 0.874, corresponding to the minimum entropy value of the initial condition. TheSSPRK33 time-integration scheme [24] is used, and the convergence rate is given in Table 2. Although theEBDG-scheme remains stable, it can be seen that the solutions do not reach the optimal rates for DGP3and DGP4. The reason for this is that the stability region is not sufficiently large for both cases. Todemonstrate this, we switch the time-integration scheme to a standard RK45. As can be seen from Table 3,the optimal convergence rates for all cases are achieved, demonstrating that the optimal convergence forsmooth solutions is preserved by the EBDG-scheme. In the following, the standard RK45 is used for allother cases. h DGP1 DGP2 DGP3 DGP4 L -error rate L -error rate L -error rate L -error rate1/10 3.074e-3 - 1.274e-4 - 4.716e-6 - 2.036e-7 -1/20 6.508e-4 2.240 1.513e-5 3.073 3.073e-7 3.940 1.980e-8 3.3621/40 1.535e-4 2.084 1.891e-6 3.000 2.182e-8 3.816 2.454e-9 3.0131/80 3.775e-5 2.024 2.364e-7 3.000 1.880e-9 3.537 3.130e-10 2.9711/160 9.398e-6 2.006 2.955e-8 3.000 2.001e-10 3.232 3.924e-11 2.9951/320 2.347e-6 2.002 3.694e-9 3.000 2.401e-11 3.059 4.922e-12 2.995 Table 2: Convergence test of 1D advection with SSPRK33, showing degradation of convergence order for DGP3 and DGP4(here we use density to evaluate the error). h DGP1 DGP2 DGP3 DGP4 L -error rate L -error rate L -error rate L -error rate1/10 3.494e-3 - 2.140e-4 - 4.650e-6 - 1.438e-7 -1/20 7.231e-4 2.273 1.513e-5 3.823 2.920e-7 3.993 4.517e-9 4.9921/40 1.630e-4 2.150 1.891e-6 3.000 1.826e-8 3.999 1.419e-10 4.9921/80 3.790e-5 2.105 2.364e-7 3.000 1.141e-9 4.000 4.444e-12 4.9971/160 9.398e-6 2.012 2.955e-8 3.000 7.134e-11 4.000 1.497e-13 4.8921/320 2.347e-6 2.002 3.694e-9 3.000 4.463e-12 3.999 8.930e-14 7.453e-1 Table 3: Convergence test of 1D advection with standard RK45 (here we use density to evaluate the error).
A moving shock-wave in a one-dimensional domain is considered as a test-case for evaluating the robust-ness and performance of EBDG for shock-capturing. A domain with x ∈ [ − . , .
1] is considered, in which24he initial shock front is located at x = 0. The domain is initialized in x < ρ = 1 . ,u = 0 ,p = 1 . Shocks are specified with different Mach numbers (Ma = u s /c ), and Ma = { , , } are considered in thiscase. For all cases considered, the initial value for the entropy, s , is set to a value of 0.620, correspondingof the minimum value in the initial condition. The simulation ends when the exact solution of the shockfront reaches the location at x = 1. Results are illustrated in Fig. 3, showing that the entropy bounded-ness guarantees the robustness and consistent performance over a wide range of shock strengths. Entropybounding ( ε (cid:54) = 0) is only activated in elements that are occupied by flow discontinuities. Compared tothe positivity-preserving method, entropy bounding entirely avoids unphysical undershoots in pressure, andprovides an improved suppression of oscillations in the post-shock region. Compared to limiting, the entropybounding shows better robustness in describing shocks at different conditions, introducing lower dissipationin the vicinity of discontinuities. In this section, we verify the convergence order of the EBDG-scheme for high-order curved elements byconsidering a two-dimensional flow over a round cylinder. The radius of the cylinder is R = 1 and thefar-field boundary is a concentric circle with R = 20. The condition in the free stream is given as: ρ ∞ = 1 . ,u ∞ = 5 . ,v ∞ = 0 . ,p ∞ = 1 . The corresponding Mach number is 0.38 and characteristic boundary conditions are imposed at the far-field. The entire domain is initialized with free-stream conditions and s = 0 . EB number from Table 2 for the corresponding shape and polynomial order, multiplied by a factor of 0.8.A main issue in these simulations is the occurrence of numerical instabilities that are initiated at theleading edge of the cylinder. As a result of this instability, DGP2 and DGP3 without any entropy-boundingdiverge (the code blows up) after few iterations. Previously, limiters have been used in this case for sta-bilizing the transient solutions [8]. However, for high-order polynomials, it is difficult to develop limiters25 p ExactEBPPlimiter+PP ε x (a) Ma = 2 , h = 1 / p ε x (b) Ma = 5 , h = 1 / p ε x (c) Ma = 100 , h = 1 / p ε x (d) Ma = 2 , h = 1 / p ε x (e) Ma = 5 , h = 1 / p ε x (f) Ma = 100 , h = 1 / − entropybounding; PP − positivity preserving [9]; Limiter+PP − WENO limiter [7]) with positivity preserving [9]. to achieve the optimal convergence rate without a nontrivial implementation. In contrast, EBDG pro-vides a considerably simpler implementation for enabling high-order simulations for such complex geometricconfigurations.Comparisons of the computational meshes, simulation results, and convergence properties are presentedin Figs. 4 and 5. It is evident that the solution is improved by increasing the mesh resolution. Theconvergence history of the residual, provided in the last row of both figures, shows that entropy bounding ismostly activated during the start-up phase of the simulation to suppress numerical oscillations and ensurestability. It is interesting to note that the number of elements that require bounding is restricted to theregion near the stagnation point upstream of the cylinder, and is confined to less than 8% of the totalnumber of elements. As the solution converges to the steady-state condition, entropy-bounding remainsdeactivated, retaining the high-order accuracy. Since the solution is smooth the physical entropy productionis zero. Therefore, the convergence rates are measured in terms of entropy error using the discrete L -norm.26esh DGP1 DGP2 DGP3 L -error rate L -error rate L -error rateQuadrilateral ElementsLevel 1 7.272e-2 - 1.694e-2 - 3.816e-3 -Level 2 1.318e-2 2.464 7.219e-4 4.552 1.827e-4 4.384Level 3 2.441e-3 2.433 6.029e-5 3.582 1.036e-5 4.141Triangular ElementsLevel 1 1.137e-1 - 2.590e-2 - 4.086e-3 -Level 2 1.865e-2 2.608 8.899e-4 4.863 1.291e-4 4.984Level 3 3.391e-3 2.459 7.222e-5 3.623 6.939e-6 4.217 Table 4: Comparisons of convergence rate for 2D flow over a cylinder (here we use entropy to evaluate the error).
A comparison of the convergence rates are presented in Table 4, confirming that the optimal convergencerate is preserved even for complex geometries with curved elements.
This test case is designed to assess the performance of EBDG for simulation of flows with strong shocksand wave structures. The numerical setup follows this of Woodward and Colella [25], representing a Mach10 shock over a 30 ◦ -wedge. All quantities are non-dimensionalized, and the computational domain is [0 , × [0 , L e = h = 0 .
02) and a mesh with triangular elements ( L e ≈ . s is set to 0.620. The CFL number is prescribed from Table2 using a safety factor of 0.8.Simulation results for density contours at time t = 0 .
25 are shown in Fig. 6. The proposed EBDG-methodcaptures all wave-features, and it is found that without enforcing the entropy constraint the solution divergesin the first iteration for these strong shock conditions. For comparison, a reference solution obtained usinga fifth-order WENO-scheme is shown in Fig. 6(e), and results from a DGP2-simulation using a WENO-limiter [7] are presented in Fig. 6(f). Comparisons between EBDGP1 and EBDGP2 results show the benefitof the high-order scheme in providing improved representations of the shock-wave structure. At the samedegrees of freedom, the DGP2-solution provides comparable predictions to that of the fifth-order WENOscheme, except for the small oscillations that cannot be removed by the linear scaling procedure. Comparedto the DG-simulation with WENO-limiter (Fig. 6(f)), EBDG effectively avoids introducing excessive nu-merical dissipation since the solution is only entropy-constrained in regions in which the entropy conditionis violated. 27esh Level 1 Mesh Level 2 Mesh Level 3 M e s h x y
10 5 0 5 10 10 50510 x y
10 5 0 5 10 10 50510 x y
10 5 0 5 10 10 50510 D G P s o l u t i o n s x y Mach: 0.05 0.2 0.35 0.5 0.65 0.8 x y x y C o n v e r g e n ce −8 −6 −4 −2 R e s i du a l I terations
DGP1DGP2DGP3 N o . o f ce ll s w i t h ε = −8 −6 −4 −2 R e s i du a l I terations N o . o f ce ll s w i t h ε = −8 −6 −4 −2 R e s i du a l I terations N o . o f ce ll s w i t h ε = Figure 4: EBDG-solution of flow over a cylinder on curved quadrilateral meshes with three different refinement levels; top:computational mesh in near-field of the cylinder; middle: Mach number; bottom: convergence history and activation of entropybounding as a function of iteration.
This test case extends the evaluation of the EBDG-method to three-dimensional configurations withcomplex geometries. Currently, robust approaches for capturing strong shocks in three-dimensional curvedelements are still subject to investigation. This test case considers a flow at a Mach number of 6.8 overa sphere. The radius of the sphere is R = 1. Due to the geometric symmetry, the computational domainconsiders only an eighth section of the domain, and it extends to 3 R in radial direction. Symmetry boundaryconditions are imposed at the planes y = 0 and z = 0, and outflow boundary conditions are prescribed at28esh Level 1 Mesh Level 2 Mesh Level 3 M e s h x y
10 5 0 5 10 10 50510 x y
10 5 0 5 10 10 50510 x y
10 5 0 5 10 10 50510 D G P s o l u t i o n s x y Mach: 0.05 0.2 0.35 0.5 0.65 0.8 x y x y C o n v e r g e n ce −8 −6 −4 −2 R e s i du a l I terations
DGP1DGP2DGP3 N o . o f ce ll s w i t h ε = −8 −6 −4 −2 R e s i du a l I terations N o . o f ce ll s w i t h ε = −8 −6 −4 −2 R e s i du a l I terations N o . o f ce ll s w i t h ε = Figure 5: DG-solution of flow over a cylinder on curved triangular meshes with three different refinement levels; top: compu-tational mesh in the near-field of the cylinder; middle: Mach number; bottom: convergence history and activation of entropybounding as a function of iteration. x = 0 . Normal velocity inflow is prescribed at the outer shell with the following specification: ρ ∞ = 1 . ,u ∞ = − . ,v ∞ = 0 . ,w ∞ = 0 . ,p ∞ = 1 . y (a) EBDGP1, quadrilateral mesh with h = 0 . . x y (b) EBDGP1, triangular mesh with h = 0 . . (c) EBDGP2, quadrilateral mesh with h = 0 . . (d) EBDGP2, triangular mesh with h = 0 . . x y (e) WENO5, quadrilateral mesh with h = 0 . . x y (f) DGP2 with WENO limiter [7], quadrilateral mesh with h = 0 . x y : (g) Instantaneous snapshot of ε for DGP2 on quadrilateral mesh.Figure 6: (Color online) Simulation results of double Mach reflection over a 30 o -wedge. x = 0is partitioned using 12 elements. DGP2 is applied for this case and the CFL number is 0.8CFL EB withCFL EB = 0 . ε can be utilized as a indicator for local mesh refinement. We sample theelements with non-zero ε values over few iterations, and then locally refine these elements. Results usingone and two levels of refinement are shown in Figs. 7(b) and 7(c), respectively. This direct comparisonshows that the shock profiles become sharper with increasing resolution. Since the bounding parameter issharp, the mesh-refinement is confined to a narrow region in the vicinity of the shock. The flow-field solutionbehind the shock is smooth, and no entropy bounding is applied in this region. To provide a quantitativeanalysis, simulation results from the EBDG-method are compared against measurements by Billig [26] inFig. 7(d), showing good agreement between experiments and computations.
9. Conclusions
A regularization technique for the discontinuous Galerkin scheme was developed using the entropy prin-ciple. Motivated by the FV entropy solution, the high-order DG-scheme is stabilized by constraining thesolution to obey the entropy condition. The implementation of the resulting entropy-bounding discontin-uous Galerkin scheme relies on two key components, namely a limiting operator and a CFL-constraint.These essential components were derived by considering first a one-dimensional setting and the subsequentextension to multi-dimensional configurations with arbitrary and curved elements. Specifically, utilizing theinterpolation basis we were able to extend the entropy bounding (also including positivity preserving) toarbitrarily shaped elements independent of specific quadrature rules. The bounding procedure is obtainedfrom algebraic operations, resulting in a computationally efficient and simple implementation. A sufficientCFL-condition was rigorously derived and proofed to ensure that the entropy constraint can be enforcedon different types and orders of elements. By considering different configurations, numerical tests wereconducted to examine accuracy and stability of the entropy-bounding DG-scheme. These test cases confirmthe efficacy in regularizing solutions in the vicinity of discontinuities, generated either by true flow physicsor during the transient solution update. The added benefit of the entropy bounding method is its utilizationas a refinement indicator.Since the herein proposed entropy bounding scheme relies on a linear scaling operator, it is not capableto remove shock-triggered oscillations of smaller magnitude, although it stabilizes the solution and preventsthe solver from diverging. As a final remark, the derivation that was presented in this study is general31 YZ (a) Baseline mesh. X YZ (b) Local refinement at level one.
X YZ (c) Local refinement at level two.
Mach76543210 (d) Comparison with measurement in [26].Figure 7: (Color online) Simulations of Mach 6.8 flow over a sphere showing the (cid:15) profile ( y = 0 plane) and Mach-numberdistribution ( z = 0 plane) on (a) baseline mesh, and simulation results with local refinement with (b) one refinement level and(c) two refinement levels. Comparisons of the shock location with measurements by Billig [26] are shown in (d). and extendable to other discontinuous schemes with sub-cell solution representations, such as spectral finitevolume schemes [27] and the flux reconstruction scheme [28]. Therefore, entropy-bounding, as an idea, hasthe potential to improve the robustness of shock-capturing for these emerging high-order numerical methods.32 cknowledgment Financial support through NSF with Award No. CBET-0844587 is gratefully acknowledged. This workused the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NationalScience Foundation grant number ASC130004. Helpful discussions with Yee Chee See on the mathematicalanalysis are appreciated.
A. Combination Rule
In this section, we derive an estimate for the maximum characteristic speed for convex state solution.For this, we consider a state vector U of Eq. (2 a ), which is written in the following form: U = (cid:88) k β k U k , (69)where β k > (cid:80) k β k = 1. The maximum characteristic speed of U is: ν ( U ) = | u ( U ) | + c ( U ) = | u ( U ) | + (cid:115) γ ( γ − (cid:18) e ( U ) − | u ( U ) | (cid:19) , in which u and E can be calculated according to Eq. (69) as u ( U ) = (cid:88) k α k u ( U k ) , e ( U ) = (cid:88) k α k e ( U k ) , and α k = β k ρ ( U k ) (cid:80) k β k ρ ( U k )is a new set of coefficients that is introduced to convert from conservative to primitive variables. Furthermore,because of γ ( γ − e ( U ) = (cid:88) k α k (cid:18) c ( U k ) + γ ( γ − | u ( U k ) | (cid:19) , | u ( U ) | = (cid:115) | (cid:88) k α k u ( U k ) | , ≤ (cid:115)(cid:88) k α k | u ( U k ) | ,
33e obtain ν ( U ) = | u ( U ) | + (cid:118)(cid:117)(cid:117)(cid:116)(cid:88) k α k c ( U k ) + γ ( γ − (cid:32)(cid:88) k α k | u ( U k ) | − | u ( U ) | (cid:33) , ≤ (cid:115)(cid:88) k α k | u ( U k ) | + (cid:115)(cid:88) k α k c ( U k ) + γ ( γ − (cid:88) k α k | u ( U k ) | , ≤ (cid:115) (cid:88) k α k c ( U k ) + (2 + γ ( γ − (cid:88) k α k | u ( U k ) | , ≤ (cid:112) γ ( γ − (cid:115)(cid:88) k α k ( c ( U k ) + | u ( U k ) | ) , ≤ (cid:112) γ ( γ − (cid:115)(cid:88) k α k ( c ( U k ) + | u ( U k ) | ) , ≤ (cid:112) γ ( γ −
1) max k { c ( U k ) + | u ( U k ) |} . We used this estimation for preselecting the dissipation coefficient λ ∗ used in Lemma 4. B. Formulation of J ∂k For a two-dimensional configuration, the curve of an edge is parameterized by g k ∈ R , and the surface Ja-cobian is written as J ∂k = (cid:104) ∂x ∂g k , ∂x ∂g k (cid:105) ; for a three-dimensional configuration, an edge surface is parameterizedby g k = (cid:104) g (1) k , g (2) k (cid:105) T ∈ R , and the Jacobian can be written as J ∂k = (cid:20) ∂x ∂g (1) k , ∂x ∂g (1) k , ∂x ∂g (1) k (cid:21) × (cid:20) ∂x ∂g (2) k , ∂x ∂g (2) k , ∂x ∂g (2) k (cid:21) . References [1] B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method forconservation laws III: One dimensional systems.
J. Comp. Phys. , 84:90–113, 1989.[2] B. Cockburn and C.-W. Shu. The Runge-Kutta discontinuous Galerkin method for conservation laws V: Multidimensionalsystems.
J. Comp. Phys. , 141:199–224, 1998.[3] R. Biswas, K. Devine, and J. E. Flaherty. Parallel adaptive finite element methods for conservation laws.
Appl. Numer.Math. , 14:255–284, 1994.[4] L. Krivodonova. Limiters for high-order discontinuous Galerkin methods.
J. Comp. Phys. , 226:879–896, 2007.[5] J. Qiu and C.-W. Shu. Runge-Kutta discontinuous Galerkin method using WENO limiters.
SIAM J. Sci. Comput. ,26:907–929, 2005.[6] J. Zhu, J. Qiu, C.-W. Shu, and M. Dumbser. Runge-Kutta discontinuous Galerkin method using WENO limiters II:Unstructured meshes.
J. Comp. Phys. , 227:4330–4353, 2008.[7] X. Zhong and C.-W. Shu. A simple weighted essentially nonoscillatory limiter for Runge-Kutta discontinuous Galerkinmethods.
J. Comp. Phys. , 232:397–415, 2013.[8] H. Luo, J. D. Baum, and R. L¨ohner. A Hermite WENO-based limiter for discontinuous Galerkin method on unstructuredgrids.
J. Comp. Phys. , 225:686–713, 2007.[9] X. Zhang and C.-W. Shu. On positivity-preserving high order discontinuous Galerkin schemes for compressible Eulerequations on rectangular meshes.
J. Comp. Phys. , 229:8918–8934, 2010.
10] X. Zhang and C.-W. Shu. Positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equationswith source terms.
J. Comp. Phys. , 230:1238–1248, 2011.[11] X. Zhang and C.-W. Shu. A minimum entropy principle of high order schemes for gas dynamics equations.
Numer. Math. ,121:545–563, 2012.[12] Y. Lv and M. Ihme. Computational analysis of re-ignition and re-initiation mechanisms of quenched detonation wavesbehind a backward facing step.
Proc. Combust. Inst. , 2014. available online.[13] Y. Lv and M. Ihme. Discontinuous Galerkin method for multicomponent chemically reacting flows and combustion.
J.Comp. Phys. , 270:105–137, 2014.[14] P. D. Lax. Shock waves and entropy. In E.H. Zarantonello, editor,
Contributions to Nonlinear Functional Analysis , pages603–634. Academic Press, New York, London, 1971.[15] E. Tadmor. A minimum entropy principle in the gas dynamics equations.
Appl. Numer. Math. , 2:211–219, 1986.[16] S. L. Sobolev and V. Vaskevich.
The theory of cubature formulas . Springer, 1997.[17] B. Perthame and C.-W. Shu. On positivity preserving finite volume schemes for Euler equations.
Numer. Math. , 73:119–130, 1996.[18] E. Tadmor. Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependant problem.
Acta Numer. , 12:451–512, 2003.[19] B. Khobalatte and B. Perthame. Maximum principle on the entropy and second-order kinetic schemes.
Math. Comp. ,62:119–131, 1994.[20] Z. J. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert, H. T.Huynh, N. Kroll, G. May, P. -O. Persson, B. van Leer, and M. Visbal. High-order CFD methods: current status andperspective.
Int. J. Numer. Meth. Eng. , 72:811–845, 2013.[21] B. Cockburn and C.-W. Shu. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems.
J. Sci.Comput. , 16(3):173–261, 2001.[22] D. Dunavant. High degree efficient symmetrical Gaussian quadrature rules for the triangle.
Int. J. Numer. Meth. Eng. ,21:1129–1148, 1985.[23] L. Zhang, T. Cui, and H Liu. A set of symmetric quadrature rules on triangles and tetrahedra.
J. Comput. Math. ,21:89–96, 2009.[24] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong stability-preserving high-order time discretization methods.
SIAM Review ,43:89–112, 2001.[25] P. R. Woodward and P. Colella. The numerical simulation of two-dimensional fluid flow with strong shocks.
J. Comp.Phys. , 54:115–173, 1984.[26] F. S. Billig. Shock-wave shapes around spherical-and cylindrical-nosed bodies.
J. Spacecraft Rockets , 4:822–823, 1967.[27] Z. J. Wang and Y. Liu. Spectral (finite) volume method for conservation laws on unstructured grids V: Extension tothree-dimensional systems.
J. Comp. Phys. , 212:454–472, 2006.[28] H. T. Huynh. A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods. In , Miami, FL, 2007. AIAA 2007-4079., Miami, FL, 2007. AIAA 2007-4079.