Efficient Entropy-Stable Discontinuous Spectral-Element Methods Using Tensor-Product Summation-by-Parts Operators on Triangles and Tetrahedra

We present a new class of efficient and robust discontinuous spectral-element methods of arbitrary order for nonlinear hyperbolic systems of conservation laws on curved triangular and tetrahedral unstructured grids. Such discretizations employ a recently introduced family of sparse tensor-product summation-by-parts (SBP) operators in collapsed coordinates within an entropy-stable modal formulation. The proposed algorithms exploit the structure of such SBP operators alongside that of the Proriol-Koornwinder-Dubiner polynomial basis, and a weight-adjusted approximation is used to efficiently invert the local mass matrix for curvilinear elements. Using such techniques, the number of required entropy-conservative two-point flux evaluations between pairs of quadrature nodes is significantly reduced relative to existing entropy-stable formulations using (non-tensor-product) multidimensional SBP operators, particularly for high polynomial degrees, with an improvement in time complexity from $\mathcal{O}(p^{2d})$ to $\mathcal{O}(p^{d+1})$, where $p$ is the polynomial degree of the approximation and $d$ is the number of spatial dimensions. In numerical experiments involving smooth solutions to the compressible Euler equations, the proposed tensor-product schemes demonstrate similar levels of accuracy for a given mesh and polynomial degree to those using multidimensional SBP operators based on symmetric quadrature rules. Furthermore, both operator families are shown to give rise to entropy-stable methods which exhibit excellent robustness for under-resolved problems. Such results suggest that the algorithmic advantages resulting from the use of tensor-product operators are obtained without compromising accuracy or robustness, enabling the efficient extension of the benefits of entropy stability to higher polynomial degrees than previously considered for triangular and tetrahedral elements.


Introduction
Hyperbolic or advection-dominated nonlinear systems of conservation laws constitute a class of partial differential equations (PDEs) of considerable importance in numerous scientific and engineering disciplines, with applications ranging from aircraft design to climate modeling.As a result of their complex multiscale behaviour and propensity to develop discontinuities, such PDEs present a significant challenge in the design of efficient, automated, and robust numerical methods.High-order discontinuous spectral-element methods (DSEMs) 1 have emerged as an attractive approach for this class of problem, having received considerable attention in recent years due to their performance on modern hardware (see, for example, Klöckner et al. [1], Abdi et al. [2], and Vermiere et al. [3]) resulting from their relatively high arithmetic intensity (i.e. the ratio of floating-point operations to memory accesses) and data locality.Moreover, such schemes are highly flexible in their support for local variation in element size as well as polynomial degree, facilitating the adaptive solution of complex problems in an efficient and automated manner, as exemplified in recent papers by Parsani et al. [4] and Mossier et al. [5].
High-order methods such as DSEMs notoriously lack robustness when applied to strongly nonlinear problems, particularly in the presence of under-resolved scales or discontinuous solutions.While such issues are traditionally addressed through ad hoc stabilization techniques such as over-integration, modal filtering, slope limiting, or numerical dissipation, modern entropy-stable formulations enable rigorous proofs of stability (under certain physical admissibility criteria) by guaranteeing that a strictly convex function of the numerical solution remains bounded for all time.For example, in the context of fluid dynamics, an entropy-stable discretization can be constructed so as to satisfy the Second Law of Thermodynamics in a discrete sense, a property which bounds the growth of the numerical solution, provided that the pressure and density remain positive.Entropy-stable schemes were introduced for nonlinear hyperbolic systems of conservation laws by Tadmor [6], who devised first-order and second-order methods using specially designed two-point fluxes which conserve a particular mathematical entropy function.LeFloch et al. [7] extended Tadmor's approach to high-order accuracy on periodic domains in one dimension.The modern era of high-order entropystable discretizations, however, began when Fisher [8] combined entropy-conservative two-point flux functions with summation-by-parts (SBP) operators, which are discrete differential operators mimetic of integration by parts on bounded domains (see, for example, the review papers by Svärd and Nordström [9] and Del Rey Fernández et al. [10]).Fisher's approach (see also Fisher and Carpenter [11] and Fisher et al. [12]) enabled the construction of entropy-stable high-order finite-difference methods for the compressible Euler and Navier-Stokes equations on curvilinear block-structured grids using affordable entropy-conservative two-point flux functions proposed by Ismail and Roe [13].Such a combination of SBP operators with two-point flux functions came to be known in the entropy stability community as flux differencing (not to be confused with the similarly named flux-difference splitting technique introduced by Roe [14] decades earlier for approximate Riemann solvers).
The extension of entropy stability to DSEM formulations was made possible through the work of Gassner [15], who recognized that the matrix operators employed within discontinuous Galerkin (DG) methods based on collocated Legendre-Gauss-Lobatto quadrature were, in fact, SBP operators, an equivalence which he exploited in the construction of provably stable discretizations of Burgers' equation using split forms first proposed in the finite-difference community.Entropy-stable DSEMs for systems of conservation laws on tensor-product quadrilateral and hexahedral elements were then introduced by Carpenter et al. [16] and Gassner et al. [17], with the latter demonstrating that a broad class of split-form and entropy-stable DSEMs could be recovered by choosing different two-point flux functions.These schemes are distinguished from the approach taken by Barth [18] as well as Hiltebrand and Mishra [19] based on the work of Hughes et al. [20], wherein space-time DG schemes are formulated in terms of the entropy variables.Unlike flux-differencing approaches, the latter methodology results in discretizations which are only entropy stable under the assumption that all integrals in the corresponding variational formulation are evaluated exactly (which is impractical, if not impossible, for many PDEs of interest to practitioners, including the compressible Euler and Navier-Stokes equations) and, furthermore, cannot be formulated explicitly in time.
Entropy-stable DSEMs were generalized beyond tensor-product quadrature rules on quadrilaterals and hexahedra through the use of multidimensional SBP operators, which were introduced by Hicken et al. [21] as an extension of the generalized SBP framework proposed in a one-dimensional setting by Del Rey Fernández et al. [22].The subsequent development of entropy-stable high-order methods on triangles and tetrahedra (which first appeared in papers by Chen and Shu [23], Crean et al. [24], and Chan [25]) is outlined in a review paper by Chen and Shu [26].While advantageous in terms of geometric flexibility, such triangular and tetrahedral SBP operators do not possess a tensor-product structure, and, as such, result in schemes which become considerably more costly than comparable tensor-product DSEMs on quadrilaterals and hexahedra as the polynomial degree is increased.Such a difference is due in part to the fact that tensor-product spectral-element operators are amenable to sum-factorization techniques, which were proposed in the context of spectral methods by Orszag [27].Sum factorization, which involves the application of tensor-product operators in a dimension-by-dimension manner, results in the number of floating-point operations required to evaluate spatial operators scaling asymptotically as O(p d+1 ) with the polynomial degree p, where d denotes the number of spatial dimensions.This compares favourably to the O(p 2d ) complexity obtained using multidimensional operators without a tensor-product structure.In the context of an entropy-stable scheme, multidimensional SBP operators are further disadvantaged due to the fact that they require the evaluation of entropy-conservative two-point flux functions (which, in the case of the Euler and Navier-Stokes equations, are relatively expensive operations involving the logarithmic mean) between all pairs of quadrature nodes, rather than simply along lines of nodes as in a tensor-product formulation.Due in part to these limitations, as well as the difficulty in constructing suitable quadrature rules, entropy-stable discretizations based on multidimensional (i.e.non-tensor-product) SBP operators rarely employ polynomial degrees greater than four or five.
In a recent paper [28], the present authors introduced tensor-product SBP operators on triangles and tetrahedra based on collapsed coordinate systems, which were used within skew-symmetric nodal and modal DSEM formulations to obtain energy-stable discretizations of the linear advection equation on curved elements.These contributions served to extend the SBP framework to tensor-product operators in collapsed coordinates, which had previously been recognized in the SEM community (see, for example, Sherwin and Karniadakis [29,30], Lomtev and Karniadakis [31], Kirby et al. [32], and Moxey et al. [33]) to be crucial for obtaining efficient discretizations on triangles and tetrahedra (as well as prisms and pyramids) at high polynomial degreees.Specifically, it was shown in [28] that by exploiting the tensor-product structure of the proposed SBP operators, and, in the case of the modal approach, making use of the "warped" tensor-product structure of the Proriol-Koornwinder-Dubiner (PKD) orthogonal polynomial basis [34][35][36] alongside a weight-adjusted approximation proposed by Chan et al. [37] to invert the curvilinear mass matrix, the time derivative could be obtained in O(p d+1 ) floating-point operations (i.e. the same asymptotic complexity as for the DSEMs on quadrilaterals and hexahedra described above) through sum factorization.Moreover, the modal formulations were shown to be similar in accuracy and in spectral radius (which, for linear problems, dictates the maximum stable time step for explicit temporal integration) to those using multidimensional operators on triangles and tetrahedra, while requiring far fewer floating-point operations at higher polynomial degrees.
Building upon the above contributions, the primary objective of this paper is to apply the tensor-product SBP operators on triangles and tetrahedra introduced in [28] to the construction of efficient entropy-stable DSEMs of arbitrary order for nonlinear hyperbolic systems of conservation laws.The proposed methods make use of a modal PKD basis and a weight-adjusted approximation of the mass matrix inverse within an entropy-stable flux-differencing formulation, in which we exploit operator sparsity as well as sum factorization in order to obtain O(p d+1 ) complexity for all local elemental operations.In the tetrahedral case, we also approximate the metric terms arising from the mapping from reference to physical coordinates using a curl formulation from Chan and Wilcox [38] to satisfy the discrete metric identities, which must hold in order to ensure free-stream preservation and entropy stability on curvilinear meshes.These contributions enable the efficient extension of entropy-stable schemes on triangles and tetrahedra to polynomial degrees beyond those typically considered for such element types.
We now outline the structure of the remainder of this paper.In Section 2, we introduce some relevant notation and review some fundamental concepts relating to hyperbolic systems of conservation laws and SBP operators as well as classical orthogonal polynomials and Gaussian quadrature rules.In Section 3, we review the relevant aspects of the approach introduced in [28] for the construction of tensor-product SBP operators on triangles and tetrahedra.Section 4 describes a general framework for the construction of entropy-stable DSEMs using tensor-product or multidimensional SBP operators on curved triangular or tetrahedral elements, while Section 5 describes algorithmic considerations for the efficient implementation of the proposed schemes.The methods are then analyzed with respect to their conservation, free-stream preservation, and entropy stability properties in Section 6. Numerical studies of accuracy as well as robustness for the compressible Euler equations are presented in Section 7, and conclusions are provided in Section 8.

Notation
The notation in this paper follows that introduced by the authors in [39], [40], and [28].Single underlines are used to denote vectors (treated as column matrices), whereas double underlines denote matrices.Symbols in bold are used specifically to denote Cartesian (i.e.spatial) vectors, for which we employ the usual dot product x • y := x 1 y 1 + • • • + x d y d and Euclidean norm ∥x∥ 2 := x • x.The symbol ∇ is used to denote componentwise differentiation with respect to either type of vector, for example, as and S d−1 denote the real numbers, the positive real numbers, the non-negative real numbers, the natural numbers (excluding zero), the natural numbers including zero, and the unit (d − 1)-sphere S d−1 := {x ∈ R d : ∥x∥ = 1}, respectively.The symbols 0 (N ) and 1 (N ) are reserved for vectors of length N ∈ N containing all zeros and all ones, respectively, where the superscript is omitted when clear from the context, and the symbol 0 (M ×N ) likewise denotes an M by N matrix of zeros.The notation {1 : N } is used as shorthand for the index set {1, 2, . . ., N }.Given an arbitrary bounded domain D ⊂ R d , we use ∂D to denote its boundary and D := D ∪ ∂D to denote its closure; the interior of a closed domain D is then given by D := D \ ∂D.The space of polynomials of maximum total degree p ∈ N 0 on D is then defined as Other relevant notational conventions and definitions are introduced as they appear.

Systems of conservation laws and the entropy inequality
We are interested in systems of conservation laws governing the evolution of subject to appropriate boundary conditions, where F m (U (x, t)) ∈ R Nc is the m th Cartesian flux component and U 0 (x) ∈ Υ represents the initial data, where Υ denotes the set of admissible solution states.For convenience, we define the flux in any direction n ∈ S d−1 for an arbitrary state U ∈ Υ as A system of conservation laws in the form of ( 1) is then called hyperbolic if the flux Jacobian ∇ U F (U , n) ∈ R Nc×Nc is diagonalizable with all real eigenvalues for all U ∈ Υ and n ∈ S d−1 .
Remark 1.In this work, we tacitly assume that numerical solutions remain within the admissible set Υ, corresponding, for example, to the requirement for thermodynamic quantities such as pressure or density to be positive.This is difficult to ensure a priori for nonlinear problems, and often requires the use of bespoke limiting procedures, which are described in the context of entropy-stable schemes, for example, in recent papers by Rueda-Ramírez et al. [41] and Lin et al. [42].
We are specifically interested in hyperbolic systems of conservation laws endowed with an entropy function and corresponding entropy flux, which are defined as follows.
Definition 1.The functions S : Υ → R and F : Υ → R d are, respectively, an entropy function and entropy flux if S is strictly convex (i.e. its Hessian is positive definite), and the relation holds for all m ∈ {1 : d}, where ∇ U S(U ), ∇ U F m (U ) ∈ R Nc are the gradients of the entropy and entropy flux components with respect to the conservative variables, respectively, and ∇ U F m (U ) ∈ R Nc×Nc denotes the Jacobian of the m th flux component with respect to the conservative variables.The entries of the vector W(U ) := ∇ U S(U ) are referred to as the entropy variables, where the mapping W has an inverse given by U due to the strict convexity of S over the admissible set Υ.As described by Friedrichs and Lax [43], the existence of an entropy-entropy flux pair in the sense of Definition 1 implies that any classical (i.e.continuously differentiable) solution to (1) satisfies an additional conservation equation of the form Integrating (4) over the spatial domain and using the divergence theorem then results in where n(x) ∈ S d−1 denotes the outward unit normal to ∂Ω, and we have defined the flux potential ).We are, however, often interested in weak solutions, which satisfy (1) in the sense of distributions and can therefore be discontinuous, describing phenomena such as shock waves.Replacing (5) with an entropy inequality of the form then provides an admissibility criterion for physically relevant weak solutions (see, for example, Kružkov [44] or Lax [45]).Provided that U (x, t) remains within Υ and that the boundary conditions are imposed correctly, it can be shown (see, for example, Dafermos [46]) that ( 6) implies a bound on the solution itself due to the strict convexity of the entropy function.We are therefore interested in constructing discretizations which respect a semi-discrete form of such an entropy bound.
Remark 2. The requirement for a strictly convex entropy function in Definition 1 is consistent with the mathematical literature, but is opposite the convention used in physics.As such, the inequality in ( 6) is a generalized statement of the Second Law of Thermodynamics, up to a change in sign.

Summation-by-parts operators
The discretizations described in this paper involve the construction of SBP operators on a canonical reference element, which are then used to discretize the PDE on an unstructured mesh through the use of a bijective mapping from the reference element to each physical element.We now require the following definition of a nodal SBP operator from [21], which extends the generalized definition proposed in [22] to the multidimensional setting.Definition 2. Let Ω ⊂ R d denote a closed, bounded, and connected reference domain on which we define a set of N q ∈ N distinct nodes {ξ (i) } i∈{1:Nq} ⊂ Ω, and let the vectors u := U (ξ (1) ), . . ., U (ξ (Nq) ) T and v := V (ξ (1) ), . . ., V (ξ (Nq) ) T (7) contain the nodal values of arbitrary functions U, V : Ω → R. A matrix D (m) ∈ R Nq×Nq approximating the partial derivative ∂/∂ξ m is then a nodal SBP operator of degree q ∈ N 0 if it satisfies and may be decomposed as , where holds for some r ≥ q, and n : ∂ Ω → S d−1 denotes the outward unit normal vector to Ω.
The SBP property is equivalent to the decomposition , where S (m) is skew-symmetric and E (m) is a symmetric matrix satisfying the conditions in (9).Moreover, SBP operators mimic integration by parts in a discrete sense, as given by We refer to any SBP operator for which the associated SPD matrix W is diagonal as a diagonal-norm SBP operator.In such a case, it was shown in [21, Theorem 3.2] that the diagonal entries of W constitute the weights {ω (i) } i∈{1:Nq} ⊂ R + for a quadrature rule satisfying which is of degree τ ≥ 2q − 1 for any diagonal-norm SBP operator of degree q.
Remark 3. The polynomial exactness condition in (8) implies that any SBP operator of degree q is, by definition, also an SBP operator of any positive integer less than q.Where such ambiguity arises, the degree q of an SBP operator is taken to refer uniquely to the maximum integer value for which (8) holds, and we adopt an analogous convention for quadrature rules such as that in (11).
Remark 4.Although a multidimensional discretization requires an SBP operator of the form to be constructed for each partial derivative ∂/∂x m , it will be assumed throughout this work, as is typical in the SBP literature, that W is the same for all coordinate indices m ∈ {1 : d}.

Decomposition of the boundary operators
Let us now assume that the reference element Ω ⊂ R d is a polytope, corresponding to a polygon in two dimensions or a polyhedron in three dimensions, and partition its boundary into N f ∈ N closed subsets { Γ(ζ) } ζ∈{1:N f } with disjoint interiors, which we denote as facets.Each facet Γ(ζ) ⊂ ∂ Ω is assumed to be flat, and we denote its (constant) outward unit normal vector by n (ζ) ∈ S d−1 .On each facet, we then introduce N (ζ) q f ∈ N quadrature nodes and weights given, respectively, by As in Del Rey Fernández et al. [47,Section 3], we restrict our attention to the class of SBP operators for which the boundary matrices in ( 9) can be constructed as where q f is a diagonal matrix with entries ×Nq is an interpolation/extrapolation operator of degree r (ζ) ≥ p, satisfying the accuracy conditions As a special case of the decomposition in ( 13), we note that E (m) can be made diagonal by constructing SBP operators for which the nodal sets in (12) all form subsets of the volume quadrature nodes, wherein R (ζ) simply selects boundary values from nodal vectors such as those in (7).Such SBP operators, denoted here as diagonal-E operators, are analogous to those on one-dimensional nodal sets including both endpoints, and were first constructed on triangles by Chen and Shu [23].

Orthogonal polynomials and Gaussian quadrature rules
The fundamental building blocks used for constructing the triangular and tetrahedral SBP operators developed by the authors in [28] are Jacobi and Legendre polynomials as well as their associated Gaussian quadrature rules and interpolants, the basic properties of which we will review here.The normalized Jacobi polynomials are denoted by Such polynomials can be constructed through recurrence relations, as shown, for example, by Hesthaven and Warburton [48,Appendix A].For a given non-negative integer q, the Gaussian quadrature rules corresponding to a Jacobi weight with exponents a and b have nodes {η 1] given by the q + 1 solutions to a polynomial equation.Such equations are given by Gauss: Gauss-Radau: (1 + η)P (a,b+1) Gauss-Lobatto: where the Gauss, Gauss-Radau, and Gauss-Lobatto families of quadrature rules include zero, one, and two endpoints of the interval, respectively. 2The Lagrange polynomials {ℓ (a,b) q,i } i∈{0:q} associated with such a nodal set constitute a basis for The corresponding Gaussian quadrature weights {ω q,i } i∈{0:q} can then be expressed as where alternative expressions for such weights can be found, for example, in references such as Karniadakis and Sherwin [49,Appendix B].The resulting quadrature rule satisfies

Tensor-product summation-by-parts operators in collapsed coordinates
In an earlier paper [28], the present authors proposed a procedure for constructing SBP operators of arbitrary degree on triangles and tetrahedra, which, unlike those developed previously, are sparse and possess a tensor-product structure, allowing for the application of sum-factorization techniques.In this section, we briefly review the essential steps involved in constructing such operators, where we refer to the above paper for details regarding the construction of SBP operators and more broadly to Karniadakis and Sherwin [49] regarding SEM formulations in collapsed coordinates.

Tensor-product summation-by-parts operators on the reference triangle
The reference element on which we construct our triangular SBP operators is taken to be where, as a convention, we number the facets (i.e.edges of the triangle) counter-clockwise as The collapsed coordinate transformation from the square to the reference triangle is then given by which is illustrated in Figure 1.Letting q m ∈ N 0 denote the degree of the approximation in the η m coordinate, tensor-product quadrature rules are obtained with nodes {ξ (i) } i∈{1:Nq} ⊂ Ω and corresponding weights {ω (i) } i∈{1:Nq} ⊂ R + given for a multi-index α ∈ {0 : q 1 } × {0 : q 2 } as where σ : {0 : q 1 } × {0 : q 2 } → {1 : N q } is a bijective mapping which defines an ordering of the N q := (q 1 + 1)(q 2 + 1) volume quadrature nodes, and we use the notation introduced in Section 2.5 for Gaussian quadrature rules.Similarly, we let q f ∈ N 0 and define the N (ζ) q f := q f + 1 facet quadrature nodes and weights in (12) as ξ (1,i)  (2,i)  (3,i) Defining the volume and facet quadrature weight matrices W ∈ R Nq×Nq and q f with entries given by respectively, derivative operators D (m) ∈ R Nq×Nq can be constructed by applying the chain rule to a tensor-product Lagrange interpolant in collapsed coordinates, resulting in Similarly, the matrices ×Nq have entries given by In this work, we construct SBP operators of degree q ∈ N 0 using LG quadrature rules in all directions, with q 1 = q 2 = q f = q.As a consequence of [28,Theorem 3.1], such operators are guaranteed to satisfy the conditions of Definition 2, with boundary operators in the form of (13).

Tensor-product summation-by-parts operators on the reference tetrahedron
For tetrahedral elements, the reference domain is given by with the facets (i.e.faces of the tetrahedron) numbered as The collapsed coordinate transformation χ : [−1, 1] 3 → Ω from the cube to the tetrahedron, as depicted in Figure 2, is constructed from three successive applications of ( 22), resulting in Similarly to the triangular case, we let q m ∈ N 0 denote the degree of the approximation in the η m coordinate and introduce N q := (q 1 + 1)(q 2 + 1)(q 3 + 1) tensor-product volume quadrature nodes, which are ordered using the bijective mapping σ : {0 : q 1 } × {0 : q 2 } × {0 : q 3 } → {1 : N q }.Taking (a, b) = (0, 0) in the η 1 and η 2 directions and (a, b) = (1, 0) in the η 3 direction, we obtain Defining a collapsed coordinate system (η f 1 , η f 2 ) on each face of the tetrahedron, we let q f 1 , q f 2 ∈ N 0 and introduce the bijective mapping σ f : {0 :

−→ χ
, the facet quadrature nodes are given by and the corresponding weights are The entries of the differentiation and interpolation/extrapolation matrices are then given by and respectively, where the diagonal matrices W and B (ζ) are defined similarly to the triangular case using the quadrature weights in ( 30) and (32), respectively.In this paper, we use Gauss quadrature rules with respect to the Jacobi weights indicated in the superscripts, and we take q 1 , q 2 , q 3 , q f 1 , and q f 2 to be all equal to q ∈ N 0 .It then follows from [28,Theorem 4.1] that such choices result in diagonal-norm SBP operators of degree q for which the boundary matrices are given as in (13).
Remark 5.The choices of ω (σ(α)) := ω q3,α3 /8 for triangles and tetrahedra, respectively, are made in [29] and [30] in order to use the Jacobi weight to subsume the factor in ( 23) or (30) resulting from the Jacobian determinant of the collapsed coordinate transformation.However, such Jacobi weights lead to nodal derivative operators in (25) and (33) which do not, in general, satisfy the SBP property on the reference triangle and tetrahedron, respectively.The operators introduced in [28], which are constructed specifically to satisfy the SBP property on the reference simplex, thus differ from existing SEM operators in collapsed coordinates.

Entropy-stable discontinuous spectral-element framework
In this section, we detail how the SBP operators described in the previous section can be used to construct entropy-stable DSEMs of any order on curvilinear unstructured grids.Due to the generality of the SBP approach, however, the formulations which we present can be used with any set of diagonal-norm SBP operators on the reference triangle or tetrahedron for which the boundary matrices can be decomposed as in (13).While we introduce the proposed schemes within the context of a generalized framework in order to enable the unified implementation and analysis of DSEMs using SBP operators and to facilitate the comparison of the proposed methods to existing schemes, efficient algorithms require consideration of the particular properties of the operators which constitute a discretization, a topic which we will address in Section 5.

Mesh and coordinate transformation
The first step in constructing a DSEM is to subdivide the spatial domain Ω into a mesh or grid, which consists of a collection {Ω (κ) } κ∈{1:Ne} of N e ∈ N closed, bounded, and connected elements Ω (κ) ⊂ Ω with nonempty interiors, satisfying Ne κ=1 We denote the characteristic element size for such a mesh as h ∈ R + .In this work, we assume that the mesh is time invariant and that each element is the image of the reference triangle or tetrahedron under a polynomial mapping X (κ) ∈ [P pg ( Ω)] d of total degree p g ∈ N.Such a mapping is given in terms of a multivariate Lagrange basis {ℓ where {x pg } i∈{1:N * pg } are the prescribed physical positions of the mapping nodes.To ensure a watertight mesh, we assume that the mapping nodes contain a subset of nodes on each facet Γ(ζ) ⊂ ∂ Ω which are unisolvent for the corresponding trace space P pg ( Γ(ζ) ) := {V | Γ(ζ) : V ∈ P pg ( Ω)}, thus resulting in continuity at element interfaces.Denoting the Jacobian of the transformation by ∇ ξ X (κ) (ξ) ∈ R d×d and defining J (κ) (ξ) := det(∇ ξ X (κ) (ξ)), we assume that the mapping is bijective and orientation preserving, satisfying The adjugate of the Jacobian matrix is given by G (κ) (ξ) := det(∇ ξ X (κ) (ξ))(∇ ξ X (κ) (ξ)) −1 with entries (often referred to as the metric terms in the literature) satisfying the metric identities, Using the metric identities, we can express (1a) in conservation form on the reference element as Further details regarding the formulation of conservation laws in curvilinear coordinates are provided, for example, in Pulliam and Zingg [50, Section 4.2] and Kopriva [51, Section 6.2].

Approximation of the metric terms
Considering a mapping using polynomials of degree p g as in (36), the metric terms are polynomials of degree p g − 1 in two dimensions and degree 2p g − 2 in three dimensions.Since operations such as differentiation using SBP operators are exact for polynomials of at most degree q, we cannot expect that a discrete analogue of the metric identities in (38) will hold unless p g ≤ q + 1 in two dimensions or p g ≤ ⌊q/2⌋ + 1 in three dimensions.To circumvent this requirement for a subparametric mapping in three-dimensional case, we use an adaptation by Chan and Wilcox [38, Section 5] of Kopriva's approximation of the metric terms in conservative curl form [52, Eq. ( 36)], which is itself based on techniques introduced by Thomas and Lombard [53].To obtain such an approximation, we introduce a Lagrange basis {ℓ and construct polynomial interpolants of the form r (κ,1) The above functions r (κ,m) q+1 ∈ [P q+1 ( Ω)] 3 are used to define the matrix of approximate metric terms which has entries of degree q and satisfies (38) by construction.Using the exact metric terms in two dimensions and the approximation (41) in three dimensions, we can compute the outward unit normal vector to the facet Γ (κ,ζ) ⊂ ∂Ω (κ) which is the image of Γ(ζ) under the mapping X (κ) as Remark 6.In general, the normals computed as in (42) are not exact when the metric terms are computed approximately using (41).However, if the analytically defined mesh is watertight and the nodes used for the interpolants in (40) define a continuous approximation space, the approximate normals remain equal and opposite at element interfaces (see, for example, [38,Theorem 5]).

Modal polynomial expansions on triangles and tetrahedra
The focus of this paper is on modal formulations, in which the degrees of freedom for the semidiscrete approximation to the e th solution variable on the k th element are taken to be the expansion coefficients with respect to a polynomial basis {ϕ (i) } i∈{1:N * p } of degree p ≤ q contained within the vector ũ(h,κ,e) (t) ∈ R N * p .Such coefficients define an approximate solution resulting in a global approximation given piecewise by U (x, t) ≈ U h (x, t) := U (h,κ) ((X (κ) ) −1 (x), t) for x ∈ Ω (κ) .Evaluating (43) at each volume quadrature node, the vector u (h,κ,e) (t) ∈ R Nq of nodal solution values can be expressed in terms of the generalized Vandermonde matrix In order to construct polynomial bases on the triangle and tetrahedron for which operations such as (44) are amenable to sum factorization when tensor-product quadrature rules in collapsed coordinates are used, we follow [49, Section 3.2] and define the one-dimensional principal functions The normalized Proriol-Koornwinder-Dubiner (PKD) orthogonal polynomials [34][35][36] are then given in terms of collapsed coordinates on the reference triangle as and on the reference tetrahedron as where we order α ∈ N (p) using the bijection π : N (p) → {1 : N * p }.The "warped" tensor-product structure of the PKD basis allows for such a matrix or its transpose to be applied in O(p d+1 ) operations through sum factorization when tensor-product quadrature rules are used with O(p) nodes in each direction, as described, for example, in [49, Sections 4.1.6.1 and 4.1.6.2].Moreover, since the principal functions in (45) have been scaled to obtain an orthonormal basis, the reference mass matrix M := V T W V is the identity matrix if the quadrature rule in ( 11) is of degree τ ≥ 2p.Remark 7. A nodal formulation is recovered by taking V to be the identity matrix and directly evolving the nodal solution vector u (h,κ,e) (t).Such collocation-based approaches using tensor-product LGL [16,17] or LG [54] quadrature rules are popular for quadrilateral and hexahedral elements due to the numerical solution being available at the volume quadrature nodes (and, in the case of LGL quadrature, the facet quadrature nodes) without the need to perform the matrix operation in (44).When using tensor-product SBP operators on triangular and tetrahedral elements, however, we are primarily interested in modal formulations, as, unlike nodal formulations based on collapsed coordinates, they are not subject to the severe explicit time step restriction resulting from the singularity of the mapping in (22) or (29).As discussed by the authors in [28,Section 7.4] in the context of the linear advection equation, the spectral radii for the modal tensor-product DSEMs on triangles and tetrahedra scale at most quadratically with the polynomial degree, similarly to DSEMs using multidimensional SBP operators based on symmetrical nodal sets.Additionally, modal formulations allow for the decoupling of the degree p of the solution polynomial from the SBP operator degree q and volume quadrature degree τ , which enables the use of over-integration while retaining a minimal number of solution degrees of freedom for a given order of accuracy.

Discontinuous spectral-element methods using summation-by-parts operators
A weak-form DG approximation of (1) can be derived either by integrating by parts against a smooth test function on the physical element and performing a change of variables within each of the resulting integrals, or by integrating the transformed conservation law in (39) by parts against a smooth test function on the reference element.In either case, we approximate each solution variable as in (43) and consider test functions belonging to the same space, resulting in where a numerical flux function F * : Υ×Υ×S d−1 → R Nc has been used to resolve the discontinuity in the global numerical solution at each facet, on which we denote the exterior solution as U (h,κ,+) (ξ, t) ∈ Υ.In order to obtain an algebraic formulation of ( 48), we begin by forming the diagonal matrices q f with entries given by Defining the physical mass matrix as M (κ) := V T W J (κ) V , we then obtain a semi-discrete formulation governing the evolution of the vector of modal coefficients for the expansion in (43) as where the form of r (h,κ,e) (t) ∈ R Nq determines the particular DSEM recovered within the present framework.Evaluating the nodal solution as u (h,κ,e) (t) := V ũ(h,κ,e) (t), we gather the solution variables at the volume and facet quadrature nodes as and form the vectors f (κ,m,e) (t) ∈ R Nq and f ( * ,κ,ζ,e) ∈ R q f with entries given by where u (t) ∈ Υ denotes the exterior solution state.With such definitions in place, the right-hand side of (48) can be discretized directly using nodal SBP operators as where it follows from a similar analysis to [40, Section 4] that the resulting scheme is conservative for an appropriate choice of numerical flux as well as energy stable for linear, constant-coefficient problems on meshes for which the mapping from reference to physical space is affine. 3o ensure energy stability for linear problems on curvilinear meshes, a property which cannot be guaranteed for discretizations in the form of ( 53), the schemes proposed by the authors in [39] and [28] employ a skew-symmetric weak formulation adapted from the work of Gassner and Kopriva [55] and Del Rey Fernández et al. [47,56], which is given by In Section 4.6, we will construct modifed formulations, which, unlike (53) or (54), are provably entropy stable for nonlinear systems of conservation laws endowed with suitable entropy functions and entropy fluxes, satisfying semi-discrete bounds analogous to (6).However, we must first address the issue of inverting the mass matrix in (50) in order to obtain the time derivative.

Weight-adjusted approximation of the inverse mass matrix
Even when using orthonormal bases on the reference element such as those presented in Section 4.3, the mass matrix M (κ) appearing on the left-hand side of ( 50) is dense when the mapping from the reference element to the physical element is not affine, and its inverse lacks a tensor-product structure amenable to sum factorization.Obtaining the time derivative for such a scheme in the context of explicit temporal integration thus requires either the storage and application of a non-tensorial factorization or inverse, or, otherwise, the solution of a dense N * p by N * p linear system, for each curved element in the mesh.To obtain a fully explicit formulation for the time derivative in (50), we use a weight-adjusted approximation given by The above approximation was initially proposed by Chan et al. [37] for the purpose of reducing storage requirements for curved elements.However, it was also demonstrated in [28] that such an approach preserves the tensor-product operator structure which would otherwise be lost by taking the inverse of the mass matrix.The time derivative can then be computed explicitly as where we can exploit sum factorization in the application of the operators V and V T in (55), as discussed in Section 4.3.While the formulation in ( 56) is not, in general, discretely conservative with respect to the quadrature rule defined by the diagonal entries of W as in (11), we restore conservation using a technique proposed by Chan and Wilcox [38, Lemma 2].In the context of a mapping in the form of (36), such a modification involves approximating the determinant of the mapping Jacobian, which is of degree 2p g − 2 in two dimensions and 3p g − 3 in three dimensions, by an interpolant of degree p g given in terms of the nodal basis used in (36) as and using such an approximation to define J (κ) in (55), noting that such a modification does not affect the (approximate) metric terms G (κ) (ξ) used to compute the right-hand side of (50).

Entropy-stable flux-differencing formulation
In this section, we will develop a formulation of the entropy-stable modal approach introduced by Chan [25] which lends itself to an efficient implementation when used with the tensor-product operators on triangles and tetrahedra described in Section 3. We begin by considering a fundamental issue in the development of entropy-stable modal DSEMs, which is the fact that the entropy variables may not lie within approximation space in which the solution is sought.As such, they cannot be taken as test functions in the discrete variational formulation, a critical step in establishing entropy stability for such schemes.To resolve this, entropy-stable modal formulations employ an entropy projection, referring to the approximation of the entropy variables as where, as in [38, Eq. ( 31)], we obtain the modal coefficients through a weight-adjusted projection as Using the generalized Vandermonde matrix to evaluate the projected entropy variables at the volume quadrature nodes as w (h,κ,e) (t) := V w(h,κ,e) (t), we then define Next, we present the following definition of an entropy-conservative two-point flux, which is an essential component of the entropy-stable methods described in this work.

Definition 3.
A continuously differentiable bivariate function F # m : Υ × Υ → R Nc is an entropyconservative two-point flux if it is symmetric with respect to its two arguments, consistent with the PDE in (1a) such that F # m (U , U ) = F m (U ) for all U ∈ Υ, and satisfies First proposed in [6], the property in ( 61) is referred to in the literature as Tadmor's condition or the shuffle condition, and enables the chain rule to be circumvented when deriving semi-discrete forms of bounds such as (6) in the context of an entropy-stable discretization.At element interfaces, we use entropy-stable or entropy-conservative directional numerical fluxes, for which the following definition is introduced (see, for example, [23, Definitions 3.1 and 3.2]).

Definition 4. A directional numerical flux function F
it satisfies the conservation and consistency conditions given by respectively, as well as the entropy condition, which is defined analogously to (61) by Such a numerical flux is entropy conservative if (63) holds as an equality for all arguments.
Although other choices are possible, we assume here that the numerical interface flux consists of an entropy-conservative two-point flux in the normal direction augmented by a local Lax-Friedrichs dissipative term (see, for example, Ranocha [57, Section 6.1]).Such a flux takes the form where the first term on the right-hand side is an entropy-conservative directional flux given by and λ(U − , U + , n) ∈ R + is an estimate of the maximum wave speed in the normal direction.
In order to construct an expression for r (h,κ,e) (t) in (56) which results in an entropy-stable scheme, we require the evaluation of the averaged volume and surface metric terms as which define the matrices {{G q f , respectively.The entropy-conservative two-point flux functions are then computed using the entropy-projected conservative variables, a term introduced in [25] in reference to the conservative variables evaluated in terms of the projected entropy variables in (59), as given by which define the matrices F (κ,m,e) (t) ∈ R Nq×Nq and F (κ,ζ,m,e) (t) ∈ R q f .Defining the exterior values of the entropy variables as w Having introduced the essential components of the scheme, flux-differencing weight-adjusted modal DSEM is now obtained by computing the time derivative as in (56), with the nodal right-hand side computed as where ⊙ denotes the Hadamard product given by [A ⊙ B] ij := A ij B ij , and we define Remark 8.The above formulation can be used with any set of SBP operators for which the boundary matrices can be decomposed as in (13), and it will be shown in Section 6 that the formulation is, in fact, mathematically equivalent to that in [38, Eq. ( 35)] when the same SBP operators are used in both schemes, a fact which we will exploit in our analysis of conservation, free-stream preservation, and entropy stability.When diagonal-E operators are used, the terms in (69) involving the correction operator C (κ,ζ,e) (t) in ( 70) cancel, and hence the scheme becomes where the first equality recovers the optimized weak formulation from Ranocha et al. [58, Section 2.1], and we have used the SBP property to obtain the second equality, recovering a strong formulation consisting of a flux-differencing term plus a penalty term (i.e. a simultaneous approximation term in SBP terminology) as in the original flux-differencing approach proposed by Fisher [8].
Remark 9. When ( 1) is a linear symmetric hyperbolic system with a constant coefficient matrix, the mean-value flux )) conserves the quadratic entropy S(U ) := 1 2 U T U , and the resulting formulation in ( 69) is then equivalent to the energy-stable scheme in (54).

Efficient implementation
In this section, we discuss and analyze several important algorithmic considerations pertaining to the implementation of the proposed schemes, particularly regarding techniques for exploiting the sparsity and tensor-product structure of the SBP operators described in Section 3 within the context of an entropy-stable flux-differencing DSEM.The algorithms described here are implemented within the open-source Julia code StableSpectralElements.jl developed by the first author. 4

Exploiting operator sparsity in flux differencing
As discussed, for example, by Ranocha et al. [58,Figure 3], the cost of an entropy-stable scheme is dominated by the flux-differencing terms, for which the primary expense is the evaluation of two-point entropy-conservative flux functions between pairs of quadrature nodes.By rewriting the volume contributions appearing in the first term on the right-hand side of ( 69) or (71) as we observe, as noted in [58, Section 2.2], that due to the symmetry of {{G (κ) lm }} ⊙ F (κ,m,e) (t) and the skew-symmetry of S (l) , it is only necessary to iterate over the indices i and j corresponding to the strictly upper-triangular parts of such matrices.Moreover, the sum need only be taken over the indices for which S (l) ij ̸ = 0, and the corresponding values of the vector 4 StableSpectralElements.jl is available at https://github.com/tristanmontoya/StableSpectralElements.jl.
within the directional two-point flux can be computed on the fly in order to avoid storing the dense matrices {{G (κ) lm }} in memory.Computing the sum on the right-hand side of (72) for all i ∈ {1 : N q } therefore requires the evaluation of the two-point flux function, which for entropy-stable schemes is typically a relatively expensive operation involving the logarithmic mean, once per nonzero entry in the strictly upper-triangular part of S (l) , with the total number of required floating-point operations being proportional to such a quantity.Similarly, expressing the averaged metric terms in (66b) as we can evaluate the contributions from C (κ,ζ,e) (t)1 ) and (C (κ,ζ,e) (t)) T 1 (Nq) on the second line of (69) simultaneously by initializing both such vectors to zero and then iterating over values of i and j such that [(R (ζ) ) T B (ζ) ] ij ̸ = 0. Within each iteration, we compute the corresponding two-point flux and multiply each component by the corresponding nonzero matrix entry to obtain Such values are then accumulated within C (κ,ζ,e) (t)1 ) and (C (κ,ζ,e) (t)) T 1 (Nq) as before proceeding to the next nonzero entry of the matrix (R (ζ) ) T B (ζ) .Letting nnz(A) denote the number of nonzero entries in a given matrix A, it follows from the above discussion that the number of element-local directional two-point flux evaluations (i.e.neglecting the numerical interface flux, which is computed in a separate routine) required to evaluate the right-hand side in ( 69) can be computed as To the authors' knowledge, the matrices S (l) are dense for all existing high-order entropy-stable discretizations on triangles or tetrahedra.Recalling from [21] and [47] that the minimum number of volume quadrature nodes for an SBP operator of degree q is the dimension N * q of the associated total-degree polynomial space, which scales as O(q d ), the number of two-point fluxes, and hence the computational work required to evaluate the flux-differencing volume terms, is therefore expected to scale as O(q 2d ).By contrast, such matrices are sparse for the operators described in Section 3, with their one-dimensional coupling along lines of nodes resulting in the same O(q d+1 ) complexity as for tensor-product DSEMs on quadrilaterals or hexahedra.
Remark 10.While the second term in (77) indeed vanishes for diagonal-E multidimensional SBP operators due to the simplifications made on the first line of ( 71), such operators require quadrature rules using a much larger number of nodes for a given degree than would otherwise be needed, and are currently only available for modest polynomial degrees. 5For more general non-tensor-product quadrature rules without collocated facet nodes, the matrices (R (ζ) ) T B (ζ) are also dense, coupling every volume quadrature node to every facet quadrature node.

Exploiting sum factorization for tensor-product operators
In addition to the reduction in computational complexity of the flux-differencing terms, the algorithmic benefits of tensor-product operators discussed in [28,Section 6] in regard to the skewsymmetric scheme in (54) extend directly to the matrix-vector products in the proposed entropy-stable formulations, which involve V and R (ζ) as well as their transposes.As such, we are able to exploit the structure of such matrices through sum factorization in the evaluation of the conservative variables at the volume and facet quadrature nodes as u (h,κ,e) (t) ← V ũ(h,κ,e) (t), (78a) as well as in the weight-adjusted projection of the entropy variables as T , (79a) where we also exploit the fact that M −1 is the identity matrix due to the PKD basis remaining orthonormal under all quadrature rules considered in this work.Furthermore, we employ sum factorization when applying (R (ζ) ) T to obtain r (h,κ,e) (t) in ( 69), and, finally, in the evaluation of the time derivative in ( 56) using the weight-adjusted inverse as Such an approach results in the entire algorithm for computing the local time derivative requiring O(p d+1 ) floating-point operations under the standard assumption that q scales as O(p) with the polynomial degree p of the modal expansion.To the authors' knowledge, this is not achieved by any prior entropy-stable method on triangles or tetrahedra, for which the required dense matrix operations are of O(p 2d ) complexity.Moreover, by avoiding the construction of physical operator matrices and averaging the metric terms on the fly, memory usage is minimized, with per-element memory requirements scaling as O(p d ) due to the only necessary storage being geometric information at each volume and facet quadrature node, and, optionally, precomputed diagonal entries of the matrices W (J (κ) ) −1 , W J (κ) , and B (ζ) J (κ,ζ) .

Comparison to multidimensional summation-by-parts operators
To make the above discussion more quantitative, we now analyze the number of two-point flux evaluations using (77) for specific SBP operators of varying polynomial degrees, where we compare our sparse tensor-product operators on triangles and tetrahedra to those constructed using symmetric quadrature rules as described by Chan [25, Lemma 1], for which the matrices S (l) and (R (ζ) ) T B (ζ)  are dense.We designate the latter class of operator as multidimensional to distinguish them from the tensor-product operators described in Section 3. On the triangle, we construct such multidimensional SBP operators using 2q Xiao-Gimbutas quadrature rules [60] for volume integration and degree 2q + 1 LG quadrature rules for facet integration.On the tetrahedron, the multidimensional SBP operators are constructed using degree 2q Jaśkowiec-Sukumar quadrature rules [61] for volume integration and degree 2q Xiao-Gimbutas triangular quadrature rules for facet integration.Examples of such quadrature nodes are shown in Figure 3 alongside those used for the tensor-product operators described in Section 3. In Figure 4, we plot the number of required two-point flux evaluations given by (77) for each class of operator over a range of polynomial degrees from q = 2 to q = 15, with the exception of the multidimensional operators on tetrahedra, for which suitable symmetric quadrature rules are only currently available, to the authors' knowledge, for SBP operators of degree q ≤ 10.Since dense matrix storage is used in our implementation of the multidimensional operators, the only zero entries considered in (77) for such operators are those along the main diagonal of the skew-symmetric matrix S (l) , although we observe numerically that a small fraction of the off-diagonal entries are, in fact, on the order of machine precision.The scaling in Figure 4 is observed to be slightly better than the asymptotic estimates for both classes of operators, and the number of two-point flux evaluations required for the proposed tensor-product approach is smaller than that required when using multidimensional SBP operators for all polynomial degrees considered.As expected, this discrepancy increases with the polynomial degree; for example, the number of two-point flux evaluations is reduced by factors of 1.56 at q = 2, 2.78 at q = 5, and 4.57 at q = 10 when using tensor-product operators on triangles, and by factors of 1.88 at q = 2, 3.44 at q = 5, and 10.99 at q = 10 when using tensor-product operators on tetrahedra.We also consider the total number of floating-point operations per variable incurred in evaluating the matrix-vector products in ( 69), ( 78), (79), and (80) on a given element, where we take p = q in all cases.The results of such an analysis, which are displayed in Figure 5, are qualitatively similar to those in Figure 4 for higher polynomial degrees, although the benefit of the tensor-product operators for low polynomial degrees is less substantial in this regard, requiring roughly the same number of floating-point operations as the multidimensional operators, for example, at p = 2.
There are several caveats which must be addressed regarding the preceding discussion and analysis.First, we note that due to their use of a larger number of volume and facet quadrature nodes than typical multidimensional operators of the same degree, the tensor-product operators require a somewhat greater number of conversions between conservative and entropy variables as well as a somewhat larger number of numerical interface flux evaluations.These operations do not, however, typically constitute the most significant contribution to the overall expense of an entropy-stable scheme and incur a cost which grows more slowly with the polynomial degree than that of the flux-differencing terms.Second, we recognize that comparisons on the basis of    floating-point operation count, while more objective than implementation-specific and hardwarespecific timing comparisons, are only truly representative of computational cost for compute-bound algorithms.For lower polynomial degrees, algorithm performance is often substantially limited by memory bandwidth, and hence the tensor-product operators may not necessarily outperform their multidimensional counterparts in practice for such regimes.However, due to the arithmetic intensity of an SEM increasing with the polynomial degree of the discretization (see, for example, the roofline analysis in [33]), floating-point operation count indeed becomes a relevant measure at higher polynomial degrees, which is precisely where the proposed tensor-product approach shows significant cost savings.The point at which such benefits become substantial is dependent on the specifics of the implementation and hardware (e.g.considering the memory access pattern, cache size, as well as the use of single-instruction-multiple-data vectorization or multithreading) as well as the PDE and choice of two-point flux, and is therefore an important topic of future investigation within the context of the high-performance implementation and evaluation of the proposed algorithms.

Conservation, free-stream preservation, and entropy stability
The previous section demonstrates that the use of sparse tensor-product operators in collapsed coordinates enables algorithmic improvements relative to comparable entropy-stable discretizations using multidimensional SBP operators, particularly for higher polynomial degrees.We will now demonstrate that the proposed schemes are conservative, free-stream preserving, and entropy stable due to their equivalence to formulations based on hybridized summation-by-parts operators. 6The construction of such hybridized operators from SBP operators satisfying the conditions of Definition 2 is demonstrated with the following lemma.

Lemma 1. Given any nodal SBP operator D
) on the reference element in the sense of Definition 2 for which E (l) takes the form of (13), we can construct a block matrix satisfying the following "SBP-like" property from [25, Eq. ( 32)]: n( 1) B (1)  . . .
Proof.Aside from the fact that our notation separates the operators on individual facets, the result is identical to [63, Theorem 1].The proof relies on the fact that when adding (81) to its transpose, the top-left block vanishes by the skew-symmetry of S (l) and the off-diagonal blocks vanish due to the blocks of the first row being the negative transposes of those in the first column, and hence the diagonal matrix on the right-hand side of ( 82) is all that remains.
The hybridized operators in (81) are of dimension Nq by Nq , where Nq := N q + N (1) . Within a flux-differencing formulation, they operate on block matrices of two-point fluxes given by where couples quadrature nodes on the facets Γ(ζ) , Γ(η) ⊂ ∂ Ω as As in [38], hybridized SBP operators on the physical element are constructed in split form as where the diagonal matrices of concatenated volume and facet metric terms are given by with the blocks G (κ,ζ,l,m) ∈ R The following lemma demonstrates the relation of such a split form to the metric-averaging approach used in Section 4.6.

Lemma 2. The split-form operator in (85) can be rewritten as
where the entries of {{ Ḡ(κ) lm }} ∈ R Nq× Nq are given by {{ Ḡ(κ) Additionally, an SBP-like property analogous to (82) is satisfied on the physical element, as given by provided that J (κ,ζ) and N (κ,ζ,m) are computed using (42) based on the volume metric terms.
Proof.Expressing (85) in indicial form and factorizing, we obtain and hence the result in (87) follows directly from the definition of the Hadamard product.The SBP-like property in (88) follows from ( 82) and ( 85), where we obtain the right-hand side using (42) and the fact that the hybridized boundary operators Ē (l) are diagonal.
Our analysis of the schemes constructed in Section 4 then requires the following assumption.
Assumption 1.The discretization given by (56) with r (h,κ,e) (t) defined as in ( 69) is constructed using a set of diagonal-norm SBP operators on the reference element sharing the same positivedefinite matrix W , with boundary operators in the form of (13).Whether computed exactly or approximately as in (57), the Jacobian determinant of the mapping from the reference element to each physical element satisfies (37).The two-point flux used to compute (67a) and (67b) is consistent and symmetric in the sense of Definition 3, while the directional numerical flux in ( 68) is consistent and conservative in the sense of Definition 4.
Remark 11.The entropy conditions in ( 61) and (63) will not be invoked until the proof of entropy stability.As such, the analysis of conservation and free-stream preservation applies more generally to a wide range of flux-differencing DSEMs, which, as shown in [17], can be constructed so as to recover various split forms in the literature through the particular choice of two-point flux.
We now have the following theorem relating the formulation proposed in Section 4.6 to one which is readily analyzed using the properties of hybridized SBP operators.Theorem 1.Under Assumption 1, the DSEM given by (56) with r (h,κ,e) (t) defined as in (69) is equivalent to the following hybridized SBP formulation proposed in [38,Eq. (35)]: where we define Proof.Substituting ( 69) into ( 56) and grouping all Hadamard products into block matrices gives where the facet correction terms in (70) have been incorporated into the hybridized operator l B (1) R (1)  . . .
Recognizing the matrix in (93) as the skew-symmetric part of Q(l) , we invoke the SBP-like property from Lemma 1 to obtain 2 S (l) = 2 Q(l) − Ē (l) .Substituting such a relation into (92) and moving the inner sum outside the Hadamard product, the flux-differencing terms can be rewritten as where the second term on the right-hand side results from the consistency of the two-point flux and the fact that Ē (l) is diagonal.Finally, we invoke Lemma 2 to express the first term on the right-hand side of (94) using the split-form operator in (85).Substituting the result into (92) then yields the scheme in (90).
As a consequence of the above equivalence, the conservation, free-stream preservation, and entropy stability of the proposed discretizations follow directly from the analysis in [25] and [38].The proofs of such results rely on the assumption of a conforming mesh in the following sense.Assumption 2. For each pair of element indices κ, ν ∈ {1 : N e } with κ ̸ = ν such that ∂Ω (κ) ∩∂Ω (ν) ̸ = ∅, there exist a unique pair of facet indices ζ, η ∈ {1 : N f } such that Γ (κ,ζ) = Γ (ν,η) .Furthermore, for every i ∈ {1 : N (ζ) q f } there exists a unique j ∈ {1 : N (η) q f } for which we have where the normals are computed as in (42) using the (exact or approximate) volume metric terms.
Remark 12. Warburton et al. [64] describe preprocessing algorithms for orienting local coordinate systems on tetrahedra so as to obtain matching node positions in physical space, with a cost scaling linearly with the number of elements.Such procedures ensure that the conditions of Assumption 2 are satisfied despite the asymmetry of the tensor-product facet quadrature nodes arising in the tetrahedral case when using collapsed coordinate systems.
Referring to [25] and [38] for further details, we summarize the critical steps of the analysis in the remainder of this section, beginning with the following theorem establishing discrete conservation.Theorem 2. Let Assumptions 1 and 2 hold and also assume that the Jacobian determinant of the mapping, whether computed exactly or approximately as in (57), satisfies J (κ) ∈ P p ( Ω) and that the quadrature rule in (11) induced by W is of degree τ ≥ 2p.The discretization given by (90), or, equivalently, by (56) with r (h,κ,e) (t) defined as in (69), is then discretely conservative, satisfying Proof.Letting the vector 1 ∈ R N * p contain the modal expansion coefficients for the constant function such that V 1 = 1 (Nq) , we use the skew-symmetry of the matrices S (l) ⊙ {{ Ḡ(κ) lm }} ⊙ F (κ,m,e) in (92) and the exactness of R (ζ) for constant functions to obtain which is a statement of local conservation with respect to the weight-adjusted mass matrix.From [38,Lemma 2], the exactness of the volume quadrature when U (h,κ) e J (κ) ∈ P 2p ( Ω) gives Using the conservation property of the numerical interface flux and invoking Assumption 2 then results in the statement of global conservation in (96).
Next, we present the following lemma proven in [38,Theorem 5], which establishes conditions under which the hybridized SBP operators on the physical element satisfy a discrete form of (38).Lemma 3. Assume that the metric terms, whether computed exactly or approximately, satisfy G (κ) mn ∈ P q ( Ω) as well as the metric identities in (38), and that the normals are computed as in (42).Then, the hybridized SBP operators in (85) satisfy the following discrete metric identities: The discrete metric identities are closely related to the free-stream preservation property, in which a uniform solution state is guaranteed to remain constant in time.This is established for the proposed schemes with the following theorem.Theorem 3. Let Assumptions 1 and 2 hold and also assume that the discrete metric identities in (99) are satisfied for all κ ∈ {1 : κ}.The scheme in (90), or, equivalently, in (56) with r (h,κ,e) (t) defined as in (69), is then free-stream preserving, such that the right-hand side of (56) vanishes for all κ ∈ {1 : N e } and e ∈ {1 : N c } for any uniform solution state satisfying the boundary conditions.
Proof.Considering the formulation in (90), the facet penalty on the second line vanishes when the solution is identical on both sides of the interface due to the consistency property of the numerical interface flux.Invoking the consistency of the two-point flux as well, we then see that the entire right-hand side of (90) vanishes when (99) holds.Since the weight-adjusted mass matrix is invertible by construction, the time derivative is zero, and the scheme is therefore free-stream preserving.
Invoking the entropy conditions in (61) and (63), we now present the following theorem, which establishes that the proposed discretizations satisfy a discrete version of the bound in (6).Theorem 4. Let Assumptions 1 and 2 hold, and also assume that the discrete metric identities in (99) are satisfied for all κ ∈ {1 : κ}, that the two-point flux used to compute (67a) and (67b) is entropy conservative in the sense of Definition 3, and that the numerical interface flux is entropy stable in the sense of Definition 4. The discretization given by (90), or, equivalently, by (56) with r (h,κ,e) (t) defined as in (69), is then discretely entropy stable, satisfying the entropy balance where the entries of the vectors s (h,κ) (t) ∈ R Nq and ψ (κ,ζ,m) (t) ∈ R are given, respectively, by Moreover, the entropy balance in (100) holds as an equality when the numerical interface flux is entropy conservative, and the right-hand side vanishes for periodic boundary conditions.
Proof.The proof of entropy conservation for a non-dissipative interface flux is identical to that of [38,Theorem 2], wherein we left-multiply (90) by (w (h,κ,e) (t)) T , sum over e ∈ {1 : N c }, and use the metric identities in (99) as well as the SBP property in (88) to obtain the local entropy balance For an entropy-conservative interface flux, summing (102) over all elements and splitting the interface contributions between adjacent elements results in a global statement of entropy conservation, corresponding to (100) being satisfied as an equality, where the boundary contributions vanish similarly to the interior interface contributions for periodic problems.The entropy inequality for an entropy-stable interface flux then follows in a straightforward manner, for example, from the analysis in [23,Theorems 3.4 and 4.3].

Numerical experiments
We now present numerical experiments in which we assess the accuracy and robustness of the entropy-conservative and entropy-stable DSEMs using tensor-product as well as multidimensional SBP operators on triangles and tetrahedra through the numerical solution of the Euler equations using StableSpectralElements.jl.The parameters used to run such simulations as well as the Jupyter notebooks used to generate the figures appearing in this section are provided in this paper's reproducibility repository, which is available at https://github.com/tristanmontoya/ReproduceEntropyStableDSEM.

Euler equations
The Euler equations constitute a system of N c = d + 2 coupled nonlinear PDEs governing the conservation of mass, momentum, and energy for a compressible, inviscid, and adiabatic fluid.Such a system takes the form of (1), with the solution variables and flux components given by where ρ(x, t) ∈ R denotes the fluid density, V (x, t) ∈ R d denotes the flow velocity, E(x, t) ∈ R denotes the total energy per unit volume, and P (x, t) ∈ R denotes the pressure.We obtain the pressure using the equation of state for an ideal gas with constant specific heat, where γ > 1 is the specific heat ratio, which we take as 1.4 for air in all numerical experiments.The Euler equations are hyperbolic for solutions belonging to the admissible set Although there exist many entropy-entropy flux pairs which satisfy the conditions of Definition 1 for the Euler equations (see, for example, Harten [65]), we restrict our attention to the choice of which is shown in [20] to be the unique pair (up to an affine transformation) of entropy function and entropy flux which also symmetrize the viscous terms of the compressible Navier-Stokes equations.

Entropy-conservative and entropy-stable flux functions
To obtain an entropy-conservative two-point flux in the sense of Definition 3 with respect to the entropy function and entropy flux in (106), we first define the arithmetic mean and logarithmic mean, respectively, as where we use Taylor-series approximations from Ranocha et al. [58,Algorithms 2 and 3] to compute the logarithmic mean and its reciprocal in cases for which a − and a + are nearly equal.The particular entropy-conservative flux used in this work was proposed by Ranocha [66,67], and is given by In addition to the entropy conservation property in (61), Ranocha's flux is kinetic energy preserving and pressure equilibrium preserving (see, for example, Ranocha and Gassner [68]).The interface flux takes the form of ( 64), where we use Davis's wave speed estimate [69], which is given by In the results which are to follow, we refer to the schemes employing local Lax-Friedrichs dissipation as entropy-stable methods.We also implement a variant without dissipation, in which the second term on the right-hand side of ( 64) is absent; such schemes are denoted as entropy-conservative methods.

Curvilinear meshes
The problems considered in this work are defined on the spatial domain Ω := (0, L) d , where L ∈ R + and d ∈ {2, 3}.The mesh is generated by beginning with a Cartesian grid with M edges in each direction and splitting each quadrilateral into two triangles or each hexahedron into six tetrahedra, resulting in N e = 2M 2 in two dimensions and N e = 6M 3 in three dimensions.Recalling Remark 12, we use the second algorithm described in [30, Section 2.3], which was originally proposed by Warburton et al. [64], to orient the local coordinate systems on each element in order to ensure that Assumption 2 is satisfied for tetrahedral meshes.For both the triangular and tetrahedral case, we use isoparametric mappings, corresponding to (36) with the choice of p g = p = q, where the mapping nodes are obtained using the interpolatory warp-and-blend procedure from Chan and Warburton [70].An affine transformation is used to obtain the positions of the mapping nodes on each element of the split Cartesian mesh.Following [54, Section 5], the mesh is then warped by perturbing the mapping node positions as  the right-hand side of (1b) as We consider the smooth density wave problem from Jiang and Shu [71, Section 7], for which the initial values of the primitive variables are given by on the domain Ω := (0, 2) d , with periodic boundary conditions applied in all directions.The solution is advanced in time for one period of wave propagation (i.e. until a final time of T = 2) using a Julia implementation [72] of the eighth-order Dormand-Prince algorithm described in [73, Section II.5], with the time step taken sufficiently small for the temporal discretization error to be negligible in comparison to that due to the spatial discretization.
Convergence is examined with respect to the nominal element size (taken here to be h := L/M ) as well as the polynomial degree p for entropy-conservative and entropy-stable DSEMs using the tensor-product and multidimensional SBP operators on triangles and tetrahedra considered in Section 5.The degree p of the solution expansion in ( 43) is taken to be equal to the degree q of the SBP operators, and we report the L 2 norm of the density error, which is computed numerically using a quadrature rule of degree 35.Similar convergence behaviour was observed for the other solution variables.Figure 7 demonstrates optimal O(h p+1 ) algebraic convergence under h-refinement for the entropy-stable discretizations (i.e.those including interface dissipation) as well as exponential convergence under p-refinement, where we recall from Section 5 that the multidimensional SBP operators on tetrahedra are only available for degrees up to ten.For a given mesh and polynomial degree, the error norms obtained for the proposed tensor-product discretizations are found to be very close to those obtained for their multidimensional counterparts, which suggests that their algorithmic advantages discussed in Section 5 with respect to the ability to exploit operator sparsity and sum factorization do not come at the expense of accuracy.Remark 13.The numerical results are consistent with the conservation property established in Theorem 2, with the time derivative on the left-hand side of (96) remaining close to machine precision for all solution variables at approximately 100 equispaced snapshots taken during each test.Furthermore, the rates of entropy dissipation given by the left-hand side of (100) are verified for the entropy-conservative and entropy-stable schemes to be zero and non-positive, respectively, at all snapshots (up to roundoff error levels close to machine precision), as is consistent with Theorem 4.

Robustness tests
High-order methods such as DSEMs are prone to numerical stability issues, particularly in the context of under-resolved nonlinear problems, in which high-frequency numerical modes are produced and potentially amplified by the discretization, resulting in non-physical blow-up or negative values of thermodynamic quantities such as pressure or density (i.e.corresponding to solutions outside the admissible set Υ).Such under-resolution commonly arises in simulations of turbulent fluid flow problems, which are characterized by the cascade of energy from larger eddies to those progressively smaller in size, eventually dissipating as heat due to the viscosity of the fluid (see, for example, Pope [74,Chapter 6]).As the Euler equations do not model physical viscosity, the only dissipation present is that inherent in the numerical scheme, which is typically small for a high-order method, thus exposing the potentially destabilizing effects of under-resolved eddies.We therefore use simulations of inviscid vortical flows to mimic worst-case scenarios with respect to under-resolved turbulence, where we are interested in assessing whether the simulations run to completion and whether the semi-discrete entropy bounds established in Section 6 are satisfied, rather than in evaluating the accuracy of such simulations.In the two-dimensional case, we consider the Kelvin-Helmholtz instability (KHI) problem described by Rueda-Ramírez and Gassner [75] and used for robustness tests by Chan et al. [76], for which we define the smoothed step function on the domain Ω := (0, 2) 2 , with periodic boundary conditions in both directions.As in [76], we integrate until a final time of T = 15.In three dimensions, we consider an inviscid Taylor-Green vortex (TGV) problem, for which the initial condition is given on the periodic domain Ω := (0, 2π) 3 as ρ 0 (x) := 1, V 0 (x) := sin(x 1 ) cos(x 2 ) cos(x 3 ), − cos(x 1 ) sin(x 2 ) cos(x 3 ), 0 where Ma ∈ R + is the nominal Mach number.We run the TGV simulations until a final time of T = 14, and, as in Pazner and Persson [77], we consider the nearly incompressible case of Ma = 0.1 in addition to the case of Ma = 0.7, where the latter is expected to pose a greater challenge to the robustness of the proposed methods, specifically with respect to positivity preservation.The Euler equations are solved for the above initial conditions using the proposed entropyconservative and entropy-stable DSEMs for polynomial degrees 4 to 8, taking M = 16 for the KHI  problem and M = 4 for the TGV problem.We integrate in time using the same explicit eighth-order Dormand-Prince method used for the accuracy tests.The time step is taken to be sufficiently small to ensure that the reported instances of instability result purely from the spatial discretization, and we did not find any of the reported instabilities to be remedied by decreasing the time step size.Without any additional stabilization beyond the interface dissipation provided by the numerical flux in (64), all simulations ran to completion for the entropy-stable schemes using multidimensional as well as tensor-product operators on triangles and tetrahedra.Figure 8 demonstrates that the entropy is nonincreasing for all time in each of such cases, as expected from Theorem 4 for periodic boundary conditions.
While the entropy-conservative simulations ran to completion for the TGV with Ma = 0.1, incurring a change in entropy close to machine precision, such computations crashed for the TGV with Ma = 0.7 as well as for the KHI due to negative densities or pressures.This reflects a well-known limitation of the entropy analysis, namely that entropy stability does not guarantee positivity of thermodynamic quantities.Although not provably positivity preserving, the entropy-stable schemes using the dissipative interface flux in (64) did not incur negative densities or pressures for any of the tests considered in this work, likely as a consequence of the dissipative term within the numerical flux serving to dampen any oscillations which would otherwise eventually lead to a violation of positivity.As a result, such schemes are highly robust even in the presence of substantially under-resolved solution features.These results are consistent with the observations in [76], where the authors demonstrated numerically that entropy-stable schemes which incorporate an entropy projection are often able to avoid negative densities or pressures for challenging under-resolved problems without the need for positivity-preserving limiters.We recognize, however, that such an approach may not be sufficient to preserve positivity for problems with discontinuities, and the extension of subcell limiting techniques such as those in [41] and [42] to the proposed tensor-product discretizations on triangles and tetrahedra is an important topic of future research.

Conclusions
We have developed discretizations of arbitrary order for nonlinear hyperbolic systems of conservation laws which combine the geometric flexibility of curved triangular and tetrahedral elements with the efficiency of tensor-product operators as well as the robustness of a provably entropy-stable formulation.The main components of our proposed methodology are outlined below.
• The numerical solution is represented using an orthonormal PKD polynomial basis on the reference triangle or tetrahedron, which supports the efficient matrix-free evaluation of polynomials at tensor-product quadrature nodes in collapsed coordinates through sum factorization.
• The flux-differencing volume terms of the DSEM are computed using sparse tensor-product SBP operators in collapsed coordinates, significantly reducing the number of required entropyconservative two-point flux evaluations relative to discretizations using multidimensional SBP operators based on symmetric quadrature rules, particularly at high polynomial degrees.
• A weight-adjusted approximation is used to invert the dense local mass matrices arising from the use of curved elements, both in the evaluation of the time derivative and in the entropy projection, allowing for the spatial residual to be obtained in O(p d+1 ) floating-point operations, comparing favourably to the O(p 2d ) complexity of existing multidimensional formulations.
• The discrete metric identities, which must hold to obtain an entropy-stable and free-stream preserving scheme, are satisfied by approximating the metric terms in conservative curl form.
This paper provides a comprehensive description of the mathematical and algorithmic aspects of the proposed tensor-product DSEMs on triangles and tetrahedra, discussing their formulation and efficient implementation as well as proving that they are conservative, free-stream preserving, and entropy stable by rewriting them in an equivalent formulation based on hybridized SBP operators.Furthermore, through the solution of the compressible Euler equations on curved triangular and tetrahedral meshes, we numerically verify that the proposed DSEMs are conservative as well as entropy conservative or dissipative (for discretizations with or without interface dissipation, respectively).We also verify the schemes' convergence properties for smooth solutions under h-refinement as well as p-refinement, providing numerical evidence that the proposed discretizations using tensorproduct operators offer very similar accuracy to those using multidimensional SBP operators for such problems.Finally, we demonstrate the robustness of our approach in a suite of challenging under-resolved problems.Future work includes the extension of the proposed methodology to prismatic elements as well as to problems with diffusive terms (e.g. the Navier-Stokes equations) and those in non-conservative form (e.g.multiphase atmospheric flows).We are also interested in the development of positivity-preserving subcell limiters for the proposed schemes as well in practical comparisons of accuracy, efficiency, and robustness relative to other entropy-stable DSEMs as well as to those using other techniques such as over-integration to achieve robustness in practice.

Figure 1 :
Figure 1: Illustration of the collapsed coordinate transformation ξ = χ(η) from the square to the reference triangle

Figure 2 :
Figure 2: Illustration of the collapsed coordinate transformation ξ = χ(η) from the cube to the reference tetrahedron (a) Multidimensional quadrature nodes on the triangle (b) Multidimensional quadrature nodes on the tetrahedron (c) Tensor-product quadrature nodes on the triangle (d) Tensor-product quadrature nodes on the tetrahedron

Figure 3 :
Figure 3: Multidimensional and tensor-product quadrature nodes for SBP operators of degree q = 4

Figure 4 :
Figure 4: Number of entropy-conservative two-point flux evaluations in local flux-differencing terms

Figure 5 :
Figure 5: Total number of floating-point operations per variable in local operator evaluation

Figure 6 :
Figure 6: Examples of warped meshes and mapping nodes for pg = 4

Figure 7 :
Figure 7: Convergence with respect to h and p for discretizations of the Euler equations; dashed and solid lines denote density error for entropy-conservative and entropy-stable schemes, respectively