Discretely entropy stable weight-adjusted discontinuous Galerkin methods on curvilinear meshes
DDISCRETELY ENTROPY STABLE WEIGHT-ADJUSTED DISCONTINUOUSGALERKIN METHODS ON CURVILINEAR MESHES
JESSE CHAN, LUCAS C. WILCOX
Abstract.
We construct entropy conservative and entropy stable high order accurate discontinuous Galerkin(DG) discretizations for time-dependent nonlinear hyperbolic conservation laws on curvilinear meshes. Theresulting schemes preserve a semi-discrete quadrature approximation of a continuous global entropy inequal-ity. The proof requires the satisfaction of a discrete geometric conservation law, which we enforce throughan appropriate polynomial approximation. We extend the construction of entropy conservative and entropystable DG schemes to the case when high order accurate curvilinear mass matrices are approximated us-ing low-storage weight-adjusted approximations, and describe how to retain global conservation propertiesunder such an approximation. The theoretical results are verified through numerical experiments for thecompressible Euler equations on triangular and tetrahedral meshes. Introduction
High order discontinuous Galerkin (DG) methods are attractive for the simulation of time-dependent wavepropagation due to their low numerical dispersion and dissipation [5, 6] and ability to handle unstructuredmeshes and complex geometries. These same properties make them attractive for the resolution of transientwaves and vortices in compressible flow [7]. However, whereas the construction of stable DG methods forwave propagation is relatively well-established, it is not possible to extend the same formulations directlyfrom linear wave problems to the nonlinear conservation laws which govern compressible fluid flow.The low numerical dissipation of high order DG methods combined with the lack of inherently stableformulations for nonlinear conservation laws has given high order discretizations the reputation of beingnon-robust and highly sensitive to instabilities and under-resolved features [7]. This instability is addressedin practice by adding additional stabilization through limiting, filtering, or artificial viscosity [8, 9, 10, 11].However, these approaches are typically ad-hoc, and do not guarantee the stability of the resulting scheme.Moreover, high order accuracy can be lost if the stabilization is too strong.The lack of inherently stable formulations was addressed for high order nodal DG methods on quadrilateraland hexahedral meshes in [12, 13]. The resulting schemes ensure that that the numerical solution satisfiesa semi-discrete version of an entropy inequality, independently of discrete effects such as under-integration.This results in significantly more robust high order simulations, where the numerical solution does notblow up even in the presence of under-resolved features such as shock discontinuities or turbulence. Highorder entropy stable schemes have since been extended to staggered grid and non-conforming [14, 15] tensorproduct elements, as well as to simplicial meshes [16, 1, 17].Entropy stable methods have largely relied on a finite difference summation-by-parts (SBP) framework.The summation-by-parts property holds for spectral element DG methods (DG-SEM), which assume apolynomial basis which collocates the solution at Gauss–Legendre–Lobatto (GLL) quadrature nodes ontensor product elements. The collocation approach requires the number of quadrature nodes to be identicalto the number of basis functions. Because analogous quadrature rules on the triangle and tetrahedron mustcontain more nodes than the dimension of the underlying polynomial space [18, 19], GLL-like collocationcannot be replicated on simplicial elements. It is still possible to construct entropy stable schemes onsimplicial elements within the SBP framework [16, 17]. However, the resulting SBP operators are not“modal”, in the sense that the matrices are not associated with an underlying basis or approximation space.Entropy stable schemes were extended to modal discretizations on affine meshes in [1], allowing for theuse of arbitrary basis functions and over-integrated quadrature rules. In this work, we show how to extendthe construction of modal high order entropy stable schemes to curved meshes. As noted in [17], this requiresthe use of a split formulation for the geometric terms involved in differentiation, as well as the satisfaction of a r X i v : . [ m a t h . NA ] J un discrete geometric conservation law (GCL) [20, 2]. We also show how to extend entropy stable schemes toaccomodate weight-adjusted mass matrices, which provide low storage approximations of inverse weightedmass matrices appearing for curved meshes [4]. These weight-adjusted approximations can be applied moreefficiently than inverse weighted mass matrices on many-core architectures such as Graphics Processing Units(GPUs) [21]. However, additional steps are required to ensure entropy stability and the conservation of mass,momentum, and energy under weight-adjusted mass matrices.The paper is organized as follows: Section 2 reviews the derivation of entropy inequalities for systemsof nonlinear conservation laws, and describes why this does not hold under a high order DG discretization.Section 3 describes the construction of entropy stable DG methods for curved meshes using weighted massmatrices, and Section 4 describes how to extend this to the case of weight-adjusted mass matrices. Section 5describes how to enforce the geometric conservation law using a modification of the approach described in[2], and Section 6 concludes by presenting two and three-dimensional numerical experiments which verifythe accuracy and stability of the presented schemes.2. Systems of nonlinear conservation laws
This work addresses high order accurate schemes for the following system of n nonlinear conservation lawsin d dimensions(1) ∂ u ∂t + d (cid:88) j =1 ∂ f j ( u ) ∂x j = 0 , u : R d × [0 , ∞ ) → R n , f j : R n → R n , where u ( x , t ) denotes the conservative variables for this system. We are interested in nonlinear conservationlaws for which an entropy function U ( u ) exists, where U ( u ) is convex with respect to the conservativevariables u . If this function exists, then it is possible to define entropy variables v ( u ) = ∂U∂ u . These functionssymmetrize the system of nonlinear conservation laws (1) [22].It can be shown (see, for example, [23]) that symmetrization is equivalent to the existence of entropy fluxfunctions F j ( u ) and entropy potentials ψ j such that v T ∂ f j ∂ u = ∂F j ( u ) ∂ u T , ψ j ( v ) = v T f j ( u ( v )) − F j ( u ( v )) , ψ (cid:48) j ( v ) = f j ( u ( v )) . Smooth solutions of (1) can be shown to satisfy a conservation of entropy by multiplying (1) by v ( u ). Usingthe definition of the entropy variables, entropy flux, and the chain rule yields(2) v T ∂ f j ( u ) ∂x j = ∂U ( u ) ∂ u T ∂ f j ( u ) ∂ u ∂ u ∂x j = ∂F j ( u ) ∂x j , ∂U ( u ) ∂t + d (cid:88) j =1 ∂F j ( u ) ∂x j = 0 . Let Ω ⊂ R d be a closed domain with boundary ∂ Ω. Integrating over Ω and using Gauss’ theorem on thespatial derivative yields(3) (cid:90) Ω ∂U ( u ) ∂t d x + (cid:90) ∂ Ω d (cid:88) j =1 (cid:16) v ( u ) T f j ( u ) − ψ j ( v ( u )) (cid:17) n j d x = 0 , where n = ( n , . . . , n d ) T denotes the unit outward normal vector on ∂ Ω.General solutions (including non-smooth solutions such as shocks) satisfy an entropy inequality (4) (cid:90) Ω ∂U ( u ) ∂t d x + (cid:90) ∂ Ω d (cid:88) j =1 (cid:16) v ( u ) T f j ( u ) − ψ j ( v ( u )) (cid:17) n j d x ≤ , which results from considering solutions of an appropriate viscous form of the equations (1) and taking thelimit as viscosity vanishes. In this work, schemes which satisfy a discrete form of (4) will be constructed byfirst enforcing a discrete version of entropy conservation (3), then adding an appropriate numerical dissipationwhich will enforce the entropy inequality (4).2.1. Standard DG formulations for nonlinear conservation laws.
We begin by reviewing the con-struction of standard high order accurate DG formulations for (1). .1.1. Mathematical notation.
Let the domain Ω ⊂ R d be decomposed into elements (subdomains) D k , andlet (cid:98) D denote a d -dimensional reference element with boundary ∂ (cid:98) D . Let (cid:98) x = { (cid:98) x , . . . , (cid:98) x d } denote coordinateson (cid:98) D , and let (cid:98) n i denote and the i th component of the unit normal vector on ∂ (cid:98) D . We assume that (cid:98) n i isconstant; i.e., that the faces of the reference element are planar (this assumption holds for all commonlyused reference elements [24]).We will assume that each physical element D k is the image of (cid:98) D under some smoothly differentiablemapping Φ k ( (cid:98) x ) such that x = Φ k ( (cid:98) x ) , x ∈ D k . This also implies that integrals over physical elements can be mapped back to the reference element as follows (cid:90) D k u d x = (cid:90) (cid:98) D uJ k d (cid:98) x , where J k denotes the determinant of the Jacobian of Φ k . Integrals over physical faces of D k can similarlybe mapped back to reference faces.We define an approximation space using degree N polynomials on the reference element. For example,on a d -dimensional reference simplex, the natural polynomial space are total degree N polynomials P N (cid:16) (cid:98) D (cid:17) = (cid:40)(cid:98) x i . . . (cid:98) x i d d , (cid:98) x ∈ (cid:98) D, ≤ d (cid:88) k =1 i k ≤ N (cid:41) . Other element types possess different natural polynomial spaces [24], but typically contain the space of totaldegree N polynomials. This work is directly applicable to other elements and spaces as well. We denote thedimension of the approximation space P N as N p = dim (cid:16) P N (cid:16) (cid:98) D (cid:17)(cid:17) . We also define trace spaces for eachface of the reference element. Let (cid:98) f be a face of the reference element (cid:98) D . The trace space over (cid:98) f is definedas the space of traces of functions in P N (cid:16) (cid:98) D (cid:17) P Nf (cid:16) (cid:98) f (cid:17) = (cid:110) u | (cid:98) f , u ∈ P N (cid:16) (cid:98) D (cid:17)(cid:111) , (cid:98) f ∈ ∂ (cid:98) D. We denote the dimension of the trace space as dim (cid:16) P Nf (cid:16) (cid:98) f (cid:17)(cid:17) = N fp .We next define the L norm and inner products over the reference element (cid:98) D and the surface of thereference element ∂ (cid:98) D as( u , v ) (cid:98) D = (cid:90) (cid:98) D u · v d (cid:98) x , (cid:107) u (cid:107) (cid:98) D = ( u , u ) (cid:98) D , (cid:104) u , v (cid:105) ∂ (cid:98) D = (cid:90) ∂ (cid:98) D u · v d (cid:98) x . We also introduce the continuous L projection operator Π N and lifting operator L . For u ∈ L (cid:16) (cid:98) D (cid:17) , the L projection Π N u is defined through(5) (cid:90) (cid:98) D Π N uv d (cid:98) x = (cid:90) (cid:98) D uv d (cid:98) x , ∀ v ∈ P N (cid:16) (cid:98) D (cid:17) . Likewise, for a boundary function u ∈ L (cid:16) ∂ (cid:98) D (cid:17) , the lifting operator L [25, 26] is defined through(6) ( Lu, v ) (cid:98) D = (cid:104) u, v (cid:105) ∂ (cid:98) D , ∀ v ∈ P N (cid:16) (cid:98) D (cid:17) . Finally, we introduce L , L ∞ Sobolev norms and spaces, which will be utilized for error estimates. The L space is defined as the space of functions with finite L norm. The Lebesgue L ∞ norm and the associated L ∞ space over a general domain Ω are (cid:107) u (cid:107) L ∞ (Ω) = inf { C ≥ | u ( x ) | ≤ C ∀ x ∈ Ω } , L ∞ (Ω) = (cid:110) u : Ω → R , (cid:107) u (cid:107) L ∞ (Ω) < ∞ (cid:111) . The L and L ∞ Sobolev norms of degree s are then defined as (cid:107) u (cid:107) W s, (Ω) = (cid:88) | α |≤ s (cid:107) D α u (cid:107) L (Ω) , (cid:107) u (cid:107) W s, ∞ (Ω) = max | α |≤ s (cid:107) D α u (cid:107) L ∞ (Ω) , espectively. Here α = { α , . . . , α d } is a multi-index of order | α | = α + · · · + α d such that D α u = ∂ α ∂x α · · · ∂ α d ∂x α d d u. The Sobolev spaces W s, and W s, ∞ are defined as the spaces of functions with finite L and L ∞ Sobolevnorms of degree s , respectively.2.1.2. Discontinuous Galerkin formulations and the L projection. Discontinuous Galerkin methods havebeen widely applied to systems of nonlinear conservation laws (1) [27, 28, 29]. The development of newdiscontinuous Galerkin methods for nonlinear conservation laws has focused heavily on the choice of numericalflux [30] or the development of slope limiters [9, 31] and artificial viscosity strategies [8, 10, 32]. However,the treatment of the underlying volume discretization remains relatively unchanged between each of theseapproaches.Ignoring terms involving filters, limiters, or artificial viscosity, a semi-discrete “weak” DG formulation for(1) can be given locally over an element D k : find u ∈ (cid:0) P N (cid:0) D k (cid:1) × [0 , ∞ ) (cid:1) n such that (cid:90) D k ∂ u ∂t · v − d (cid:88) j =1 f j ( u ) · ∂ v ∂x i d x + d (cid:88) j =1 (cid:90) ∂D k (cid:0) f ∗ j (cid:0) u + , u (cid:1)(cid:1) · v n j d x = 0 , ∀ v ∈ (cid:0) P N (cid:0) D k (cid:1)(cid:1) n , (7)where the numerical flux f ∗ is a function of the solution u on both D k and neighboring elements.Unfortunately, solutions to (7) do not (in general) obey a discrete version of the entropy inequality (4).Since (4) is a generalized statement of energy stability, the lack of a discrete entropy inequality implies thatthe discrete solution can blow up in finite time. The reason for this is due to the fact that, in practice, theintegrals in (7) are not computed exactly and are instead approximated using polynomially exact quadratures.This is compounded by the fact that the nonlinear flux function f j ( u ) is often rational and impossible tointegrate exactly using polynomial quadratures.While this paper focuses on curved meshes, the inexactness of quadrature leads to the loss of the chainrule and thus the loss of entropy conservation or entropy dissipation (4) even on affine meshes. We canrewrite (7) in a strong form using a discrete quadrature-based L projection. For polynomial approximationspaces on affine meshes, ∂ v ∂x i is polynomial. Then, mapping (7) back to the reference element (cid:98) D and usingthe L projection and (5), we have that (cid:90) D k f j ( u ) · ∂ v ∂x i d x = (cid:90) (cid:98) D Π N f j ( u ) · ∂ v ∂x i J k d x . Thus, integrating by parts (7) recovers a “strong” DG formulation involving the projection operator (cid:90) D k ∂ u ∂t − d (cid:88) j =1 ∂ Π N f j ( u ) ∂x j · v d x + d (cid:88) j =1 (cid:90) ∂D k (cid:0) f ∗ j (cid:0) u + , u (cid:1) − Π N f j ( u ) (cid:1) · v n j d x = 0 , ∀ v ∈ (cid:0) P N (cid:0) D k (cid:1)(cid:1) n . (8)From this, we see that our discrete scheme does not differentiate the nonlinear flux function f j ( u ) exactly,but instead differentiates the projection of Π N f j ( u ) onto polynomials of degree N . Because the L projectionoperator is introduced, the chain rule no longer holds at the discrete level and step (2) of the proof of entropyconservation is no longer valid. Thus, ensuring discrete entropy stability will require a discrete formulation ofthe system of nonlinear conservation laws (1) from which we can prove a discrete entropy inequality withoutrelying on the chain rule.3. Discretely entropy stable DG methods on curved meshes
We will first show how to construct discretely entropy stable high order accurate DG methods on curvi-linear meshes, but will present this using a matrix formulation as opposed to a continuous formulation. Thisis to ensure that the effects of discretization, nonlinear, and quadrature are accounted for in the proof of emi-discrete entropy stability. We first introduce quadrature-based matrices, which we will then use toconstruct discretely entropy stable DG formulations.3.1. Basis and quadrature rules.
We now introduce quadrature-based matrices for the d -dimensionalreference element (cid:98) D , which we will use to construct matrix-vector formulations of DG methods. Assuming u ( (cid:98) x ) ∈ P N (cid:16) (cid:98) D (cid:17) , it can be represented in terms of the vector of coefficients u using some polynomial basis φ i of degree N and dimension N p u ( (cid:98) x ) = N p (cid:88) j =1 u j φ j ( (cid:98) x ) , P N (cid:16) (cid:98) D (cid:17) = span { φ i ( (cid:98) x ) } N p i =1 . We construct quadrature-based matrices based on φ i and appropriate volume and surface quadraturerules. The volume and surface quadrature rules are given by points and positive weights { ( (cid:98) x i , (cid:98) w i ) } N q i =1 and (cid:110) ( (cid:98) x fi , (cid:98) w fi ) (cid:111) N fq i =1 , respectively. We make the following assumptions on the strength of these quadratures: Assumption 1 (Integration by parts under quadrature) . The volume quadrature rule { ( (cid:98) x i , (cid:98) w i ) } N q i =1 is exactfor polynomials of degree N − . Additionally, for any u, v ∈ P N (cid:16) (cid:98) D (cid:17) , integration by parts (cid:18) ∂u∂ (cid:98) x i , v (cid:19) (cid:98) D = (cid:104) u, v (cid:98) n i (cid:105) ∂ (cid:98) D − (cid:18) u, ∂v∂ (cid:98) x i (cid:19) (cid:98) D holds when the volume and surface integrals are approximated using quadrature. Assumption 1 holds, for example, for any surface quadrature rule which is exact for degree 2 N polynomialson the boundary of the reference element ∂ (cid:98) D .3.2. Reference element matrices.
Let W , W f denote diagonal matrices whose entries are volume andsurface quadrature weights, respectively. The surface quadrature weights are given by quadrature weightson reference faces, which are mapped to faces of the reference element. We define the volume and surfacequadrature interpolation matrices V q and V f such that( V q ) ij = φ j ( (cid:98) x i ) , ≤ j ≤ N p , ≤ i ≤ N q , ( V f ) ij = φ j ( (cid:98) x fi ) , ≤ j ≤ N p , ≤ i ≤ N fq , (9)which map coefficients u to evaluations of u at volume and surface quadrature points.Next, let D i denote the differentiation matrix with respect to the i th coordinate, defined implicitly throughthe relations u ( (cid:98) x ) = N p (cid:88) j =1 u j φ j ( (cid:98) x ) , ∂u∂ (cid:98) x i = N p (cid:88) j =1 ( D i u ) j φ j ( (cid:98) x ) . The matrix D i maps basis coefficients of some polynomial u ∈ P N (cid:16) (cid:98) D (cid:17) to coefficients of its i th derivativewith respect to the reference coordinate (cid:98) x , and is sometimes referred to as a “modal” differentiation matrix(with respect to a general non-nodal “modal” basis [19]).Using the volume quadrature interpolation matrix V q , we can compute a quadrature-based mass matrix M by evaluating L inner products of different basis functions using quadrature M = V Tq W V q , M ij = N q (cid:88) k =1 (cid:98) w k φ j ( (cid:98) x k ) φ i ( (cid:98) x k ) ≈ (cid:90) (cid:98) D φ j φ i d (cid:98) x = ( φ j , φ i ) (cid:98) D . The approximation in the formula for the mass matrix becomes an equality if the volume quadrature rule isexact for polynomials of degree 2 N . The mass matrix is symmetric and positive definite under Assumption 1;however, we do not make any distinctions between diagonal and dense (lumped) mass matrices in this work.The mass matrix appears in the discretization of L projection (5) and lift operators (6) using quadrature.The result are quadrature-based L projection and lift operators P q , L q ,(10) P q = M − V Tq W , L q = M − V Tf W f , hich are discretizations of the continuous L projection operator Π N and continuous lift operator L . Thematrix P q maps a function (in terms of its evaluation at quadrature points) to coefficients of the L projectionin the basis φ j ( x ), while the matrix L q “lifts” a function (evaluated at surface quadrature points) from theboundary of an element to coefficients of a basis defined in the interior of the element.Finally, we introduce quadrature-based operators D iN which will be used to construct discretizations ofour nonlinear conservation laws. This operator was introduced in [1] as a “decoupled summation-by-parts”operator(11) D iN = V q D i P q − V q L q diag (cid:16)(cid:98) n i ◦ (cid:98) J f (cid:17) V f P q V q L q diag (cid:16)(cid:98) n i ◦ (cid:98) J f (cid:17) − diag (cid:16)(cid:98) n i ◦ (cid:98) J f (cid:17) V f P q diag (cid:16)(cid:98) n i ◦ (cid:98) J f (cid:17) where (cid:98) n i is the vector containing values of the i th component of the unit normal on the surface of thereference element (cid:98) D , and (cid:98) J f is the vector containing values of the face Jacobian factor (cid:98) J f which result frommapping a face of (cid:98) D to a reference face. Here (cid:98) n i ◦ (cid:98) J f is the Hadamard product (i.e., the entrywise product)of the vectors (cid:98) n i and (cid:98) J f . When combined with projection and lifting matrices, D iN produces a high orderapproximation of non-conservative products. Let f , g denote vectors containing the evaluation of functions f ( x ) , g ( x ) at both volume and surface quadrature points (cid:2) P q L q (cid:3) diag ( f ) D iN g ≈ f ∂g∂ (cid:98) x i . It was shown in [1] that the matrix D iN satisfies several key properties. First, it can be observed that D iN = 0, where is the vector of all ones. Second, D iN satisfies a summation-by-parts property. Let Q iN be the scaling of D iN by the diagonal matrix of volume and surface quadrature weights Q iN = W N D iN , W N = (cid:18) W W f (cid:19) . Then, Q iN satisfies the following discrete analogue of integration by parts(12) Q iN + (cid:0) Q iN (cid:1) T = B iN , B N = (cid:32) W f diag (cid:16)(cid:98) n i ◦ (cid:98) J f (cid:17) (cid:33) . The matrix D iN reduces to polynomial differentiation when applied to polynomials, in the sense that D iN (cid:20) V q V f (cid:21) = (cid:20) V q D i (cid:21) . (13)3.3. Matrices on curved physical elements.
The key difference between curvilinear and affine meshesis that geometric terms now vary spatially over each element. In practice, derivatives are computed over thereference element and mapped to the physical element D k through a change of variables formula J k ∂u∂x i = d (cid:88) j =1 G kij ∂u∂ (cid:98) x j , G kij = J k ∂ (cid:98) x j ∂x i , where we have defined the elements of the matrix G k as the derivatives of the reference coordinates (cid:98) x j withrespect to the physical coordinates x i on D k times the Jacobian of the transformation from reference tophysical coordinates J k . We denote evaluations of G kij at both volume and surface quadrature points as thevector G kij .We assume in this work that the mesh is stationary. It can be shown at the continuous level that, for anydifferentiable and invertible mapping, the quantity G k satisfies a geometric conservation law (GCL) [2, 20](14) d (cid:88) j =1 ∂∂ (cid:98) x j G kij = 0 , r that (cid:98) ∇ · G k = 0. Using (14), the scaled physical derivative J k ∂u∂x i can be computed via(15) J k ∂u∂x i = 12 d (cid:88) j =1 (cid:32) G kij ∂u∂ (cid:98) x j + ∂ (cid:0) G kij u (cid:1) ∂ (cid:98) x j (cid:33) . We will require the following assumptions on the mesh, as well as the geometric terms and outward normalvectors:
Assumption 2 (Mesh assumptions) . We assume that the mesh is quasi-uniform. The mesh is also assumedto be watertight, such that normals are consistent across neighboring elements as follows: for a shared face f between D k and D k, + , the scaled outward normal vectors for each element are equal and opposite at allpoints such that n J kf = − n + J k, + f . (16) We also assume that the scaled matrix of geometric terms transforms scaled reference normal vectors toscaled physical normals, such that d (cid:88) j =1 G kij (cid:98) n j (cid:98) J f = n i J kf , d (cid:88) j =1 (cid:0) G kij (cid:1) f ◦ (cid:16)(cid:98) n j ◦ (cid:98) J f (cid:17) = (cid:0) n i ◦ J kf (cid:1) , (17) where n j and J kf are vectors containing evaluations of the physical unit normals and face Jacobian factorsfor D k at surface quadrature points, respectively. Likewise, (cid:0) G kij (cid:1) f is a vector containing evaluations of G kij at the surface quadrature points. The properties (16) and (17) hold at the continuous level for a watertight mesh [33], and thus at all pointswhere the geometric terms are computed exactly. However, we will also consider cases where the geometricterms G kij are modified to enforce a discrete form of (14); in these situations, it will be important to ensurethat (17) holds after such modifications.Similar to what is done to stabilize finite difference discretizations [34, 35], we define physical differentiationmatrices based on the approximation of (15). Define D ik as D ik = 12 d (cid:88) j =1 (cid:16) diag (cid:0) G kij (cid:1) D jN + D jN diag (cid:0) G kij (cid:1)(cid:17) . Using properties of the Hadamard product [36], we can rewrite D ik as(18) D ik = d (cid:88) j =1 (cid:16) D jN ◦ (cid:8) { G kij } (cid:9)(cid:17) , (cid:8) { G kij } (cid:9) mn = 12 (cid:16)(cid:0) G kij (cid:1) m + (cid:0) G kij (cid:1) n (cid:17) , where (cid:8) { G kij } (cid:9) denotes the matrix of averages between each of the entries of G kij . From Assumption 2and (18), it is straightforward to show that (because (cid:8) { G kij } (cid:9) is symmetric) Q ik = W N D ik also satisfies asummation-by-parts property(19) Q ik + (cid:0) Q ik (cid:1) T = B ik , B ik = (cid:32) W f diag (cid:16) n i ◦ J kf (cid:17) (cid:33) . Curvilinear mappings also imply that integrals over each physical element D k are no longer simple scalingsof integrals over (cid:98) D . The L projection of u ∈ L (cid:0) D k (cid:1) over a curvilinear element D k is defined through(20) (cid:0) Π kN u, v (cid:1) D k = ( u, v ) D k , ∀ v ∈ P N (cid:16) (cid:98) D (cid:17) . Mapping integrals to the reference element (cid:98) D yields(21) (cid:0) Π kN u, vJ k (cid:1) (cid:98) D = (cid:0) u, vJ k (cid:1) (cid:98) D , ∀ v ∈ P N (cid:16) (cid:98) D (cid:17) . For affine elements, J k is constant and can be cancelled. Thus, the L projection over affine elements isequivalent to simply taking the L projection of a function over the reference element. However, for curvedelements, J acts as a spatially varying weight within the L inner product. iscretizing (21) requires a weighted mass matrix. We define a curved mass matrix over an element D k by weighting the discrete L norm with values of J at quadrature points(22) M k = V Tq W diag (cid:0) J k (cid:1) V q , where J k is a vector containing evaluation of the physical Jacobian factors for D k at volume quadraturepoints. Then, curvilinear L projection and lift matrices can be defined in a manner analogous to (10)(23) P kq = (cid:0) M k (cid:1) − V Tq W diag (cid:0) J k (cid:1) , L kq = (cid:0) M k (cid:1) − V Tf W f diag (cid:0) J kf (cid:1) . These matrices are distinct from element to element, reflecting the fact that problem (21) is distinct fromelement to element.3.4.
A discretely entropy stable DG formulation on curved meshes.
Given the matrices in Sec-tion 3.3, we can now define a local entropy stable DG formulation on an element D k . Here, we seek anapproximation solution u N ( x , t ) to (1), which is represented using vector-valued coefficients u h ( t ) such that u N ( x , t ) = N p (cid:88) j =1 ( u h ( t )) j φ j ( x ) , ( u h ( t )) j ∈ R n . Since the coefficients are vector valued, we assume that all matrices act component-wise on u h in theKronecker product sense.We first define the numerical fluxes f i,S ( u L , u R ) as the bivariate function of “left” and “right” conservativevariable states u L , u R . Definition 1.
The numerical flux f i,S ( u L , u R ) is entropy conservative (or entropy stable) if it satisfies thefollowing conditions:(1) f i,S ( u L , u R ) = f i,S ( u R , u L ) (symmetry).(2) f i,S ( u , u ) = f i ( u ) (consistency).(3) f i,S is referred to as entropy conservative if it satisfies conditions 1, 2, and( v L − v R ) T f i,S ( u L , u R ) = ψ i ( u L ) − ψ i ( u R ) . We now introduce the L projection of the entropy variables v h and the entropy-projected conservativevariables (cid:101) uu q = V q u h , v h = P kq v ( u q ) , (cid:101) v = (cid:20) (cid:101) v q (cid:101) v f (cid:21) = (cid:20) V q V f (cid:21) v h , (cid:101) u = (cid:20) (cid:101) u q (cid:101) u f (cid:21) = u ( (cid:101) v ) . (24)In (24), the entropy-projected conservative variables (cid:101) u denote the evaluation of the conservative variablesin terms of the projected entropy variables at volume and face quadrature points. We note that, underan appropriate choice of quadrature on quadrilaterals and hexahedra, this approach is equivalent to theapproach taken in [14], where the entropy variables are evaluated at Gauss nodes, then interpolated to adifferent set of nodes and used to compute the nonlinear fluxes.We now introduce a semi-discrete DG formulation for u h d u h dt + (cid:2) P kq L kq (cid:3) d (cid:88) j =1 (cid:16) D jk ◦ F j,S (cid:17) + d (cid:88) j =1 L kq diag ( n j ) (cid:0) f ∗ j − f j ( (cid:101) u f ) (cid:1) = 0 , (25) ( F j,S ) mn = f j,S (( (cid:101) u ) m , ( (cid:101) u ) n ) , ≤ m, n ≤ N q + N fq , f ∗ j = f j,S ( (cid:101) u + f , (cid:101) u f ) on interior faces , where (cid:101) u + denotes the values of the entropy-projected conservative variables on the neighboring elementacross each face of D k , and f ∗ j on the boundary denotes the j th component of some numerical flux throughwhich boundary conditions are imposed. Note that the face/surface Jacobian factors J kf are incorporatedinto the definition of L kq .Define the diagonal boundary quadrature matrix W ∂ Ω such that( W ∂ Ω ) ii = (cid:40) W f , if (cid:98) x fi is on the ∂ Ω0 , otherwise . e have the following semi-discrete statement of entropy conservation: Theorem 1.
Let f i,S be an entropy conservative flux from Definition 1 and assume that Q jk = 0 for j = 1 , . . . , d over each element D k . Then, (25) is entropy conservative in the sense that (cid:88) k T diag (cid:0) J k (cid:1) W d U ( u q )dt = (cid:88) k d (cid:88) j =1 T diag (cid:0) n j ◦ J kf (cid:1) W ∂ Ω (cid:0) ψ j ( (cid:101) u f ) − (cid:101) v Tf f ∗ j (cid:1) . Proof.
Under the assumption that Q jk = 0 for i = 1 , . . . , d over each element and (19), the proof of entropyconservation is identical to that of [1]. (cid:3) An entropy stable scheme can be constructed by adding an entropy-dissipating penalty term, such asa Lax-Friedrichs penalization or the matrix dissipation terms introduced in [37, 38]. For example, Lax-Friedrichs penalization can be incorporated by replacing the flux term with L kq diag ( n j ) (cid:0) f ∗ j − f ( u ) (cid:1) = ⇒ L kq (cid:18) diag ( n j ) (cid:0) f ∗ j − f ( u ) (cid:1) − λ (cid:74) (cid:101) u f (cid:75) (cid:19) , where λ is an estimate of the maximum eigenvalue of ∂ f ( u ) ∂ u [16, 1].4. Discretely stable and low storage DG methods on curved meshes
A disadvantage of the formulation (25) is high storage costs, especially at high orders of approximation.While the matrices Q ik can be applied to a vector without needing to explicitly store the matrix, the projectionand lifting matrices (23) differ from element to element, necessitating either explicit pre-computation andstorage or the assembly and inversion of a weighted mass matrix for each right hand side evaluation. Thelatter option is computationally expensive, while the former option increases storage costs. This increasein storage can result in suboptimal performance on modern computational architectures [21], due to theincreasing cost of memory operations and data movement compared to arithmetic operations.In this section, we present a discretely entropy stable scheme which avoids this high storage cost throughthe use of a low-storage weight-adjusted approximation to the inverse of a weighted mass matrix. To ensurea discrete entropy conservation or a discrete entropy inequality, we also modify the formulation (25) to takeinto account the use of a weight-adjusted mass matrix.4.1. A weight-adjusted approximation to the curvilinear mass matrix.
The presence of the weighted L inner product (cid:0) J k Π kN u, v (cid:1) (cid:98) D in (21) results in the presence of a weighted mass matrix. Because the weight J k varies spatially over each element, the inverse of a weighted mass matrix is no longer a scaling of theinverse reference mass matrix. The motivation for the weight-adjusted mass matrix is to replace the inversionof weighted mass matrices over each element with the application of inverse reference mass matrices andquadrature-based operations involving the spatially varying weights J k [3, 4].To define a weight-adjusted approximation to the curvilinear L inner product, we first define the operator T − w : L → P N as follows(26) (cid:0) wT − w u, v (cid:1) (cid:98) D = ( u, v ) (cid:98) D , ∀ v ∈ P N (cid:16) (cid:98) D (cid:17) . Roughly speaking, T w u approximates u/w . Thus, taking w = 1 /J k provides an approximation of thecurvilinear L inner product (cid:0) J k u, v (cid:1) (cid:98) D ≈ (cid:16) T − /J k u, v (cid:17) (cid:98) D . Computing T − /J k u requires solving (26). Let u ∈ P N (cid:16) (cid:98) D (cid:17) , and let u J denote coefficients for the polyno-mial T − /J k u . This results in the following matrix system M /J k u J = M u , M /J k = V Tq W diag (cid:0) / J k (cid:1) V q , hich implies that, when restricted to polynomials, the matrix form of T − /J k is M − /J k M . Then, the weight-adjusted mass matrix is the Gram matrix with respect to the weight-adjusted inner product (cid:16) T − /J k u, v (cid:17) (cid:98) D ,such that M k ≈ M M − /J k M , (cid:0) M k (cid:1) − ≈ M − M /J k M − . The inverse of the weight-adjusted mass matrix can be applied in a matrix-free fashion by using quadrature toform M /J k . This requires storage of the inverse reference mass matrix and the values of J k at quadraturepoints. Assuming that the number of quadrature points scales as O ( N d ) in d dimensions, this yields astorage cost of O ( N d ) per-element compared to an O ( N d ) per element storage cost required for the storageof inverse weighted mass matrices (cid:0) M k (cid:1) − . This application of the weight-adjusted mass matrix is typicallyapplied using the L projection matrix P q as follows M − M /J k M − = P q diag (cid:0) / J k (cid:1) V q M − . When evaluating the right hand side of a semi-discrete formulation such as (25), the inverse mass matrix istypically merged into operations on the right hand side, such that the main work in applying the weight-adjusted mass matrix consists of applying the interpolation matrix V q , scaling by pointwise values of 1 / J k at quadrature points, and multiplying by the L projection matrix P q .4.2. A discretely entropy stable low storage DG formulation on curved meshes.
Given the weight-adjusted inverse mass matrix, we can also define a weight-adjusted version of the L projection over a curvedelement D k . We refer to this operator as (cid:101) Π kN : L → P N , which satisfies (cid:16) T − /J k (cid:101) Π kN u, v (cid:17) (cid:98) D = (cid:0) uJ k , v (cid:1) (cid:98) D , ∀ v ∈ P N (cid:16) (cid:98) D (cid:17) . It was shown in [4] that (cid:101) Π kN is given explicitly by(27) (cid:101) Π kN u = Π N (cid:18) J k Π N (cid:0) uJ k (cid:1)(cid:19) , where Π N is the L projection operator on the reference element (cid:98) D . We can discretize (cid:101) Π kN using quadratureto yield a weight-adjusted projection matrix (cid:101) P kq (cid:101) P kq = M − M / J k M − V Tq W diag (cid:0) J k (cid:1) = M − V Tq W diag (cid:0) / J k (cid:1) V q P q diag (cid:0) J k (cid:1) = P q diag (cid:0) / J k (cid:1) V q P q diag (cid:0) J k (cid:1) . (28)We can similarly define a weight-adjusted lifting matrix (cid:101) L q by replacing the weighted mass matrix in (23)with the weight-adjusted mass matrix (cid:101) L kq = M − M / J k M − V Tf W f diag (cid:0) J kf (cid:1) = M − V Tq W diag (cid:0) / J k (cid:1) V q L q diag (cid:0) J kf (cid:1) = P q diag (cid:0) / J k (cid:1) V q L q diag (cid:0) J kf (cid:1) . (29)We can now introduce the weight-adjusted projection of the entropy variables v h and the correspondingentropy-projected conservative variables (cid:101) uu q = V q u h , v h = (cid:101) P kq v ( u q ) , (cid:101) v = (cid:20) V q V f (cid:21) v h , (cid:101) u = (cid:20) (cid:101) u q (cid:101) u f (cid:21) = u ( (cid:101) v ) . (30)A semi-discrete DG formulation for u h can be constructed using the variables defined in (30)d u h dt + (cid:104) (cid:101) P kq (cid:101) L kq (cid:105) d (cid:88) j =1 (cid:16) D jk ◦ F j,S (cid:17) + d (cid:88) j =1 (cid:101) L kq diag ( n j ) (cid:0) f ∗ j − f j ( (cid:101) u f ) (cid:1) = 0 , (31) ( F j,S ) mn = f j,S (( (cid:101) u ) m , ( (cid:101) u ) n ) , ≤ m, n ≤ N q + N fq , f ∗ j = f j,S ( (cid:101) u + f , (cid:101) u f ) on interior faces . Since the weight-adjusted mass matrix inverse is low-storage, and since the matrices D jk in (18) can beassembled from reference matrices D iN and the values of geometric terms at quadrature points, the overall cheme requires only O ( N d ) storage per element. We can additionally show that formulation (31) is entropyconservative in the same sense as Theorem 1: Theorem 2.
Let f i,S be an entropy conservative flux from Definition 1 and assume that Q jk = 0 for j = 1 , . . . , d over each element D k . Then, (31) is entropy conservative in the sense that (cid:88) k T diag (cid:0) J k (cid:1) W d U ( u q )dt = (cid:88) k d (cid:88) j =1 T diag (cid:0) n j ◦ J kf (cid:1) W ∂ Ω (cid:0) ψ j ( (cid:101) u f ) − (cid:101) v Tf f ∗ j (cid:1) . Proof.
The proof is similar to that of Theorem 1. First, recall from the definitions of the weight-adjustedprojection matrix (28) and weight-adjusted lift matrix (29) that (cid:101) P kq = M − M /J k M − V Tq W diag (cid:0) J k (cid:1) , (32) (cid:101) L kq = M − M /J k M − V Tf W f diag (cid:0) J kf (cid:1) . (33)Then, multiplying by the weight-adjusted mass matrix M M − /J k M on both sides of (31) and using (32),(33) yields the weak form of (31)(34) M M − /J k M d u h dt + (cid:2) V Tq V Tf (cid:3) d (cid:88) j =1 (cid:16) Q jk ◦ F j,S (cid:17) + d (cid:88) j =1 V Tf W f diag ( n j ) (cid:0) f ∗ j − f j ( (cid:101) u f ) (cid:1) = 0 . Testing with the weight-adjusted projection of the entropy variables v h = (cid:101) P kq v ( u q ) and using (32) thenyields for the time term (cid:16) (cid:101) P kq v ( u q ) (cid:17) T M M − /J k M d u h dt = v ( u q ) T W diag (cid:0) J k (cid:1) V q M − M /J k M − M M − /J k M d u h dt= v ( u q ) T W diag (cid:0) J k (cid:1) d V q u h dt = T W diag (cid:0) J k (cid:1) (cid:18) diag ( v ( u q )) d u q dt (cid:19) = T W diag (cid:0) J k (cid:1) d U ( u q )dt . The remainder of the proof is identical to that of Theorem 1 and [1]. (cid:3)
Analysis of weight-adjusted projection.
The construction of the discretely entropy stable weight-adjusted DG formulation (31) replaces the L projection operator Π kN with the weight-adjusted projectionoperator (cid:101) Π kN . While this preserves entropy stability, it is unclear whether (cid:101) Π kN is high order accurate. In thissection, we prove that the weight-adjusted projection is high order accurate due to the fact that, for a fixedgeometric mapping and sufficiently regular u , the difference between the L and weight-adjusted projectionis (cid:13)(cid:13)(cid:13) Π kN u − (cid:101) Π kN u (cid:13)(cid:13)(cid:13) L (Ω) = O ( h N +2 ). Because the approximation error for the L projection is O ( h N +1 ) forsufficiently regular u , the difference between the L and weight-adjusted projection converges faster thanthe L best approximation error. Consequentially, solutions computed using the L and weight-adjustedprojection are typically indistinguishable for a fixed geometric mapping [39].We first note that (cid:101) Π kN = Π N (cid:0) J k Π N (cid:0) uJ k (cid:1)(cid:1) is self-adjoint with respect to the J -weighted L inner product(35) (cid:16) J k (cid:101) Π kN u, v (cid:17) (cid:98) D = (cid:18) Π N (cid:18) J k Π N (cid:0) uJ k (cid:1)(cid:19) , vJ k (cid:19) (cid:98) D = (cid:18) uJ k , Π N (cid:18) J k Π N (cid:0) vJ k (cid:1)(cid:19)(cid:19) (cid:98) D = (cid:16) uJ k , (cid:101) Π kN v (cid:17) (cid:98) D . Furthermore, using that the operator T − /J k is self-adjoint for v ∈ P N (cid:16) (cid:98) D (cid:17) with respect to the L innerproduct [3], we find that a projection-like property holds for the weight-adjusted L inner product(36) (cid:16) T − /J k (cid:101) Π kN u, v (cid:17) (cid:98) D = (cid:18) J k Π N ( uJ k ) , T − /J k v (cid:19) (cid:98) D = (cid:0) Π N ( uJ k ) , v (cid:1) (cid:98) D = (cid:0) uJ k , v (cid:1) (cid:98) D , ∀ v ∈ P N (cid:98) D. To prove (cid:13)(cid:13)(cid:13) Π kN u − (cid:101) Π kN u (cid:13)(cid:13)(cid:13) L ( (cid:98) D ) = O ( h N +2 ), we use a generalized inverse inequality and results from [3, 4].We first introduce a modification of Theorem 3.1 in [40, 3] heorem 3. Let D k be a quasi-regular element with representative size h = diam (cid:0) D k (cid:1) , and let (cid:98) D be thereference element. For N ≥ , w ∈ W N +1 , ∞ (cid:0) D k (cid:1) , and u ∈ W r, (cid:0) D k (cid:1) , (cid:13)(cid:13)(cid:13)(cid:13) u − w Π N ( uw ) (cid:13)(cid:13)(cid:13)(cid:13) L ( (cid:98) D ) ≤ Ch min( r,N +1) (cid:13)(cid:13)(cid:13)(cid:13) √ J k (cid:13)(cid:13)(cid:13)(cid:13) L ∞ ( D k ) (cid:13)(cid:13)(cid:13)(cid:13) w (cid:13)(cid:13)(cid:13)(cid:13) L ∞ ( D k ) (cid:107) w (cid:107) W N +1 , ∞ ( D k ) (cid:107) u (cid:107) W r, ( D k ) . The proof is a straightforward modification of the proofs presented in [40, 3] accounting for reducedregularity of u when r < ( N + 1). The next result we need is a generalized inverse inequality. Lemma 1.
Let v ∈ P N (cid:16) (cid:98) D (cid:17) , and let h = diam (cid:0) D k (cid:1) . Then, (cid:107) v (cid:107) W N +1 , ( D k ) ≤ C N h − N (cid:13)(cid:13)(cid:13) √ J k (cid:13)(cid:13)(cid:13) L ∞ (cid:13)(cid:13)(cid:13)(cid:13) √ J k (cid:13)(cid:13)(cid:13)(cid:13) L ∞ (cid:13)(cid:13)(cid:13)(cid:13) J k G k (cid:13)(cid:13)(cid:13)(cid:13) W N, ∞ ( D k ) (cid:107) v (cid:107) L ( D k ) . where C N depends on N , but is independent of h .Proof. By applying Fa`a di Bruno’s formula, we can express the degree ( N + 1) Sobolev norm of v on D k interms of derivatives of v on the reference element (cid:98) D . Noting that all ( N + 1) derivatives of v disappear for v ∈ P N (cid:16) (cid:98) D (cid:17) allows us to bound the degree ( N + 1) Sobolev norm of v by its degree N Sobolev norm (cid:107) v (cid:107) W N +1 , ( D k ) ≤ C N (cid:13)(cid:13)(cid:13)(cid:13) J k G k (cid:13)(cid:13)(cid:13)(cid:13) W N, ∞ ( D k ) (cid:107) v (cid:107) W N, ( D k ) , where G k is the matrix of scaled geometric terms for D k . Then, a scaling argument [33, 41] yields (cid:107) v (cid:107) W N, ( D k ) ≤ C h − N (cid:13)(cid:13)(cid:13) √ J k (cid:13)(cid:13)(cid:13) L ∞ (cid:107) v (cid:107) W N, ( (cid:98) D ) . (37)The quantity (cid:107) v (cid:107) W N, ( (cid:98) D ) can be bounded by noting that v ∈ P N (cid:16) (cid:98) D (cid:17) . Since P N (cid:16) (cid:98) D (cid:17) is finite-dimensional,the Sobolev norm can be bounded from above by the L norm of (cid:98) D with a constant C depending on N, d .By another scaling argument, we have (cid:107) v (cid:107) W N +1 , ( (cid:98) D ) ≤ C (cid:107) v (cid:107) L ( (cid:98) D ) ≤ (cid:13)(cid:13)(cid:13)(cid:13) √ J k (cid:13)(cid:13)(cid:13)(cid:13) L ∞ C (cid:107) v (cid:107) L ( D k ) . (38) (cid:3) We can now prove that (cid:101) Π kN u is superconvergent to the curvilinear L projection Π kN u : Theorem 4.
Let u ∈ W r, (cid:0) D k (cid:1) . The difference between the L projection Π kN u and the weight-adjustedprojection (cid:101) Π kN u is (cid:13)(cid:13)(cid:13) Π kN u − (cid:101) Π kN u (cid:13)(cid:13)(cid:13) L ( D k ) ≤ C N C J h min( r,N +1)+1 (cid:107) u (cid:107) W N +1 , ( D k ) , where C N is a mesh-independent constant which depends on N, d and C J is C J = (cid:13)(cid:13) J k (cid:13)(cid:13) . L ∞ ( D k ) (cid:13)(cid:13)(cid:13)(cid:13) J k (cid:13)(cid:13)(cid:13)(cid:13) . L ∞ ( D k ) (cid:13)(cid:13)(cid:13)(cid:13) J k (cid:13)(cid:13)(cid:13)(cid:13) W N +1 , ∞ ( D k ) (cid:13)(cid:13) J k (cid:13)(cid:13) W N +1 , ∞ ( D k ) (cid:13)(cid:13)(cid:13)(cid:13) J k G k (cid:13)(cid:13)(cid:13)(cid:13) W N, ∞ ( D k ) . Proof.
We can rewrite the norm of the difference between the weight-adjusted and L projections (cid:13)(cid:13)(cid:13) Π kN u − (cid:101) Π kN u (cid:13)(cid:13)(cid:13) L ( D k ) = (cid:16) Π kN u − (cid:101) Π kN u, vJ k (cid:17) (cid:98) D , v = Π kN u − (cid:101) Π kN u. Because v ∈ P N (cid:16) (cid:98) D (cid:17) , we can also evaluate the squared error as (cid:13)(cid:13)(cid:13) Π kN u − (cid:101) Π kN u (cid:13)(cid:13)(cid:13) L ( D k ) = (cid:12)(cid:12)(cid:12)(cid:0) Π kN u, vJ k (cid:1) (cid:98) D − (cid:16)(cid:101) Π kN u, vJ k (cid:17) (cid:98) D (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:0) u, vJ k (cid:1) (cid:98) D − (cid:16)(cid:101) Π kN u, vJ k (cid:17) (cid:98) D (cid:12)(cid:12)(cid:12) ≤ (cid:13)(cid:13) J k (cid:13)(cid:13) L ∞ ( D k ) (cid:12)(cid:12)(cid:12)(cid:16) u − (cid:101) Π kN u, v (cid:17) (cid:98) D (cid:12)(cid:12)(cid:12) . e can then note that (cid:101) Π kN u = Π N (cid:0) J k Π N (cid:0) uJ k (cid:1)(cid:1) to show that (cid:16) u − (cid:101) Π kN u, v (cid:17) (cid:98) D = (cid:16) uJ k , vJ k (cid:17) (cid:98) D − (cid:18) Π N (cid:18) J k Π N (cid:0) uJ k (cid:1)(cid:19) , v (cid:19) (cid:98) D = (cid:16) uJ k , vJ k (cid:17) (cid:98) D − (cid:16) Π N (cid:0) uJ k (cid:1) , vJ k (cid:17) (cid:98) D . Adding and subtracting (cid:0) Π N (cid:0) uJ k (cid:1) , Π N (cid:0) vJ k (cid:1)(cid:1) (cid:98) D and using Theorem 3 (noting that v ∈ W N +1 , (cid:0) D k (cid:1) sinceit is polynomial) gives (cid:13)(cid:13)(cid:13) Π kN u − (cid:101) Π kN u (cid:13)(cid:13)(cid:13) L ( D k ) ≤ (cid:13)(cid:13) J k (cid:13)(cid:13) L ∞ ( D k ) (cid:12)(cid:12)(cid:12)(cid:16) u − (cid:101) Π kN u, v (cid:17) (cid:98) D (cid:12)(cid:12)(cid:12) = (cid:13)(cid:13) J k (cid:13)(cid:13) L ∞ ( D k ) (cid:12)(cid:12)(cid:12)(cid:16) uJ k − Π N (cid:0) uJ k (cid:1) , vJ k − Π N (cid:16) vJ k (cid:17)(cid:17) (cid:98) D (cid:12)(cid:12)(cid:12) ≤ (cid:13)(cid:13) J k (cid:13)(cid:13) L ∞ ( D k ) (cid:13)(cid:13) uJ k − Π N (cid:0) uJ k (cid:1)(cid:13)(cid:13) (cid:98) D (cid:13)(cid:13)(cid:13) vJ k − Π N (cid:16) vJ k (cid:17)(cid:13)(cid:13)(cid:13) (cid:98) D ≤ Ch min( r,N +1)+ N +1 (cid:101) C J (cid:107) u (cid:107) W r, ( D k ) (cid:107) v (cid:107) W N +1 , ( D k ) , where (cid:101) C J = (cid:13)(cid:13) J k (cid:13)(cid:13) L ∞ ( D k ) (cid:13)(cid:13)(cid:13)(cid:13) J k (cid:13)(cid:13)(cid:13)(cid:13) L ∞ ( D k ) (cid:13)(cid:13) J k (cid:13)(cid:13) W N +1 , ∞ ( D k ) (cid:13)(cid:13)(cid:13)(cid:13) J k (cid:13)(cid:13)(cid:13)(cid:13) W N +1 , ∞ ( D k ) . Applying Lemma 1 then yields (cid:13)(cid:13)(cid:13) Π kN u − (cid:101) Π kN u (cid:13)(cid:13)(cid:13) L ( D k ) ≤ C N h min( r,N +1)+1 C J (cid:107) u (cid:107) W N +1 , ( D k ) (cid:107) v (cid:107) L ( D k ) . Dividing through by (cid:107) v (cid:107) L ( D k ) = (cid:13)(cid:13)(cid:13) Π kN u − (cid:101) Π kN u (cid:13)(cid:13)(cid:13) L ( D k ) gives the desired result. (cid:3) Theorem 3 can be used to show that the L error between Π kN u, (cid:101) Π kN u and u ∈ W r, (cid:0) D k (cid:1) is O ( h r ).Theorem 4 demonstrates that the L difference between Π kN u, (cid:101) Π kN u is O ( h r +1 ), or at least one order higherthan the approximation error. We note that optimal convergence of the weight-adjusted projection requiresthat the geometric mapping Φ k is asymptotically affine (i.e., the Sobolev norm of J k , G k does not growunder mesh refinement), which is ensured under nested mesh refinement and appropriate curvilinear blendingstrategies [42, 40, 4].4.3.1. Local and global conservation.
We next address local conservation of the weight-adjusted scheme (31),which is also referred to as primary conservation [43, 12, 13, 15]. We begin by noting that (31) is locallyconservative with respect to the weight-adjusted inner product. Testing (34) with yields T M M − /J k M d u h dt + d (cid:88) j =1 T (cid:16) Q jk ◦ F j,S (cid:17) + d (cid:88) j =1 T W f diag ( n j ) (cid:0) f ∗ j − f j ( (cid:101) u f ) (cid:1) = 0 . Local conservation can be shown by applying the SBP property (19) and noting that (cid:0) B ik ◦ F j,S (cid:1) = W f diag (cid:16) n j ◦ J kf (cid:17) f j ( u ) (using the consistency of f S and diagonal nature of B ik ) yields T M M − /J k M d u h dt + d (cid:88) j =1 T (cid:18)(cid:18) Q jk − (cid:16) Q jk (cid:17) T (cid:19) ◦ F j,S (cid:19) + d (cid:88) j =1 T W f diag ( n j ) (cid:0) f ∗ j (cid:1) = 0 . Since (cid:18) Q jk − (cid:16) Q jk (cid:17) T (cid:19) is skew-symmetric and F j,S is symmetric, the Hadamard product of these two matricesis skew-symmetric. As a result, T (cid:18)(cid:18) Q jk − (cid:16) Q jk (cid:17) T (cid:19) ◦ F j,S (cid:19) = 0 and(39) T M M − /J k M d u h dt + d (cid:88) j =1 T W f diag ( n j ) (cid:0) f ∗ j (cid:1) = 0 . lobal conservation is shown by summing (39) over all elements D k . Because the flux f ∗ j is single-valued on each face and the outward normal n j changes sign on adjacent elements, the contributions T W f diag ( n j ) (cid:0) f ∗ j (cid:1) cancel on all interior interfaces. On periodic meshes, this yields conservation withrespect to the weight-adjusted mass matrix (cid:88) k T M M − /J k M d u h dt = 0 . However, while the weight-adjusted approximation of the mass matrix is high order accurate and efficient,it does not preserve the average over a physical element, which is equivalent to the J -weighted average overthe reference element. This is due to the fact that, in general,(40) (cid:90) (cid:98) D uJ k d (cid:98) x − (cid:90) (cid:98) D T − /J k u d (cid:98) x ≈ T M J k u − T M M − /J k M u (cid:54) = 0 . Results in [3] show that the difference between the true mean and weight-adjusted mean in (40) convergesextremely fast at a rate of O ( h N +2 ). However, for systems of conservation laws, it is often desired that thelocal element average is preserved exactly up to machine precision. We present two simpler approaches toensuring local conservation in this section.The simplest way to ensure local conservation is to approximate J using a degree N polynomial andutilize a sufficiently accurate quadrature. Then, we have the following lemma: Lemma 2.
Let J k ∈ P N , and let integrals be computed using quadrature which is exact for degree N polynomials. Then, T M M − /J k M u = (cid:90) (cid:98) D T − /J k u d (cid:98) x = (cid:90) (cid:98) D uJ k d (cid:98) x . Proof.
The proof relies on (26), which states that (cid:16) J k T − /J k u, v (cid:17) (cid:98) D = ( u, v ) (cid:98) D for all v ∈ P N . If J k is apolynomial of degree N , then taking v = 1 yields (cid:16) T − /J k u, (cid:17) (cid:98) D = (cid:18) J k T − /J k u, J k (cid:19) (cid:98) D = (cid:0) u, J k (cid:1) (cid:98) D = (cid:90) (cid:98) D uJ k d (cid:98) x . Additionally, the proof still holds if integrals are approximated using a quadrature rule which exactly inte-grates uJ k ∈ P N . (cid:3) We note that, for isoparametric curved elements, J k (cid:54)∈ P N in general (in 2D, J k ∈ P N − , while in 3D, J k ∈ P N − [44]). Thus, to ensure local conservation, we will approximate J k using a degree N polynomial(for example, the interpolant or L projection onto P N ). We note that this approximation is required onlyin the weight-adjusted mass matrix, and does not modify the scaled geometric terms G kij .The second approach relies on a simple correction which restores exact conservation of the true mean.In [4], it was shown that a rank one correction of the weight-adjusted mass matrix inverse preserves themean exactly. However, this requires the use of the Sherman–Morrison formula to compute the inverse of arank one matrix update, which can be cumbersome to incorporate. We present a simpler explicit correctionformula which does not involve matrices. Let u J and u WADG be defined through( u J J, v ) (cid:98) D = ( f, v ) (cid:98) D , ∀ v ∈ P N (cid:16) (cid:98) D (cid:17) , (cid:16) T − /J u WADG , v (cid:17) (cid:98) D = ( f, v ) (cid:98) D , ∀ v ∈ P N (cid:16) (cid:98) D (cid:17) , where u J corresponds to the inversion of the weighted mass matrix and u WADG corresponds to the inversionof a weight-adjusted mass matrix. For example, if f = uJ for some function u ( x ), then u J = Π kN u and u WADG = (cid:101) Π kN u . To ensure both local and global conservation, we require that the weighted average of Conservation still holds when incorporating entropy dissipation through penalty terms involving jumps. This is because thedefinition of the jump changes sign on adjacent elements, such that jump contributions cancel when summing over all elements. WADG is the same as the weighted average of u J . Let the conservative weight-adjusted (cid:101) u WADG be definedas(41) (cid:101) u WADG = u WADG + (cid:82) (cid:98) D ( f − Ju WADG ) d (cid:98) x (cid:82) (cid:98) D J d (cid:98) x . Taking the weighted integral of (cid:101) u WADG yields (cid:90) (cid:98) D (cid:101) u WADG J d (cid:98) x = (cid:90) (cid:98) D f d (cid:98) x . In other words, (41) ensures local conservation by correcting the weighted average of u WADG to match thatof u J . Applying this correction to the right hand side of (31) then yields a scheme which locally and globallyconserves mean values of the conservative variables.This approach is applicable to an arbitrary weight, and can be generalized to matrix-valued weights aswell [21]. Moreover, using Theorem 6 in [3], one can show that the L norm of the difference (cid:101) u WADG − u WADG is O ( h N +1 ), and does not affect high order accuracy.5. Enforcing the discrete geometric conservation law
An important aspect of Theorem 1 is the assumption that Q jk = 0 for j = 1 , . . . , d . However, this is notalways guaranteed to hold for Q jk as defined through (18). In this section, we discuss methods of constructingthe geometric terms G kij for curvilinear meshes in a way that ensures Q jk = 0.From (18), the condition Q jk = 0 is equivalent to Q jk = W N d (cid:88) j =1 D jN ◦ (cid:8) { G kij } (cid:9) = 12 W N d (cid:88) j =1 (cid:16) diag (cid:0) G kij (cid:1) D jN + D jN diag (cid:0) G kij (cid:1) (cid:17) = 12 W N d (cid:88) j =1 D jN G kij = 0 , where we have used that D jN = 0 to eliminate the first term. Since W N is a diagonal matrix with positiveentries, Q jk = 0 is equivalent to ensuring that a discrete version of the GCL (14) holds(42) d (cid:88) j =1 D jN G kij = 0 . This condition is required to ensure that free-stream preservation holds at the discrete level. In other words,we wish to ensure that the semi-discrete scheme preserves (for u constant) ∂ u ∂t + ∇ · f ( u ) = ∂ u ∂t = 0 . For isoparametric geometric mappings (where the degree of the mapping matches the degree of the polynomialapproximation) in two dimensions, the GCL is naturally enforced by the “cross-product” form, noting thatthe scaled metric terms G kij are exactly polynomials of degree N . As a result, computing the metric termsexactly automatically enforces that both the continuous GCL (14) and the discrete GCL (42) are satisfied.However, the discrete GCL is not always maintained at the discrete level in 3D.In three dimensions, geometric terms are typically computed in “cross-product” form G k i G k i G k i = ∂ x ∂ (cid:98) x j × ∂ x ∂ (cid:98) x k , ( i, j, k ) = (1 , , , cyclic . (43)Note the abuse of notation here and in the sequel, the superscript k refers to the element number and thesubscript k to the cyclic index. This formula can be used to compute the geometric terms exactly at volumeand surface quadrature points. However, because ∂ x ∂ (cid:98) x j , ∂ x ∂ (cid:98) x k ∈ P N − , the geometric terms G kij are polynomialsof degree P N − . The discrete GCL condition holds only if G kij ∈ P N − are differentiated exactly; however, ecause applying D jN involves the L projection, and because G kij and its L projection onto degree N polynomials can differ, the discrete GCL (42) does not hold in general [2].This can be remedied by using an alternative form of the geometric terms, which ensures that (42) issatisfied a-priori [20, 45, 2]. The geometric terms G kij can also be computed using a “conservative curl” form,where G kij are the image of the curl applied to some quantity as follows: G k i G k i G k i = (cid:16) − (cid:98) ∇ × (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i (cid:16) (cid:98) ∇ × (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i (cid:16) (cid:98) ∇ × (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i , (44)where the subscript i denotes the i th component of the vector quantity. From (44), it can be observed that,because the divergence of a curl vanishes, the continuous GCL condition (14) holds. The central idea of[45, 2] is to use (44), but to interpolate before applying the curl G k i G k i G k i = (cid:16) − (cid:98) ∇ × I N (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i (cid:16) (cid:98) ∇ × I N (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i (cid:16) (cid:98) ∇ × I N (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i , (45)where I N denotes the degree N polynomial interpolation operator. Since the geometric terms are stillcomputed by applying a curl, the continuous GCL condition (14) is still satisfied. We shall also show thatthis approximation also satisfies the discrete GCL condition.We adopt a slight modification of (45) in this work which is tailored towards triangular and tetrahedralelements. Because the geometric terms are computed by applying the curl, the geometric terms are approxi-mated as degree ( N −
1) polynomials rather than degree N polynomials, which can reduce accuracy. Instead,we approximate geometric terms by using the interpolation operator I N +1 onto degree ( N + 1) polynomials,then interpolating back to degree N polynomials G k i G k i G k i = I N (cid:16) − (cid:98) ∇ × I N +1 (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i I N (cid:16) (cid:98) ∇ × I N +1 (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i I N (cid:16) (cid:98) ∇ × I N +1 (cid:16) x (cid:98) ∇ x (cid:17)(cid:17) i . (46)For any u ∈ P N (cid:16) (cid:98) D (cid:17) , ∂u∂ (cid:98) x i ∈ P N − (cid:16) (cid:98) D (cid:17) for i = 1 , ,
3. Thus, the interpolation to degree N polynomials isexact, since the derivatives of a degree ( N + 1) polynomial are degree N on triangles and tetrahedra. Remark.
The accuracy of (46) depends on the choice of interpolation points. It is well known that interpo-lating at equispaced points can result in inaccurate polynomial approximations. One can determine good inter-polation point sets by optimizing over some measure of interpolation quality (such as the Lebesgue constant),and in practice, sets of interpolation points are pre-computed for some polynomial degrees N = 1 , . . . , N max on the reference element and stored [46, 47] , or explicitly computed as the image of equispaced points underan appropriately defined mapping [48, 49, 50] . To prove that the construction (46) satisfies Assumption 2, we must assume that the interpolation pointsfor a degree N element include an appropriate number of points on each face. We note that these assumptionsexclude interpolation points which lie purely in the interior of an element, such as those introduced in [51, 52].We can now show that the geometric terms satisfy all conditions necessary to guarantee entropy stability: Theorem 5.
Let the mesh consist of triangles or tetrahedral elements which satisfy Assumption 2, and letthe interpolation points which define the degree N interpolation operator be distributed such that N fp points lieon each face. Then, the approximate geometric terms (cid:94) J k G kij and approximate scaled normals (cid:93) n i J kf computedusing (46) and (17)) satisfy both the discrete GCL condition (42) and Assumption 2. Additionally, the error n the approximation satisfies (cid:13)(cid:13)(cid:13) G kij − (cid:103) G kij (cid:13)(cid:13)(cid:13) L (Ω) ≤ C N | Ω | h N +2 (cid:118)(cid:117)(cid:117)(cid:116) d (cid:88) i =1 (cid:88) k (cid:107) r i (cid:107) W N +2 , ( D k ) . Proof.
The satisfaction of (42) relies on the fact that (cid:103) G kij is a degree N polynomial and is equal to its own L projection. Let (cid:103) G kij denote the polynomial coefficients of (cid:103) G kij . Then, applying D jN to evaluations of (cid:103) G kij at volume and surface quadrature points and using (13), we have d (cid:88) j =1 D jN (cid:20) V q V f (cid:21) (cid:103) G kij = d (cid:88) j =1 (cid:34) V q D i (cid:103) G kij (cid:35) . The entries of V q D i (cid:103) G kij correspond to values of the derivatives of (cid:103) G kij evaluated at quadrature points. Since d (cid:88) j =1 ∂∂ (cid:98) x j (cid:103) G kij = 0by construction using (46), (cid:80) j =1 ,...,d V q D i (cid:103) G kij = as well.Equation (17) of Assumption 2 is satisfied by directly constructing the scaled normals (cid:101) n J kf using values of (cid:103) G kij at quadrature points. We must now prove that the construction of (cid:101) n J kf implies that Equation (16) holds.This is not immediately clear; since the normals are constructed from the approximate geometric terms andthe formula (17), it is not guaranteed that (cid:101) n + J k, + f = − (cid:101) n J kf will hold across a shared face. However, thescaled normal vectors involve only nodal values on the shared face because the normals are computed in termsof the tangential reference derivatives [2]. Thus, assuming a watertight mesh, the interpolation nodes on twoneighboring elements will coincide for a shared face f , such that the trace of (cid:103) G kij from either neighboringelement will be the same lower-dimensional polynomial on f . This is sufficient to ensure that the scalednormal vectors (cid:101) n + J k, + f , (cid:101) n J kf will be equal and opposite.The local L error (cid:13)(cid:13)(cid:13) G kij − (cid:103) G kij (cid:13)(cid:13)(cid:13) L ( D k ) can be bounded by noting that, since the error G kij − (cid:103) G kij consistsof linear combinations of derivatives of the interpolation error r i − I N +1 r i , it can be bounded by the H -seminorm of the latter quantity (cid:13)(cid:13)(cid:13) G kij − (cid:103) G kij (cid:13)(cid:13)(cid:13) L ( D k ) ≤ d (cid:88) i =1 C (cid:13)(cid:13)(cid:13) √ J k (cid:13)(cid:13)(cid:13) L ( (cid:98) D ) | ( r i − I N +1 r i ) | H ( (cid:98) D ) ≤ d (cid:88) i =1 (cid:101) C N (cid:13)(cid:13)(cid:13) √ J k (cid:13)(cid:13)(cid:13) L ( (cid:98) D ) | r i | W N +2 , ( (cid:98) D ) , where we have used the Bramble-Hilbert lemma [41] on the reference element in the last step. Since it isapplied on the reference element (cid:98) D rather than the physical element D k , the constant (cid:101) C N depends on thereference element and order of approximation, but not the mesh size h . A scaling argument for quasi-uniformmeshes then yields that | r i | W N +2 , ( (cid:98) D ) ≤ C h N +2 (cid:107) r i (cid:107) W N +2 , ( D k ) . (47)The global estimate results from squaring (47), summing over all elements and using (cid:13)(cid:13)(cid:13) √ J k (cid:13)(cid:13)(cid:13) L ( (cid:98) D ) = (cid:12)(cid:12) D k (cid:12)(cid:12) . (cid:3) Remark.
It should be pointed out that this approach does not work on hexahedral elements. This is due tothe fact that the natural polynomial space on hexahedral elements is the tensor product space Q N (cid:16) (cid:98) D (cid:17) Q N (cid:16) (cid:98) D (cid:17) = (cid:110)(cid:98) x i . . . (cid:98) x i d d , (cid:98) x ∈ (cid:98) D, ≤ i k ≤ N (cid:111) . or u ∈ Q N +1 (cid:16) (cid:98) D (cid:17) , ∂u∂ (cid:98) x i is at most degree N in the coordinate (cid:98) x i , but can remain a polynomial of degree N + 1 in all other coordinates. Thus, interpolating from degree N + 1 to degree N polynomials in (46)introduces aliasing errors and is no longer exact. The result of (46) is no longer the image of a curl, andthus does not satisfy the discrete GCL by construction. We briefly outline how to compute (cid:103) G kij in three dimensions. Let (cid:8)(cid:98) x Ni (cid:9) N p i =1 denote the set of degree N interpolation points, and let (cid:96) Ni ( (cid:98) x ) denote the i th degree N Lagrange basis function on the reference element.We define interpolation matrices V N +1 N and V NN +1 between degree N and N + 1 polynomials such that (cid:0) V N +1 N (cid:1) ij = (cid:96) Nj ( (cid:98) x N +1 i ) , ≤ i ≤ N p , ≤ i ≤ ( N + 1) p (48) (cid:0) V NN +1 (cid:1) ij = (cid:96) N +1 j ( (cid:98) x Ni ) , ≤ i ≤ ( N + 1) p , ≤ i ≤ N p , where N p , ( N + 1) p denotes the number of interpolation points for degree N and N + 1 polynomials, respec-tively. Next, let x , x , x denote vectors containing x , x , x coordinates of degree N interpolation pointson a curved physical element D k , and let (cid:101) x , (cid:101) x , (cid:101) x denote their evaluation at degree ( N + 1) interpolationpoints (cid:101) x = V N +1 N x , (cid:101) x = V N +1 N x , (cid:101) x = V N +1 N x . (49)Let (cid:101) D N +1 i denote the nodal differentiation matrix of degree N +1 with respect to the i th coordinate direction.The geometric factors are computed as follows: G k = V NN +1 (cid:16) (cid:101) D N +13 (cid:16)(cid:16) (cid:101) D N +12 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +12 (cid:16)(cid:16) (cid:101) D N +13 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) (50) G k = V NN +1 (cid:16) (cid:101) D N +11 (cid:16)(cid:16) (cid:101) D N +13 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +13 (cid:16)(cid:16) (cid:101) D N +11 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) G k = V NN +1 (cid:16) (cid:101) D N +12 (cid:16)(cid:16) (cid:101) D N +11 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +11 (cid:16)(cid:16) (cid:101) D N +12 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) G k = − V NN +1 (cid:16) (cid:101) D N +13 (cid:16)(cid:16) (cid:101) D N +12 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +12 (cid:16)(cid:16) (cid:101) D N +13 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) G k = − V NN +1 (cid:16) (cid:101) D N +11 (cid:16)(cid:16) (cid:101) D N +13 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +13 (cid:16)(cid:16) (cid:101) D N +11 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) G k = − V NN +1 (cid:16) (cid:101) D N +12 (cid:16)(cid:16) (cid:101) D N +11 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +11 (cid:16)(cid:16) (cid:101) D N +12 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) G k = − V NN +1 (cid:16) (cid:101) D N +13 (cid:16)(cid:16) (cid:101) D N +12 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +12 (cid:16)(cid:16) (cid:101) D N +13 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) G k = − V NN +1 (cid:16) (cid:101) D N +11 (cid:16)(cid:16) (cid:101) D N +13 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +13 (cid:16)(cid:16) (cid:101) D N +11 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) G k = − V NN +1 (cid:16) (cid:101) D N +12 (cid:16)(cid:16) (cid:101) D N +11 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17) − (cid:101) D N +11 (cid:16)(cid:16) (cid:101) D N +12 (cid:101) x (cid:17) ◦ (cid:101) x (cid:17)(cid:17) . Remark.
We note that the discrete GCL (42) can also be enforced directly through a local constrainedminimization problem [53, 17] , which yields a solution in terms of a pseudo-inverse. However, we have notfound a straightforward way to simultaneously enforce both the discrete GCL condition (42) and Assumption 2using this approach. Numerical experiments
In this section, we present numerical results which verify the theoretical results in this work. We firstverify that the weight-adjusted L projection (cid:101) Π kN and the GCL-satisfying geometric factors G kij obey theerror estimates in Theorem 4 and Theorem 5. Next, we verify the semi-discrete entropy conservation,primary conservation, and accuracy of the proposed high order accurate methods for the compressible Eulerequations on curved meshes in two and three dimensions. For all curved meshes, we utilize the low storageweight-adjusted formulation (31).In choosing the timestep dt , we follow [24] and set(51) dt = min k CFL (cid:107) J (cid:107) L ∞ ( D k ) C N (cid:107) J f (cid:107) L ∞ ( ∂D k ) , a) Warped curvilinear mesh − . − − . − − − Mesh size h L e rr o r L projectionWeight-adjustedDifference (b) Smooth function − . − − . − − − .52.5 Mesh size h (c) Discontinuous function
Figure 1. L errors in approximating both smooth and discontinuous functions using L and weight-adjusted projections on a curved mesh. The approximation order is N = 4, anda degree 2 N quadrature rule is used to compute integrals over the reference triangle.where C N is the O ( N ) order-dependent constant in the surface trace inequality [24] and CFL is a user-chosenconstant. All numerical experiments in this work utilize the five-stage fourth order low storage Runge–Kutta(LSRK-45) time-stepper [54].6.1. Accuracy of weight-adjusted projection and geometric terms.
In this section, we verify The-orems 4 and 5 concerning the accuracy of the weight-adjusted projection and modified construction ofgeometric terms satisfying the discrete geometric conservation law. Figure 1 shows L errors for both thestandard L projection (20) and the weight-adjusted projection (28) for a series of warped meshes of degree N = 4. Errors are estimated using degree 2 N + 1 quadratures for triangles and tetrahedra [55]. We compute L errors for both smooth and discontinuous functions f ( x ) = e x + x sin ( πx ) sin ( πx ) , g ( x ) = f ( x ) + H ( x + x − sin( πx )) , where H is the Heaviside function. For the smooth function f ( x ), we observe that the difference betweenthe L and weight-adjusted projections indeed converges at a rate of O (cid:0) h N +2 (cid:1) as predicted by Theorem 1,such that the L errors for each projection appear identical. The L errors for the L and weight-adjustedprojections of the discontinuous function g ( x ) are also virtually identical. However, the difference beteweenthe L and weight-adjusted projection converges faster than estimated by Theorem 1, with the L errorconverging as O ( h / ) and the difference converging as O ( h / ).We next compare the approximation of the geometric factors on a curved three-dimensional mesh. Wegenerate a sequence of quasi-uniform unstructured tetrahedral meshes using GMSH [56] and construct acurvilinear mesh from the distorted coordinates (cid:101) x = x + cos (cid:0) π x (cid:1) cos (cid:0) π x (cid:1) cos (cid:0) π x (cid:1) . We then computethe L error in approximating geometric terms for each element D k by computing G kij − (cid:103) G kij at quadraturepoints. We estimate the mesh size as h = max k (cid:13)(cid:13)(cid:13) J k /J kf (cid:13)(cid:13)(cid:13) L ∞ , since J k = O ( h d ) and J kf = O ( h d − ) in d dimensions [24]. Figure 2 shows errors for an N = 3 and N = 4 mesh. We refer to the construction ofapproximate geometric terms (cid:103) G kij introduced in [2, 57] as “Geo- N ”, since the interpolation is performedusing a degree N interpolation operator. We refer to the construction of (cid:103) G kij in (46) and Theorem 5 as“Geo-( N + 1)”, since the main interpolation step is performed on degree ( N + 1) polynomials instead. It canbe observed that the Geo-N scheme converges at a rate of O ( h N +1 ), while the Geo-(N+1) scheme convergesat a rate of O ( h N +2 ). We note that the error in both the Geo- N and Geo-( N + 1) approximations of (cid:103) G kij converge at the same rate or faster than the best approximation error. − − . − − − Mesh size h L e rr o r Geo-( N + 1) h N +2 Geo- Nh N +1 (a) N = 3 − − . − − − − Mesh size h Geo-( N + 1) h N +2 Geo- Nh N +1 (b) N = 4 Figure 2. L errors in the approximation of metric terms G kij with metric terms (cid:103) G kij satisfying the discrete GCL condition (42) and Assumption 2.6.2. Two-dimensional compressible Euler equations.
The compressible Euler equations in two dimen-sions are given as follows: ∂ρ∂t + ∂ ( ρu ) ∂x + ∂ ( ρv ) ∂x = 0 , (52) ∂ρu∂t + ∂ (cid:0) ρu + p (cid:1) ∂x + ∂ ( ρuv ) ∂x = 0 ,∂ρv∂t + ∂ ( ρuv ) ∂x + ∂ (cid:0) ρv + p (cid:1) ∂x = 0 ,∂E∂t + ∂ ( u ( E + p )) ∂x + ∂ ( v ( E + p )) ∂x = 0 . In two dimensions, the pressure is p = ( γ − (cid:0) E − ρ ( u + v ) (cid:1) , and the specific internal energy is ρe = E − ρ ( u + v ).The choice of convex entropy for the Euler equations is non-unique [58]. However, a unique entropy canbe chosen by restricting to choices of entropy variables which symmetrize the viscous heat conduction termin the compressible Navier-Stokes equations [22]. This leads to U ( u ) of the form(53) U ( u ) = − ρsγ − , where s = log (cid:16) pρ γ (cid:17) is the physical specific entropy. The entropy variables in two dimensions are v = ρe ( γ + 1 − s ) − Eρe , v = ρuρe , v = ρvρe , v = − ρρe . (54)The conservation variables in terms of the entropy variables are given by(55) ρ = − ( ρe ) v , ρu = ( ρe ) v , ρv = ( ρe ) v , E = ( ρe ) (cid:18) − v + v v (cid:19) , where ρe and s in terms of the entropy variables are(56) ρe = (cid:18) ( γ − − v ) γ (cid:19) / ( γ − e − sγ − , s = γ − v + v + v v . a) Uniform mesh (b)
Warped mesh
Figure 3.
2D curved meshes used for testing entropy conservation and primary conservation.The entropy conservative numerical fluxes for the two-dimensional compressible Euler equations are givenby Chandrashekar [37] f ,S ( u L , u R ) = log , f ,S ( u L , u R ) = log , (57) f ,S ( u L , u R ) = f ,S + p avg , f ,S ( u L , u R ) = f ,S ,f ,S ( u L , u R ) = f ,S , f ,S ( u L , u R ) = f ,S + p avg ,f ,S ( u L , u R ) = ( E avg + p avg ) , f ,S ( u L , u R ) = ( E avg + p avg ) , where we have defined the auxiliary quantities p avg = , E avg = log log ( γ −
1) + (cid:107) u (cid:107) , (58) (cid:107) u (cid:107) = 2( + ) − (cid:0)(cid:8) { u } (cid:9) + (cid:8) { v } (cid:9)(cid:1) . Entropy conservation.
We begin by testing the propagation of a shock on a two-dimensional curvedmesh using a discontinuous profile on the domain Ω = [0 , × [ − , ρ ( x , t ) = (cid:40) | x | < / | x | < /
22 otherwise , u ( x , t ) = v ( x , t ) = 0 , p ( x , t ) = ρ γ . To test the scheme (31), we utilize an entropy conservative flux and run it on a uniform triangular mesh witha curvilinear warping shown in Figure 3. Theorem 2 ensures that, under an entropy conservative flux, (31)is semi-discretely entropy conservative. This does not hold at the fully discrete level; however, it is possibleto verify that (31) is entropy stable using other approaches.First, can examine the entropy RHS, which we define as the right hand side of (31) tested with v h (60) entropy RHS = − d (cid:88) j =1 (cid:16)(cid:101) v T (cid:16) Q jk ◦ F j,S (cid:17) + (cid:101) v Tf W kf diag ( n j ) (cid:0) f ∗ j − f j ( (cid:101) u f ) (cid:1)(cid:17) . For positive density and pressure, (60) should be zero to machine precision. We can also track the changein entropy ∆ U = | U ( u ( x , t )) − U ( u ( x , | , which should converge to 0 as the timestep dt approaches zero.Furthermore, the rate of convergence should match the order of the time-stepper used [59, 1].Figure 4 shows the evolution of entropy U ( u ) over time [0 , T ], where the final time T = 2. We comparewhen the entropy-projected conservative variables (cid:101) u are computed using the standard L projection on thereference element, and when (cid:101) u are computed using the weight-adjusted projection. For the standard L projection, the change in entropy does not decrease as dt decreases (for a fixed mesh and order N ). When (cid:101) u is defined using the weight-adjusted projection, the entropy decreases as dt decreases. Moreover, the rateof convergence is approximately O ( dt . ), which is slightly higher than the expected rate of O ( dt ) whenusing LSRK-45. Regardless of whether the standard and weight-adjusted projection was used, the entropyRHS (60) is O (10 − ), indicating that the proposed scheme is implemented correctly. . . − − − − − − Time t C h a n g e i n e n t r o p y ∆ U ( u ) CFL = . . . (a) With weight-adjusted projection . . − − − − − − Time t CFL = . . . (b) Without weight-adjusted projection
Figure 4.
Change in entropy under an entropy conservative formulation with N = 4. Inboth cases, the magnitude of the entropy RHS (60) is O (cid:0) − (cid:1) .These results suggest that the weight-adjusted projection is necessary to produce an entropy conservativescheme on curvilinear meshes. However, entropy conservative schemes result in spurious oscillations andlower convergence rates [1]. In practice, dissipative interface terms are added to produce entropy stableschemes. In the presence of interface dissipation, the evolution of entropy over time, L errors for smoothsolutions, and the qualitative behavior of the solution are very similar with or without the weight-adjustedprojection. This may reflect the fact that aliasing-driven instabilities arise from the spatial discretization,and the entropy RHS (60) is machine precision zero with or without the weight-adjusted projection. Incontrast, the presence of the weight-adjusted projection affects only the time-derivative on the left-handside, and may not play as large a role in suppressing aliasing instabilities.6.2.2. Local and global conservation.
Next, we check that primary conservation is maintained numerically oncurved meshes. We follow [15] and examine the semi-discrete evolution of the average of u over the domainΩ ∂∂t (cid:90) Ω u d x = (cid:88) k (cid:90) D k ∂ u ∂t d x = (cid:88) k (cid:90) (cid:98) D ∂ u ∂t J k d (cid:98) x = (cid:88) k T W diag (cid:0) J k (cid:1) V q d u h dt . (61)The quantity (61) is computed using quadrature by multiplying the semi-discrete system (31) by T W diag (cid:0) J k (cid:1) V q (cid:88) k T W diag (cid:0) J k (cid:1) V q d u h dt = T W diag (cid:0) J k (cid:1) V q (cid:104) (cid:101) P kq (cid:101) L kq (cid:105) d (cid:88) j =1 (cid:16) D jk ◦ F j,S (cid:17) + d (cid:88) j =1 (cid:101) L kq diag ( n j ) (cid:0) f ∗ j − f j ( (cid:101) u f ) (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) conservation residual . (62)We consider three different cases:(1) Unmodified: the weight-adjusted DG scheme (31) is used without any special modifications,(2)
Polynomial approximation: the weight-adjusted DG scheme (31) utilizes the L projection of J k onto P N ( (cid:98) D ) in the weight-adjusted mass matrix (such that Lemma 2 holds),(3) Conservative correction: the conservation correction (41) is applied to the right hand side of thescheme (31).All results use a Lax-Friedrichs flux, a CFL of 1 / N = 4, and the curved mesh warping shown in Figure 3.We use volume quadrature rules from [55] which are exact for degree 2 N +1 polynomials, and use ( N +1)-node1D Gauss–Legendre quadrature rules on the faces. − − − − Time C o n s e r v a t i o n r e s i du a l Unmodified (a)
Conservation residual over time(unmodified) − − Time C o n s e r v a t i o n r e s i du a l Polynomial approx.Conservative correction (b)
Conservation residual over time (modified)
Figure 5.
Conservation residuals for a discontinuous initial condition using the unmodifiedWADG formulation (31) and two different techniques for restoring local conservation. Bothexperiments use a Lax-Friedrichs flux, a CFL of 1 / N = 4, and the curved mesh warpingshown in Figure 3. (a) Affine mesh (b)
Curved mesh
Figure 6.
Example of affine and warped meshes ( h = 1 shown) used in convergence studies.Figure 5a shows computed conservation residuals from (62) on the curved mesh in Figure 3 for the discon-tinuous pulse initial condition (59). The conservation residual oscillates around O (10 − ) for the unmodifiedscheme (31). Modifying the scheme using either the “Polynomial approximation” or “Conservative correc-tion” approaches reduces the conservation error to between O (10 − ) and O (10 − ). We have also studiedthe accuracy of each case in by comparing L errors on a curved mesh for a smooth vortex solution at time T = 5. We observe that in all cases, the L errors are virtually identical. For N = 4 on an 8 × L errors for each case differ only in the 5th digit, while for the 16 ×
16 and 32 ×
32 meshes, they differ onlyin the 8th digit. This confirms that both approaches outlined in Section 4.3.1 restore primary conservationwith no perceivable effect on accuracy.6.2.3.
Accuracy and convergence.
Finally, we test the accuracy of the proposed schemes for smooth solutionson curved meshes in two dimensions. We use the isentropic vortex problem [60], which has an analytical − − − − − − N = 2 3 N = 3 4 N = 4 5 Mesh size h L e rr o r AffineCurved
Figure 7.
Convergence of L errors for the 2D isentropic vortex problem on affine andcurved meshes.solution ρ ( x , t ) = (cid:32) − ( γ − βe − r ( x ,t ) ) γπ (cid:33) γ − , p = ρ γ , (63) u ( x , t ) = 1 − β π e − r ( x ,t ) ( x − c ) , v ( x , t ) = β π e − r ( x ,t ) ( x − c ) , where u, v are the x and x velocity and r ( x , t ) = (cid:112) ( x − c − t ) + ( x − c ) . Here, we take c = 5 , c = 0and β = 5. The solution is computed on a periodic rectangular domain [0 , × [ − ,
5] at final time T = 5.Quasi-uniform triangular meshes are generated using GMSH [56], and a curvilinear warping is applied tothe mesh to test the effect of non-affine mappings. This warping is shown in Figure 6 and is defined bymapping nodal positions on each triangle to warped nodal positions ( (cid:102) x , (cid:102) x ) via (cid:102) x = x + sin ( πx /
20) sin (2 π ( x + 5) / (cid:102) x = x −
12 sin (2 πx /
20) sin ( π ( x + 5) / . To ensure primary conservation, we use the “Polynomial Approximation” strategy described in Section 4.3.1and compute the inverse of the weight-adjusted mass matrix using the degree
N L projection (cid:98) Π N J k insteadof J k . The computed L errors are shown in Figure 7. We observe optimal O ( h N +1 ) rates of convergencefor N = 2 , N = 4. For degree N = 3, the rate of convergence is slightly higher than O ( h ) and may indicatethat the mesh is not yet sufficiently fine for the L error to show the asymptotic convergence rate.6.3. Three dimensional compressible Euler equations.
In three dimensions, the compressible Eulerequations are given by ∂∂t ρρuρvρwE + ∂∂x ρuρu + pρuvρuwu ( E + p ) + ∂∂x ρvρuvρv + pρvwv ( E + p ) + ∂∂x ρwρuwρvwρw + pw ( E + p ) = 0 , (64)where the pressure p and specific internal energy ρe are defined p = ( γ − (cid:18) E − ρ ( u + v + w ) (cid:19) , ρe = E − ρ ( u + v + w ) . (65) he formula for the entropy U ( u ) in three dimensions is the same as the two-dimensional formula (53). Theentropy variables in three dimensions are v = ρe ( γ + 1 − s ) − Eρe , v = ρuρe , v = ρvρe , v = ρwρe , v = − ρρe . (66)The conservation variables in terms of the entropy variables are given by(67) ρ = − ( ρe ) v , ρu = ( ρe ) v , ρv = ( ρe ) v , ρw = ( ρe ) v , E = ( ρe ) (cid:18) − v + v + v v (cid:19) , where ρe and s in terms of the entropy variables are(68) ρe = (cid:18) ( γ − − v ) γ (cid:19) / ( γ − e − sγ − , s = γ − v + v + v + v v . A set of entropy conservative numerical fluxes for the three-dimensional compressible Euler equations canbe written as f ,S = log log + p avg log log ( E avg + p avg ) , f ,S = log log log + p avg log ( E avg + p avg ) , (69) f ,S = log log log log + p avg ( E avg + p avg ) . where we have defined the auxiliary quantities p avg = , E avg = log γ − log + 12 log (cid:107) u (cid:107) (70) (cid:107) u (cid:107) = 2( + + ) − (cid:0)(cid:8) { u } (cid:9) + (cid:8) { v } (cid:9) + (cid:8) { w } (cid:9)(cid:1) . Accuracy and convergence.
As before, we test the accuracy of the proposed scheme using an isentropicvortex solution adapted to three dimensions. We take the solution to be the extruded 2D vortex propagatingin the x direction, whose analytic expression is derived from [61] ρ ( x , t ) = (cid:18) − ( γ − (cid:19) γ − u ( x , t ) = Π r ,E ( x , t ) = p γ − (cid:18) − γ −
12 Π (cid:19) γγ − + ρ | u | . where u = ( u, v, w ) T is the velocity vector andΠ = Π max e − r T r , r = − ( x − c − t ) x − c . In this problem, we take c = c = 5, p = 1 /γ , and Π max = 0 .
4. The problem is solved on the domain[0 , × [0 , × [0 , GMSH to construct three unstructured meshes consisting of 1354, 9543, a) Mesh for h = 1 / − . − . − . . − − − − N = 23 N = 34 N = 4 5 Mesh size h L e rr o r AffineCurved (b) L errors Figure 8.
Convergence of L errors for the 3D isentropic vortex problem on unstructuredaffine and curved meshes, with optimal O ( h N +1 ) rates of convergence shown for reference.and 72923 affine tetrahedra, corresponding to h = 2, h = 1, and h = 1 / (cid:102) x , (cid:102) x , (cid:102) x ) via (cid:102) x = x + 12 sin (cid:16) π x (cid:17) sin (cid:16) π x (cid:17) sin (cid:16) π x (cid:17) , (cid:102) x = x − sin (cid:16) π x (cid:17) sin (cid:16) π x (cid:17) sin (cid:16) π x (cid:17) , (cid:102) x = x + 12 sin (cid:16) π x (cid:17) sin (cid:16) π x (cid:17) sin (cid:16) π x (cid:17) . The geometric terms are approximated using (46) to ensure the satisfaction of the discrete GCL. Figure 8shows L errors at final time T = 5 for N = 2 , , O ( h N +1 ) asymptotic rates of convergence for N = 2 ,
3, while for N = 4 we observe a ratewhich is slightly higher than O ( h N +1 / ) but not quite O ( h N +1 ). This may be due to the fact that thetime-stepper is 4th order while the spatial discretization is 5th order.We also compared approximations of geometric factors using degree N polynomials via (46) and degree N − N − Inviscid Taylor–Green vortex.
Our last numerical experiment investigates the behavior of entropystable DG schemes for the inviscid Taylor–Green vortex [63, 35, 17]. The domain is the periodic box[ − π, π ] , and the initial conditions are given as ρ = 1 u = sin( x ) cos( x ) cos( x ) ,v = − cos( x ) sin( x ) cos( x ) ,w = 0 ,p = 100 γ + 116 (cos(2 x ) + cos(2 x )) (2 + cos(2 x )) . t − ∂ κ ∂ t AffineCurved (a)
KE dissipation rate for N = 3, h = π/ − . − . − − . − . − − − − CFL A v e r ag ee n t r o p y AffineCurved (b)
Entropy conservative convergence of (cid:82) Ω U ( u ) Figure 9.
Evolution of the kinetic energy dissipation rate over time on affine and curvi-linear meshes, as well as dependence of average entropy over the domain (cid:82) Ω U ( u ) at time T = 20 for an entropy conservative formulation.The Taylor–Green vortex is used to study the transition and decay of turbulence [64]. In the absence ofviscosity, the Taylor–Green vortex develops smaller and smaller scales, implying that for sufficiently largetimes, the solution contains under-resolved features. We study the evolution of the kinetic energy κ ( t ) κ ( t ) = 1 | Ω | (cid:90) Ω ρ u · u d x , as well as the kinetic energy dissipation rate − ∂κ∂t , which we approximate by differencing κ ( t ). Figure 9shows the evolution of κ ( t ) over time for both affine and curvilinear meshes with h = π/ N = 3.We also use a curved coarse mesh with element size h = π test the convergence of the average entropy fora non-dissipative entropy conservative formulation as the timestep decreases. In both cases, the mesh isdefined by constructing nodal positions (cid:101) u from warpings of affine nodal positions x as follows (cid:101) x = x + .
125 sin( x ) sin( x ) sin( x ) . The kinetic energy is plotted at 100 equally spaced times between [0 , t = 4 andpeaking with a value of roughly 0 .
014 before t = 9 [64, 35]. We also observe that, for an entropy conservativeformulation, the average entropy appears to converge at a rate of O ( dt ) for a 4th order time-steppingmethod. The same phenomena was also observed in [1, 17].7. Conclusions
This work describes how to extend entropy conservative and entropy stable DG “modal” discretizationsto curvilinear meshes using both weighted and weight-adjusted mass matrices. Assuming that the geometricterms satisfy a discrete geometric conservation law, the presented schemes allow for the use of over-integrationwhile satisfying a semi-discrete entropy equality or inequality on curved meshes. Numerical results showthe presented schemes achieve optimal rates of convergence for smooth solutions on both two and threedimensional affine and curvilinear meshes while remaining robust in the presence of under-resolved solutionssuch as shocks and turbulence.Several outstanding computational questions remain to be answered. First, the use of weight-adjustedmass matrices is motivated by the low storage requirements and their efficient application on GPUs. However,to guarantee discrete entropy stability, the computation of a weight-adjusted projection is necessary, which dds additional computational cost compared to the affine case. Furthermore, numerical results suggestthat, while computing the entropy-projected conservative variables (cid:101) u using the weight-adjusted projection isnecessary to ensure a semi-discrete conservation of entropy, the difference between using a weight-adjustedprojection and regular projection may be negligible in practice. The necessity of the weight-adjusted pro-jection in computing (cid:101) u should be examined further. Secondly, on triangular and tetrahedral meshes, theproposed schemes involve more computational work than under-integrated collocation-style SBP schemes[19, 16, 17]. A careful computational comparison of the presented schemes with existing methods should bedone to weigh the benefits of improved accuracy with additional computational costs.Finally, we note that while this work has focused on triangular and tetrahedral elements, the approachesoutlined here can be extended to more arbitrary pairings of (possibly non-polynomial) approximation spacesand quadratures. For example, similar techniques can be used to construct entropy stable B-spline orGalerkin difference discretizations on curved meshes [65, 39].8. Acknowledgements
The authors thank Mark H. Carpenter and David C. Del Rey Fernandez for informative discussions. JesseChan is supported by NSF DMS-1719818 and DMS-1712639.
References [1] Jesse Chan. On discretely entropy conservative and entropy stable discontinuous Galerkin methods.
Journal of Computa-tional Physics , 362:346–374, June 2018. doi:10.1016/j.jcp.2018.02.033 .[2] David A Kopriva. Metric identities and the discontinuous spectral element method on curvilinear meshes.
Journal ofScientific Computing , 26(3):301–327, March 2006. doi:10.1007/s10915-005-9070-8 .[3] Jesse Chan, Russell J Hewett, and T Warburton. Weight-adjusted discontinuous Galerkin methods: Wave propagation inheterogeneous media.
SIAM Journal on Scientific Computing , 39(6):A2935–A2961, 2017. doi:10.1137/16m1089186 .[4] Jesse Chan, Russell J Hewett, and T Warburton. Weight-adjusted discontinuous Galerkin methods: Curvilinear meshes.
SIAM Journal on Scientific Computing , 39(6):A2395–A2421, 2017. doi:10.1137/16m1089198 .[5] Fang Q Hu, MY Hussaini, and Patrick Rasetarinera. An analysis of the discontinuous Galerkin method for wave propagationproblems.
Journal of Computational Physics , 151(2):921–946, May 1999. doi:10.1006/jcph.1999.6227 .[6] Mark Ainsworth. Dispersive and dissipative behaviour of high order discontinuous Galerkin finite element methods.
Journalof Computational Physics , 198(1):106–130, July 2004. doi:10.1016/j.jcp.2004.01.004 .[7] Zhijian J. Wang, Krzysztof Fidkowski, R´emi Abgrall, Francesco Bassi, Doru Caraeni, Andrew Cary, Herman Deconinck,Ralf Hartmann, Koen Hillewaert, Hung T Huynh, Norbert Kroll, Georg May, Per-Olof Persson, Bram van Leer, and MiguelVisbal. High-order CFD methods: Current status and perspective.
International Journal for Numerical Methods in Fluids ,72(8):811–845, January 2013. doi:10.1002/fld.3767 .[8] Per-Olof Persson and Jaime Peraire. Sub-cell shock capturing for discontinuous Galerkin methods. In . American Institute of Aeronautics and Astronautics, January 2006. doi:10.2514/6.2006-112 .[9] Lilia Krivodonova. Limiters for high-order discontinuous Galerkin methods.
Journal of Computational Physics , 226(1):879–896, September 2007. doi:10.1016/j.jcp.2007.05.011 .[10] Garrett E Barter and David L Darmofal. Shock capturing with PDE-based artificial viscosity for DGFEM: Part I. Formu-lation.
Journal of Computational Physics , 229(5):1810–1827, March 2010. doi:10.1016/j.jcp.2009.11.010 .[11] Jean-Luc Guermond, Richard Pasquetti, and Bojan Popov. Entropy viscosity method for nonlinear conservation laws.
Journal of Computational Physics , 230(11):4248–4267, May 2011. doi:10.1016/j.jcp.2010.11.043 .[12] Travis C Fisher and Mark H Carpenter. High-order entropy stable finite difference schemes for nonlinear conservation laws:Finite domains.
Journal of Computational Physics , 252:518–557, November 2013. doi:10.1016/j.jcp.2013.06.014 .[13] Mark H Carpenter, Travis C Fisher, Eric J Nielsen, and Steven H Frankel. Entropy stable spectral collocation schemes forthe Navier–Stokes equations: Discontinuous interfaces.
SIAM Journal on Scientific Computing , 36(5):B835–B867, 2014. doi:10.1137/130932193 .[14] Matteo Parsani, Mark H Carpenter, Travis C Fisher, and Eric J Nielsen. Entropy stable staggered grid discontinuousspectral collocation methods of any order for the compressible Navier–Stokes equations.
SIAM Journal on ScientificComputing , 38(5):A3129–A3162, 2016. doi:10.1137/15m1043510 .[15] Lucas Friedrich, Andrew R Winters, David C Del Rey Fern´andez, Gregor J Gassner, Matteo Parsani, and Mark H Carpen-ter. An entropy stable h/p non-conforming discontinuous Galerkin method with the summation-by-parts property. arXivpreprint arXiv:1712.10234 , 2017. URL: https://arxiv.org/abs/1712.10234 .[16] Tianheng Chen and Chi-Wang Shu. Entropy stable high order discontinuous Galerkin methods with suitable quadraturerules for hyperbolic conservation laws.
Journal of Computational Physics , 345:427–461, September 2017. doi:10.1016/j.jcp.2017.05.025 .
17] Jared Crean, Jason E Hicken, David C Del Rey Fern´andez, David W Zingg, and Mark H Carpenter. Entropy-stablesummation-by-parts discretization of the Euler equations on general curved elements.
Journal of Computational Physics ,356:410–438, March 2018. doi:10.1016/j.jcp.2017.12.015 .[18] Brian T Helenbrook. On the existence of explicit hp -finite element methods using Gauss–Lobatto integration on the triangle. SIAM Journal on Numerical Analysis , 47(2):1304–1318, 2009. doi:10.1137/070685439 .[19] Jason E Hicken, David C Del Rey Fern´andez, and David W Zingg. Multidimensional summation-by-parts operators:General theory and application to simplex elements.
SIAM Journal on Scientific Computing , 38(4):A1935–A1958, 2016. doi:10.1137/15m1038360 .[20] PD Thomas and CK Lombard. Geometric conservation law and its application to flow computations on moving grids.
AIAA Journal , 17(10):1030–1037, October 1979. doi:10.2514/3.61273 .[21] Jesse Chan. Weight-adjusted discontinuous Galerkin methods: Matrix-valued weights and elastic wave propagation inheterogeneous media.
International Journal for Numerical Methods in Engineering , 113(12):1779–1809, March 2018. doi:10.1002/nme.5720 .[22] Thomas JR Hughes, LP Franca, and M Mallet. 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.
ComputerMethods in Applied Mechanics and Engineering , 54(2):223–234, February 1986. doi:10.1016/0045-7825(86)90127-1 .[23] Michael S Mock. Systems of conservation laws of mixed type.
Journal of Differential Equations , 37(1):70–88, July 1980. doi:10.1016/0022-0396(80)90089-3 .[24] Jesse Chan, Zheng Wang, Axel Modave, Jean-Francois Remacle, and T Warburton. GPU-accelerated discontinuousGalerkin methods on hybrid meshes.
Journal of Computational Physics , 318:142–168, August 2016. doi:10.1016/j.jcp.2016.04.003 .[25] Jan S Hesthaven and Tim Warburton.
Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications .Springer New York, 2008. doi:10.1007/978-0-387-72067-8 .[26] Daniele Antonio Di Pietro and Alexandre Ern.
Mathematical Aspects of Discontinuous Galerkin Methods . Springer BerlinHeidelberg, 2012. doi:10.1007/978-3-642-22980-0 .[27] Bernardo Cockburn and Chi-Wang Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element methodfor conservation laws. II. General framework.
Mathematics of Computation , 52(186):411–435, April 1989. doi:10.1090/s0025-5718-1989-0983311-4 .[28] Bernardo Cockburn and Chi-Wang Shu. The Runge–Kutta discontinuous Galerkin method for conservation laws V: Mul-tidimensional systems.
Journal of Computational Physics , 141(2):199–224, April 1998. doi:10.1006/jcph.1998.5892 .[29] Bernardo Cockburn. Devising discontinuous Galerkin methods for non-linear hyperbolic conservation laws.
Journal ofComputational and Applied Mathematics , 128(1-2):187–204, March 2001. doi:10.1016/s0377-0427(00)00512-4 .[30] Jianxian Qiu, Boo Cheong Khoo, and Chi-Wang Shu. A numerical study for the performance of the Runge–Kutta discon-tinuous Galerkin method based on different numerical fluxes.
Journal of Computational Physics , 212(2):540–565, March2006. doi:10.1016/j.jcp.2005.07.011 .[31] Xiangxiong Zhang, Yinhua Xia, and Chi-Wang Shu. Maximum-principle-satisfying and positivity-preserving high orderdiscontinuous Galerkin schemes for conservation laws on triangular meshes.
Journal of Scientific Computing , 50(1):29–62,January 2012. doi:10.1007/s10915-011-9472-8 .[32] A Kl¨ockner, T Warburton, and JS Hesthaven. Viscous shock capturing in a time-explicit discontinuous Galerkin method.
Mathematical Modelling of Natural Phenomena , 6(3):57–83, 2011. doi:10.1051/mmnp/20116303 .[33] Philippe G Ciarlet.
The finite element method for elliptic problems . North-Holland Publishing Company, 1978.[34] Jan Nordstr¨om. Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation.
Journal of Scientific Computing , 29(3):375–404, December 2006. doi:10.1007/s10915-005-9013-4 .[35] Gregor J Gassner, Andrew R Winters, and David A Kopriva. Split form nodal discontinuous Galerkin schemes withsummation-by-parts property for the compressible Euler equations.
Journal of Computational Physics , 327:39–66, Decem-ber 2016. doi:10.1016/j.jcp.2016.09.013 .[36] Roger A Horn and Charles R Johnson.
Matrix Analysis . Cambridge University Press, 2nd edition, 2012.[37] Praveen Chandrashekar. Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler andNavier-Stokes equations.
Communications in Computational Physics , 14(5):1252–1286, November 2013. doi:10.4208/cicp.170712.010313a .[38] Andrew R Winters, Dominik Derigs, Gregor J Gassner, and Stefanie Walch. A uniquely defined entropy stable matrixdissipation operator for high Mach number ideal MHD and compressible Euler simulations.
Journal of ComputationalPhysics , 332:274–289, March 2017. doi:10.1016/j.jcp.2016.12.006 .[39] Jesse Chan and John A Evans. Multi-patch discontinuous Galerkin isogeometric analysis for wave propagation: Explicittime-stepping and efficient mass matrix inversion.
Computer Methods in Applied Mechanics and Engineering , 333:22–54,May 2018. doi:10.1016/j.cma.2018.01.022 .[40] T Warburton. A low-storage curvilinear discontinuous Galerkin method for wave problems.
SIAM Journal on ScientificComputing , 35(4):A1987–A2012, 2013. doi:10.1137/120899662 .[41] Susanne C Brenner and L Ridgway Scott.
The Mathematical Theory of Finite Element Methods . Springer New York, 3rdedition, 2008. doi:10.1007/978-0-387-75934-0 .[42] Mo Lenoir. Optimal isoparametric finite elements and error estimates for domains involving curved boundaries.
SIAMJournal on Numerical Analysis , 23(3):562–580, 1986. doi:10.1137/0723036 .
43] Travis C Fisher, Mark H Carpenter, Jan Nordstr¨om, Nail K Yamaleev, and Charles Swanson. Discretely conservativefinite-difference formulations for nonlinear conservation laws in split form: Theory and boundary conditions.
Journal ofComputational Physics , 234:353–375, February 2013. doi:10.1016/j.jcp.2012.09.026 .[44] A Johnen, J-F Remacle, and C Geuzaine. Geometrical validity of curvilinear finite elements.
Journal of ComputationalPhysics , 233:359–372, January 2013. doi:10.1016/j.jcp.2012.08.051 .[45] Miguel R Visbal and Datta V Gaitonde. On the use of higher-order finite-difference schemes on curvilinear and deformingmeshes.
Journal of Computational Physics , 181(1):155–185, September 2002. doi:10.1006/jcph.2002.7117 .[46] Qi Chen and Ivo Babuˇska. The optimal symmetrical points for polynomial interpolation of real functions in the tetrahe-dron.
Computer Methods in Applied Mechanics and Engineering , 137(1):89–94, October 1996. doi:10.1016/0045-7825(96)01051-1 .[47] JS Hesthaven. From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex.
SIAM Journalon Numerical Analysis , 35(2):655–676, 1998. doi:10.1137/s003614299630587x .[48] MG Blyth and C Pozrikidis. A Lobatto interpolation grid over the triangle.
IMA Journal of Applied Mathematics ,71(1):153–169, February 2006. doi:10.1093/imamat/hxh077 .[49] T Warburton. An explicit construction of interpolation nodes on the simplex.
Journal of Engineering Mathematics ,56(3):247–262, November 2006. doi:10.1007/s10665-006-9086-6 .[50] Jesse Chan and T Warburton. A comparison of high order interpolation nodes for the pyramid.
SIAM Journal on ScientificComputing , 37(5):A2151–A2170, 2015. doi:10.1137/141000105 .[51] DM Williams, Lee Shunn, and Antony Jameson. Symmetric quadrature rules for simplexes based on sphere close packedlattice arrangements.
Journal of Computational and Applied Mathematics , 266:18–38, August 2014. doi:10.1016/j.cam.2014.01.007 .[52] Freddie D Witherden and Peter E Vincent. On the identification of symmetric quadrature rules for finite element methods.
Computers & Mathematics with Applications , 69(10):1232–1241, May 2015. doi:10.1016/j.camwa.2015.03.017 .[53] David C Del Rey Fern´andez, Jason E Hicken, and David W Zingg. Simultaneous approximation terms for multi-dimensional summation-by-parts operators.
Journal of Scientific Computing , 75(1):83–110, April 2018. doi:10.1007/s10915-017-0523-7 .[54] Mark H Carpenter and Christopher A Kennedy. Fourth-order 2 N -storage Runge-Kutta schemes. Technical Report NASA-TM-109112, NAS 1.15:109112, NASA Langley Research Center, 1994. doi:2060/19940028444 .[55] Hong Xiao and Zydrunas Gimbutas. A numerical algorithm for the construction of efficient quadrature rules in two andhigher dimensions. Computers & Mathematics with Applications , 59(2):663–676, January 2010. doi:10.1016/j.camwa.2009.10.027 .[56] Christophe Geuzaine and Jean-Fran¸cois Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities.
International Journal for Numerical Methods in Engineering , 79(11):1309–1331, September 2009. doi:10.1002/nme.2579 .[57] Florian Hindenlang, Gregor J Gassner, Christoph Altmann, Andrea Beck, Marc Staudenmaier, and Claus-Dieter Munz.Explicit discontinuous Galerkin methods for unsteady problems.
Computers & Fluids , 61:86–93, May 2012. doi:10.1016/j.compfluid.2012.03.006 .[58] Amiram Harten. On the symmetric form of systems of conservation laws with entropy.
Journal of Computational Physics ,49(1):151–164, January 1983. doi:10.1016/0021-9991(83)90118-3 .[59] Gregor J Gassner, Andrew R Winters, and David A Kopriva. A well balanced and entropy conservative discontinuousGalerkin spectral element method for the shallow water equations.
Applied Mathematics and Computation , 272:291–308,January 2016. doi:10.1016/j.amc.2015.07.014 .[60] Chi-Wang Shu. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservationlaws. In
Advanced Numerical Approximation of Nonlinear Hyperbolic Equations , pages 325–432. Springer Berlin Heidelberg,1998. doi:10.1007/bfb0096355 .[61] David M Williams and Antony Jameson. Nodal points and the nonlinear stability of high-order methods for unsteadyflow problems on tetrahedral meshes. In . American Institute ofAeronautics and Astronautics, June 2013. doi:10.2514/6.2013-2830 .[62] Thomas Toulorge, Jonathan Lambrechts, and Jean-Fran¸cois Remacle. Optimizing the geometrical accuracy of curvilinearmeshes.
Journal of Computational Physics , 310:361–380, April 2016. doi:10.1016/j.jcp.2016.01.023 .[63] Geoffrey I Taylor and Albert E Green. Mechanism of the production of small eddies from large ones.
Proceedings of theRoyal Society A: Mathematical, Physical and Engineering Sciences , 158(895):499–521, February 1937. doi:10.1098/rspa.1937.0036 .[64] James DeBonis. Solutions of the Taylor-Green vortex problem using high-resolution explicit finite difference methods. In . American Instituteof Aeronautics and Astronautics, January 2013. doi:10.2514/6.2013-382 .[65] Jeffrey W Banks and Thomas Hagstrom. On Galerkin difference methods.
Journal of Computational Physics , 313:310–327,2016., 313:310–327,2016.