An Extended Range of Energy Stable Flux Reconstruction Methods on Triangles

We present an extended range of stable flux reconstruction (FR) methods on triangles through the development and application of the summation-by-parts framework in two-dimensions. This extended range of stable schemes is then shown to contain the single parameter schemes of Castonguay et al. (J Sci Comput 51:224–256, 2011) on triangles, and our definition enables wider stability bounds to be developed for those single parameter families. Stable upwinded spectral difference (SD) schemes on triangular elements have previously been found using Fourier analysis. We used our extended range of FR schemes to investigate the linear stability of SD methods on triangles, and it was found that a only first order SD scheme could be recovered within this set of FR methods.


Introduction.
The flux reconstruction (FR) [18] method is a versatile numerical approach to approximating the solution of advection-diffusion equations and can be generalised to arbitrary order accuracy.A compelling advantage of the FR method is the dominance of locally structured computation in the algorithm, making it highly efficient on modern GPU hardware [34].Furthermore, the FR method can be thought of as a generalisation of the discontinuous Galerkin (DG) [1,25] method over the set of different lifting functions.Within the FR literature, these lifting functions are realised through what are called a correction functions.In the seminal work of Huynh [18], it was shown that changes to the correction function could result in significantly different numerical characteristics.
Families of correction functions can be formed by considering different norms in order to prove linear stability.The proofs of linear stability generally look to show that for a system such as: (1. 1) ∂u(x, t) ∂t + ∂f (u) ∂x = 0, for u(x, t = 0) = u 0 (x), x ∈ K ⊂ R, and f = au.
The objective of stability proofs for FR was to find some correction function where the following is true: for some norm A, such as in the work of Vincent et al. [31].It was later shown by Allaneau and Jameson [2] and Zwanenburg and Nadarajah [39] that there is an equivalence between FR and linearly filtered DG.The work of Vincent et al. [33] defined a wider -multi-parameter -set of stable FR schemes which can be interpreted as a wider set of filters, with the filter implicitly defined by A. Concurrently, Ranocha et al. [23] showed how FR could be cast into a summation-by-parts (SBP) formulation to achieve the same results.
In the work of Castonguay et al. [7], and later that of Williams et al. [36], a stable set of FR schemes on triangles was defined that is analogous to those defined in 1D by Vincent et al. [31].For simulations of real flows, triangular elements are advantageous due to the ease of mesh generation for complex geometries, such as via Delaunay triangulation.Concurrently to these works on FR, Balan et al. [3] defined a spectral difference (SD) method on triangles.The FR and SD methods are closely related [19]; however, previous works like that of Veilleux et al. [30] have used Fourier analysis with upwinding to find stable SD schemes.In general, SD schemes on triangles are weakly unstable, and we would like to answer if stable SD methods on triangles can be found as a subset of linearly stable FR methods.
The objective of this work is twofold, firstly we investigate the FR method on triangles further through the use of the SBP framework.We seek to extend the current definition of linearly stable FR methods on triangles from that of Castonguay et al. [7] to a definition analogous to that of Vincent et al. [33].Furthermore, we will study the conditions needed for stability, symmetry, and conservation.Secondly, using this extended definition of stable FR on triangles we investigate the connection between FR and SD, and whether linearly stable SD methods can be found within this set of FR methods.With this in mind, this work is structured as follows.In section 2 we introduce the FR method, previous correction functions, and the SBP framework.The main results of this work are presented in section 3, where linear stability of the FR method on triangles is explored, and conditions are set out for conservation and symmetry.In section 4, we define the new set of stable FR methods for several orders and prescribe there stability conditions.Then in section 5 we use this extended range of stable FR of schemes to investigate SD on triangles.In section 6 we present some brief numerical results from using the newly found FR schemes and finally conclusions are drawn in section 7.

Preliminaries.
2.1.Flux Reconstruction.The flux reconstruction (FR) method of Huynh [18] can be applied to advection and advection-diffusion systems as a method of spatial discretisation.In this work we solely focus on advection systems, and we now give a brief introduction to the FR method in one dimension.For a more detailed description of the method, readers can see references Huynh [18], Witherden and Vincent [38] and references therein.
Beginning with (3.1), the first stage of FR is to subdivide the domain K into N compatible sub-domains, {K} i≤N , such that: A reference sub-domain K = [−1, 1] is defined with the transformation T i : K i → K, such that interpolation and differentiation operators can be more efficiently applied.Here we will only consider affine elements, i.e.where T i is a linear functional.For a given polynomial order, k, k + 1 solution points are placed in the sub-domains at the physical locations {x ij } 0≤j≤k , and in the reference sub-domain at {ξ j } 0≤j≤k .
Polynomials of the solution and flux functions from element K i can be fitted in K using a Lagrange finite-element as To the flux we have added superscript δ to symbolise its correspondence to a discontinuous solution.An approximation to the continuous flux gradient is then formed via the correction procedure as The last and penultimate terms on the right-hand side are used to correct the discontinuous flux to continuous; required for the method to be conservative.Subscripts L and R are used to denote a quantity at the left or right interface and in the case of f δ these values are interpolated.For f num these are interface values that are common to the all the elements that share that interface point.Later for linear equations we will define how this is set, but for alternative equation sets approximate Riemann solver offer a suitable means of setting f num .In (2.3) we introduced the functions h L and h R , these are the correction functions with the boundary conditions h Although it is not necessary, it is typical to set h L and h R as degree k + 1 polynomials.The primary aim of this paper is then how to set these correction functions such that methods are linearly stable.
Lastly, the corrected flux gradient can be coupled to an explicit time integration method, such as SSP-RK3, or coupled to some more complex implicit time-integration system, such as diagonally implicit Runge-Kutta [35].

Correction Functions.
In the earliest paper on the subject of FR, Huynh [18] noted the different numerical properties realised by changing the correction function.In the later works of Vincent et al .[31,32], these correction functions were extended to form continuous classes of functions in one-dimension.With all but one of the correction functions put forward by Huynh [18] being found to be in that class.In the more recent work by Trojak and Witherden [28], a weighted norm was used in the continuous analysis framework of Vincent et al. [31] to produce yet another one parameter family of correction functions with Jacobi orthogonal polynomials as the basis.This idea can be taken to the natural conclusion for any weight function that is positive almost every by using the three-term recurrence relation to generate sets of monic orthogonal polynomials [16].
An alternative approach that was taken by Trojak [29] was to extend the norm used to define stability.Previously, a limited Sobolev norm was used that was sufficient to define the topology of the approximation space, but does not fully capture it.The high order terms at the interfaces that occur in the stability analysis can not be reconciled with the analysis of Vincent et al. [32], it nonetheless showed that it was possible to construct vast sets of correction functions.
In the work of Castonguay et al. [7], an analogous family of correction functions to Vincent et al. [31] was defined on equilateral triangles.To define this, first consider the (m, k)differentiation operator in two dimensions as The Dubiner basis [11] can then be defined as Definition 2.1 (Dubiner Basis).The set is orthogonal over a reference equilateral triangle, where The orders w and v, are then integer solutions to: where ψ is a normalised Jacobi polynomial and ψ (0,0) i = ψ i is a normalised Legendre polynomial.With this definition we can define the set of basis polynomial as We can now define the correction function family of Castonguay et al. [7] in the following lemma Lemma 2.2 (Castonguay et al. correction functions).For a flux point j at the x j and with surface quadrature wight w j , then defining the reconstructed divergence of the correction corresponding to point j as then if the modal coefficients are found from a sufficient condition for linear stability of the resulting FR method is that c > 0.
Remark 2.3.Both the Dubiner basis and the Castonguay correction functions are defined on a triangle with a total order basis [27].This is the most commonly used basis for simplex elements due to the trace space being a polynomial space, and so we will restrict our investigation to these elements.
The work discussed until now used a continuous approach to find correction functions.A different insight may be gained if a discrete approach is used.This was the approach used by Vincent et al. [33] to produce an extended range of stable corrections, which was later encompassed in the work of Trojak [29].This discrete approach was further formalised within the summation-by-parts framework in the work of Ranocha et al. [23].In that work it was also shown how a skew-symmetric split form with a lumped mass matrix could be used to prove stability for Burgers' equation, but not in the general case.First let us define the basic discrete operators that will be used throughout this work.If u i is an approximation in element K i to the function u ∈ C 1 (K), where for our domain we have K ⊂ R d .Then for some set of N s solution points x i = {x i,j } j≤Ns we have the vector u i = u i (x i ), which is the evaluation of u i at the solution points.If l j (x) is the Lagrange polynomial corresponding to l j (x i,k ) = δ jk , then we can define a mass matrix as If we call the cardinal axes x 1 , . . ., x d , then we can define the differentiation operators (2.12) In the following we will drop the subscript i for clarity except where it is explicitly needed.Moving on to define SBP in higher-dimensions, we start with the analogy of integrationby-parts in higher dimensions Definition 2.4 (Divergence integration-by-parts).For a scalar field, v ∈ C 1 (K) and a vector field, W ∈ (C 1 (K)) d , in the closed domain K ⊂ R d with boundary ∂K, then With this we may then define the generalised SBP relation as Definition 2.5 (Generalised Summation-by-parts).
For solutions u ∈ C 1 (K) and U ∈ (C 1 (K)) d , let u i and U i be approximations in element K i such that for some nodal point set x i {x i,j } j≤Ns we have u i = u i (x i ) and U i = U i (x i ), then a set of operators is said to satisfy the generalised SBP property if: , where the divergence and gradient operators are defined as Then defining the interpolation operator L ∂K : K → ∂K, and surface mass matrix, W ∂K , such that where n is a vector of outwards facing normals at the interface.Finally, the Kronecker product of a matrix with the identity matrix is denoted by Remark 2.6.From this definition of SBP we see that the restriction on the mass matrix is that it should accurately integrate all functions in at least Q 2k−1 .From the definition of the mass matrix in ?? this is true, however, in many applications this may not be true if using a quadrature.This case is explicitly handled in section A.
SBP is simply a discrete restatement of integration-by-parts.The advantage is it permits the development of discrete analogues of continuous properties of the physical system.The earliest discussion of SBP in the context of finite element methods, to the authors' knowledge, is that of Fisher and Carpenter [15].This was an adaptation of ideas previously used throughout the finite difference community.There are many works studying SBP in a finite difference context, with some of the earliest works including Carpenter and Otto [5] and Olsson [22].An important work in concreting the utility of these approaches is that of Carpenter et al. [6].There it was shown that for finite differences applied to hyperbolic systems, a scheme with energy bounded via SBP leads to the solution being bounded in the continuous problem.This is important as it shows consistency of the discrete stability analysis and the continuous problem.
In general the exploration of SBP operators has largely focused on one-dimension, but some recent works have move beyond this.For example, on tensor-product elements [14].In the work of Hicken et al. [17] they were able to extend the theory to simplex elements using the generalisation of SBP by Fernández et al. [13].The form given in [13] is analogous to that shown in Definition 2.5.
The operators set out in Definition 2.5 can then be used to construct the FR scheme.First consider a linear advection equation such as where f = au.The FR discretisation of the spatial derivatives can be written within the SBP framework as: where C is the correction function matrix.
In the previous work of Castonguay et al. [7] and Vincent et al. [33], the modal form was used in the presentation of the stable correction functions as the forms are far sparser.Transformation from a nodal to modal representation is defined via where ũ is a vector of modal coefficients, and V is the Vandermonde matrix.Throughout the rest of this work we will use a tilde to denote a matrix of vector in the modal representation.

Linear stability analysis.
To study linear stability, we prescribe the system being solved as a generalised linear advection equation with the form: First we consider the known correction function where FR corresponds to DG, here stability is found in the L 2 norm induced by M. We do this to demonstrate the use of the SBP framework in higher dimensions and to more clearly set out the interface treatment.Lemma 3.1 (Linear Stability).Setting the nodal solution as u i = u(x i ) and linear nodal flux as F = [a x 1 , a x 2 , . . .] T ⊗ u, then for the FR scheme applied to (3.1), the following conditions Proof.The FR method applied to (3.1) can be written as Then multiplying this by u T M we get we can then apply (2.14) to obtain a second equation (3.6) and (3.7) can then be combined to give here we have used the symmetry of M which leads to u T G T MF − u T MDF = 0. Then applying (3.2) we recover Now considering a mesh of multiple elements and focusing on a single point on the boundary of an element, say point j.Using the condition that W ∂ is diagonal, then the global contribution to the right-hand-side of (3.9) from point j is: Here − and + are the contributions from either side of the interface, and from (3.3) w j is the positive surface quadrature weight at j. Then setting the numerical flux from (3.4) we obtain , where we have used n + j = −n − j by definition for a conformal mesh.Applying (3.15a) we recover (3.12) ≤ 0, where we further assume that w − j = w + j i.e. the flux points used here have some degree of rotational symmetry.And hence summing over the domain we recover the required result of As described in section 2, Vincent et al. [33] and later Ranocha et al. [23] were able to derive a multi-parameter extended-range set of FR methods in one dimension that were linearly stable.These methods were found to be stable in a modified norm such that: We now generalise this set of methods to higher dimensions in the following lemma Lemma 3.2 (Extended-range linear stability).For the conditions set out in Lemma 3.1, with the additional constraints that and the modified condition that Proof.Following the same steps as in the proof of Lemma 3.1 and using the modified condition in (3.16) we obtain Then applying (3.15a) and (3.15b), the first term of the right-hand-side is found to be zero, and so proceeding with the proof of Lemma 3.1, we achieve the required results of The final condition (3.15c) is used to ensure that the norm induced by M + Q is a valid norm.
Remark 3.3.A stricter condition on Q is QD = −G T Q; however, when looking for a Q that satisfies this the solution Q = 0 is typically recovered.Alternatively, if the less strict condition (3.15b) is used, a wider range of valid Q matrices can be found that still guarantee linear stability.
Remark 3.4.By finding a norm u M +Q where the solution norm monotonically decays in time, we can use the equivalence of norms to infer stability.Therefore, as c u M ≤ u| M +Q ≤ C u M , the norm u M may not decay monotonically in time, but its rate of decay must remain bounded.
It is often convenient to consider methods in the modal form rather than the nodal form, but to be confident that a scheme found to be stable in modal form is stable in nodal form consider the following: Corollary 3.5 (Nodal-modal equivalence).The stability of a scheme that satisfies conditions (3.3), (3.4), (3.15), and (3.16) is independent of modal or nodal representation, provided V is invertible.
Proof.To prove this it is sufficient to show that the conditions (3.15), if satisfied in one frame, are satisfied in the other.First consider the transform of Q as Next considering the skew symmetry property we have Considering this property we have that: which holds as V is full rank.This completes the proof.

Conservation.
Lemma 3.6 (Linear conservation).Consider the d-dimensional FR method with linear flux function such that, for Banach space V , u ∈ V and F = F (u) ∈ (V , V ) d , then sufficient conditions for conservation are that the gradient operator is are least first order accurate, i.e.
(3.24) G1 = 0, and that the lifting operator is such that ) and F = F (u), then the FR method can be written as Then multiplying by 1 T M we obtain Then applying (2.14) we obtain If (3.24) holds, then we obtain and proceeding to apply (3.25) we get The term on the right-hand side is discrete statement of divergence theorem, which in 1D would give f R − f L .Therefore, the scheme is conservative.
Remark 3.7.A similar lemma can be defined for non-linear flux functions, if an intermediate set of quadrature points is used, see Chan [8].
Here we set Q = 0, but Vincent et al. [33] showed how changing Q could lead some methods to be non-conservative in an integral with unit measure.Conservation of the extended range of stable schemes is then considered in the following lemma Lemma 3.8 (Conservation of extended schemes).For an FR scheme that satisfies (3.3), (3.4), (3.15), and (3.16), with a linear flux function, then if the following condition is also satisfied This can be straightforwardly seen from Lemma 3.6 and (3.16).
Remark 3.9.What can be seen from Lemma 3.6 and Lemma 3.2 is that any FR method satisfying (3.3), (3.4), (3.15), and (3.16) for a linear flux, is automatically conservative in terms of the norm induced by M + Q.However, this is not physical and could lead to schemes that are not consistent.

Symmetry Conditions.
It is taken as axiomatic that the numerical method should be independent of node ordering, or element orientation.As the correction function can change the numerical characteristics of the FR method; therefore, the correction function is required to have some degree of symmetry.
Defining the cardinal axes for different face reference frames, as in Figure 1a, the transformation of reference coordinates can be made via: (3.29) x (θ) = x cos θ + y sin θ, and y (θ) = −x sin θ + y cos θ.
A transformation matrix, Tmn , can then be defined which transforms the basis from the face reference space with θ m to θ n .This allows rotational symmetry conditions to be imposed on Q to give This condition ensures that a function such as φ i (x, y) and the same function rotated to the new reference, φ i (x , y ), then have the same value in the norm induced by M + Q.A further symmetry condition is that, given a pair of flux points on a face that are symmetric about some axis, the corresponding correction functions should be symmetric.This gives the condition that where x is the axis of symmetry and Sx is a matrix that reflects the modes about the axis x.A comparable axial symmetry condition was proposed by Ranocha et al. [23] for use in one-dimension.Finally, when applying symmetry conditions care must be taken to not over constrain the system.This occurs when one symmetry in the reference frame of a face is a linear combination of other symmetry conditions, and is often indicated by the erroneous recovery of Q := 0. For the reference triangle this means that applying two rotational symmetries and one axial symmetry is over constrained as one rational symmetry can be expressed using the other rotation and the axial symmetry.
4. Extended-Range Scheme on Triangle.Vincent et al. [33] and Ranocha et al. [23] derived an extended range of energy stable 1D correction functions and the analysis presented in section 3 derived the conditions required to extend this family to triangles.We now set out this generalised extended range of stable correction functions for several orders using the reference triangle shown in Figure 1b.Furthermore, schemes will be defined in the modal form due to the sparsity of the matrices, and is supported by Corollary 3.5.
In this section we also look to recover the single parameter schemes of Castonguay et al. [7].This set can be cast into the SBP framework with the following definition: Definition 4.1 (Castonguay et al. simplex method).Given the reference triangular element of Figure 1b and a total order basis, FR is found to be stable in the broken Sobolev norm: where A = | K|.This can then be used to define QC required to recover this set of methods in the extended range of stable schemes defined here.Therefore: In what follows we will then look if and how this QC matrix be recovered in the new set of schemes defined, and what the constraints on c are for stability.
As an example of how the correction functions are effected by Q, consider Figure 2 which shows the divergence of the DG correction field and Figure 3 which shows the divergence of the correction field for Q(q 0 = 1, q 1 = 1, q 2 = 1).

Table 1
Stability limits for single parameter Castonguay et al. [7] correction functions on triangles.
5. Spectral Difference Methods.The spectral difference (SD) method [21] is a highorder method similar to FR but with the flux evaluated at staggered set of points, analogous to the method of Kopriva and Kolias [20].The nodal basis of the flux is then chosen such that it lies in the Raviart-Thomas [24] space of the approximation space * .This has the effect of elevating the flux function order, which has been found to give rise to better aliasing properties [9].
In the 1D linear case, SD can be found to be a member of the one parameter class of FR methods [28,31].Using the generalisation of Trojak and Witherden [28], the SD correction functions can be expressed as Jameson [19] has previously shown that in 1D the only linearly stable SD scheme is that corresponding to (α, β) = (0, 0), i.e. the interior flux points are located at the Gauss-Legendre nodes.Trojak and Witherden [28] showed that Fourier analysis can be effectively used to find stable SD schemes with alternative point layouts for which linear stability proofs could not be constructed.A long-standing difficulty has existed in defining linearly stable SD schemes on triangles.Schemes can be constructed for tensor product elements on a maximal order basis, such as quadrilaterals and hexahedrals.However, simplex elements have proven to be more difficult, with some schemes found via Fourier analysis that are stable under interface upwinding.The broad set of stable schemes outlined in section 4 offers a promising route to find generally stable SD method.

One-dimension.
As an initial test of a procedure to find SD correction functions, we look to confirm the SD stability theorem of Jameson [19] in 1D.Here we assume that the interior flux points are placed symmetrically within an element, and when there is an odd number of interior points a single point is placed at the centre.A numerical version of this study has previously been performed by den Abeele et al. [10].* For this reason SD is sometimes referred to as Raviart-Thomas SD.

−z
Using a point layout similar to the example shown in Figure 4, the nodal representation of the 1D correction functions can be written as (5.2) , for m = p/2 and n = p mod 2. This can then be differentiated and transformed into a modal representation allowing for Q to be found via: Here the interpolation operators have been factored out by using the DG correction matrix to give system that is more straightforwardly solved.For k = 3, we have m = 1 and n = 1, and we find Q as However, from (3.15b) we have the additional constraint that q 1 = −5q 2 /3, which can only be satisfied if q 1 = q 2 = 0.This occurs when z 1 = ± 3/5, which when combined with the flux point at x = 0, gives the interior flux points as the nodes of the Gauss-Legendre quadrature.This confirms the result of Jameson [19] and when z 1 is substituted into q 0 we find q 0 = 3/14, as reported by Vincent et al. [33].
Repeating this for k = 4, we find that (5.5) Again using the condition of (3.15b), we find that q i = 0 ∀i ∈ {1, . . ., 5} which is achieved when or vice versa.This again corresponds to the Gauss-Legendre quadrature and gives q 0 = 8/45, corresponding to Vincent et al. [33].
Remark 5.1.The point symmetry imposed and the irrelevance of the ordering of the zeros is why there are multiple solutions that give a valid Q.The correction functions recovered from each is the same.This symmetry can be further identified from the form of the q terms and their lack odd powers of z.

Triangular elements.
Next we extend this procedure to triangles.From the work of Balan et al. [4], we start by defining the Raviart-Thomas (RT) space in two-dimensions as (5.6)In the FR method, the correction functions are within an RT space, and likewise the analogy of corrections in SD are within an RT space.For SD, this two-dimensional basis is then defined via a staggered or flux point set, {σ σ σ}, and requires normals to be associated with each point, n s .An example of these flux points and there normal can be seen in Figure 5.The Lagrange basis can then be defined via the Vandermonde as (5.8) where V RT,x and V RT,x are the Vandermonde matrices over RT k • [1, 0] T and RT k • [0, 1] T respectively.The Lagrange basis is then found from V −1 RT .Finally, the corrections are set using this basis where, from the definition of the SD method, the trace of the SD flux points are located at the FR flux points.5.2.1.k = 1.The most straightforward SD method to define on simplex element is for k = 1.In this case a single interior flux point is required at the element centroid, with normals in x and y.This case was not considered in section 4, but the extended range of stable schemes can be found to be: The SD correction function is then found to be recovered when q 0 = 1/3.5.2.2.k = 2.We now consider k = 2, here the number of flux points is 15, but this can be reduced to 12 degrees of freedom by repeating the interior flux point with orthogonal normals [4].Similar to the method used in one-dimension, we parameterise the interior flux point locations by z 1 , which can be placed using Barycentric coordinates in a manner ensuring rotational symmetry, see Figure 5a.Using this construction, the following matrix can be formed: (5.10) where a Q compatible with subsection 4.1 is sought such that Ã = 0.It was shown by Balan et al. [4] that the stability of the method is independent of the boundary flux point locations, at least for linear equations, and so to reduce the complexity of the resulting matrices we place these points in an equispaced configuration.Focusing on the value of Ã0,0 we find that (5.11) Ã0,0 = 2 − 3 √ 2 18 4  √ 3z 1 (2z 1 − 1) .
The assumption of collocated interior flux points can be relaxed and a second parameter z 2 can be introduced.Repeating the procedure above now with two variables, likewise no pair of variables, (z 1 , z 2 ), can be found for a norm in this class for which the energy monotonically decays in time.For brevity, a full display of the contradictions encountered is not given.5.3.k = 3.For k = 3, assuming collocated interior flux points, there are six degrees of freedom.For symmetric placement of these points there are two possible choices of orbits (0, 2, 0) and (0, 0, 1), based on the work of Witherden and Vincent [38], i.e. two three-point orbits (parameterised by z 1 and z 2 ) or one six-point orbit (also parameterised by z 1 and z 2 ).Examples of these orbits are shown in Figures 5b and 5c.
Starting with the (0, 2, 0) configuration, we again use the result of Balan et al. [3] that stability is independent of the boundary flux point location and use equispaced boundary flux points.Forming the Raviart-Thomas space and finding Ã, we find from the second column of Ã that (5.12) , q 0 = 3/5, and q 2 = 0. and a second solution with z 1 and z 2 swapped.Substituting these into Ã and studying the first column of Ã, we see this leads to the contradiction of (5.13) 45q 1 − 7 = 0 and 15q 1 + 11 = 0.
Therefore, there is no k = 3 SD scheme with the interior flux points in the configuration (0, 2, 0) that is a form of filtered DG.
Repeating this for interior flux points in the orbit (0, 0, 1), we find in the first column of Ã that Clearly this does not satisfy the condition that Ã = 0, from which we can draw the conclusion that there is no stable k = 3 SD scheme, in either (0, 2, 0) or (0, 0, 1), that is a form of filtered DG.
Remark 5.2.This instability is similar to a finding presented by den Abeele et al. [10], however, in that work only Fourier analysis was used to explore stability.Stable schemes where found by Veilleux et al. [30] and Balan et al. [4] using Fourier analysis, however, in that analysis interface upwinding was required to find stable schemes.Therefore, they are not strictly linearly stable.Summarising, we found a linearly stable SD scheme for k = 1, but show that for k = 2 and k = 3 none exist in this set of stable FR methods.It is unlikely that at yet higher orders stable SD schemes will be found in this set of FR methods, and taken with previous results, such as those of Balan et al. [4], den Abeele et al. [10], and Veilleux et al. [30], it is unlikely that linearly stable SD methods on triangles can be found at all without upwind stabilisation.
6. Numerical Experiments.To perform a numerical evaluation of the schemes defined here we considered the Euler vortex case [26], a two-dimensional test case for the Euler equations.A periodic domain Ω = [−10, −10] 2 subdivided into 2(n x − 1) 2 regular right-angles triangles was used, see Figure 6.The system of equations was then (6.1) for pressure P and energy E, and the initial condition was set as: where M is the Mach number, β is the vortex strength, and R is the vortex width, set as 0.4, 13.5, and 1.5 respectively.The error with time can then be calculated for a series of meshes, specifically we used the definition of L 1 and L 2 error of (6.3) where the integrals are approximated with a degree 23 quadrature.In these tests, solution points were positioned at the quadrature points defined by Williams et al. [37].The common interface flux was calculated using a Rusanov flux and Einfeldt wavespeed predictions [12] at flux points located with the Gauss-Legendre quadrature.For time integration, a standard explicit RK4 method was used.Results for k = 3 are presented in Table 2 for Q = 0, Q 1 (q 0 = 0.1, q 1 = 0.1, q 2 = 0.01), and Q 2 (q 0 = 0, q 1 = 0, q 2 = 0.1).Here a constant time step of ∆t = 5 × 10 −3 was used and the L 1 and L 2 error is calculated at t = 100.Table 3 shows the results for the test repeated for k = 4, with Q = 0, Q 1 (q 0 = 0.01, q 1 = 0.01, q 2 = 0.01, q 3 = 0.01), and Q 2 (q 0 = 0.1, q 1 = 0, q 2 = 0, q 3 = 0).At k = 4 a constant time step of ∆t = 2 × 10 −3 was used, again with error calculated at t = 100.These data show that the correction functions tested were stable for t ∈ [0, 100] and furthermore the expected order of accuracy was recovered.The variation in error is evidence of the changing numerical properties caused by varying the correction function.Table 2 Error and order of the Euler vortex for k = 3 FR with DG at t = 100, Q1(q0 = 0.1, q1 = 0.1, q2 = 0.01), and Q2(q0 = 0, q1 = 0, q2 = 0.1).Table 3 Error and order of the Euler vortex with k = 4 FR for DG at t = 100, Q1(q0 = 0.01, q1 = 0.01, q2 = 0.01, q3 = 0.01), and Q2(q0 = 0.1, q1 = 0, q2 = 0, q3 = 0).

Conclusions.
A new multi-parameter set of stable flux reconstruction (FR) methods on triangles was constructed by using the summation-by-parts framework.The correction functions of Castonguay et al. [7] were found to be a subset of this new stable set of FR methods, moreover we were able to successfully expand the stability region of Castonguay et al. [7].Using this new set of FR methods, we investigated if stable SD methods could be defined within it.We found that a stable SD scheme could be produced for k = 1 and that none can be produced in this set of FR methods for k = 2 and k = 3. Numerical experiments were performed for a number of the correction functions outlined in this work and it was shown that the desired order of accuracy was recovered.The approaches outlined here can be used to find similar sets of methods on other element topologies which will be the subject of future work.
Data Availability.The data that support the findings of this study are available from the corresponding author upon reasonable request.

Declaration.
for stability.By definition, we have P T qk M q P qk R qk M −1 q R T qk = I and hence we obtain Finally, the definition of Q q = (V −1 R qk ) T QV −1 R qk follows naturally.This concludes the proof.

2. 3 .
Summation-by-parts.Many advances have been made in the theory of DG and FR methods by considering the discrete problem.A foundation of these analytic techniques is the definition of summation-by-parts (SBP) operators.
) are sufficient for energy stability in the norm induced by the mass matrix, M, i.e.

Figure 2 .
Figure 2. Divergence of the correction field for k = 3 with q0 = q1 = q2 = 0, i.e.DG, for two flux points shown in red.

Figure 3 .
Figure 3. Divergence of the correction field for k = 3 with q0 = q1 = q2 = 1 for two flux points shown in red.

4. 3 .
k = 4. Finally repeating this analysis for k = 4, we find the following definition of the Q matrix.

Figure 5 .
Figure 5. Interior and boundary flux point locations and normals for SD on triangles.

Figure 6 .
Figure 6.Density contour for the Euler vortex with mesh shown for nx = 20.