Efficient Entropy Stable Gauss Collocation Methods

EFFICIENT ENTROPY STABLE GAUSS COLLOCATION METHODS arXiv:1809.01178v2 [math.NA] 3 Aug 2019 JESSE CHAN, DAVID C. DEL REY FERNANDEZ, MARK H. CARPENTER Abstract. The construction of high order entropy stable collocation schemes on quadrilateral and hexahedral elements has relied on the use of Gauss-Legendre-Lobatto collocation points [1, 2, 3] and their equivalence with summation-by-parts (SBP) finite difference operators [4]. In this work, we show how to efficiently generalize the construction of semi-discretely entropy stable schemes on tensor product elements to Gauss points and generalized SBP operators. Numerical experiments suggest that the use of Gauss points significantly improves accuracy on curved meshes. 1. Introduction. Time dependent nonlinear conservation laws are ubiquitous in computational fluid dynamics, for which high order methods are increasingly of interest. Such methods are more accurate per degree of freedom than low order methods, while also possessing much smaller numerical dispersion and dissipation errors. This makes high order methods especially well suited to time-dependent simulations. In this work, we focus specifically on discontinuous Galerkin methods on unstructured quadrilateral and hexahedral meshes. These methods combine properties of high order approximations with the geometric flexibility of unstructured meshing. However, high order methods are notorious for being more prone to instability compared to low order methods [5]. This instability is addressed through various stabilization techniques (e.g. artificial viscosity, filtering, slope limiting). However, these techniques often reduce accuracy to first or second order, and can prevent solvers from realizing the advantages of high order approximations. Moreover, it is often not possible to prove that a high order scheme does not blow up even in the presence of stabilization. This ambiguity can necessitate the re-tuning of stabilization parameters, as parameters which are both stable and accurate for one problem or discretization setting may provide either too little or too much numerical dissipation for another. The instability of high order methods is rooted in the fact that discretizations of nonlinear conservation laws do not typically satisfy a discrete analogue of the conservation or dissipation of energy (entropy). For low order methods, the lack of discrete stability can be offset by the presence of numerical dissipation, which serves as a stabilization mechanism. Because high order methods possess low numerical dissipation, the absence of a discrete stability property becomes more noticeable, manifesting itself through increased sensitivity and instability. The dissipation of entropy serves as an energetic principle for nonlinear conservation laws [6], and requires the use of the chain rule in its proof. Discrete instability is typically tied to the fact that, when discretizing systems of nonlinear PDEs, the chain rule does not typically hold at the discrete level. The lack of a chain rule was circumvented by using a non-standard “flux differencing” formulation [1, 2, 3, 7], which is key to constructing semi-discretely entropy stable high order schemes on unstructured quadrilateral and hexahedral meshes. Flux differencing replaces the derivative of the nonlinear flux with the discrete differentiation of an auxiliary quantity. This auxiliary quantity is computed through the evaluation of a two-point entropy conservative flux [8] using pairs of solution values at quadrature points. These entropy stable schemes were later extended to non-tensor product elements using GLL-like quadrature points on triangles and tetrahedra [9, 10]. More recently, the construction of efficient entropy stable schemes was extended to more arbitrary choices of basis and quadrature [11, 12]. While entropy stable collocation schemes have been constructed on Gauss-like 1 quadrature points without boundary nodes [13], the inter-element coupling terms for such schemes introduce an “all-to-all” coupling between degrees of freedom on two neighboring elements in one dimension (on tensor product elements in higher dimensions, these coupling terms couple together lines of nodes). These coupling terms require evaluating two-point fluxes between solution states at all collocation nodes on two neighboring elements, resulting in significantly more communication and computation compared to collocation schemes based on point sets containing boundary nodes. This work introduces efficient and entropy stable inter-element coupling terms for Gauss collocation schemes which require only communication of face values between neighboring elements. The construction of these terms follows the framework introduced in [11, 12] for triangles and tetrahedra. The main motivation for exploring tensor product (quadrilateral and hexahedral) elements is the significant reduction in the number of operations required compared to high order entropy stable schemes on simplicial meshes [11, 12]. Entropy stability schemes on simplicial elements require evaluating two-point fluxes between solution states at all quadrature points on an element. For a degree N approximation, the number of quadrature points on a simplex scales as O(N d ) in d dimensions, and results in O(N 2d ) two-point flux evaluations per element. In contrast, entropy stable schemes on quadrilateral and hexahedral elements require only the evaluation of two-point fluxes along lines of nodes due to a tensor product structure, resulting in O(N d+1 ) evaluations in d dimensions. In Section 2, we briefly review the derivation of continuous entropy inequalities for systems of nonlinear conservation laws. In Section 3, we describe how to construct entropy stable discretizations of nonlinear conservation laws using different quadrature points on affine tensor product elements. In Section 4, we describe how to extend this construction to curvilinear elements, and Section 5 presents numerical results which confirm the high order accuracy and stability of the proposed method for smooth, discontinuous, and under-resolved (turbulent) solutions of the compressible Euler equations in two and three dimensions. 2. A brief review of entropy stability theory. We are interested in methods for the numerical solution of systems of conservation laws in d dimensions d ∂u X ∂fi (u) = 0, + ∂t ∂xi i=1 (1) where u denotes the conservative variables, fi (u) are nonlinear fluxes, and xi denotes the ith coordinate. Many physically motivated conservation laws admit a statement of stability involving a convex scalar entropy S(u). We first define the entropy variables v(u) to be the gradient of the entropy S(u) with respect to the conservative variables ∂S(u) . ∂u For a convex entropy, v(u) defines an invertible mapping from conservative to entropy variables. We denote the inverse of this mapping (from entropy to conservative variables) by u(v). At the continuous level, it can be shown (for example, in [6]) that vanishing viscosity solutions to (1) satisfy the strong form of an entropy inequality def v= d (2) ∂S(u) X ∂Fi (u) ≤ 0, + ∂t ∂xi i=1 2 Fi (u) = v T ∂fi , ∂u where Fi denotes the ith scalar entropy flux function. Integrating (2) over a domain Ω and applying the divergence theorem yields an integrated entropy inequality (3) Z Ω ∂S(u) + ∂t Z d X ∂Ω i=1
ni v T fi (u) − ψi (u) ≤ 0, where ψi (u) = v T fi (u) − Fi (u) denotes the ith entropy potential, ∂Ω denotes the boundary of Ω and ni denotes the ith component of the outward normal on ∂Ω. Roughly speaking, this implies that the time rate of change of entropy is less than or equal to the flux of entropy through the boundary. 3. Entropy stable Gauss and Gauss-Legendre-Lobatto collocation methods. The focus of this paper is on entropy stable high order collocation methods which satisfy a semi-discrete version of the entropy inequality (3). These methods collocate the solution at some choice of collocation nodes, and are applicable to tensor product meshes consisting of quadrilateral and hexahedral elements. Entropy stable collocation methods have largely utilized Gauss-Legendre-Lobatto (GLL) nodes [1, 2, 3, 7], which contain points on the boundary. The popularity of GLL nodes can be attributed in part to a connection made in [4], where it was shown by Gassner that collocation DG discretizations based on GLL nodes could be recast in terms of summation-by-parts (SBP) operators. This equivalence allowed Gassner to leverage existing finite difference formulations to produce stable high order discretizations of the nonlinear Burgers’ equation. GLL quadratures contain boundary points, which greatly simplifies the construction of inter-element coupling terms for entropy stable collocation schemes. However, it is also known that the use of GLL quadrature within DG methods under-integrates the mass matrix, which can lead to solution “aliasing” and lower accuracy [14]. In this work, we explore entropy stable collocation schemes based on Gauss quadrature points instead of GLL points. This comparison is motivated by the accuracy of each respective quadrature rule. While (N + 1)-point GLL quadrature rules are exact for polynomial integrands of degree (2N − 1), (N + 1)-point Gauss quadrature rules are exact for polynomials of degree (2N + 1). This higher accuracy of Gauss quadrature has been shown to translate to lower errors and slightly improved rates of convergence in simulations of wave propagation and fluid flow [15, 16, 17]. However, Gauss points have not been widely used to construct entropy stable discretizations due to the lack of efficient, stable, and high order accurate inter-element coupling terms, known as simultaneous approximation terms (SAT) in the finite difference literature [18, 13, 19]. SATs for Gauss points are non-compact, in the sense that they introduce all-to-all coupling between degrees of freedom on neighboring elements in one dimension. This results in greater communication between elements, as well as a significantly larger number of two-point flux evaluations and floating point operations. It is possible to realize the improved accuracy of Gauss points while avoiding non-compact SATs through a staggered grid formulation, where the solution is stored at Gauss nodes but interpolated to a set of higher degree (N + 2) GLL “flux” points for computation [14]. Because GLL nodes include boundary points, compact and high order accurate SAT terms can be constructed for the flux points. After performing computations on the flux points, the results are interpolated back to Gauss points and used to evolve the solution forward in time. Figure 1 shows an illustration of GLL, staggered grid, and Gauss point sets for a 2D quadrilateral element. 3 (a) GLL nodes (b) Staggered grid nodes (c) Gauss nodes Fig. 1: Examples of nodal sets under which efficient entropy stable schemes can be constructed. This work focuses on the construction of efficient and accurate SAT terms for Gauss nodal sets. The following sections will describe how to construct efficient high order entropy stable schemes using Gauss points. These schemes are based on “decoupled” SBP operators introduced in [11, 12], which are applicable to general choices of basis and quadrature. By choosing a tensor product Lagrange polynomial basis and (N + 1) point Gauss quadrature rules, we recover a Gauss collocation scheme. The high order accuracy and entropy stability of this scheme are direct results of theorems presented in [11, 12]. However, we will also present a different proof of entropy stability in one dimension for completeness. 3.1. Gauss nodes and generalized summation by parts operators. We assume the solution is collocated at (N + 1) quadrature points xi with associated quadrature weights wi . We do not make any assumptions on the points, in order to accommodate both GLL and Gauss nodes using this notation. The collocation assumption is equivalent to approximating the solution using a degree N Lagrange basis ℓj (x) defined over the (N + 1) quadrature points. Let D denote the nodal differentiation matrix, and let Vf denote the 2 × (N + 1) matrix which interpolates polynomials at Gauss nodes to values at endpoints. These two matrices are defined entrywise as Dij = ∂ℓj ∂x , x=xi (Vf )1i = ℓi (−1), (Vf )2i = ℓi (1). We also introduce the diagonal matrix of quadrature weights Wij = δij wi , as well as the one-dimensional mass matrix M whose entries are L2 inner products of basis functions. We assume that these inner products are computed using the quadrature rule (xi , wi ) at which the solution is collocated. Under such an assumption, the mass matrix is diagonal with entries equal to the quadrature weights Mij = Z 1 −1 ℓi (x)ℓj (x) ≈ N +1 X ℓi (xk )ℓj (xk )wk = δij wi = Wij . k=1 Since M = W under a collocation assumption, we utilize W for the remainder of this work to emphasize that the mass matrix is diagonal and related to the quadrature weights wi . The treatment of non-diagonal mass matrices is covered in [11, 12]. 4 It can be shown that the mass and differentiation matrices for Gauss nodes fall under the class of generalized SBP (GSBP) operators [20]. Lemma 1. Q = W D satisfies the generalized summation by parts property   −1 T T Q = Vf BVf − Q , B= . 1 The proof is a direct restatement of integration by parts, and can be found in [20, 21, 22, 23]. Lemma 1 holds for both GLL and Gauss nodes, and switching between these two nodal sets simply results in a redefinition of the matrices D, Vf . For example, because GLL nodes include boundary points, the interpolation matrix Vf reduces to a generalized permutation matrix which extracts the nodal values associated with the left and right endpoints. 3.2. Existing entropy stable SATs for generalized SBP operators. In this section, we will review the construction of semi-discretely entropy stable discretizations. Entropy stable discretizations can be constructed by first introducing an entropy conservative scheme, then adding appropriate interface dissipation to produce an entropy inequality. The construction of entropy conservative schemes relies on the existence of an two-point (dyadic) entropy conservative flux [8]. Definition 2. Let fS (uL , uR ) be a bivariate function which is symmetric and consistent with the flux function f (u) fS (uL , uR ) = fS (uR , uL ), fS (u, u) = f (u) The numerical flux fS (uL , uR ) is entropy conservative if, for entropy variables vL = v(uL ), vR = v(uR ), the Tadmor “shuffle” condition is satisfied T (vL − vR ) fS (uL , uR ) = (ψL − ψR ), ψL = ψ(v(uL )), ψR = ψ(v(uR )). For illustrative purposes, we will prove a semi-discrete entropy inequality on a onedimensional mesh consisting of two elements of degree N . We assume both meshes are translations of a reference element [−1, 1], such that derivatives with respect to physical coordinates are identical to derivatives with respect to reference coordinates. The extension to multiple elements and variable mesh sizes is straightforward. The construction of entropy conservative schemes relies on appropriate SATs for Gauss collocation schemes [18, 13, 19]. Let the rows of Vf be denoted by column vectors tL , tR   (tL )1 , . . . , (tL )N +1 Vf = , (tL )j = ℓj (−1), (tR )j = ℓj (1). (tR )1 , . . . , (tR )N +1 The inter-element coupling terms in [18, 13, 19] utilize a decomposition of the surface matrix VfT BVf as (4) VfT BVf = tR tTR − tL tTL . The construction of entropy conservative schemes on multiple elements requires appropriate inter-element coupling terms (SATs) involving tL , tR . We consider a two element mesh, and show how when coupled with SATs, the resulting discretization matrices can be interpreted as constructing a global SBP operator. Let u1N , u2N denote nodal degrees of freedom of the vector valued solution u(x) on the first and second element, respectively. To simplify notation, we assume that 5 all following operators are defined in terms of Kronecker products [9], such that they are applied to each component of u1N , u2N . We first define the matrix 1 def S = Q − VfT BVf . 2 (5) It is straightforward to show (using Lemma 1) that S is skew-symmetric. We can now define an SBP operator Dh = Wh−1 Qh over two elements    1    1 T T S def def W 2 tR tL + − 2 tL tL (6) Qh = = , W . h 1 T W S − 12 tL tTR 2 t R tR | {z } | {z } Sh 1 2 Bh It can be shown that Dh is high order accurate such that, if uh is a polynomial of degree N , it is differentiated exactly. Straightforward computations show that Qh also satisfies an SBP property Qh + QTh = Bh . Ignoring boundary conditions, an entropy conservative scheme for (1) on two elements can then be given as  1 d u Wh uh + 2 (Qh ◦ FS ) 1 = 0, uh = N (7) u2N dt   (FS )ij = fS (uh )i , (uh )j , 1 ≤ i, j ≤ 2(N + 1), where ◦ denotes the Hadamard product [24]. It should be emphasized that here, (uh )i , (uh )j denote vectors containing solution components at nodes i, j, and that (because fS is a vector-valued flux) the term (FS )ij should be interpreted as a diagonal matrix whose diagonal entries consist of the components of fS ((uh )i , (uh )j ). T Multiplying (7) by vhT = v (uh ) will yield a semi-discrete version of the conservation of entropy (mimicking (3) with the inequality replaced by an equality) (8) d Wh S(uh ) + vhT (Bh ◦ FS ) 1 − 1T Bh ψ (uh ) = 0. dt We refer to [13, 10] for the proof of (8). The drawback of the SATs introduced in this section lies in the nature of the offdiagonal matrices tR tL and −tL tR . For Gauss nodes, these blocks are dense, which implies that inter-element coupling terms produce a non-compact stencil. Evaluating (7) requires computing two-point fluxes fS between all nodes on two neighboring elements, which significantly increases both the computational work, as well as communication between neighboring elements. This leads to all-to-all coupling between degrees of freedom in 1D, and to coupling along one-dimensional lines of nodes in higher dimensions due to the tensor product structure. The main goal of this work is to circumvent this tighter coupling of degrees of freedom introduced by the SATs described in this section, which can be done through the use of “decoupled” SBP operators. 3.3. Decoupled SBP operators. Decoupled SBP operators were first introduced in [11] and used to construct entropy stable schemes on simplicial elements. These operators (and simplifications under a collocation assumption) are presented in a more general setting in [11, 12] and in Appendix A. In this section, decoupled SBP operators are introduced in one dimension for GLL and Gauss nodal sets. 6 Decoupled SBP operators build upon the GSBP matrices W , Q, interpolation matrix Vf , and boundary matrix B introduced in Section 3.1. The decoupled SBP operator QN is defined as the block matrix   1 1 T T def Q − 2 Vf BVf 2 Vf B . (9) QN = 1 − 12 BVf 2B Lemma 1 and straightforward computations show that QN also satisfies the following SBP property Lemma 3. Let QN be defined through (9). Then,   0 T QN + QN = . B We note that the matrix QN acts not only on volume nodes, but on both volume and surface nodes. Thus, it is not immediately clear how to apply this operator to GSBP discretizations of nonlinear conservation laws. It is straightforward to evaluate the nonlinear flux at volume nodes since the solution is collocated at these points; however, evaluating the nonlinear flux at surface nodes is less straightforward. Moreover, QN does not directly define a difference operator, and must be combined with another operation to produce an approximation to the derivative. We will discuss how to apply QN in two steps. First, we will show how to approximate the derivative of an arbitrary function using QN given function values at both volume and surface nodes. Then, we will describe how to apply this approximation to compute derivatives of nonlinear flux functions given collocated solution values at volume nodes. Let f (x), g(x) denote two functions, and let f , g denote the values of f, g at interior nodal points. We also define vectors fN , gN denoting the values of f, g at both interior and boundary points     g (x1 ) f (x1 )         .. ..     . . f g     , gN = g (xN +1 ) = (10) fN = f (xN +1 ) = . f g     f f  g(−1)   f (−1)  g(1) f (1) ∂g can be computed using QN . Let u denote Then, a polynomial approximation to f ∂x ∂g the nodal values of the polynomial u(x) ≈ f ∂x . These coefficients are computed via  T I (11) Wu = diag (fN ) QN gN . Vf The approximation (11) can be rewritten in “strong” form as follows  T I u = W −1 diag (fN ) QN gN Vf 1 = diag (f ) Dg + diag (f ) W −1 VfT B (gf − Vf g) 2 1 −1 T + W Vf Bdiag (ff ) (gf − Vf g) , 2 where we have used the fact that diagonal matrices commute to simplify expressions. The decoupled SBP operator QN can thus be interpreted as adding boundary corrections to the GSBP operator D in a skew-symmetric fashion. More specifically, the 7 GSBP GSBP 2 100 Decoupled Decoupled SBP L2 errors Exact 0 10−1 10−2 10−3 10−4 −2 −1 0 −0.5 0.5 0 1 5 x coordinate 10 Degree N 15 (b) L2 errors (a) N = 5 approximation 2 Fig. 2: Approximations of derivatives of a Gaussian e−4x using the generalized SBP operator D and the decoupled SBP operator QN via (11). In Figure 2a, the colored circles and squares denote values at Gauss points for a degree N = 5 approximation. Figure 2b shows the convergence of L2 errors as N increases. expression (11) corresponds to a quadrature approximation of the following variational approximation of the derivative [11]: find u ∈ P N ([−1, 1]) such that, Z 1 u(x)v(x) = −1 Z 1 IN f −1 (f v + IN (f v)) ∂IN g v + (g − IN g) ∂x 2 1 , −1 ∀v ∈ P N ([−1, 1]), where IN denotes the degree N polynomial approximation at the (N + 1) Gauss points. The approximation (11) can also be applied to Gauss collocation schemes for nonlinear conservation laws. Let u ∈ P N be represented by the vector u of values at Gauss points, and let f (x), g(x) denote two nonlinear functions. The operator QN can be used to approximate the quantity f (u) ∂g(u) ∂x using (11) if f , g are defined as     f (u1 ) g (u1 )       .. ..       . . f (u) g(u)     , gN = g (uN +1 ) = fN = f (uN +1 ) = , uf = Vf u. f (uf ) g(uf )  

  f tT u   g tT u  L
L
f tTR u g tTR u It was shown in [11] that u is a high order accurate approximation to the quantity ∂g . Both the generalized SBP operator D and the expression in (11) involving the f ∂x decoupled SBP operator recover exact derivatives of high order polynomials. However, when applied to non-polynomial functions, the decoupled SBP operator QN improves accuracy near the boundaries. Figure 2 illustrates this by using both oper2 ators to approximate the derivative of a Gaussian e−4x on [−1, 1]. The decoupled SBP operator results in an improved approximation at all orders of approximation. 3.4. An entropy stable Gauss collocation scheme based on decoupled SBP operators. We can now construct an entropy conservative Gauss collocation 8 scheme with compact SATs using decoupled SBP operators. As in Section 3.2, we will construct a Gauss collocation scheme and provide a proof of semi-discrete entropy conservation for a two-element mesh. We first note that B can be trivially decomposed into the sum of two outer products as in (4)       −1 0 1 0 T T = eR eR − eL eL , eL = , eR = . (12) B= 0 1 0 1 The vectors eL , eR are related to tL , tR through the interpolation matrix Vf . Because tL , tR are rows of Vf , tL 1 = tR 1 = 1. This can be used to show, for example, that
(13) BVf 1 = eR eTR − eL eTL 1 = eR − eL . We can define a decoupled SBP matrix Qh over two  1 T S 2 Vf B  − 1 BVf − 1 eL eT def  L 2 Qh =  2 (14)  S − 21 eL eTR − 12 BVf elements as follows  1 T 2 eR eL 1 T 2 Vf B 1 T 2 eR eR   ,  where we have abused notation and redefined Qh . We can also show that Qh 1 = 0. Using (13), we have that        Q − 12 VfT BVf + 21 VfT BVf 1 S + 21 VfT B 1     1 1  − BVf 1 + 1 (eR − eL )  2  2 2 B (−1 + 1)  =     = 0. (15) Qh 1 =      S + 21 VfT B 1   Q − 21 VfT BVf + 21 VfT BVf 1  1 − 12 BVf 1 + 12 (eR − eL ) 2 B (−1 + 1) Here, we have used the definition of B in (12), the fact that Vf 1 = 1 and [20], and the definition of S in (5). This property (15) will be used in the entropy conservation. It can be helpful to split up Qh into two matrices  1 T S 2 Vf B 1 1 T  − BVf 1 def  2 eR eL Sh =  2 Qh = S h + Bh , 1 T 2  S 2 Vf B 1 1 T − 2 eL eR − 2 BVf    def Bh =   −eL eTL eR eTR  .  Q1 = 0 proof of      The matrix Sh is skew-symmetric, while the matrix Bh functions as a boundary operator which extracts boundary values over the two-element domain (i.e. Bh extracts a left boundary value from element 1 and a right boundary value from element 2). Note that here, Bh is diagonal, unlike the boundary operator defined in Section 3.2. It should also be noted that the two-element decoupled SBP operator Qh in (14) is related to the GSBP operator in (6) through a block interpolation matrix   I   1   1 T T  Vf S def T 2 tR tL + − 2 tL tL , I Q I = Ih =  h h 1 T . h  I S − 21 tL tTR 2 tR tR Vf 9 We will now utilize the two-element operators described in this section to construct an entropy stable Gauss collocation scheme on two elements. This scheme will differ from that of Section 3.2 in that neighboring elements will only be coupled together through face values. Let u1N , u2N denote the values of the conservative variables at Gauss points on elements 1 and 2, respectively, and let uh denote their concatenation as defined in (7). Let v(uiN ) denote the evaluation of entropy variables at Gauss points on element i, and define the “entropy-projected conservative variables” e 1f , u e 2f by evaluating the conservative variables in terms of the interpolated values of u the entropy variables at element boundaries e if = u(vfi ), u def vfi = Vf v(uiN ), def i = 1, 2. T  We now introduce uh = u1N u2N as the concatenated vector of solution values at e, and u e as follows: Gauss points, and define vh , v  1 uN  1  u e 1f  uN def def def  e = Ih v h , e = u (e vh = v (uh ) = v , v u v) =  2 u2N  . uN e 2f u e corresponds The term vh corresponds to the vector of entropy variables at uh , while v to the concatenated vector of Gauss point values and interpolated boundary values e denotes the evaluation of the conservative vfi of the entropy variables. The term u e. variables in terms of v Finally, we define FS as the matrix of evaluations of the two-point flux fS at e combinations of values of u def ej ) , (FS )ij = fS (e ui , u 1 ≤ i, j ≤ 2(N + 3). Note that, due to the consistency of fS , the diagonal of FS reduces to flux evaluations e i ) = f (e (FS )ii = fS (e ui , u ui ) . (16) We can now construct a semi-discretely entropy conservative formulation based on decoupled SBP operators: Theorem 4. Let Qh be defined by (14), and let uh denote the two-element solution of the following formulation: (17) Wh d uh + 2IhT (Qh ◦ FS ) 1 = 0. dt Then, uh satisfies a semi-discrete conservation of entropy
dS(uh ) + 1T B vfT f (e uf ) − ψ (e uf ) = 0 dt
 tT v u1N
def def e f = u(vf ). , u vf = TL tR v u2N 1T Wh Proof. The proof results from testing with vhT . Since v(u) = diagonal, the time term yields vhT ∂S(u) ∂u d d d Wh uh = 1T Wh diag (vh ) uh = 1T Wh S(uh ). dt dt dt 10 and Wh is The spatial term can be manipulated as follows   1 T 2vhT IhT (Qh ◦ FS ) 1 = 2 (Ih vh ) (Sh ◦ FS ) 1 + (Bh ◦ FS ) 1 2 eT (Bh ◦ FS ) 1 + v eT (Sh ◦ FS ) 1 − 1T (Sh ◦ FS ) v e, =v where we have used the skew-symmetry of Sh in the last step. The boundary term reduces to eT (Bh ◦ FS ) 1 = 1T Bh v eT f (e efT f (u ff ) , v u) = 1T B v eT f (e efT f (e u), v uf ) where we have used (16) and the fact that Bh is diagonal. Here, v e and denote vectors whose entries are the vector inner products of components of v f (e u) at volume and face points, respectively. The volume terms can be manipulated using the definition of FS and the Tadmor shuffle condition in Definition 2 X eT (Sh ◦ FS ) 1 − 1T (Sh ◦ FS ) v e= ej )T fS (e ej ) v (Sh )ij (e vi − v ui , u ij = X ij (Sh )ij (ψ (e ui ) − ψ (e uj ))   T T = ψ (e u) Sh 1 − 1T Sh ψ (e u) = 2 ψ (e u) S h 1 , where we have again used the skew-symmetry of Sh in the last step. Substituting Sh = Qh − 21 Bh and using (15) yields   T T 2 ψ (e u) Sh 1 = ψ (e u) (2Qh − Bh ) 1 T = −ψ (e u) Bh 1 = −1T Bh ψ (e u) = −1T Bψ (e uf ) , where we have used the symmetry of Bh in the second to last step. Remark. The proof of Theorem 4 follows directly from choosing either GLL or Gauss quadratures in Theorem 4 of [11]. The proof is reproduced here for clarity, as the two-element case illuminates the skew-symmetric nature and structure of the inter-element coupling more explicitly. Extending the two-element case to multiple elements can be done by defining analogous “global” differentiation matrices. Suppose that there are now three elements. Then, the matrices Sh , Bh can be defined as   1 T S 2 Vf B 1 T   − 1 BVf 2 eR eL   2   1 T def   V B S 2 f Sh =  1 T   − 12 eL eTR − 12 BVf e e R L  2   1 T   S 2 Vf B 1 1 T − 2 eL eR − 2 BVf     I  Vf     −eL eTL     I def . ,  0 I = Bh =  h     Vf      I eR eTR Vf 11 These global operators again satisfy Qh 1 = 0 and the SBP property Qh + QTh = Bh , where Qh = Sh + 21 Bh . Moreover, Qh can be used to produce high order accurate approximations of derivatives, and can be used to construct entropy conservative schemes as in Theorem 4. The extension to a general number of elements is done in a similar fashion. It is possible to convert the semi-discrete entropy equality in Theorem 4 to a semi-discrete entropy inequality by adding appropriate interface dissipation terms, such as Lax-Friedrichs or matrix dissipation [25]. We note that these terms must be e f in order to ensure a discrete dissipation of entropy [9, 11]. computed in terms of u Boundary conditions can also be incorporated into the formulation (17) in a weak fashion. Let uL , uR denote the values of the solution at the left and right domain boundaries of a 1D domain. Using the SBP property, we can modify (17) to   0 f (uL )   duh   T T Wh + 2Ih (Sh ◦ FS ) 1 + Ih Bh  ...  = 0,   dt  0  f (uR ) where we have used the fact that Bh is diagonal and the consistency of fS (implying that the diagonal of FS reduces to the evaluation of the flux f (u)) to simplify the boundary term (Bh ◦ FS ) 1. Boundary conditions are incorporated by replacing f (uL ), f (uR ) with numerical fluxes fL∗ , fR∗ as follows:   0 fL∗    duh   (18) Wh + 2IhT (Sh ◦ FS ) 1 + IhT Bh  ...  = 0.   dt 0 fR∗ If the boundary numerical flux satisfies the entropy stability conditions T ∗ ψL − v L fL ≤ 0, T ∗ ψR − vR fR ≤ 0 then the resulting scheme satisfies a global entropy inequality [9]. Assuming that the ODE system (17) exactly integrated in time and that entropy dissipative numerical fluxes and boundary conditions are used, the solution will satisfy a discrete entropy inequality (3). In practice, the system (17) is solved using an ODE time-stepper. For an explicit time-stepper, then all that is necessary is to invert the diagonal matrix Wh and evaluate the spatial terms in (18). We note that the conservation or dissipation of entropy is guaranteed up to timestepper accuracy. In practice, we observe that entropy is dissipated for all problems considered, despite the fact that the proof does not hold at the fully discrete level. A fully discrete entropy inequality can be guaranteed using implicit or space-time discretizations [26, 27]. 4. Extension to higher dimensions and non-affine meshes. The formulation in Theorem 4 can be naturally extended to Cartesian meshes in higher dimensions through a tensor product construction. We first consider the construction of higher b asdimensional differentiation matrices on a two-dimensional reference element Ω, suming a two dimensional tensor product grid of quadrature nodes (the construction 12 of decoupled SBP operators in three dimensions is straightforward and similar to the two-dimensional case). We then construct physical differentiation matrices on mapped elements Ωk , through which we construct an entropy conservative scheme. Let D1D , W1D denote the 1D differentiation and mass matrices, respectively, on the reference interval [−1, 1]. Let W denote the 2D reference mass matrix, and let D i denote the differentiation matrices with respect to the ith reference coordinate. These matrices can be expressed in terms of Kronecker products D 1 = D1D ⊗ IN +1 , D 2 = IN +1 ⊗ D1D , def def def W = W1D ⊗ W1D , where IN +1 denotes the (N + 1) × (N + 1) identity matrix. We also construct higher dimensional face interpolation matrices. Let Vf,1D denote the one-dimensional interpolation matrix, and let B1D denote the boundary matrix defined in Lemma 1. For an appropriate ordering of face quadrature points, the two-dimensional face interpolation matrix Vf and reference boundary matrices B 1 , B 2 can be expressed as the concatenation of Kronecker product matrices       def Vf,1D ⊗ I2 2 def 0 1 def B1D ⊗ I2 . , B = , B = Vf = 0 I2 ⊗ Vf,1D I2 ⊗ B1D In three dimensions, the face interpolation matrix would be expressed as the concatenation of three Kronecker products involving Vf,1D . The higher dimensional differentiation and interpolation matrices D i , Vf can now be used to construct higher dimensional decoupled SBP operators. Let QiN denote the decoupled SBP operator for the ith coordinate on the reference element, where QiN is defined as   i 1 T i 1 T i def def Q − 2 Vf B Vf 2 Vf B , Qi = W D i . QiN = 1 i B − 12 B i Vf 2 Let the domain now be decomposed into non-overlapping elements Ωk , such that b under a degree N polynomial mapping Φk . We define geometric Ω is the image of Ω xj k ∂b k b with respect to terms Gij = J ∂xi as scaled derivatives of reference coordinates x physical coordinates x. We also introduce the scaled normals ni Jfk , which can be computed on quadrilateral and tensor product elements via k ni Jfk = ± d X j=1 Gkij = ± d X j=1 Jk ∂b xj , ∂xi where the sign of ni Jfk is negative for a “left” face and positive for a “right” face. These geometric terms introduce scalings J k , Jfk , where J k is the determinant of the Jacobian of Φk and Jfk denotes the determinant of the Jacobian of the mapping from a physical face to a reference face. The quantities Gkij , ni Jfk can now be used to define matrices over each physical element Ωk

def def Bki = B i diag nki , Wk = W diag J k , d (19) Qik = def


1X diag Gkij QiN + QiN diag Gkij , 2 j=1 where we have discretized the curved differentiation matrix in split form [28, 29]. Here, J k denotes the vector of values of J k at volume quadrature points, Gkij denotes 13 the vector of values of Gkij at volume and face quadrature points, and nki denotes the ith scaled outward normal ni Jfk at face quadrature points. A 2D Gauss collocation scheme can now be given in terms of Wk , Qik , Bki . Let k uN denote the vector of solution values on Ωk , and define the quantities FSi
en ) , = fSi (e um , u 1 ≤ m, n ≤ (N + 1)2 + 4(N + 1)  k

def u def def e= N e kf = u vfk , u vfk = Vf v ukN , , u k ef u def mn where (N + 1)2 + 4(N + 1) is the total number of volume and face points. We then have the following theorem on the semi-discrete conservation of entropy: Theorem 5. Assume that Qk 1 = 0, and that Ωk are mapped curvilinear elements. Let ukN satisfy the local formulation on Ωk    
dukN X  e+ e kf − f i (e + uf ) = 0, 2 Qik ◦ FSi 1 + VfT Bki fSi u f ,u dt i=1 d (20) Wk e + denotes boundary values of the entropy-projected conservative variables on where u neighboring elements. Then, ukN satisfies the discrete conservation of entropy !
d  X X

 dS ukN T k T i k k T i ef 1 Wk + = 0. uf ) − ψ u 1 Bk vf f (e dt i=1 k The proof is a special case of the proof of Theorem 1 in [12], with the quadrature rule taken to be a tensor product rule with (N + 1) Gauss (or GLL) points in each coordinate direction. Remark. The operator Qik in (19) can be applied without needing to explicitly store geometric terms Gkij at face quadrature points. By using the structure of the boundary matrix B i , expressions involving geometric terms on faces can be replaced by expressions involving components of scaled outward normals ni Jfk on Ωk . 4.1. Discrete geometric conservation law. We note that Theorem 5 relies on the assumption that Qk 1 = 0. This condition is equivalent to ensuring that the scaled geometric terms Gkij satisfy a discrete geometric conservation law (GCL) [2, 7, 10, 12]. For two-dimensional degree N (isoparametric) mappings, the GCL is automatically satisfied. However, in three dimensions, the GCL is not guaranteed to be satisfied at the discrete level, due to the fact that geometric terms for isoparametric mappings are polynomials of degree higher than N . It is possible to ensure the satisfaction of a discrete GCL by using a sub-parametric polynomial geometric mapping. Let Ngeo denote the degree of a polynomial geometric mapping. In three dimensions, the exact geometric terms for a degree Ngeo polynomial   mapping are polynomials of degree 2Ngeo [30, 16, 10]. If 2Ngeo ≤ N , or if Ngeo ≤ N2 , 1 then the discrete  NGCL  is automatically satisfied. For Ngeo ≥ 2 , modifications to the computation of geometric terms are required to ensure that the GCL is satisfied at the discrete level. For general SBP operators, the discrete GCL can be enforced through the solution of a local least squares problem [10]. 1 We note that this condition is specific to three-dimensional hexahedral elements, and does not necessarily hold for other element types. For example, the discrete GCL is satisfied if 2Ngeo − 2 ≤ N for isoparametric tetrahedral elements [12]. 14 We take an alternative approach and construct geometric terms using the approach of Kopriva [30]. This construction takes advantage of the fact that GLL and Gauss collocation methods correspond to polynomial discretizations. Kopriva’s construction is based on rewriting the geometric terms as the reference curl of an interpolated auxiliary quantity     b × IN z ∇y b −∇  k  G1j   j    b × IN z ∇x b , Gk2j  =  ∇ (21) j = 1, 2, 3.   j    Gk3j b × IN x∇y b ∇ j Here, IN denotes the polynomial interpolation operator using GLL nodes. By interpolating the auxiliary quantity in (21) using polynomial interpolation prior to applying the curl, the geometric terms are approximated by degree N polynomials which satisfy the discrete GCL by construction. These GCL-satisfying geometric terms can then be used to compute normal vectors. For a watertight mesh, the constructed normal vectors are guaranteed to be continuous across faces [12]. To summarize, extending higher dimensional entropy stable Gauss collocation schemes to curved meshes requires the following steps: 1. Construct polynomial approximations of the geometric terms using equation (21) and interpolation at GLL nodes [30]. 2. Evaluate approximate geometric terms at volume and surface points, and compute normal vectors in terms of the approximate geometric terms. 3. Compute physical derivatives using the split form (19). Apart from evaluating the GCL-satisfying geometric terms at separate volume and surface points, entropy stable Gauss collocation schemes are extended to curved meshes in the same manner as GLL and staggered-grid collocation schemes [2, 14, 12]. 5. Numerical results. In this section, we present numerical examples for the compressible Euler equations, which are given in d dimensions as follows: d ∂ρ X ∂ (ρuj ) = 0, + ∂t j=1 ∂xj d ∂ρui X ∂ (ρui uj + pδij ) = 0, + ∂t ∂xj j=1 i = 1, . . . , d d ∂E X ∂ (uj (E + p)) = 0. + ∂t ∂xj j=1 Here, ρ is density, ui denotes the ith component of velocity, and E is the total energy. The pressure p and specific internal energy ρe are given by      d d X X 1 1 p = (γ − 1) E − ρ  ρe = E − ρ  u2j  , u2j  . 2 2 j=1 j=1 There exists an infinite family of suitable convex entropies for the compressible Euler equations [31]. However, there is only a single unique entropy which appropriately treats the viscous heat conduction term in the compressible Navier-Stokes 15 equations [32]. This entropy S(u) is given by S(u) = − ρs , γ−1   where s = log ρpγ is the physical specific entropy, and the dimension d = 1, 2, 3. The entropy variables in d dimensions are given by ρe(γ + 1 − s) − E , ρe ρui = , i = 1, . . . , d, ρe ρ =− , ρe v1 = v1+i vd+2 while the conservation variables in terms of the entropy variables are given by ρ = −(ρe)vd+2 , ρui = (ρe)v1+i , i = 1, . . . , d ! Pd 2 j=1 v1+j , E = (ρe) 1 − 2vd+2 where the quantities ρe and s in terms of the entropy variables are ρe =  (γ − 1) γ (−vd+2 ) 1/(γ−1) e −s γ−1 , s = γ − v1 + Pd j=1 2 v1+j 2vd+2 . Let f denote some piecewise polynomial function, and let f + denote the exterior value of f across an element face. We define the average and logarithmic averages as follows: f+ + f f+ − f log {{f }} = , {{f }} = . 2 log (f + ) − log (f ) The average and logarithmic average are assumed to act component-wise on vectorvalued functions. We evaluate the logarithmic average using the numerically stable algorithm of [33]. Explicit expressions for entropy conservative numerical fluxes in two dimensions are given by Chandrashekar [34] 1 f1,S (uL , uR ) = {{ρ}} 2 f1,S (uL , uR ) 3 f1,S (uL , uR ) 4 f1,S (uL , uR ) = = log {{u1 }} , 1 f1,S {{u1 }} 2 f2,S , + pavg , = (Eavg + pavg ) {{u1 }} , 1 f2,S (uL , uR ) = {{ρ}} 2 f2,S (uL , uR ) 3 f2,S (uL , uR ) 4 f2,S (uL , uR ) = = 1 f2,S 1 f2,S log {{u2 }} , {{u1 }} , {{u2 }} + pavg , = (Eavg + pavg ) {{u2 }} , where we have introduced the auxiliary quantities (22) β= ρ , 2p pavg = u2avg {{ρ}} , 2 {{β}} Eavg = {{ρ}} log log 2 {{β}} (γ − 1)  2 
= 2({{u1 }} + {{u1 }} ) − {u1 } + {u22 } . 2 2 16 + u2avg , 2 Expressions for entropy conservative numerical fluxes for the three-dimensional compressible Euler equations can also be explicitly written as f1,S     =   log {{ρ}} {{u1 }} log 2 {{ρ}} {{u1 }} + pavg log {{ρ}} {{u1 }} {{u2 }} log {{ρ}} {{u1 }} {{u3 }} (Eavg + pavg ) {{u1 }}  f3,S     ,   f2,S log  log {{ρ}} {{u2 }} log {{ρ}} {{u1 }} {{u2 }} log 2 {{ρ}} {{u2 }} + pavg log {{ρ}} {{u2 }} {{u3 }} (Eavg + pavg ) {{u2 }}     =   {{ρ}} {{u3 }}  log {{ρ}} {{u1 }} {{u3 }}   =  {{ρ}}log {{u2 }} {{u3 }}   {{ρ}}log {{u3 }}2 + pavg (Eavg + pavg ) {{u3 }}     ,      ,   where the auxiliary quantities are defined as β= ρ , 2p u2avg log 1 log {{ρ}} u2avg , 2 2(γ − 1) {{β}}
   2 2 2 = 2({{u1 }} + {{u2 }} + {{u3 }} ) − {u21 } + {u22 } + {u23 } . pavg = {{ρ}} , 2 {{β}} Eavg = {{ρ}} log + In all problems, we use an explicit low storage RK-45 time-stepper [35] and estimate the timestep size dt using J, Jfk , and degree-dependent L2 trace constants CN dt = CCFL h , aCN h= 1 kJ −1 kL∞ Jfk , L∞ where a is an estimate of the maximum wave speed, h estimates the mesh size, and CCFL is some user-defined CFL constant. For isotropic elements, the ratio of J k to Jfk scales as the mesh size h, while CN captures the dependence of the maximum stable timestep on the polynomial degree N . For hexahedral elements, CN varies depending on the choice of quadrature. It was shown in [17] that CN = ( d N (N2+1) +2) d (N +1)(N 2 for GLL nodes . for Gauss nodes Thus, based on this rough estimate of the maximum stable timestep, GLL collocation schemes should be able to take a timestep which is roughly (1+2/N ) times larger than the maximum stable timestep for Gauss collocation schemes. We do not account for this discrepancy in this work, and instead set the timestep for both GLL and Gauss collocation schemes based on the more conservative Gauss collocation estimate of dt. Numerical results in 1D are similar to those presented in [11]. Thus, we focus on two and three dimensional problems and comparisons of entropy stable GLL and Gauss collocation schemes. 5.1. 2D isentropic vortex problem. We begin by examining high order convergence of the proposed methods in two dimensions using the isentropic vortex prob17 (a) Lightly warped mesh (b) Moderately warped mesh (c) Heavily warped mesh Fig. 3: Lightly, moderately, and heavily warped meshes for N = 4, K = 16. lem [36, 13]. The analytical solution is (23) ρ(x, t) = 1− u1 (x, t) = 1 − 1 2 (γ 2 − 1)(βe1−r(x,t) )2 8γπ 2 β 1−r(x,t)2 e (y − y0 ), 2π 1 ! γ−1 , p = ργ , u2 (x, t) = β 1−r(x,t)2 e (x − x0 − t), 2π p where u1 , u2 are the x and y velocity and r(x, t) = (x − x0 − t)2 + (y − y0 )2 . Here, we take x0 = 5, y0 = 0 and β = 5. We solve on a periodic rectangular domain [0, 20] × [−5, 5] until final time T = 5, and compute errors over all all solution fields. For a degree N approximation, we approximate the L2 error using an (N + 2) point Gauss quadrature rule. We also examine the influence of element curvature for both GLL and Gauss collocation schemes by examining L2 errors on a sequence of moderately and heavily warped curvilinear meshes (see Figure 3). These warpings are constructed by modifying nodal positions according to the following mapping      3π Lx π cos y , x− x e = x + Lx α cos Lx 2 Ly      4π π Lx ye = y + Ly α sin cos y , x e− Lx 2 Ly where Lx = 20, Ly = 10 denote the lengths of the domain in the x and y directions, respectively. The lightly warped mesh corresponds to α = 1/64, the moderately warped mesh corresponds to α = 1/16, and the heavily warped mesh corresponds to α = 1/8. All results are computed using CCFL = 1/2 and Lax-Friedrichs interface dissipation. 18 L2 errors 10−1 10−1 10−4 10−4 10−7 10−7 GLL 10 GLL Gauss 10 −10 10−1.5 Mesh size h Gauss −10 10−1 10−1.5 Mesh size h L2 errors (a) Affine mesh (b) Lightly warped mesh 10−1 10−1 10−4 10−4 10−7 10−7 GLL 10 GLL Gauss 10 −10 10−1.5 Mesh size h 10−1 10−1 Gauss −10 10−1.5 Mesh size h (c) Moderately warped mesh 10−1 (d) Heavily warped mesh Fig. 4: L2 errors for the 2D isentropic vortex at time T = 5 for degree N = 2, . . . , 7 GLL and Gauss collocation schemes. Figure 4 shows the L2 errors for affine, moderately warped, and heavily warped meshes. For affine meshes, Gauss collocation results in a lower errors than GLL collocation at all orders. However, the difference between both schemes decreases as N increases. This is not too surprising: on a Cartesian domain, the discrete L2 inner product resulting from GLL quadrature converges to exact L2 inner product over the space of polynomials as N increases [37]. However, GLL and Gauss collocation differ more significantly on curved meshes. For both moderately and heavily warped meshes, the errors for a degree N Gauss collocation scheme are nearly identical to errors for a higher order GLL collocation scheme of degree (N + 1). These results are in line with numerical experiments in [14], which show that GLL collocation schemes lose one order of convergence in the L2 norm on unstructured non-uniform meshes. Both results show that increasing quadrature accuracy significantly reduces the effect of polynomial aliasing due to curved meshes and spatially varying geometric terms. We note that L2 approximation estimates on curved meshes [38, 39] assume that the mesh size is small enough to be in the asymptotic regime (such that asymptotic error estimates hold). On curved meshes, this requires that the mesh is sufficiently 19 fine to resolve both the solution and geometric mapping. The results in Figure 4 suggest that, compared to Gauss collocation schemes, under-integrated GLL collocation schemes require a finer mesh resolution to reach the asymptotic regime. 5.2. 3D isentropic vortex problem. As in two dimensions, we test the accuracy of the proposed scheme using an isentropic vortex solution adapted to three dimensions. The solution is the extruded 2D vortex propagating in the y direction, with an analytic expression given in [40] ρ(x, t) =  (γ − 1) 2 1− Π 2 ui (x, t) = Πri , p0 E(x, t) = γ−1  1  γ−1 γ−1 2 1− Π 2 γ  γ−1 d + ρX 2 u . 2 i=1 i where (u1 , u2 , u3 )T is the three-dimensional velocity vector and     P −(x2 − c2 − t) r1 2 r 1− d i=1 i . r2  =  2 x1 − c 1 Π = Πmax e , 0 r3 We take c1 = c2 = 7.5, p0 = 1/γ, and Πmax = 0.4, and solve on the domain [0, 15] × [0, 20] × [0, 5] until final time T = 5. We decompose the domain into uniform hexahedral elements with edge length h. As in the 2D case, we also examine the effect of curvilinear mesh warping. We construct a curved warping of the initial Cartesian mesh by mapping nodes on each hexahedron to warped nodal positions (e x, ye, ze) through the transformation       1 (y − 10) (z − 2.5) (x − 7.5) cos π cos π ye = y + Ly cos 3π 8 15 20 5       (e y − 10) (z − 2.5) (x − 7.5) 1 sin 4π cos π x e = x + Lx cos π 8 15 20 5       (e y − 10) (z − 2.5) (e x − 7.5) 1 cos 2π cos π . ze = z + Lz cos π 8 15 20 5 where Lx = 15, Ly = 20, Lz = 5. On curved meshes, the geometric terms are constructed using the approach of Kopriva described in Section 4. A CCFL = .75 is used for all experiments. Figure 5 shows the L2 errors for degrees N = 2, . . . , 5.2 As in the 2D case, Gauss collocation schemes produce smaller errors than GLL collocation. The difference between the two schemes on affine meshes is slightly more pronounced than in 2D, while the difference between the two schemes on curved meshes is less significant than observed in 2D experiments. We expect that this difference is due to the specific curved mapping. The warped 2D mesh used in Figure 4 was generated to mimic a severe “vorticular” warping. The 3D mapping is less severely warped, due to different domain dimensions and difficulties ensuring invertibility of the map from the reference to physical element. 2 For the N = 2 GLL collocation scheme on the coarsest 3D mesh, the error is not shown because the solution diverged. 20 100 10−2 10−2 L2 errors L2 errors 100 10−4 GLL 10−6 GLL 10−6 Gauss 10−0.5 10−4 100 Mesh size h Gauss 10−0.5 (a) Affine mesh 100 Mesh size h (b) Curved mesh Fig. 5: L2 errors for the 3D isentropic vortex for N = 2, . . . , 5 on sequences of Cartesian and curved meshes. We also compared L2 errors for both isoparametric and sub-parametric geometric mappings. In both cases, the discrete GCL is satisfied. For sub-parametric mappings,   we chose the degree of approximation of geometry Ngeo = N2 + 1, such that the geometric terms are computed exactly and the GCL is satisfied by default. This test is intended to address the fact that the GCL-preserving interpolation of Kopriva introduces a small approximation error, as the geometric terms are no longer exact. For these sub-parametric mappings, the gap between GLL and Gauss collocation widens slightly at N = 2. However, the results for sub-parametric mappings are nearly identical to the results in the isoparametric case for higher polynomial degrees. 5.3. Shock-vortex interaction. The next problem considered is the shockvortex interaction described in [36]. The domain is taken to be [0, 2] × [0, 1], and is triangulated with uniform quadrilateral elements. Wall boundary conditions are imposed on the top and bottom boundaries, and inflow boundary conditions are typically imposed on the left and right boundaries. We modify the problem setup such that periodic boundary conditions are imposed at the left and right boundaries. Wall boundary conditions are imposed using a mirror state for the normal velocity, which was shown to be entropy stable in [41, 9]. The initial condition is taken to be the superposition of a stationary shock and a vortex propagating towards the right. The stationary Mach Ms = 1.1 shock is posi
√ tioned at x = .5 normal to the x axis, with left state (ρL , uL , vL , pL ) = 1, γ, 0, 1 , where uL , vL denote x, y components of the left-state velocity. The right state is a scaling of the left state computed using the Rankine-Hugoniot conditions, such that the ratio of upstream and downstream states is ρL uL 2 + (γ − 1)Ms2 = = , ρR uR (γ + 1)Ms2
2γ pL =1+ Ms2 − 1 , pR γ+1 vR = 0. The isentropic vortex is centered at (xc , yc ) = (.25, .5) and given in terms of velocity fluctuations δu and δv, which are functions of the tangential velocity vθ δu = vθ sin(θ), δv = −vθ cos(θ), 21 2 vθ = ǫτ eα(1−τ ) , (a) Entropy conservative flux, T = .3 (b) Entropy conservative flux, T = .7 (c) Lax-Friedrichs flux, T = .3 (d) Lax-Friedrichs flux, T = .7 (e) Matrix dissipation flux, T = .3 (f) Matrix dissipation flux, T = .7 Fig. 6: Shock vortex solution at time T = .7 using entropy stable Gauss collocation schemes with N = 4, h = 1/100. p where r =  (x −xc )2 + (y − yc )2 is the radius from the vortex center, τ = r/rc , and y−yc . We follow [36] and take ǫ = .3, α = .204, and rc = .05. The vortex θ = tan−1 x−x c temperature is computed as a fluctuation δT of the upstream state TL = pL /ρL 2 (γ − 1)ǫ2 e2α(1−τ ) δT = − . 4αγ The vortex density and pressure are computed using an isentropic assumption. To summarize, the initial condition for the shock-vortex interaction problem is ρ = ρs  Tvor TL 1  γ−1 , u1 = us + δu, u2 = vs + δv, p = ps  Tvor TL 1  γ−1 , where (ρs , us , vs , ps ) denote the discontinuous stationary shock solution given by the left and right states (ρL , uL , vL , pL ) , (ρR , uR , vR , pR ). We compare three different entropy stable Gauss collocation schemes. All three utilize the entropy conservative flux of Chandrashekar [34]. For the first scheme, 22 we do not introduce any additional interface dissipation, which produces an entropy conservative scheme. The second scheme introduces an entropy-dissipative interface term using Lax-Friedrichs penalization, while the third scheme utilizes the matrix dissipation flux introduced in [25]. This flux adds a dissipation of the form RDRT Jvfk K, where Jvfk K denotes the jump in the entropy variables. In two dimensions, the matrices R, D are   1 1 0 1 {{u1 }} − ānx {{u1 }} ny {{u1 }} + ānx   R= {{u2 }} − āny {{u2 }} −nx {{u2 }} + āny  1 2 h − āūn ūn ny − {{u2 }} nx h + āūn 2 uavg   log log {{ρ}} (γ−1) {{ρ}}log D = diag |ūn − a| {{ρ}} , , |ū | , |ū | p , |ū + a| n n avg n 2γ γ 2γ where ā, ūn are defined in 2D as ūn = {{u1 }} nx + {{u2 }} ny , ā = s γpavg {{ρ}} log , h= γ 2(γ − 1) {{β}} log 1 + ū, 2 and pavg , u2avg are defined as in (22). Figure 6 shows density solutions for N = 4 and h = 1/100 Gauss collocation schemes using a non-dissipative entropy conservative flux, a dissipative Lax-Friedrichs flux, and a matrix dissipation flux. In all cases, the vortex passes through the shock stably without the use of additional slope limiting, filtering, or artificial viscosity. However, the entropy conservative scheme produces a large number of spurious oscillations in the solution. These are reduced away from the shock under Lax-Friedrichs dissipation, though oscillations still persist around a large neighborhood of the discontinuity. The Gibbs-type oscillations are most localized under the matrix dissipation flux. We note that this experiment also verifies that entropy stable decoupled SBP schemes (including the over-integrated case [11]) are compatible with entropy stable wall boundary conditions. As far as the authors know, the stable and high order accurate imposition of such boundary conditions for existing GSBP couplings described in [13] and Section 3.2 remains an open problem. 5.4. Inviscid Taylor-Green vortex. We conclude by investigating the behavior of entropy stable Gauss collocation schemes for the inviscid Taylor–Green vortex [42, 3, 10]. This problem is posed on the periodic box [−π, π]3 , with initial conditions 100 1 + (cos(2x1 ) + cos(2x2 )) (2 + cos(2x3 )) , γ 16 u1 = sin(x1 ) cos(x2 ) cos(x3 ), u2 = − cos(x1 ) sin(x2 ) cos(x3 ), u3 = 0. ρ = 1, p= The Taylor–Green vortex is used to study the transition and decay of turbulence. In the absence of viscosity, the Taylor–Green vortex develops successively smaller scales as time increases. As a result, the solution is guaranteed to contain under-resolved features after a sufficiently large time. We study the evolution of kinetic energy κ(t) 1 κ(t) = |Ω| Z ρ Ω d X i=1 23 u2i ! dx, ·10−2 ·10−2 1.5 1.5 Affine Affine Curved Curved 1 0.5 0.5 0 0 − ∂κ ∂t 1 0 5 10 Time t 15 20 (a) GLL collocation 0 5 10 Time t 15 20 (b) Gauss collocation Fig. 7: Kinetic energy dissipation rate for entropy stable GLL and Gauss collocation schemes with N = 7 and h = π/8. as well as the kinetic energy dissipation rate − ∂κ ∂t , which is approximated by differencing κ(t). For both GLL and Gauss collocation schemes, integrals in the kinetic energy formulaare evaluated using an (N + 1)-point Gauss quadrature rule. Figure 7 shows the evolution of the kinetic energy dissipation rate from t ∈ [0, 20] for Gauss and GLL collocation schemes on affine and curved meshes. The curved meshes used here are constructed by modifying nodal positions through the mapping e =x+ x 1 sin(x) sin(y) sin(z). 2 All cases utilize N = 7 and h = π/8 (corresponding to 8 elements per side), as well as a CFL of .25. Lax-Friedrichs dissipation is used for all simulations. For both affine and curved meshes, the presented Gauss collocation schemes remain stable in the presence of highly under-resolved solution features. Kinetic energy dissipation rates for both GLL and Gauss collocation are qualitatively similar and are consistent with existing results in the literature for the inviscid Taylor-Green vortex [3, 12]. 6. A theoretical cost comparison. While the numerical experiments presented in previous sections demonstrate several advantages of Gauss collocation methods over GLL collocation schemes, these do not account for additional costs associated with Gauss collocation schemes. While a detailed time-to-solution comparison is outside of the scope of this work, we can compare computational costs associated with Gauss, staggered-grid, and GLL collocation schemes. The main computational costs associated with entropy stable schemes are volume operations, which include evaluations of two-point fluxes and applications of one-dimensional differentiation and interpolation matrices. We do not count interelement communication or flux computations, as these are typically less expensive than volume operations (especially for higher polynomial degrees). The total number of flux computations for each scheme can also be reduced using the symmetry of fS ; however, this does not affect relative differences in the number of two-point flux evaluations between GLL, Gauss, and staggered grid schemes. However, exploiting 24 ·104 ·104 2 GLL Staggered Gauss GLL Staggered Gauss 4 1 2 0 0 1 2 3 4 5 Degree N 6 7 1 (a) Two-point flux evaluations per dimension 2 3 4 5 Degree N 6 7 (b) Est. matrix computation operations Fig. 8: Comparison of GLL, staggered grid, and Gauss collocation schemes in terms of number of two-point flux evaluations and operations associated with matrix computations in three dimensions. symmetry in the Gauss scheme is slightly more straightforward (compared to GLL and staggered grid schemes) due to the block structure of the decoupled SBP operator. In 3D, a degree N GLL collocation scheme contains (N + 1)3 nodes. Two-point fluxes are computed between states at one node and states at (N + 1) additional nodes, resulting in a total of 3(N + 1)4 two-point flux evaluations. For a staggered grid scheme, two-point fluxes are computed on a degree (N + 1) GLL grid consisting of (N + 2)3 nodes, resulting in 3(N + 2)4 evaluations. Gauss collocation schemes compute two-point fluxes on a grid of Gauss nodes, resulting in 3(N + 1)4 flux evaluations in 3D. One must also evaluate two-point fluxes between face nodes and volume nodes, and vice versa. As a result, Gauss schemes evaluate two-point fluxes twice between each face node and the line of (N + 1) volume nodes normal to that face, resulting in an extra 12(N + 1)3 flux evaluations in 3D. We also consider costs associated with applying operator matrices. GLL and Gauss collocation schemes require applying a one-dimensional differentiation matrix
4 to each line of nodes, resulting in O 3(N + 1) operations in 3D. Staggered-grid
schemes require O 3(N + 2)4 operations (corresponding to differentiation in each
3 coordinate on a degree (N + 1) GLL grid) as well as O 6(N + 2)(N + 1) operations required for interpolation to and from a (N + 1) point Gauss grid to a (N + 2) point GLL grid in three dimensions. Gauss collocation schemes require only interpolation
to and from points on 6 faces, resulting in an additional O 12(N + 1)3 operations per dimension per element. Figure 8 shows the estimated number of two-point flux evaluations and matrix operations for GLL, staggered grid, and Gauss collocation schemes in three dimensions. We observe that a straightforward implementation of Gauss collocation does not significantly reduce the number of flux evaluations compared to staggered grid schemes, though Gauss collocation schemes result in a significantly smaller number of operations from matrix computations. 25 We note that the number of two-point flux evaluations and matrix operations will impact runtime differently depending on the implementation and computational architecture. For example, while flux evaluations typically dominate runtimes for serial CPU implementations at all orders of approximation, they do not contribute significantly to runtimes at polynomial degrees N = 1, . . . , 8 for implementations on Graphics Processing Units (GPUs) [43]. 7. Conclusion. This work shows how to construct efficient entropy stable high order Gauss collocation DG schemes on quadrilateral and hexahedral meshes. Key to the construction of efficient methods are decoupled SBP operators, which deliver entropy stability and high order accuracy while maintaining compact inter-element coupling terms. These operators are also compatible with existing entropy stable methods for applying interface dissipation [25] or imposing boundary conditions. Numerical experiments demonstrate both the stability and high order accuracy of the proposed Gauss collocation schemes on both affine and curvilinear meshes. We note that, while the numerical experiments presented here consider only mapped Cartesian domains, the method is also applicable to complex geometries, and future work will focus on studying the performance of such methods on curvilinear quadrilateral and hexahedral unstructured meshes. We note that results for Gauss collocation are similar to those attained by entropy stable staggered-grid schemes [14], and require a similar number of two-point flux evaluations. However, Gauss collocation schemes result in a lower number floating point operations from matrix computations compared to staggered-grid methods. Finally, while a rigorous computational comparison between GLL and Gauss collocation schemes remains to be done, Gauss collocation schemes show significant improvements in accuracy compared to GLL collocation schemes on non-Cartesian meshes. In particular, for sufficiently warped curvilinear mappings, degree N Gauss collocation schemes achieve an accuracy comparable to degree (N + 1) GLL collocation schemes in two and three dimensions. 8. Acknowledgments. Jesse Chan is supported by NSF DMS-1719818 and DMS-1712639. The authors thank Florian Hindenlang and Jeremy Kozdon for helpful discussions and suggestions. Appendix A. Decoupled SBP operators for general choices of quadrature and basis. For general choices of quadrature and basis, decoupled projection operators involve a volume quadrature interpolation matrix Vq , a face quadrature interpolation Np matrix Vf , and a quadrature-based L2 projection matrix Pq . Let {φj }j=1 denote N q a set of Np basis functions, and let {xi , wi }i=1 denote a set of Nq volume quadrature points and weights in d dimensions. We also introduce the set of Nqf surface n oNqf quadrature points and weights xfi , wif . Then, Vq , Vf are given as i=1 (Vq )ij = φj (xi ),   (Vf )ij = φj xfi , 1 ≤ i ≤ Nq , 1 ≤ i ≤ Nqf , 1 ≤ j ≤ Np , 1 ≤ j ≤ Np . These matrices can be used to define the quadrature-based L2 projection matrix Pq . Let W , Wf denote the diagonal matrix of volume and surface quadrature weights, respectively. Then, M = VqT W Vq , Pq = M −1 VqT W . 26 Let D i now denote a modal differentiation matrix with respect to the ith coordinate, which maps coefficients in the basis φj to coefficients of the ith derivative. By composing this matrix with interpolation and projection matrices, one can define differencing operators Dqi = Vq D i Pq which map values at quadrature points to values of approximate derivatives at quadrature points. Moreover, Qi = W Dqi satisfies a generalized SBP property involving the face interpolation and projection matrices Vf , Pq [11]. The decoupled SBP operator QiN is then given as (24) QiN =  Qi − 1 2 T (Vf Pq ) Wf diag (ni ) Vf Pq − 12 Wf diag (ni ) Vf Pq 1 2  T (Vf Pq ) Wf diag (ni ) . 1 2 Wf diag (ni ) A straightforward computation shows that QiN satisfies an SBP property [11]. It is worth noting that the form of QiN does not depend on the choice of basis. So long as the approximation space spanned by the basis φj does not change, the domain and range of QiN depend solely on the choice of volume and surface quadrature points. A collocation scheme assumes that the number of quadrature points is identical to the number of basis functions. If the solution is represented using degree N Lagrange polynomials at quadrature points, the matrices Vq , Pq simplify to (Vq )ij = δij , M = W, Pq = M −1 VqT W = I. Plugging these simplifications into (24) and restricting to one spatial dimension recovers the decoupled SBP operator (9). REFERENCES [1] 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, 2013. [2] Mark H Carpenter, Travis C Fisher, Eric J Nielsen, and Steven H Frankel. Entropy Stable Spectral Collocation Schemes for the Navier–Stokes Equations: Discontinuous Interfaces. SIAM Journal on Scientific Computing, 36(5):B835–B867, 2014. [3] Gregor J Gassner, Andrew R Winters, and David A Kopriva. Split form nodal discontinuous Galerkin schemes with summation-by-parts property for the compressible Euler equations. Journal of Computational Physics, 327:39–66, 2016. [4] Gregor J Gassner. A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods. SIAM Journal on Scientific Computing, 35(3):A1233–A1253, 2013. [5] Zhijian J Wang, Krzysztof Fidkowski, Rémi Abgrall, Francesco Bassi, Doru Caraeni, Andrew Cary, Herman Deconinck, Ralf Hartmann, Koen Hillewaert, Hung T Huynh, et al. Highorder CFD methods: current status and perspective. International Journal for Numerical Methods in Fluids, 72(8):811–845, 2013. [6] Constantine M Dafermos. Hyperbolic conservation laws in continuum physics. Springer, 2005. [7] Gregor J Gassner, Andrew R Winters, Florian J Hindenlang, and David A Kopriva. The BR1 scheme is stable for the compressible Navier–Stokes equations. Journal of Scientific Computing, pages 1–47, 2017. [8] Eitan Tadmor. The numerical viscosity of entropy stable schemes for systems of conservation laws. I. Mathematics of Computation, 49(179):91–103, 1987. [9] Tianheng Chen and Chi-Wang Shu. Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws. Journal of Computational Physics, 345:427–461, 2017. [10] Jared Crean, Jason E Hicken, David C Del Rey Fernández, David W Zingg, and Mark H Carpenter. Entropy-stable summation-by-parts discretization of the Euler equations on general curved elements. Journal of Computational Physics, 356:410–438, 2018. [11] Jesse Chan. On discretely entropy conservative and entropy stable discontinuous Galerkin methods. Journal of Computational Physics, 362:346 – 374, 2018. 27 [12] Jesse Chan and Lucas C Wilcox. Discretely entropy stable weight-adjusted discontinuous Galerkin methods on curvilinear meshes. Journal of Computational Physics, 378:366 – 393, 2019. [13] Jared Crean, Jason E Hicken, David C Del Rey Fernández, David W Zingg, and Mark H Carpenter. High-Order, Entropy-Stable Discretizations of the Euler Equations for Complex Geometries. In 23rd AIAA Computational Fluid Dynamics Conference. American Institute of Aeronautics and Astronautics, 2017. [14] Matteo Parsani, Mark H Carpenter, Travis C Fisher, and Eric J Nielsen. Entropy Stable Staggered Grid Discontinuous Spectral Collocation Methods of any Order for the Compressible Navier–Stokes Equations. SIAM Journal on Scientific Computing, 38(5):A3129–A3162, 2016. [15] David A Kopriva and Gregor Gassner. On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods. Journal of Scientific Computing, 44(2):136–155, 2010. [16] 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, 2012. [17] Jesse Chan, Zheng Wang, Axel Modave, Jean-Francois Remacle, and T Warburton. GPUaccelerated discontinuous Galerkin methods on hybrid meshes. Journal of Computational Physics, 318:142–168, 2016. [18] David C Del Rey Fernández, Jason E Hicken, and David W Zingg. Review of summationby-parts operators with simultaneous approximation terms for the numerical solution of partial differential equations. Computers & Fluids, 95:171–196, 2014. [19] David C Del Rey Fernández, 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, 2018. [20] David C Del Rey Fernández, Pieter D Boom, and David W Zingg. A generalized framework for nodal first derivative summation-by-parts operators. Journal of Computational Physics, 266:214–239, 2014. [21] Sigrun Ortleb. Kinetic energy preserving DG schemes based on summation-by-parts operators on interior node distributions. PAMM, 16(1):857–858, 2016. [22] Sigrun Ortleb. A Kinetic Energy Preserving DG Scheme Based on Gauss–Legendre Points. Journal of Scientific Computing, 71(3):1135–1168, 2017. [23] Hendrik Ranocha. Generalised summation-by-parts operators and variable coefficients. Journal of Computational Physics, 362:20 – 48, 2018. [24] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012. [25] Andrew R Winters, Dominik Derigs, Gregor J Gassner, and Stefanie Walch. A uniquely defined entropy stable matrix dissipation operator for high Mach number ideal MHD and compressible Euler simulations. Journal of Computational Physics, 332:274–289, 2017. [26] Eitan Tadmor. Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. Acta Numerica, 12:451–512, 2003. [27] Lucas Friedrich, Gero Schnücke, Andrew R Winters, David C Fernández, Gregor J Gassner, and Mark H Carpenter. Entropy Stable Space-Time Discontinuous Galerkin Schemes with Summation-by-Parts Property for Hyperbolic Conservation Laws. arXiv preprint arXiv:1808.08218, 2018. [28] Jan Nordström. Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation. Journal of Scientific Computing, 29(3):375–404, 2006. [29] David A Kopriva and Gregor J Gassner. Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable. Applied Mathematics and Computation, 272:274–290, 2016. [30] David A Kopriva. Metric identities and the discontinuous spectral element method on curvilinear meshes. Journal of Scientific Computing, 26(3):301–327, 2006. [31] Amiram Harten. On the symmetric form of systems of conservation laws with entropy. Journal of computational physics, 49(1):151–164, 1983. [32] 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. Computer Methods in Applied Mechanics and Engineering, 54(2):223–234, 1986. [33] Farzad Ismail and Philip L Roe. Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks. Journal of Computational Physics, 228(15):5410–5436, 2009. [34] Praveen Chandrashekar. Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier-Stokes equations. Communications in Computational 28 Physics, 14(5):1252–1286, 2013. [35] Mark H Carpenter and Christopher A Kennedy. Fourth-order 2n-storage Runge-Kutta schemes. Technical Report NASA-TM-109112, NAS 1.15:109112, NASA Langley Research Center, 1994. [36] Chi-Wang Shu. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In Advanced numerical approximation of nonlinear hyperbolic equations, pages 325–432. Springer, 1998. [37] Alfio Quarteroni and Alberto Valli. Numerical Approximation of Partial Differential Equations. Springer, 1994. [38] Mo Lenoir. Optimal isoparametric finite elements and error estimates for domains involving curved boundaries. SIAM Journal on Numerical Analysis, 23(3):562–580, 1986. [39] T Warburton. A low-storage curvilinear discontinuous Galerkin method for wave problems. SIAM Journal on Scientific Computing, 35(4):A1987–A2012, 2013. [40] David M Williams and Antony Jameson. Nodal points and the nonlinear stability of high-order methods for unsteady flow problems on tetrahedral meshes. In 21st AIAA Computational Fluid Dynamics Conference. American Institute of Aeronautics and Astronautics, June 2013. [41] Magnus Svärd and Hatice Özcan. Entropy-stable schemes for the Euler equations with far-field and wall boundary conditions. Journal of Scientific Computing, 58(1):61–89, 2014. [42] Geoffrey I Taylor and Albert E Green. Mechanism of the production of small eddies from large ones. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 158(895):499–521, February 1937. [43] Niklas Wintermeyer, Andrew R Winters, Gregor J Gassner, and Timothy Warburton. An entropy stable discontinuous Galerkin method for the shallow water equations on curvilinear meshes with wet/dry fronts accelerated by GPUs. Journal of Computational Physics, 375:447–480, 2018. 29