{'fraction_non_alphanumeric': 0.05363104835426133, 'fraction_numerical': 0.016677252813472092, 'mean_word_length': 4.228371934604905, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A unifying algebraic framework for discontinuous Galerkin and flux reconstruction methods based on the summation-by-parts property

We propose a unifying framework for the matrix-based formulation and analysis of discontinuous Galerkin (DG) and flux reconstruction (FR) methods for conservation laws on general unstructured grids. Within such an algebraic framework, the multidimensional summation-by-parts (SBP) property is used to establish the discrete equivalence of strong and weak formulations, as well as the conservation and energy stability properties of a broad class of DG and FR schemes. Specifically, the analysis enables the extension of the equivalence between the strong and weak forms of the discontinuous Galerkin collocation spectral-element method demonstrated by Kopriva and Gassner (J Sci Comput 44:136-155, 2010) to more general nodal and modal DG formulations, as well as to the Vincent-Castonguay-Jameson-Huynh (VCJH) family of FR methods. Moreover, new algebraic proofs of conservation and energy stability for DG and VCJH schemes with respect to suitable quadrature rules and discrete norms are presented, in which the SBP property serves as a unifying mechanism for establishing such results. Numerical experiments are provided for the two-dimensional linear advection and Euler equations, highlighting the design choices afforded for methods within the proposed framework and corroborating the theoretical analysis.

scientific and engineering disciplines relying heavily on efficient and robust numerical methods for the solution of such partial differential equations. High-order methods (i.e. those of third order or higher) have shown significant promise for problems involving the propagation of waves over long distances (an early example being the work of Kreiss and Oliger [1]) and the fine-scale resolution of turbulent flow structures (as surveyed by Wang et al. [2]), both of which incur prohibitive computational expense when conventional second-order spatial discretizations are used. The most popular high-order methods in current use by practitioners for such problems are arguably the well-established discontinuous Galerkin (DG) methods and the more recent flux reconstruction (FR) schemes. DG methods were originally introduced for the steady neutron transport equation by Reed and Hill [3], and have since been applied successfully to a wide range of problems following extensive development by Cockburn, Shu, and collaborators (see, for example, [4][5][6][7]). The FR approach was first proposed by Huynh [8] for one-dimensional conservation laws, extended to triangular elements by Wang and Gao [9] through the so-called lifting collocation penalty (LCP) formulation, and further developed by several other authors, including Vincent et al. [10], Castonguay et al. [11], and Williams and Jameson [12], who identified energy-stable families of FR schemes on one-dimensional, triangular, and tetrahedral elements, respectively.
Although DG methods are derived from a weak (i.e. variational) formulation whereas FR methods are derived from a strong (i.e. differential) formulation, the two approaches are closely related, as discussed in Huynh's seminal paper [8], as well as in further investigations by Allaneau and Jameson [13], De Grazia et al. [14], Mengaldo et al. [15], and Zwanenburg and Nadarajah [16]. Moreover, certain DG and FR schemes have been recast in terms of summation-by-parts (SBP) operators, which are discrete differential operators equipped with compatible inner products such that integration by parts is mimicked algebraically (see, for example, the review papers of Del Rey Fernández et al. [17] and Svärd and Nordström [18]). Beginning with the work of Kreiss and Scherer [19], this mimetic property, referred to as the SBP property, has been instrumental in the development and analysis of high-order finite-difference methods with discrete stability and conservation properties, and there have been significant advances in recent years towards extending the SBP methodology to encompass a wider range of new and existing methods, facilitated by the one-dimensional and multidimensional generalizations of Del Rey Fernández et al. [20] and Hicken et al. [21], respectively.
Due to its reliance on discretization-agnostic matrix properties, the SBP approach allows for existing theoretical results established for particular numerical methods to be reinterpreted within a more general and arguably simpler setting, and facilitates the cross-pollination of techniques between different schemes sharing common algebraic properties. This cross-pollination has proved fruitful over the past decade, with Gassner [22] and Ranocha et al. [23] employing the SBP property in order to develop nonlinearly stable DG and FR methods, respectively, for Burgers' equation in one space dimension using skew-symmetric split forms originally introduced for finite-difference discretizations. Motivated by these contributions, as well as subsequent extensions to entropy-stable and kinetic-energy-preserving DG-type schemes exploiting the SBP property (see, for example, Gassner et al. [24] and Chan [25]), the goal of this work is to develop a comprehensive framework enabling the unification of existing algebraic techniques for the analysis of semi-discrete DG and FR methods applied to timedependent conservation laws and the extension of the existing theory based on the SBP property to more general DG and FR formulations. To this end, the analysis has resulted in the following novel contributions.

The discrete equivalence of strong and weak formulations, demonstrated by
Kopriva and Gassner [26] in the context of DG methods employing collocated tensor-product Legendre-Gauss (LG) and Legendre-Gauss-Lobatto (LGL) quadrature rules (i.e. the LG and LGL variants discontinuous Galerkin collocation spectral-element method, which we refer to as the DGSEM-LG and DGSEM-LGL, respectively), is extended to more general quadrature-based and collocation-based DG formulations by way of the SBP property. 2. The Vincent-Castonguay-Jameson-Huynh (VCJH) family of FR schemes (as described in [10][11][12]), is recast in a generalized matrix form applicable to general element types, allowing for the use of modal bases and over-integration 1 of the volume and facet terms. Through the SBP property, these formulations are shown to be equivalent to filtered DG schemes in strong or weak form, in the latter case leading to methods which are algebraically equivalent to conventional strong-form FR methods while offering the potential for improved computational efficiency. 3. Proofs of discrete conservation with respect to suitable quadrature rules are presented for DG and VCJH methods in strong and weak form, differing from the existing conservation proofs for the FR approach due to their reliance on the discrete divergence theorem (which holds as a consequence of the SBP property) rather than the continuity of the corrected flux, which is not explicitly reconstructed for multidimensional FR schemes on non-tensor-product elements. 4. Energy stability (with respect to suitable discrete norms) is established as a consequence of the SBP property for DG and FR schemes applied to constantcoefficient linear advection problems on affine polytopal meshes with periodic or inflow-outflow boundary conditions, recovering the existing stability results for VCJH schemes (i.e. those in [10][11][12]) as special cases of a more general theory.
We now outline the content of the remainder of this paper. In §2, we describe the systems of conservation laws to which we seek numerical solutions, as well as the meshes and approximation spaces employed for the discretizations considered in this work; we then provide a brief review of the DG and FR approaches. In §3, we develop a unifying matrix formulation for a broad class of DG and FR methods, which we employ in §4 for the analysis of such schemes based on the SBP property. In §5, we present numerical experiments for the two-dimensional linear advection and Euler equations in support of the theory. Concluding remarks are provided in §6.

Notation
Where possible, the symbols in this paper follow those used by Ranocha et al. [23] and Chan [25]. As in the former, we use a single underline to denote vectors (treated as column matrices), and a double underline to denote matrices. Entries of matrices and vectors are indexed using subscripts, while elements of sequences are typically indexed using superscripts in parentheses. The symbols 0, 1, 0, and I are reserved for vectors of zeros, vectors of ones, matrices of zeros, and identity matrices, respectively, with dimensions inferred from the contexts in which they appear.

Problem formulation
We consider a first-order system of time-dependent conservation laws governing the evolution of the N Eq ∈ N solution variables U (x, t) ∈ Υ ⊆ R NEq , where Υ denotes the set of admissible solution states. Such a system of partial differential equations takes the form For future reference, we also define the outward unit normal vector to the domain boundary as n : denotes the unit (d − 1)-sphere and x 2 := x 2 1 + · · · + x 2 d is the standard Euclidean norm. While the present framework applies to a very broad class of problems which may be formulated as in (2.1), we highlight two representative examples, which serve as test cases for the numerical experiments in §5. We first consider the linear advection equation, which governs the transport of a scalar quantity U (x, t) ∈ R at a velocity a ∈ R d , assumed here to be constant, and is given in conservation form by We also consider the Euler equations, which govern the conservation of mass, momentum, and energy for a compressible, inviscid, adiabatic fluid, given as a system of N Eq = d + 2 equations in the form of (2.1) by In the above, ρ(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, which is related to the other variables through the equation of state, which is given by 2 2 for an ideal gas with constant specific heat and specific heat ratio r. The set of admissible solution states is therefore given by Υ : Remark 2.1 Throughout this paper, we tacitly assume that the values of numerical solutions to (2.1) remain within Υ. This is difficult to ensure a priori for nonlinear problems, and often requires the use of specially designed limiting procedures (see, for example, Zhang and Shu [29]), which are not considered in this work.

Mesh and coordinate transformation
As is standard for finite-element methods (see, for example, Ern and Guermond [30, §1.3]), we consider a mesh T h := {Ω (κ) } K κ=1 consisting of K ∈ N compact elements of characteristic size h ∈ R + with disjoint interiors, satisfying ∪ K κ=1 Ω (κ) =Ω. It is assumed that the boundary of each element is composed of a finite number of smooth facets, 2 which are denoted by Γ (κ,γ) ⊂ ∂Ω (κ) and indexed as γ ∈ {1, . . . , N f }, each with an outward unit normal vector given by n (κ,γ) : Γ (κ,γ) → S d−1 . We further assume that each element is the image of a polytopal reference elementΩ ⊂ R d under a smooth, time-invariant mapping X (κ) :Ω → Ω (κ) , where the Jacobian of the transformation is given by ) > 0, thereby ensuring that the mapping is bijective and orientation preserving. Each facetΓ (γ) ⊂ ∂Ω of the reference element is then assumed to be flat, mapping onto the (potentially curved) physical facet Γ (κ,γ) under X (κ) , with the corresponding outward unit normal vectorn (γ) ∈ S d−1 transforming as where we note that such a transformation is well known in the continuum mechanics literature (see, for example, Gurtin et al. With such definitions and assumptions in place, we recall from Vinokur [32] that systems in the form of (2.1) retain a conservative form when mapped to the reference element, given in the present context by Since periodic boundaries may be treated identically to interior interfaces for systems formulated as in (2.5), we hereinafter use the symbol ∂Ω to refer only to non-periodic portions of the boundary. For the purposes of our analysis, we may then define the sets of element-facet index pairs K := {1, . . . , K} × {1, . . . , N f }, K ∂Ω := {(κ, γ) ∈ K : Γ (κ,γ) ⊂ ∂Ω}, and K Ω := K \ K ∂Ω , with the subsets K ∂Ω and K Ω corresponding to facets lying on the non-periodic portion of the domain boundary, and those forming interior or periodic interfaces, respectively.

Polynomial approximation space
Considering the transformed system in (2.5), we seek a semi-discrete approximation U (h,κ) (x, t) of U (X (κ) (x), t) such that for all κ ∈ {1, . . . , K}, e ∈ {1, . . . , N Eq }, and t ∈ (0, T ), each component U (h,κ) e (·, t) of the numerical solution lies within a finite-dimensional polynomial space of the form 6) where N ⊇ {α 1 ∈ N d 0 : |α| ≤ p} is a finite set of multi-indices 3 with cardinality N p ∈ N. The standard total-degree and tensor-product polynomial spaces P p (Ω) and Q p (Ω) are then recovered with N := {α ∈ N d 0 : |α| ≤ p} and N := {α ∈ N d 0 : max d m=1 α m ≤ p}, respectively, and it is hereinafter assumed that, as is the case for P p (Ω) and Q p (Ω) with any p ∈ N 0 , the space V is closed under partial differentiation with respect to each coordinate and contains the space of constant functions on the reference element, which we denote as P 0 (Ω). The global semi-discrete solution is therefore , . . . , N Eq } and t ∈ (0, T ), where we define the global approximation space with L 2 (Ω) denoting the space of square-integrable functions on the domain Ω.

Remark 2.2 Rather than approximating
This modification of the approximation space has no effect for affine elements, but simplifies the treatment of curvilinear coordinates for time-dependent problems. While the analysis in §4 extends in a straightforward manner to such modified formulations, we note that, as discussed by Yu et al. [33, §3], methods employing such approximation spaces on curved elements may be less accurate than those employing (2.7) as described in this paper.

Discontinuous Galerkin method
Integrating each equation in (2.5) by parts against a test function belonging to the space V, we obtain a local weak formulation of the conservation law, where the resulting integrals may be expressed in terms of the standard L 2 inner products U, After applying a consistent and conservative numerical flux function F * : Υ × Υ × S d−1 → R NEq , typically corresponding to an approximate Riemann solver developed in the context of 3 As is customary, we define the order of the multi-index α ∈ N d 0 as |α| := α 1 + · · · + α d .
finite-volume methods (see, for example, Toro [34]), we obtain the standard weak-form DG approximation on the reference element, which is given by for all κ ∈ {1, . . . , K}, e ∈ {1, . . . , N Eq }, and t ∈ (0, T ). In the above, the volumeweighted contravariant flux components appearing in parentheses in (2.5) are evaluated in terms of the numerical solution as (2.9) and the appropriately scaled numerical flux components are given as where U (+,κ,γ) (x, t) denotes the external data corresponding to a weakly imposed boundary or interface condition.

Discrete operators and the summation-by-parts property
Before presenting our algebraic formulation of the DG and FR methods, we must first define certain operators associated with the approximation space V and their associated matrix representations, the properties of which are discussed in this section.

Discrete inner products
In order to construct approximations of the L 2 inner products U, V Ω and U, V Γ(γ) based on the values of U and V on nodal sets S and S (γ) defined as in §2.6, we relax the assumption of unisolvency to allow for N ≥ N p and define discrete inner products given by in terms of the matrices W ∈ R N ×N and W (γ) ∈ R Nγ ×Nγ , respectively. These bilinear forms approximate each term of the multidimensional integration-by-parts (IBP) relation on the reference element as where we make the following assumption regarding such approximations.
Introducing an arbitrary nodal (i.e. Lagrange) or modal basis B := {φ (i) } Np i=1 for the space V, we may define the generalized Vandermonde matrices V ∈ R N ×Np and V (γ) ∈ R Nγ ×Np with entries given by   The mass matrix with entries consisting of the discrete inner products ·, · W of all pairs of basis functions in B is then characterized by the following lemma. Proof The symmetry and positivity of the corresponding quadratic form are clear from the structure of M and the fact that W is SPD, while the definiteness property follows from the fact that the only vector in the nullspace of V is the zero vector due to it being of full column rank.
the discrete inner product of any two functions U, V ∈ V may be expressed as (U, V W = u T M v, where M is SPD by Lemma 3.1, and since the mapping between any function in V and its expansion coefficients is an isomorphism, the discrete inner product satisfies V , V W > 0 for all nonzero V ∈ V. The discrete inner product is therefore a true inner product for the space V, although definiteness is violated for more general arguments outside of the space V. We also define the mass matrix M (κ) :

Projection operators
Associated with the discrete inner product ·, · W is a corresponding orthogonal projection operator Π onto the space V, which we introduce with the following lemma.
, and corresponds to an interpolation of U on the nodes in S when N = N p , satisfying (ΠU )(x (k) ) = U (x (k) ) for all k ∈ {1, . . . , N }.
Proof As with Lemma 3.3, the result follows from expanding V and ΠU in terms of B, as shown, for example, in Karniadakis and Sherwin [36, Eq. (4.1.33)], and noting that the invertibility of M follows from its positive-definiteness, which was established in Lemma 3.1. If S is unisolvent for V, as is implied by Assumption 3.2 when N = N p , we use the invertibility of V to obtain P = V −1 W −1 V −T V T W = V −1 , and hence evaluating the projection at each node we obtain (ΠU )( While Lemma 3.2 describes an approximation of the L 2 (i.e. Galerkin) projection on the reference element, we may also define a volume-weighted projection corresponding to the approximation of the L 2 projection on the physical element Ω (κ) ∈ T h by replacing P with P (κ) := (M (κ) ) −1 V T W J (κ) , for which we require the following assumption. Since the assumptions of §2.3 imply that the matrix J (κ) is SPD, W J (κ) and M (κ) are also SPD (and therefore invertible) under Assumption 3.2 when W is diagonal, as is the case for the quadrature-based approximations in §B.1. For the collocation-based approach in §B.2, W is not necessarily diagonal, but since M (κ) is equal to W J (κ) , which is the product of two SPD matrices, Assumption 3.3 remains satisfied.

Remark 3.3
If the mapping is affine (i.e. J (κ) is constant) or if N = N p , it follows that P (κ) = P , and the weighting therefore has no influence on the projection.

Summation-by-parts property
As the space V is assumed in §2.4 to be closed under partial differentiation, the operator ∂/∂x m : V → V may be represented algebraically for any m ∈ {1, . . . , d} in terms of the basis B through the matrix D (m) ∈ R Np×Np satisfying . The relation in (3.4) may then be represented algebraically, as established with the following lemma.

Lemma 3.3 Under Assumption 3.1, the matrix D (m) defined implicitly by (3.5) satisfies the SBP property
(3.6) Proof As noted by Chan [25,Eq. (26)], the result follows directly from expanding V and W and their partial derivatives in terms of the basis B and expressing the discrete inner products in matrix form as given in (3.1) and (3.2).

Remark 3.4
As the operators in (3.6) act on expansion coefficients corresponding to a given basis B, which do not necessarily correspond to nodal values, the relation is sometimes referred to as a "modal" SBP property (see, for example, Chen and Shu [37] is itself a nodal SBP operator, recovering a "staggered" SBP formulation as described by Parsani et al. [38] or Del Rey Fernández et al. [39] whenS = S.

Algebraic formulation of the discontinuous Galerkin method
Given a basis B for the space V, the components of the transformed numerical solution for a DG or FR method may be expressed as U In the particular case of a DG method, replacing the L 2 inner products in (2.8) with the corresponding discrete inner products defined in (3.1) and (3.2) results in a coupled system of ordinary differential equations given by for all κ ∈ {1, . . . , K}, e ∈ {1, . . . , N Eq }, and t ∈ (0, T ), where the nodal values of the flux components are given as

Remark 3.5
In order to obtain f (h,κ,m,e) (t) and f ( * ,κ,γ,e) (t), the fluxes in (2.9) and (2.10) may be evaluated pointwise following the pre-multiplication of u (h,κ,e) by V and V (γ) , respectively, in order to obtain the corresponding nodal values of the numerical solution. If B is taken to be a nodal basis defined onS as in §A.2, the pre-multiplication by V may be avoided if collocation is exploited with S =S, while the pre-multiplication by V (γ) may be avoided if S (γ) ⊂S.

Algebraic formulation of the flux reconstruction method
Separating the divergence of the projected flux in (2.11) from that of the correction term, defining the correction field H (γ,j) ∈ V as the divergence of the corresponding correction function G (γ,j) as in Castonguay et al. [11] or Williams and Jameson [12], and replacing the interpolation operator I introduced in §2.6 with the more general projection operator Π described in Lemma 3.2, including within the normal flux difference defined in (2.12), the FR scheme may be expressed as for all κ ∈ {1, . . . , K}, e ∈ {1, . . . , N Eq }, and t ∈ (0, T ). An algebraic formulation of the FR method may then be obtained analogously to that described for the DG method in §3.2 by expressing all functions in V appearing in (3.8) in terms of the basis B, from which we obtain for all κ ∈ {1, . . . , K}, e ∈ {1, . . . , N Eq }, and t ∈ (0, T ), with the initial condition applied as described for the DG method in §3.2. The correction fields appearing in (3.8) are encoded within the lifting matrices L (γ) ∈ R Np×Nγ , which are defined such that These operators allow for more general correction procedures than those in [25,Eq. (22)] and are analogous to the correction matrices introduced in [23] for one-dimensional collocated formulations.

Remark 3.6
The matrix P J (κ) V pre-multiplying the time derivative in (3.9) may be expressed equivalently as M −1 M (κ) and is therefore invertible under Assumption 3.3.

Vincent-Castonguay-Jameson-Huynh correction functions
To complete our algebraic formulation of the FR method in (3.9), the entries of the matrices L (γ) , or, equivalently, the expansion coefficients of the correction fields in terms of the basis B, must be specified. The following theorem demonstrates how such a matrix may be obtained, beginning with the fundamental assumptions associated with the VCJH family of correction functions.

Theorem 3.1 Suppose that there exist correction functions
for all V ∈ V, where C : M → R defines the coefficients for each term in the sum over the finite set of multi-indices M ⊂ N (with N defined in §2.4), and |Ω| is the volume (i.e. Lebesgue measure) of the reference element. Given an arbitrary basis B for V and discrete inner products of sufficient accuracy such that . . , N γ } may then be expanded as in (3.10) with coefficients given by

under the assumption that the coefficients C(α) have been chosen such that M + K is SPD.
Proof Applying IBP to the first term of the integrand in (3.11) and noting that the resulting integrals are L 2 inner products of functions in V and their traces on each facet, we can use the exactness of the discrete inner products and the definition of the correction field H (γ,j) as the divergence of the corresponding correction function (3.12) Using the fact that all terms in (3.12) are linear with respect to the test function V and the fact that the basis B spans the space V, we can test with each basis function and apply the expansion to H (γ,j) ∈ V in order to obtain Recognizing the entries of M and K on the left-hand side and expanding the discrete inner product on the right-hand side, (3.13) may be expressed as (3.14) Considering the above system of equations and using the property G (γ,j) (x (η,i) )·n (η) = δ γη δ ij , we see that the j th column of the matrix L (γ) , containing the expansion coefficients for the corresponding correction field H (γ,j) ∈ V in terms of the basis B, may be obtained by solving the linear system of equations given by ij for all k ∈ {1, . . . , N p }, where M + K is invertible under the assumption of positive-definiteness. As such a system is obtained from (3.11) for all j ∈ {1, . . . , N γ } and γ ∈ {1, . . . , N f }, the VCJH correction fields may be fully specified by defining L (γ) : Remark 3.7 While the conditions on the correction functions in Theorem 3.1 imply the existence of correction fields given by (3.10) in terms of a particular form of L (γ) , the converse of such a statement does not hold in the general case. We therefore cannot necessarily associate a given correction field H (γ,j) obtained from L (γ) as in (3.10) with a unique correction function G (γ,j) satisfying the conditions of Theorem 3.1. In fact, correction functions satisfying such conditions have not, to the authors' knowledge, been explicitly constructed for simplicial elements in two or more space dimensions. 4 The analysis of FR schemes based on properties of the lifting matrix, which we present in §4, is therefore more general than that explicitly relying on properties of the correction functions, and, moreover, ensures that all mathematical objects involved in the proofs are well defined.
To facilitate the analysis in §4, the choices of M and C are constrained by the following assumption.
for all V ∈ V, where we have abused notation to treat the second term on the left-hand side as a constant.

Connection between FR methods and strong-form filtered DG schemes
Considering an FR method for which the correction fields are given as in Theorem 3.1 with K = 0, as is obtained when (3.11) is satisfied with C(α) = 0 for all α ∈ M, we may pre-multiply the resulting scheme in the form of (3.9) by M in order to obtain a strong-form DG method given by for all κ ∈ {1, . . . , K}, e ∈ {1, . . . , N Eq }, and t ∈ (0, T ). In the above, we note that the projection is applied prior to differentiating each transformed flux component as well as prior to evaluating the normal component of the transformed flux on each facet of the reference element.
Remark 3.8 It is also possible to compute the volume terms of a strong-form DG or FR method by differentiating the flux components in (2.9) exactly using the chain rule and projecting afterwards, and to compute the facet terms by evaluating the normal component of the flux on each facet directly in terms of the numerical solution on S (γ) without first projecting. 5 However, the theory in §4 does not apply to such alternative formulations, which, as discussed by Wang and Gao [9, §3.3] and Kopriva and Gassner [26, §4], generally differ from those in (3.9) and (3.16) when applied to nonlinear or variable-coefficient problems and may therefore result in nonconservative strong-form discretizations which are not equivalent to their weak-form counterparts.
For more general choices of K satisfying Assumption 3.4, it was recognized in the one-dimensional case by Allaneau and Jameson [13] that a strong-form linearly filtered DG method is recovered in the sense that the FR residual is obtained through pre-multiplication of the local DG residual resulting from (3.16) by a constant filter matrix F (κ) . Through transformation to a modal basis, it can be shown that such a filter acts only on polynomial modes of the highest degree contained within the space V (see, for example, Allaneau and Jameson [13, §3.2] and Williams and Jameson [12, Appendix B]). The following lemma illustrates such an equivalence in the general multidimensional case with arbitrary choices of basis and projection.

Lemma 3.4 Provided that Assumptions 3.2, 3.3, and 3.4 hold, the FR method in
Proof Pre-multiplying (3.9) by M +K and simplifying the first term by using the formulation may be expressed equivalently as Noting that the third term of (3.17) vanishes when KD (m) = 0 for all m ∈ {1, . . . , d}, as is the case for methods satisfying Assumption 3.4, we therefore recover a filtered strong-form DG method differing from the unfiltered scheme in (3.16) only by the mass matrix.

Remark 3.9
The above result is distinguished from that in Zwanenburg and Nadarajah [16] due to the application of the filter to the entire semi-discrete residual, as opposed to only the facet terms.

Equivalence between strong and weak forms
Just as IBP is essential to obtaining the weak form of a partial differential equation from the strong form, the SBP property in (3.6) may be used analogously to transform a strong-form discretization into a weak-form discretization, and vice versa (as discussed in [17, §8.1]). Although not relying explicitly on the matrix form of the SBP property in (3.6), Kopriva and Gassner proved in [26] that for the DGSEM-LG and DGSEM-LGL on curvilinear hexahedral elements, the strong form in (3.16) is equivalent to the weak form in (3.7). Such an equivalence is extended to more general DG formulations with the following theorem.
Recognizing that the terms containing the normal component of the projected flux at each facet node cancel from the second and third terms of (4.1) and noting that M P = V T W , we therefore recover the weak form in (3.7).
Exploiting the connection between FR methods and filtered DG schemes established in Lemma 3.4, Theorem 4.1 also implies that an algebraically equivalent weak form may be recovered from the FR method in (3.9), as is demonstrated with the following theorem.
Proof From Lemma 3.4, the FR method in (3.9) may be expressed as a strong-form filtered DG scheme given by  are evaluated prior to computing the volume terms and that V (γ) is brought outside of the sum, the use of the weak form leads to an elimination of K · N Eq · N f dense matrix-vector multiplications of size N γ × N p unless B is taken to be a Lagrange basis with solution nodesS containing S (γ) , in which case the interpolation/extrapolation fromS to S (γ) through pre-multiplication by V (γ) is not required. 6 which is algebraically equivalent to the strong formulation in (3.9) under the assumptions of Theorem 4.2, but has a higher computational cost than the weak-form DG method in (3.7) when K = 0. This does not, however, imply that VCJH schemes with K = 0 are inherently more computationally expensive than DG methods, since by instead using the weak formulation in (4.2), we are, in fact, able to attain equivalent computational expense to (3.7) with general choices of K satisfying Assumption 3.4.

Conservation
The proofs of conservation in this section and the proof of energy stability in §4.3 make use of the following assumption, which ensures that the mesh is conforming, with facet nodes aligning in physical space such that the numerical flux may be evaluated pointwise without the need for additional interpolation operators or mortar elements (see, for example, Laughton et al. [40] for an overview and comparison of such techniques).
is satisfied for all i ∈ {1, . . . , N γ }, with an additional translation in the periodic case.
Proof Considering the filtered DG formulation in (4.2), which recovers the standard DG method in (3.7) when K = 0, we first consider a single element and pre-multiply by (P 1) T . Recalling from Lemma 3.4
which sum to zero as a result of the conservation property of the numerical flux. Each term of the second sum in (4.10) then vanishes, resulting in the statement of global conservation in (4.8).
The relation in (4.7) represents a quadrature-based approximation of the integral balance on the physical element Ω (κ) ∈ T h , which may be obtained from the differential form of the conservation law in (2.1) using the divergence theorem. On the reference element, the divergence theorem may be expressed under the discrete inner products for vector fields with components belonging to the approximation space as which is satisfied with Assumption 3.1 and P 0 (Ω) ⊆ V as sufficient but not necessary conditions. 7 We now show that this property allows for the extension of Theorem 4.3 to strong formulations. (3.9) with L (γ) := (M + K) −1 (V (γ) ) T W (γ) and the strong-form DG scheme in (3.16), provided that the discrete divergence theorem in (4.14) is satisfied.

Theorem 4.4 Theorem 4.3 is valid for the strong-form FR scheme in
Proof Expressing (4.14) in matrix form by expanding in terms of the basis B to obtain , we may premultiply the FR method in (3.9) by (P 1) T (M + K), which for K = 0 is equivalent to pre-multiplying the strong-form DG method in (3.16) by (P 1) T , and use the above relation to arrive at (4.9). The remainder of the proof is therefore identical to that of Theorem 4.3, and the result follows.
The invariance of such a quantity therefore provides a useful test for global conservation, which we employ for the numerical experiments in §5.

Remark 4.4
Following the discussion in Remark 3.7, we note that the proofs of conservation in the FR literature (for example, those of Huynh [8] and Castonguay et al. [11, §4.2]) rely explicitly on properties of the correction functions G (γ,j) , notably that G (γ,j) (x (η,i) ) ·n (η) = δ γη δ ij , in order to show that the normal trace of the corrected flux in physical space is continuous at all interfaces. Our proof, in contrast, is instead based only on well-defined algebraic properties of the correction fields encoded in the lifting matrices.

Energy stability
While the theory in §4.1 and §4.2 applies to DG and FR methods for systems in the form of (2.1) on general curvilinear elements, the analysis in this section is restricted to the particular case of the constant-coefficient linear advection equation in (2.2) on meshes consisting of affinely mapped polytopes. We therefore require the following assumption.
Under Assumption 4.2, the matrix P J (κ) V pre-multiplying the time derivative in (3.9) corresponds simply to a scaling by a constant factor of J (κ) . This ensures that the modified mass matrix M (κ) (F (κ) ) −1 corresponding to a VCJH method or the equivalent filtered DG scheme is SPD under Assumption 3.4, allowing for a norm to be defined on the global solution space V h , as shown with the following lemma.

as in (2.7), is a norm under Assumptions 3.4 and 4.2.
Proof As B h := {φ (i) • (X (κ) ) −1 : (κ, i) ∈ {1, . . . , K} × {1, . . . , N p }} is a basis for the space V h , expanding V ∈ V h as above defines an isomorphism between the vector spaces V h and R K·Np . Using such an isomorphism along with the fact that the matrix M (κ) (F (κ) ) −1 = J (κ) (M + K) is SPD for every κ ∈ {1, . . . , K} under the present assumptions allows us to conclude, as outlined in Remark 3.2 for ·, · W , that the bilinear form defined by U, is an inner product for the space V h . As every inner product induces a corresponding norm on the same space, the desired result then follows from expresseing which is the broken Sobolev-type norm introduced by Jameson [41] in the context of one-dimensional spectral-difference methods and adopted in several stability proofs for FR schemes (e.g. those in Vincent et al. [10] and Castonguay et al. [11]). The stability theory for VCJH schemes, which is recovered within the present framework as a consequence of Theorem 3.1, is therefore obtained as a special case of that developed in this section.
In order for advection problems on bounded domains to be well posed in the continuous case, boundary conditions must be prescribed appropriately on portions of ∂Ω such that only incoming waves are specified. In the discrete case, we follow a similar approach and impose boundary and interface conditions weakly in a manner consistent with the physics of the problem through an appropriate choice of numerical flux. This is analogous to the use of simultaneous approximation terms, which are penalty terms added to strong-form discretizations at boundaries and interfaces in order to obtain stable approximations of initial-boundary-value problems with SBP operators, 8 and were introduced in the finite-difference context by Carpenter et al. [42], extending the penalty approach developed for spectral methods by Funaro and Gottlieb [43].
where the parameter λ ∈ R ≥0 takes the same value on either side of an interior interface, and is taken to be unity for any facet lying on ∂Ω.
The numerical flux in (4.16) recovers an upwind flux for λ = 1 and a central flux for λ = 0, corresponding to upwind and symmetric SATs, respectively, in SBP terminology.
We are now equipped to present the following generalized energy stability result for DG and FR schemes in strong and weak form applied to the linear advection equation, in which (as discussed, for example, in Gustafsson (3.7) and the FR method in (3.9) with L (γ) := (M + K) −1 (V (γ) ) T W (γ) applied to the linear advection equation satisfy the discrete energy balance ,1) , t)), . . . , B(X (κ) (x (γ,Nγ ) , t))] T for (κ, γ) ∈ K ∂Ω − , implying that the methods are energy stable in the norm · V h .
Proof Energy stability will be demonstrated for the strong-form FR method in (3.9), which is equivalent to the filtered weak-form DG scheme in (4.2) under the present assumptions as a result of Theorem 4.2, recovering the standard DG method in (3.7) with K = 0. Beginning with a single element and pre-multiplying (3.9) by (u (h,κ) (t)) T (M + K), we may use the definition of the filter matrix and apply the chain rule in time to obtain d dt , (4.18) where the contributions E is satisfied for all m ∈ {1, . . . , d}. Noting the skew-symmetry of the first term on the second line of (4.19), we then obtain (4.20) Since the mapping is assumed to be affine and the flux is assumed to be linear with constant coefficients, the first two terms of the energy balance in (4.18) may be expressed as and respectively, where the second equality in (4.22) results from (2.4). We may therefore apply (4.20) to (4.21) in order to obtain E (κ) (t)/2. Eliminating the factor of 1/2 from both sides, the local energy balance in (4.18) is then given by where we note that the right-hand side contains only facet contributions. As with the proof of global conservation, summing (4.23) over all elements and splitting the net contribution E (κ,γ) III (t)) from each shared interface equally between the index pairs (κ, γ), (ν, η) ∈ K Ω corresponding to the coincident or periodically connected facets Γ (κ,γ) and Γ (ν,η) , the global energy balance is given by As we will now demonstrate, the contributions to the right-hand side of (4.24) arising from interior or periodic interfaces and those arising from facets on the inflow and outflow portions of the boundary may be analyzed separately in order to establish a bound for the quantity on the left-hand side.

Interior or periodic interfaces
Similarly to the proof of Theorem 4.3, it follows from Assumption 4.1 that for any interior or periodic interface corresponding to (κ, γ) ∈ K Ω , the contributions arising from a numerical flux function taking the form of (4.16) may be expressed as E (κ,γ) and E (ν,η) (4.26) Noting that the first terms of (4.25) and (4.26) cancel out the contributions of E net (t), respectively, and that the matrix W (γ) J (κ,γ) is SPD under Assumption 4.2, the net contribution to the energy balance in (4.24) corresponding to any (κ, γ) ∈ K Ω is then given by where for λ = 0, the interior or periodic interfaces do not contribute to the energy balance, whereas for λ > 0, the upwind bias of the numerical flux in (4.16) dissipates energy at a rate depending quadratically on the size of the inter-element jump. We have therefore demonstrated that the interior/periodic interface coupling procedure is energy stable, and now turn our attention to the contributions to (4.24) resulting from the inflow and outflow portions of the boundary.

Outflow boundary
Using the fact that Γ (κ,γ) ⊂ ∂Ω + implies that a · n (κ,γ) = |a · n (κ,γ) |, it follows from taking λ = 1 in (4.25) that E (κ,γ) The resulting contribution to the second sum on the right-hand side of (4.24) for any (κ, γ) ∈ K ∂Ω + is therefore simply −E (κ,γ) II (t), which is non-positive as a result of (4.22) and the positive-definiteness of W (γ) J (κ,γ) , and thus does not contribute to growth in the solution energy. The final energy balance is therefore given by where the second term on the right-hand side corresponds to the rate of energy exiting the domain, while the final two terms correspond to dissipation resulting from the facet jumps. We therefore obtain the inequality in (4.17), bounding the time rate of change in the norm U h (·, t) V h .

Remark 4.6
For homogeneous or periodic boundary conditions (i.e. B(x, t) = 0 for all x ∈ ∂Ω − or K ∂Ω = ∅, respectively), the energy balance in (4.17) implies that d U h (·, t) 2 V h /dt ≤ 0, resulting in energy conservation for periodic problems when a central flux (i.e. λ = 0) is chosen for all interfaces.

Design choices and implementation
The methods described in §3 were implemented within a freely available 9 Python code offering substantial flexibility in the design choices afforded for DG and FR methods, which are detailed in this section.

Discrete inner products
Three options for discrete inner products are employed for the numerical experiments in this section, which are referred to as Quadrature I, Quadrature II, and Collocation schemes. The first two make use of the quadrature-based approach described in §B.1, while the third employs the collocation-based approach in §B.2. For the quadraturebased schemes, all volume terms are computed using positive quadrature rules of degree 2p due to Xiao and Gimbutas [45]. To compute the facet terms, we use LG quadrature rules with p + 1 nodes on each facet (i.e. degree 2p + 1) for the Quadrature I schemes and LGL quadrature rules with p + 1 nodes on each facet (i.e. degree 2p − 1) for the Quadrature II schemes. Although both approaches result in Assumption 3.2 being satisfied (where, as discussed in Remark B.2, the rank of V was verified numerically), the insufficient accuracy of the facet integration for the Quadrature II schemes precludes the use of Lemma B.1. As a result, Assumption 3.1 is satisfied for the Quadrature I schemes, but not for the Quadrature II schemes, leading to the loss of the SBP property in the latter case. The Collocation schemes interpolate the volume terms on the "warp & blend" nodes introduced by Warburton [46], which form unisolvent sets for the corresponding polynomial spaces P p (Ω), and include p + 1 LGL nodes on each facet, which are used to interpolate the facet terms. Therefore, as discussed in §B.2, the Collocation schemes satisfy Assumptions 3.1 and 3.2 by construction. The nodal sets corresponding to each choice of discrete inner product for polynomial degrees 2, 3, and 4 are shown in Figures 1, 2, 3, respectively, where the symbols e and represent nodes in S and S (γ) , respectively.

Basis
For the quadrature-based schemes, orthonormal bases with respect to the L 2 inner product are used, constructed through normalization of the modal basis functions defined in §A.1. For the collocation-based schemes, nodal bases defined as in §A.2 are used withS = S, in this case corresponding to the "warp & blend" nodes in [46].

Spatial residual formulation
Either the strong formulation in (3.9) or the weak formulation in (4.2) may be used, in the case of K = 0 recovering the standard unfiltered strong-form and weak-form DG schemes in (3.16) and (3.7), respectively.

Lifting and filter matrices
The lifting operators employed for the FR schemes take the form L (γ) := (M + K) −1 (V (γ) ) T W (γ) as in Theorem 3.1, with K = K 2D as defined as in §C, and the corresponding filter operator given by F := (I + M −1 K) −1 as in Lemma 3.4. The coefficient c ∈ R scaling the matrix K is chosen either as c = c DG = 0, recovering a DG scheme in strong or weak form, or as c = c + , corresponding to the VCJH scheme found by Castonguay et al. [11] to allow for the largest stable time step for a given polynomial degree when solving the linear advection equation. The latter choice is merely used as a reference value in order to demonstrate that the results in §4 hold for c > 0, and we emphasize that the optimality of such a choice in terms of time step size is highly dependent on the time-marching method, the mesh, and the particular problem being solved. The numerical experiments in Castonguay et al. [11] employ values of c + given by 4.3 × 10 −2 , 6.0 × 10 −4 , and 5.6 × 10 −6 for polynomial degrees 2, 3, and 4, respectively, and we make use of the same values in order to recover discretizations analogous to those described in their work.

Numerical flux
The numerical flux for the linear advection equation is given as in (4.16), where we report results for both λ = 1 and λ = 0, corresponding to upwind and central fluxes, respectively. In the case of the Euler equations, we employ Roe's approximate Riemann solver [47].

Mesh and coordinate transformation
As the test problems considered in this work are posed on the square domain Ω := (0, L) 2 , we begin with a regular Cartesian grid with M ∈ N equal intervals in each direction, dividing each square into two triangles of equal size to obtain a total of K = 2M 2 elements. Using an affine transformation to map the interpolation nodes for a Lagrange basis B map := {φ of degree p map ∈ N (chosen here as the "warp & blend" nodes from [46]) on the reference element onto each physical element ⊂ Ω (κ) , we then interpolate the warping function considered by Del Rey Fernández et al. [48] given by on such a Lagrange basis, resulting in . Examples of such meshes obtained through warping of a split Cartesian mesh with M = 5 are shown in Figure 4, where the symbol • is used to represent the image of each interpolation node for the basis B map under the mapping.

Temporal discretization
The standard explicit fourth-order Runge-Kutta method is used to advance the solution in time, where the time step is chosen as ∆t ≈ C t h/a, 10 with a characteristic element size given by h = L/M and a characteristic wave speed a ∈ R + taken in §5.2 to be the magnitude of the advection velocity and in §5.3 to be the magnitude of the background velocity for the isentropic vortex. Motivated by the Courant-Friedrichs-Lewy (CFL) condition for the DG method in [49, §2.2], we choose C t = β/(2p + 1) with β = 2.5×10 −3 in order to ensure that the error due to the temporal discretization is dominated by the error due to the spatial discretization.

Quantities of interest
In order to provide numerical evidence in support of the theoretical results established in §4, we define the following quantities of interest. Gimbubtas [45], and summing over all κ ∈ {1, . . . , K}. This quantity is expected to be zero for the Quadrature I and Collocation schemes (but not necessarily for the Quadrature II schemes) as a result of Theorems 4.1 and 4.2.

Conservation.
Considering the functional C as defined in Remark 4.3, we compute C(U h e (·, T )) − C(U h e (·, 0)), corresponding to the net production or destruction of each conserved quantity between t = 0 and t = T , which is expected to be zero as a result of Theorems 4.3 and 4.4, the conditions of which are satisfied for all discretizations described in this section. 3. Energy difference. In the case of the linear advection equation, we also evaluate the energy difference U h (·, T ) 2 V h − U h (·, 0) 2 V h , measured in terms of the square of the discrete norm defined in Lemma 4.1, which, as discussed in Remark 4.6, is expected to be nonincreasing with λ = 1 and invariant with λ = 0 for the Quadrature I and Collocation schemes (but not necessarily for the Quadrature II schemes) applied the periodic problem considered here.

Remark 5.1
In the results tabulated in §5.2 and §5.3, values corresponding to properties established theoretically in §4 that are violated (i.e. due to design choices leading to the stated assumptions not being satisfied) are shown in italics. Furthermore, a horizontal line denotes a computation in which the scheme is observed to be unstable.
Remark 5.2 Due to roundoff error and the influence of the temporal discretization (which is not accounted for in the theoretical analysis), we consider a quantity to be zero if it is on the order of 10 −12 or less. This is many orders of magnitude smaller than the solution error for the grids considered in this work, which were deliberately chosen to be coarse in order to unambiguously demonstrate that the properties proven in §4 hold discretely regardless of the grid resolution, rather than asymptotically.

Linear advection equation
We solve the constant-coefficient linear advection equation in (2.2) on a square domain with L = 1 and periodic boundary conditions, with the initial data given by U 0 (x) := sin(2πx 1 /L) sin(2πx 2 /L) and an advection velocity given by a := a[cos(θ), sin(θ)] T with a = √ 2 and θ = π/4, stopping after one advection period (i.e. with T = 1). The mesh is constructed as described in §5.1 with M = 8, where we use p map = 1 such that Assumption 4.2 is satisfied. The results for polynomial degrees 2, 3, and 4 are shown in Tables 1, 2, and 3, respectively. In all cases, the results for the Quadrature I and Collocation schemes are consistent with the theory, satisfying the equivalence, conservation, and stability (i.e. energy conservation for λ = 0 and energy dissipation for λ = 1) properties which were proven in §4. As discussed in §5.1, the Quadrature II schemes violate Assumption 3.1, and thus the theory in §4.1 and §4.3 does not apply; accordingly, the strong-form and weak-form numerical solutions are different, and both the strong and weak forms are found to be unstable when a central flux (i.e. λ = 0) is used. Although the Quadrature II schemes are stable in practice for an upwind flux (i.e. λ = 1), no a priori energy estimate exists due to the violation of the SBP property. The use of an upwind flux in such a case can therefore be interpreted as the addition of numerical dissipation to stabilize an unstable baseline discretization, where the required amount of dissipation is dependent on the particular problem and grid resolution.

Euler equations
The propagation of an isentropic vortex is commonly used as a test case for numerical methods applied to unsteady nonlinear hyperbolic systems. Numerous versions of this problem have been posed in the literature, with the form considered in this work being a modification of that presented by Shu [50, §5.1] following the general formulation described by Spiegel et al. [51]. Considering the Euler equations in (2.3) and denoting the Mach number and direction of the background flow as Ma ∞ ∈ R ≥0 and θ ∈ [0, 2π), respectively, the initial velocity field is given for a vortex of strength ε ∈ R + centred at x 0 ∈ Ω as the initial temperature field is prescribed as and constant entropy is enforced throughout the spatial domain. Using (5.2) and (5.3), as well as the equation of state and isentropic relations for an ideal gas with constant specific heat, the initial condition is then given by , and x 0 = [5,5] T to recover a subsonic vortex which is otherwise similar to the supersonic problem described by Shu [50]. The solution is evolved in time for one period (i.e. with T = √ 2L/Ma ∞ ) on a mesh constructed as described in §5.1 with M = 16 and p map = p.
The results for polynomial degrees 2, 3, and 4 are shown in Tables 4, 5, and 6, respectively, where we see that the strong and weak forms are equivalent for the Quadrature I and Collocation schemes, and the discrete integrals of all four solution variables are invariant for all three approaches, as expected from the theoretical analysis in §4.1 and §4.2. For the Quadrature II schemes, the strong and weak forms again result in different numerical solutions due to Assumption 3.1 not being satisfied. We do not investigate energy balances for this problem, as the theory in §4.3 is not applicable as a result of the system being nonlinear and the mapping not being affine. Hence, the numerical results presented in this section for the Euler equations, as with those presented in §5.2 for the linear advection equation, are consistent with the theory in §4.

Conclusions
We have proposed a unifying framework enabling the algebraic formulation of a broad class of high-order DG and FR methods applied to systems of conservation laws, facilitating the unified theoretical analysis of such schemes through matrix-based techniques. The role of the multidimensional SBP property is examined in detail for methods within the proposed framework, revealing new insights and generalizing existing results pertaining to the equivalence, conservation, and stability properties of DG and FR methods. Potential avenues for future research include the application of similar unifying principles to nonlinearly stable schemes and problems with diffusive terms as well as the comparative evaluation of methods within and outside of the present framework. Specifically, it would be interesting to examine the difference in computational expense between strong-form and weak-form implementations of methods which we have shown in this paper to be mathematically equivalent, and to compare the DG and FR methods considered in this paper with collocated multidimensional SBP schemes without an analytical basis (e.g. those described in [21] and [52]). Additionally, the theory could be extended to enable proofs of convergence, establishing a priori error estimates for linear problems with smooth solutions as consequences of the SBP property, which was discussed by Svärd and Nordström [53] in the context of finite-difference schemes in one space dimension. polynomial basis functions, interpolation nodes, and quadrature rules were generated using modepy, developed by Andreas Klöckner. The initial uniform triangular grids and the LGL quadrature rules were generated using meshzoo and quadpy, respectively, both developed by Nico Schlömer. All plots appearing in this paper were produced using Matplotlib [55].

Code availability
The code developed for this work is publicly available on GitHub under the GNU General Public License at https://github.com/tristanmontoya/GHOST/tree/algebraic_framework_v2.

Availability of data and material
The simulation data analyzed in §5 are available on GitHub within the above repository.

Funding
The first author acknowledges the financial support provided by the University of Toronto, the Natural Sciences and Engineering Research Council of Canada, and the Government of Ontario. SciNet is funded by the Canada Foundation for Innovation, the Government of Ontario, the Ontario Research Fund -Research Excellence, and the University of Toronto.

Conflicts of interest/Competing interests
The authors declare that there are no conflicts of interest or competing interests associated with this work.

A.1 Modal basis
for the space V is orthogonal with respect to the standard L 2 inner product if φ (i) 0 , φ (j) 0 Ω = 0 for all i = j. While such basis functions may be obtained numerically through the Gram-Schmidt orthogonalization of a monomial basis, analytical expressions exist for certain spaces and element types. A two-variable analogue of the Legendre polynomial system was first proposed for triangular elements by Proriol [56] (see also Koornwinder [57,§3.3] and Dubiner [58, §5]), and is given as a basis for Pp(Ω) withΩ := {x ∈ R 2 : where I : N → {1, . . . , Np} reduces the multi-index set N := {α ∈ N d 0 : |α| ≤ p} to a single index, and we define P prismatic, and pyramidal elements, where we refer to the textbook of Karniadakis and Sherwin [36,Ch. 3] for further details regarding multivariate orthogonal polynomials and their use as basis functions for high-order methods.

A.2 Nodal basis
While it is possible to directly employ a modal basis to represent the numerical solution such that B = B 0 , we may alternatively choose B to be a nodal (i.e. Lagrange) basis associated withS := {x (i) } Np i=1 ⊂Ω such that φ (i) (x (j) ) = δ ij , which is satisfied for d = 1 with φ (i) (x) := p+1 j=1,j =i (x −x (j) )/(x (i) −x (j) ) when the nodes inS are distinct. When d ≥ 2, such analytical expressions for Lagrange bases do not exist for arbitrary non-tensorial nodal sets, and we may then take a more general approach by constructing a Vandermonde matrix V ∈ R Np×Np for the modal basis B 0 described in §A.1 such thatṼ ij = φ , we may approximate the L 2 inner products on the reference element and each of its facets using the discrete inner products U, V W := N i=1 ω (i) U (x (i) )V (x (i) ) and U, V W (γ) := Nγ i=1 ω (γ,i) U (x (γ,i) )V (x (γ,i) ), corresponding to (3.1) and (3.2), respectively, with W := diag([ω (1) , . . . , ω (N ) ] T ) and W (γ) := diag([ω (γ,1) , . . . , ω (γ,Nγ ) ] T ). The following lemmas provide sufficient conditions on such quadrature rules for Assumption 3.1, and consequently, the SBP property in (3.6), to be satisfied. Lemma B.1 Assumption 3.1 is satisfied for any reference element type with V = Pp(Ω) when the discrete inner products computed using quadrature rules on S and S (γ) of total degree greater than or equal to 2p − 1 and 2p, respectively.
Proof For any U, V ∈ Pp(Ω), the volume integrals in (3.3) contain functions belonging to the space P 2p−1 (Ω), while the surface integrals contain functions in P 2p (Γ (γ) ). For volume and facet quadrature rules of total degree 2p − 1 or greater and 2p or greater, respectively, all terms in Proof The proof is similar to that in Kopriva and Gassner [26, §3] and involves expanding U, V ∈ Qp(Ω) in terms of the tensor product of B 1D := {φ 1D (x d ), where we recall the definition of N for a tensor-product polynomial space in §2.4. Defining the abscissae and weights of the one-dimensional quadrature rule as S 1D := {x i=1 , respectively, we may factorize the discrete inner product as in [26,Eq. (32)], and the relation in (3.4) therefore follows from the fact that is satisfied for all i, j ∈ {1, . . . , p + 1} when the quadrature is of degree 2p − 1 or higher.
Remark B.1 Unlike Lemma B.1, Lemma B.2 does not presuppose any accuracy conditions on the facet quadrature, and therefore ensures that Assumption 3.1 is satisfied, for example, in the case of tensor-product LGL quadrature rules with p + 1 nodes in each direction (which were considered for the DGSEM-LGL in [26]), for which the corresponding facet quadrature is only of degree 2p − 1.

Remark B.2
The matrices W and W (γ) are SPD if and only if the quadrature weights are strictly positive, which is easily checked. However, Assumption 3.2 also requires V to be of rank Np, which is difficult to prove when d ≥ 2. In practice, the rank of V is verified numerically, and we note that Assumption 3.2 was confirmed to hold for all schemes described in §5.1.

B.2 Collocation-based approximation
Assuming, as in §2.6, that S is unisolvent for the space V, we may take B to be a nodal basis (as , where the integrals may be evaluated in a preprocessing stage, either analytically or using quadrature rules of sufficient degree. Since functions belonging to V and V (γ) are invariant under the interpolation operators I and I (γ) , respectively, the discrete inner products are exactly equal to the corresponding L 2 inner products for such functions. Moreover, the matrices W and W (γ) are generally dense, and are ensured to be SPD due to the linear independence of the functions in B and B (γ) , respectively. Assumptions 3.1 and 3.2 are therefore both satisfied by construction, enabling the formulation of linearly stable discretizations on unisolvent nodal sets for which the quadrature accuracy requirements of Lemmas B.1 and B.2 may not be satisfied.

Remark B.3
If the interpolatory quadrature rule with abscissae S and corresponding weights given by ω (i) = Np j=1 W ij = Ω φ (i) (x) dx is exact for all products of two functions in V, it can be shown that W reduces to a diagonal matrix of quadrature weights, recovering an approximation identical to that described in §B.1. An analogous equivalence holds for W (γ) .

Remark B.4
Collocation-based discrete inner products may also be used with modal bases and over-integrated formulations (i.e. those with N > Np or Nγ > N (γ) p ) by defining I using an auxiliary nodal basisB for the spaceṼ, distinct from the solution basis B, and similarly defining the operators I (γ) using nodal bases for spacesṼ (γ) which may differ from V (γ) . Assumption 3.1 is then satisfied for such schemes whenṼ ⊇ V andṼ (γ) ⊇ V (γ) for all γ ∈ {1, . . . , N f }. In the case ofṼ = V, we may express the orthogonality condition in Lemma 3.2 as V, ΠU − IU Ω = 0 for all V ∈ V, and it therefore follows (see, for example, [58, §5]) that the discrete projection Π corresponds an interpolation I onto the spaceṼ, followed by an exact L 2 projection onto V.