Spectral Analysis of High Order Continuous FEM for Hyperbolic PDEs on Triangular Meshes: Influence of Approximation, Stabilization, and Time-Stepping

In this work we study various continuous finite element discretization for two dimensional hyperbolic partial differential equations, varying the polynomial space (Lagrangian on equispaced, Lagrangian on quadrature points (Cubature) and Bernstein), the stabilization techniques (streamline-upwind Petrov–Galerkin, continuous interior penalty, orthogonal subscale stabilization) and the time discretization (Runge–Kutta (RK), strong stability preserving RK and deferred correction). This is an extension of the one dimensional study by Michel et al. (J Sci Comput 89(2):31, 2021. 10.1007/s10915-021-01632-7), whose results do not hold in multi-dimensional frameworks. The study ranks these schemes based on efficiency (most of them are mass-matrix free), stability and dispersion error, providing the best CFL and stabilization coefficients. The challenges in two-dimensions are related to the Fourier analysis. Here, we perform it on two types of periodic triangular meshes varying the angle of the advection, and we combine all the results for a general stability analysis. Furthermore, we introduce additional high order viscosity to stabilize the discontinuities, in order to show how to use these methods for tests of practical interest. All the theoretical results are thoroughly validated numerically both on linear and non-linear problems, and error-CPU time curves are provided. Our final conclusions suggest that Cubature elements combined with SSPRK and OSS stabilization is the most promising combination.


Introduction
We study several continuous finite element formulations to approximate the solution of the two dimensional hyperbolic conservation laws where ⊂ R 2 is the domain, f : R D → R 2×D is the flux function and u : → R D is the unknown of the system of equations. The largest part of the paper is dedicated to the two-dimensional spectral analysis of different stabilized approaches applied to the scalar (D = 1) transport equations obtained for One of the main objectives of this paper is to identify strategies to build (linearly) stable fully explicit high order continuous finite element schemes to discretize (1) on triangulations of the spatial domain . To this end we will vary the basis functions, the stabilization technique and the time discretization. In general, the standard Finite Element Method (FEM) derived by this approach require the inversion of a large sparse mass matrix. This procedure can be expensive as either inverting the mass matrix and, hence, the matrix multiplications must be repeated for every time step or the linear solver must be applied at each time step. Various techniques have been introduced to overcome the mass matrix inversion while keeping the high order accuracy of the scheme. The first strategy we study is the one proposed in [1]. In the reference it is suggested to combine mass lumping with a deferred correction (DeC) iterative time integration method allowing to introduce appropriate corrections in the right-hand side in order to recover the original order of accuracy. This approach can only be used in combination with finite elements whose basis functions have positive integrals. Another approach is based on a careful choice of approximation points defining sufficiently accurate quadrature formulas with all positive weights. If the variational form is evaluated with this underlying quadrature, as in spectral element methods, we obtain a diagonal mass matrix without loosing the order of accuracy. We refer to this case as cubature elements [40]. For this choice, the classical use of Runge-Kutta methods will provide the high order accuracy also for the time discretization.
Secondly, we will study the influence of the stabilization strategy. When solving (1) with continuous finite elements some additional stabilization operator is necessary to enforce the L 2 stability. Several stabilization techniques can be devised to introduce a level of dissipation comparable to that of discontinuous Galerkin methods with upwind fluxes [46,47]. Three approaches will be studied: the streamline upwind Petrov-Galerkin (SUPG) stabilization [12,18], which is strongly consistent, but it is also introduces new terms in the mass matrix; the continuous interior penalty (CIP) method [14,16,19], consisting in adding edge penalty terms proportional to the jump of the first derivative of the solution; the orthogonal subscale stabilization (OSS) [23], a term that penalizes the L 2 projection of the gradient of the error within the elements. As the CIP stabilization, this technique does not affect the mass matrix, but it requires the solution of another linear system for the L 2 projection. In this respect, the choice of the finite element space and of the quadrature have enormous impact on the cost of the method. Note that the strategy to impose boundary conditions also plays a major role in ensuring stability [4,5], but this will not be considered here.
Our objective is to perform a fully discrete spectral analysis on triangulations of the spatial domain to characterize the stability and accuracy of different combinations of approximation, quadrature, stabilization, and time stepping. In the linear case, this allows to propose optimal values of the CFL and stabilization parameters. Moreover, we analyze a further non-linear high order diffusion operator that can be used to stabilize discontinuities and to provide extra stability to the schemes that show to be unstable with the previous techniques. Numerical simulations for both linear and non-linear scalar problems, and for the shallow water system confirm the theoretical results, and allow to further investigate the impact of the discretization choices on the performance of the schemes and on their cost.
The paper is organized as follows. In Sect. 2 we describe the continuous Galerkin discretization, the stabilization techniques, the basis functions and the time integration techniques. In Sect. 3 we introduce the Fourier analysis space definitions that lead to von Neumann analysis, we discuss some technical details on the passage from physical functions to Fourier modes for different meshes and we find the parameters for which the schemes are stable for some mesh configurations. In Sect. 3.9, we also propose to introduce a viscosity term in order to enforce stability when the previous von Neumann analysis reveals instabilities. In Sect. 4 and Sect. 5 we test the found parameters on some linear and nonlinear problems, checking the order of accuracy and the computational times. Finally, in Sect. 6 we derive some conclusions on the presented schemes and possible applications of the found results.

Numerical Discretization
In this section we describe the discretization of the hyperbolic conservation law (1). We consider a tessellation of the spatial domain consisting of non overlapping (triangular) cells, which we denote by h ⊂ R 2 . The generic element of the tessellation h will be denoted by K , so that h = K . We denote the set of internal element boundaries (edges) of h by F h , using f for a general element. h denotes the characteristic mesh size of h . Despite of the fact that most of the discussion is performed for the scalar case, most of it generalizes readily to systems. If a significant difference arises in this generalization, this will be explicitly discussed.
The discrete solution is sought in a continuous finite element space V p h = {v h ∈ C 0 ( h ) : v h |K ∈ P p (K ), ∀K ∈ h }. We will use nodal and modal finite elements, and we will denote by ϕ j the basis functions associated to the degree of freedom j, so that V p h = span ϕ j j∈ h and we can write where, with an abuse of notation, with j ∈ h we mean the set of degrees of freedom with support in h . With a similar meaning, we will also use the notation j ∈ K to mean the degrees of freedom with support on the cell K . The unstabilized CG approximation of (1) reads: find u h ∈ V p h such that for any where n is the outward facing normal to the boundary of the domain. The choice of W h will be based on V h , but it may take different forms for different stabilizations.
As already mentioned, we will consider several stabilized variants of (3) which can be all formulated in the form: (4) where the flux term is written before the integration by part as we will consider only continuous piecewise polynomials approximations, whose derivatives are integrable. Here, S denotes a bilinear stabilization operator defined on V p h × V p h . Several different choices for S exist, and are discussed in detail in the following sections.

Streamline-Upwind/Petrov-Galerkin: SUPG
The SUPG method was introduced in [31] (see also [12,32] and references therein) and is strongly consistent in the sense that it vanishes when replacing the discrete solution with the exact one. It can be written as a Petrov-Galerkin method replacing v h in (3) with a test function belonging to the space Here, ∇ u f (u h ) ∈ R D×D×2 is the Jacobian of the flux, D the dimensions of the system, τ K denotes a positive definite stabilization parameter with the dimensions of D × D that we will assume to be constant for every element. Although other definitions are possible, here we will evaluate this parameter as where h K is the cell diameter and J K represents the norm of the flux Jacobian on a reference value of the element K . In the scalar case, J K = ||∇ u f (u)|| K . The final stabilized variational formulation of (4) reads The main problem of this stabilization method is that it depends on the time derivative of u and, hence, it does not maintain the structure of the mass matrix in most cases.
To characterize the accuracy of the method, we can use the consistency analysis discussed inter alia in [7, §3.1.1 and §3.2]. In particular, of a finite element polynomial approximation of degree p we can easily show that given a smooth exact solution u e (t, x), replacing formally u h by the projection of u e on the finite element space, we can write with C a constant independent of h, for all functions ψ of class at least C 1 ( ), of which ψ h denotes the finite element projection. A key point in this estimate is the strong consistency of the method allowing to subtract its formal application to the exact solution (thus subtracting zero), and obtaining the above expression featuring differences between the exact solution/flux and its evaluation on the finite element space. Preserving this error estimate precludes the possibility of lumping the mass matrix, and in particular the entries associated to the stabilization term. This makes the scheme relatively inefficient when using standard explicit time stepping. As a final note, for a linear flux (2), exact integration, with τ K = τ and in the time continuous case, a classical result is obtained for homogeneous boundary conditions by testing with v h = u h + τ ∂ t u h [12]: For periodic, or homogeneous boundary conditions, this shows that the norm The interested reader can refer to [12] for the analysis of some (implicit) fully discrete schemes.

Note on the SUPG Technique Applied to Non Scalar Problems
The extension of the SUPG method to a non scalar problem is not straightforward. Here we used the following formulation. First, we define the following system of dimension D: with U ∈ R D , F (U ) ∈ R 2×D and S(U ) ∈ R D . For example, in the results section we will consider the shallow water equations with D = 3 which read where S(U ) is the source term given by a topography term. Equation (10) can also be written in its quasi-linear form where Following the definition of the SUPG method and [52, sec. 5], we define a positive-definite stabilization matrix ø K ∈ R D×D constant for every element K . Here this matrix is evaluated as [52] with S K the set of vertices of K , and n j the outward normal of the edge opposite to the vertex j ∈ S K . h K is the cell diameter and ∇ u F (Ū K ) represents the flux Jacobian of the the average value of U h on the element K . The SUPG stabilized formulation reads, for each equation where (V ) i denotes the i-th component of a vector V ∈ R D .

Continuous Interior Penalty -CIP
Another stabilization technique, which maintains sparsity and symmetry of the Galerkin matrix, is the continuous interior penalty (CIP) method. It was developed by Burman and Hansbo originally in [15] and then in a series of works [14,16,19]. It can also be seen as a variation of the method proposed by Douglas and Dupont [26]. The method stabilizes the Galerkin formulation by adding edge penalty terms proportional to the jump of the gradient of the derivatives of the solution across the cell interfaces. The CIP introduces high order viscosity to the formulation, allowing the solution to tend to the vanishing viscosity limit. This term does not affect the structure of the mass matrix. The method reads where [[·]] denotes the jump of a quantity across a face f , n f is a normal to the face f and where F h is the collection of internal boundaries, and f are its elements. Although other definitions are possible, we evaluate the scaling parameter in the stabilization as where ∇ u f f a reference value of the norm of the flux Jacobian on f and h f a characteristic size of the mesh neighboring f . As stated above, a clear advantage of CIP is that it does not modify the mass matrix, resulting in efficient schemes if a mass lumping strategy can be devised. On the other hand, the stencil of the scheme increases as the jump of a degree of freedom interacts with cells which are not next to the degree of freedom itself (up to 2 cells distance). Note that for higher order approximations [17,38] suggest the use of jumps in higher derivatives to improve the stability of the method. However, here we consider the jump in the first derivatives in order to be able to apply the stability analysis and to study the influence of δ on the stability of the method. We note that the results presented herein might be improved by adding stabilization of higher derivatives.
The accuracy of CIP can be assessed with a consistency analysis as discussed in [7, §3.1.1 and §3.2]. This consists in, formally substituting u h by the projection onto the finite element polynomial of degree p space of u e , a given smooth exact solution u e (t, x), we can show that for all functions ψ of class at least C 1 ( ), of which ψ h denotes the finite element projection, we have the truncation error estimate with C a constant independent of h. The estimate can be derived from standard approximation results applied to u e h − u e and to its derivatives, noting that τ f is an O(h 2 ), leading to the aimed order of accuracy.
The symmetry of the stabilization allows to easily derive an energy stability estimate for the space discretized scheme only. In particular, for periodic boundary conditions and a linear flux we can easily show that which gives a bound in time on the L 2 norm of the solution.
Note that for higher than second order it may be relevant to consider additional penalty terms based on higher derivatives (see e.g. [3,13,17]). We did not do this in this work.

Orthogonal Subscale Stabilization -OSS
Another symmetric stabilization approach is the Orthogonal Subscale Stabilization (OSS) method. Originally introduced as Pressure Gradient Projection (GPS) in [24] for Stokes equations, it was extended to the OSS method in [11,23] for different problems with numerical instabilities, such as convection-diffusion-reaction problems. This stabilization penalizes the fluctuations of the gradient of the solution with a projection of the gradient onto the finite element space. The method applied to (3) reads: find For this method, the stabilization parameter is evaluated as The drawback of this method, with respect to CIP, is the requirement of a matrix inversion to project the gradient of the solution in the second equation of (18). This cost can be alleviated by the choice of elements and quadrature rules if they result in a diagonal mass matrix, as is the case for Cubature elements described below.
As before we can easily characterize the accuracy of this method. The truncation error estimate for a polynomial approximation of degree p reads in this case where the last term is readily estimated using the projection error and the boundness of ψ h as Finally, for a linear flux, periodic boundaries and taking τ K = τ constant along the mesh, we can test with v h = u h in the first equation of (18), and with v h = τ w h in the second one and sum up the result to get which can be integrated in time to obtain a bound on the L 2 norm of the solution.
The truncation consistency error analysis presented above for the three stabilization terms regards only consistency error, but it does not prove stability and convergence for these schemes. These estimations tell us that the stabilization terms that we introduced are of the wanted order of accuracy and that they are usable to aim at the prescribed order of accuracy. This type of analysis has been already done for multidimensional problems inter alia in [2]. More rigorous proof of error bounds with h p+ 1 2 estimates can be found in [13] for the CIP. We did not consider in this work projection stabilizations involving higher derivatives.

Finite Element Spaces and Quadrature Rules
In this section we describe three finite element polynomial approximation strategies used in the paper. In particular, on a triangular element K of h , we define in this section the restriction of the basis functions of V p h on each element K , which are polynomials of degree at most p. We denote by {ϕ 1 , . . . , ϕ N } the basis functions and they will have degree at most p, and their definitions amounts to describe the degrees of freedom, i.e., the dual basis.

Basic Lagrangian Equispaced Elements
On triangles, we consider Lagrange polynomials with degrees at most p: We define the barycentric coordinates λ i (x, y) which are affine functions on R 2 satisfying the following relations where v j = (x j , y j ) are the vertexes of the triangle and, with an abuse of notation, they can be written in barycentric coordinates as v j = (δ 1 j , δ 2 j , δ 3 j ). Using these coordinates, we can define the Lagrangian polynomials on equispaced points on triangles. The equispaced points are defined on the intersection of the lines λ j = k p for k = 0, . . . , p. A way to define the basis functions corresponding to the point (x α , y α ) = (α 1 / p, α 2 / p, α 3 / p) in barycentric coordinates, with α i ∈ {0, . . . , p} and i α i = 1, is in Algorithm 1.

Algorithm 1 Lagrangian basis function in barycentric coordinates
The polynomials so defined in a triangle form a partition of unity, but they have also negative values. This leads to negative or zero values of their integrals. This is problematic for some time discretization and we will see why. We will use these polynomials in combination with exact Gauss-Lobatto quadrature formulae for such polynomials and we will refer to them as Basic elements.

Bernstein Polynomials
Bernstein polynomials are as well a basis of P p but they are not Lagrangian polynomials, hence, there is not a unique correspondence between point values and coefficients of the polynomials. Anyway, there exist a geometrical identification with the Greville points (x α , y α ) = (α 1 / p, α 2 / p, α 3 / p). Given a triplet α ∈ N 3 with α i ∈ 0, . . . , p and i α i = p, the Bernstein polynomials are defined as Bernstein polynomials satisfy additional properties besides the one already cited for Lagrangian points. As before, they form a partition of unity, the basis functions are nonnegative in any point of the triangle, and so their integrals are strictly positive. More precisely These properties lead also to the fact that the value at each point is a convex combination of the coefficients of the polynomials, so that it is easy to bound minimum and maximum of the function by the minimum and maximum of the coefficients. This has been used in different techniques to preserve positivity of the solution [10,37]. We will use these polynomials with corresponding high order accurate quadrature formulae. We will denote these elements with the symbol B p and we refer to them as Bernstein elements.

Cubature Elements
Contrary to the work done in 1D [42], the extension of Legendre-Gauss-Lobatto points which minimize the interpolation error do not exist for the triangle. They have to be computed numerically such as Fekete points [34,55,57]. The problem of this approach is that it requires as classical finite elements the inversion of a sparse global mass matrix. Cubature elements were introduced by G. Cohen and P. Joly in 2001 [25] for the wave equation (second order hyperbolic equation), and are an extension of Lagrange polynomials with the goal of optimizing the underlying quadrature formula error. We will denote the with the symbolP p and they will be contained in another larger space of Lagrange elements, i.e., P p ⊆P p ⊆ P p , with p the smallest possible integer. Similar techniques have been used to minimize the interpolation error [34,55,57]. The objective of these polynomials is to use the points of the Lagrangian interpolation of the polynomials as quadrature points. This means that the obtained quadrature is K f (x, y) = α ω α f (x α , y α ), where K ϕ α = ω α and ϕ α (x β , y β ) = δ αβ . This approach can be considered an extension of the Gauss-Lobatto quadrature in 1D for non Cartesian meshes. The biggest advantage of this approach is to obtain a diagonal mass matrix. The drawback is that one needs to increase the number of basis function inside one element to obtain an accurate enough quadrature rule. In our work, we propose to extend this approach to first order hyperbolic equations. A successful extension to elliptic problem is proposed in [51]. A comparison between the equispace repartition and the Cubature repartition for elements of degree p = 3 is shown in Fig. 1.
For completeness we detail further the construction of the basis functions. The challenges of this approach are the following: -Obtain a quadrature which is highly accurate, at least p + p − 2 order accurate [22]; -Obtain positive quadrature weights ω α > 0 for stability reasons [58]; -Minimize the number of basis functions ofP p ; -The set of quadrature points has to beP p -unisolvent, so that the DoFs coincide with the quadrature points without ambiguity [33]; -The number of quadrature points of edges has to be sufficient ensure the conformity of the finite element.
The optimization procedure that lead to these elements consists of several steps where the different goals are optimized one by one. The optimization strategy exploits heavily the symmetry properties that the quadrature point must have. For p = 1 the Cubature elements do not differ from the Basic elements but in the quadrature formula. For p = 2 the Cubature elements introduce an other degree of freedom at the center of the triangle, leading to 7 quadrature points and basis functions per element. For p = 3 the additional degree of freedom in the triangle are 3, leading to 13 basis functions per triangle. All the details of such elements can be found in [25,33]. We provide in Sect. 1 the detailed expressions of the polynomials used in this work. We will use the symbolP p and the name Cubature elements to refer to them.
Other elements such as Fekete-Gauss points [29,50] exist in the literature. They are optimized to interpolate and integrate with high accuracy. However, it is shown that they require more computing time to achieve similar results than cubature points for high order of accuracy.

Time Integration
The spatial discretization leads to a coupled system of ordinary differential equation which can be written as where U is the vector of all the degrees of freedom on all the domain, M and r are the global mass matrix and right-hand side terms obtained through the discretization of the previous section with some finite elements and stabilization terms. We remark that M is diagonal only in the case of the Cubature elements without the SUPG stabilization, while, for all other choices, it is a sparse non-diagonal matrix.
In the following, we describe two different time integration method: explicit Runge-Kutta (RK) methods and their strong stability preserving (SSP) variants; and the Deferred Correction (DeC) algorithm, which avoids the mass matrix inversion through the correction iterations.

Explicit Runge-Kutta and Strong Stability Preserving Runge-Kutta Schemes
Runge-Kutta time integration methods are one step methods consisting in S stages defined by Here, we use for the solution the superscript n to indicate the timestep and the superscript in brackets (s) to denote the stage of the method. The coefficients α s j and β s j can be defined in many different ways. In particular, we will refer to Heun's method with RK2, to Kutta's method with RK3 and the original Runge-Kutta fourth order method as RK4. The respective Butcher tables can be found in Sect. 2 in Table 12, see [20]. A subset of the RK methods are the SSPRK introduced in [56]. They consist in convex combinations of forward Euler steps, and can be rewritten as follows with γ s j , μ s j ≥ 0 for all j, s = 1, . . . , S. We will consider here the second order 3 stages SSPRK(3,2) presented by Shu and Osher in [56], the third order SSPRK(4,3) presented in [54,Page 189], and the fourth order SSPRK(5,4) defined in [54, Table 3]. For complete reproducibility of the results, we put all their Butcher's tableaux in Sect. 2 in Table 13.

The Deferred Correction Scheme
Deferred Correction methods were introduced in [27] as explicit time integration methods for ODEs, but soon implicit [45], linearly implicit positivity preserving [48] versions and extensions to PDE solvers [1] were studied. In particular, in [1,3,6,8] the DeC is used in a different formulation for finite element methods and it introduces two operator through which it is possible to use a diagonal mass matrix without losing the order of accuracy. This is only achievable when the lumped matrix (defined as the sum on the rows of the full mass matrix) has only positive values on its diagonal. Hence, the use of Bernstein polynomials is recommended in [1], but also Cubature elements can serve the purpose.
Consider a discretization of each timestep into M subtimesteps as in Fig. 2. For each subtimestep we define a high order approximation of the integral form of the ODE (24) from t n,0 to t n,m , i.e., with U = U n,0 , . . . , U n,M . Moreover, the quadrature rule in time uses the subtimesteps t n,m as quadrature points. The corresponding weights ρ m z for every different subinterval are defined by Lagrangian basis functions in these subtimesteps (see [1,3,8] for details). The algebraic system L 2 (U * ) = 0 is in general implicit and nonlinear and, in order not to recast to nonlinear solvers, the DeC procedure approximates the solution of L 2 (U * ) = 0 by successive iterations relying on a low order easy-to-invert operator L 1 . This operator is typically a first order forward Euler approximation with a lumped mass matrix, i.e., Here, D denotes a diagonal matrix obtained from the lumping of M, i.e., D ii := j M i j , and The values of the coefficients β m and ρ m z for equispaced subtimesteps can be found in Sect. 2. Denoting with the superscript (k) index the iteration step, we describe the DeC algorithm as It has been proven [1] that if L 1 is coercive, L 1 − L 2 is Lipschitz with a constant α 1 t > 0 and the solution of L 2 (U * ) = 0 exists and is unique, then, the method converges with an error of O( t K ). Hence, choosing K = M + 1 we obtain a K -th order accurate scheme.
Relying only on the inversion of the low order operator, the method gets rid of the computational costs of the solution of the linear systems, leaving in the right hand side the mass matrix of the L 2 operator, that should not be inverted. The only requirement that is necessary for the DeC approach is the invertibility of the lumped mass matrix D, which limits its application to spatial elements which possess this property. Beyond degree one, basic Lagrange polynomials are not guaranteed to satisfy this property. Hence, only other polynomials as Bernstein and Cubature can be used in combination with DeC.
Finally, for the following analysis we note that the DeC method can be cast in a form similar to a Runge-Kutta method by rewriting (29c) as Comparing with the system of equations (26), we can immediately define the SSPRK coefficients associated to DeC as γ = ρ m r for m, r = 0, . . . , M and k = 0, . . . , K − 1 and instead of the mass matrix, we use the diagonal one.

Remark 1 (DeC with SUPG)
The iterative procedure of the DeC method overcomes the difficulties presented by some implicit stabilization methods such as SUPG. Indeed, the SUPG stabilization term can be added only to the L 2 operator, keeping the high order accuracy of this operator. Since the L 2 operator is applied to the previously computed iteration, all the terms of the SUPG, included the time derivative of u in (7), can be explicitly computed on U (k−1) , keeping then the diagonal mass matrix for the whole scheme.

Preliminaries and Time Continuous Analysis
In order to study the stability and the dispersion properties of the previously presented numerical schemes, we will perform a dispersion analysis on the linear advection problem with periodic boundary conditions: We then introduce the ansatz Here, denotes the damping rate, while the wavenumbers are denoted by k = (k x , k y ), with k x = 2π/L x and k y = 2π/L y with L x and L y the wavelengths in x and y directions respectively. The phase velocity c can be defined from and represents the celerity with which waves propagate in space. It is in general a function of the wavenumber. Substituting (32) in the advection equation (31) for an exact solution we obtain that ω = k · a , c = a and = 0.
In other words The objective of the next sections is to provide the semi-and fully-discrete equivalents of the above relations for the finite element methods introduced earlier. We will consider polynomial degrees up to 3, for all combinations of stabilization methods and time integration techniques. This will also allow to investigate the parametric stability with respect to the time step (through the CFL number) and stabilization parameter δ. In practice, for each choice we will evaluate the accuracy of the discrete approximation of ω and , and we will provide conditions for the non-positivity of the damping , i.e., the von Neumann stability of the method.

The Eigenvalue System
The Fourier analysis for numerical schemes on the periodic domain is based on a discrete Parseval theorem. Thanks to this theorem, we can study the amplification and the dispersion of the basis functions of the Fourier space. The key ingredient of this study is the repetition of the stencil of the scheme from one cell to another one. In particular, using the ansatz (32) we can write local equations coupling degrees of freedom belonging to neighbouring cells through a multiplication by factors e iθ x and e iθ y representing the shift in space along the oscillating solution. The dimensionless coefficient are the discrete reduced wave numbers which naturally appear all along the analysis. Here, x and y are defined by the size of the elementary periodic unit that is highlighted with a red square as an example in Fig. 3.
Formally replacing the ansatz in the scheme we end up with a dense algebraic problem of dimension N dof , where N dof is the number of all the degrees of freedom in the mesh. The obtained system with dimension N dof in the time continuous case reads Equations (31) and (32) with φ j being any finite element basis functions, U the array of all the degrees of freedom and S being the stabilization matrix defined through one of the stabilization techniques of Sect. 2.1. Although system (38) is in general a global eigenvalue problem, we can reduce its complexity by exploiting more explicitly the ansatz (32). The choice of the mesh is crucial in order to exploit the ansatz and to find a unit block that repeats periodically in space. Hence, we must consider structured periodic meshes and we will focus, in particular, on two types of meshes. The first one is the X -mesh that is depicted in Fig. 3 and the second one is the T -mesh depicted in Fig. 4. In those pictures also the distribution of some P 2 elements are represented as an example. More precisely, as it is done in [55] we can introduce elemental vectors of unknowns U Z i j , where Z i j is the stencil denoted by the red square in Fig. 3, which repeats periodically on the domain. So that U Z i j , for continuous finite elements, is an array of d degrees of freedom inside a periodic unitary block Z i j , excluding two boundaries (one on the top and one on the right for example). This number depends on the chosen (periodic) mesh type and on the elements. As an example, in Fig. 3 we display for the X type mesh the periodic elementary unit (in the red square) with Basic and cubature degrees of freedom with p = 2. In the X mesh for Basic elements p = 2 we have d = 8, while for Cubature p = 2 we have d = 12. Using the periodicity of the solution and the ansatz (32) and denoting by Z i±1, j±1 the neighboring elementary units, we can write the neighboring degrees of freedom by and by induction all other degrees of freedom of the mesh. This allows to show that the system (38) is equivalent to a compact system of dimension d (we drop the subscript K as they system is equivalent for all cells) where the matrices M, K x , K y and S are readily obtained from the elemental discretization matrices by using Equations (40). For the discrete Parseval theorem, we know that the norm or the reduced variable U is equivalent to the norm of the discrete vector U. Hence, studying the amplification factor of the two is equivalent.
We apply the same analysis to stabilized methods. The interested reader can access all 2D dispersion plots online [43]. From the plot we can see that the increase in polynomial degree provides the expected large reduction in dispersion error, while retaining a small amount of numerical dissipation, which permits the damping of parasite modes.
An example of dispersion curves is given in Fig. 5. The method used CubatureP 2 elements, the CIP stabilization technique, and a wave angle θ = 5π/4. We here show all 12 parasite modes (see Fig. 3). The principal mode of this system is represented in green. This figure also show the complexity of the analysis because of the number of modes to consider.  We summarize the number of modes for the X mesh in Table 1. A representation of each mesh is done in Sect. 1 for element of degree p = 2 and 3.

The Fully Discrete Analysis
We analyze now the fully discrete schemes obtained using the RK, SSPRK and DeC time marching methods. Let us consider as an example the SSPRK schemes. If we define as  Bern.
Expanding all the stages, we can obtain the following representation of the final stage: where coefficients ν j in (43) are obtained as combination of coefficient γ s j and μ s j in (42) and I is the identity matrix. For example, coefficients of the fourth order of accuracy scheme RK4 are ν 1 = 1, ν 2 = 1/2, ν 3 = 1/6 and ν 4 = 1/24. We can now compress the problem proceeding as in the time continuous case. In particular, using Equations (40) one easily shows that the problem can be written in terms of the local d × d matrices A := M −1 a x K x + a y K y + δ S and in particular that where G ∈ R d×d is the amplification matrix depending on θ, δ, t, x and y. Considering each eigenvalue λ i of G, we can write the following formulae for the corresponding phase ω i and damping coefficient i For the DeC method we can proceed with the same analysis transforming also the other involved matrices into their Fourier equivalent ones. Using (30) these terms would contribute to the construction of G not only in the A matrix, but also in the coefficients ν j , which become matrices as well. At the end we just study the final matrix G and its eigenstructure, whatever process was needed to build it up.
The matrix G describes one timestep evolution of the Fourier modes for all the d different types of degrees of freedom. The damping coefficients i indicate if the modes are increasing or decreasing in amplitude and the phase coefficients ω i describe the phases of such modes.
We remark that a necessary condition for stability of the scheme is that |λ i | ≤ 1 or, equivalently, i ≤ 0 for all the eigenvalues. The goal of our study is to find the largest CFL number for which the stability condition is fulfilled and such that the dispersion error is not too large.
For our analysis, we focus on the X type triangular mesh in Fig. 3 with elements of degree 1, 2 and 3. This X type triangular mesh is also used in [39] for Fourier analysis of the acoustic wave propagation system.

Methodology
The methodology we explain in the following, will be applied to all the combination of schemes we presented above (in time: RK, SSPRK and DeC, discretisation in space: Basic, Cubature and Bernstein, stabilization techniques: CIP, OSS and SUPG), in order to find the best coefficients (CFL, δ), as in [42].
It must be remarked that the dispersion analysis must satisfy the Nyquist stability criterion, i.e., x max ≤ L 2 with x max the maximal distance between two nodes on edges. In other words, k max = 2π L min = 2π 2 x max = π x max . This tells us where k should vary, i.e., k ∈ [0, π/ x max ].
The goal of this section is to minimize the dispersion error and guarantee stability, varying the stabilization parameter and the CFL number. Hence, we look for an algorithm that provides these optimal values. With the notation of [42], we will set for the different stabilizations One of our objectives is to explore the space of parameters (CFL,δ), and to propose criteria allowing to set these parameters to provide the most stable, least dispersive and least expensive methods. A clear and natural criterion is to exclude all parameter values for which there exists at least a wavenumber θ or an angle ∈ [0, 2π] such that we obtain an amplification of the mode, i.e., (θ) > 10 −12 (taking into account the machine precision errors that might occur). Doing so, we obtain what we will denote as stable area in (CFL, θ) space. For all the other points we propose 3 strategies to minimize a combination of dispersion error and computational cost.
In the following we describe the strategy we adopt to find the best parameters couple (CFL,δ) that minimizes a global solution error, denoted by η u , while maximizing the CFL in the stable area. In particular, we start from the relative square error of u Here, we denote with and ω the damping and phase of the principal mode and with ω ex = k · a the exact phase. For a small enough dispersion error |ω − ω ex | 1, we can expand the cosine in the previous formula in a truncated Taylor series as We then compute an error at the final time T = 1, over the whole phase domain, using at least 3 points per wave 0 ≤ k x p ≤ 2π 3 , with x p = x p , and p the degree of the polynomials. We obtain the following L 2 error definition, Recalling that = (k x, CFL, δ, ) and ω = ω(k, x, CFL, δ, ), we need to further set the parameter x p . We choose it to be large x p = 1, with the hope that for finer grids the error will be smaller. Moreover, we need to check that the stability condition holds for all the possible angles ∈ [0, 2π]. Finally, we seek for the couple (CFL * , δ * ) such that where the dependence on of η is highlighted with an abuse of notation. For this strategy, the parameter μ must be chosen in order to balance the requirements on stability and accuracy. After having tried different values, we have set μ to 10 providing a sufficient flexibility to obtain results of practical usefulness. Indeed, the found values will be tested in the numerical section.
To show the influence of the angle on the optimization problem we show an example for the X mesh. For a given couple of parameters (CFL,δ) = (0.4, 0.01) we compare the results for = 0 and = 3π/16. In Fig. 6 we compare the phases ω i and the damping coefficients i for the two angles. It is clear that for the angle = 0, on the left, there are some modes which are not stable i > 0, while for = 3π/16 all modes are stable.
The angle can widely influence the whole analysis as one can observe in the plot of max i i in Fig. 7, where we observe that for the only angle = 3π/16 we would obtain an optimal parameter in (CFL,δ) = (0.4, 0.01), while, using all angles, this value is not stable anymore.

Remark 2
To define the stable region, we should only consider configurations for which the damping is below machine accuracy. In practice, this cannot be done due to the fact that the   (44) is only solved approximately using the linear algebra package of numpy. This introduces some uncertainty in the definition of the stability region as machine accuracy needs to be replaced by some other finite threshold.

Results of the Fourier Analysis Using the X Type Mesh
In this section, we illustrate the result obtained with the methodology explained above. For clarity not all the results are reported in this work, however we place all the plots for all possible combination of schemes in an online repository [42]. We will provide some examples here and a summary of the main results that we obtained.
The first type of plot we introduce helps us in understanding how we can define the stability region in the (CFL, δ) plane. Thus, for every (CFL, δ) we plot the maximum of log( i ) over all modes and angles ∈ [0, 2π] (thanks to the symmetry of the mesh we can reduce this interval). An example is given in the right plot of Fig. 7, it is clear that the whole blue area is stable and the yellow/orange area is unstable. In other cases, this boundary is not so clear and setting a threshold to determine the stable area can be challenging. In Fig. 8 we compare different stabilizations for DeC with B 3 elements. In the CIP stabilization case, we clearly see that there is no clear discontinuity between unstable values and stable ones, as in SUPG, because there is a transient region where max i i varies between 10 −7 and 10 −4 .   The second type of plot combines the chosen stability region with the error η u . We plot on the (CFL, δ) plane some black crosses on the unstable region, where there exists an i and such that i > 10 −7 . The color represents log(η u ) and the best value according to the previously described method is marked with a red dot. In Figs. 9, 10, 11 and 12, we show some examples of these plots for some schemes, for different p = 1, 2, 3. In Figs. 9 and 10 we test the Basic elements with the SSPRK time discretization, while in Figs. 11, 12 we use the Cubature elements with DeC time discretization. We compare also different stabilization technique: in Figs. 9 and 11 we use the OSS, while in Figs. 10 and 12 the CIP. One can observe many differences among the schemes. For instance, for p = 3 we see a much wider stable area for SSPRK than with DeC and, in the Cubature DeC case, we see that the CIP requires a reduction in the CFL number with respect to the OSS stabilization.
We summarize the results obtained by the optimization strategy in Table 2 for all the combinations of spatial, time and stabilization discretization. The CFL and δ presented there are optimal values obtained by the process above described, which we aim to use in simulations to obtain stable and efficient schemes. Unfortunately, as already mentioned above, for some schemes the stability area is not so well defined for several reasons. One of these reasons is the "shape" of the stability area as for one-dimensional problems, see [42]. Other issues that affect this analysis are the numerical precision, see Sect. 3.6, and the mesh configuration, see Sect. 3.7. In the following we study more in details these cases and how one can find better values (Fig. 13).

Comparison with a Space-Time Split Stability Analysis
In this section, we show another stability analysis to slightly improve the results obtained above. Indeed, the solution of the eigenvalue problem (44) is only obtained within some approximation from the numpy numerical library. In some cases, the threshold used to define the stability region is defined in a somewhat heuristic manner. So to confirm the results, we use independently another criterion. To this end we treat independently the temporal and spatial discretizations as in the method of lines. We then study only the spectral properties of the spatial discretization alone, computing the eigenvalues of the corresponding matrix A (cf. (42)). With this information, we then check whether they belong to the stability area of the time discretization.
In particular, following [21], we write the time discretization for Dahlquist's equation in this example, we consider the SSPRK discretization (42). From (43) we can write the amplification coefficient (λ), i.e., "-" means that the fourier analysis shown that the scheme is unstable. * These values are not reliable, see Sect. 3.6 The stability condition for this SSPRK scheme is given by (λ) ≤ 1. Now, when we substitute the Fourier transform of the spatial semidiscretization A to the coefficient λ and we diagonalize the system (or we put it in Jordan's form), we obtain a condition on the eigenvalues of A. Then, studying the Cubature case with SUPG stabilization of order 4 with parameters (CFL,δ)=(0.234, 0.011), found in Fig. 13, see also Table 2, we plot the eigenvalues of A and the stability region of the SSPRK scheme for different θ ∈ [0, π]. We notice that for some values of θ some of the eigenvalues fall slightly outside the stable area, see Fig. 14a. There are, indeed, few eigenvalues dangerously close to the imaginary axis and some of them have actually positive real part (blue dots). As suggested before, if we decrease the CFL and increase δ, we move towards a safer region, so considering (CFL,δ)=(0.18, 0.04) with the same θ , we obtain all stable eigenvalues, as shown in Fig. 14b.  Fig. 13 Logarithm of the amplification coefficient log(max i (ε i )) for SUPG stabilization withP 3 Cubature elements and the SSPRK method. Unstable region in yellow, the red dot is the optimal parameter according to (50) Fig. 14 Eigenvalues of A using cubature discretization and the SUPG stabilization (varying k) and stability area of the SSPRK method. In red the stable eigenvalues, in blue the unstable ones. (color figure online) The summary of the optimal parameters of Table 2 updated taking into account also a larger safety region in the (CFL, δ) plane (as explained in this section) can be found in Table 15 in Appendix 2.

Different Mesh Patterns
Another important aspect about this stability analysis is the influence of the mesh structure on the results. As an example, we use the T -mesh, another regular and structured mesh type depicted in Fig. 4. In Fig.4 we plot also the degrees of freedom for elements of degree 2 and the periodic elementary unit that we take into consideration for the Fourier analysis. The number of modes in the periodic unit for this mesh type are summarized in Table 3. The elements of degree 3 can be found in Fig. 28 in Appendix 1. Cub. 1 6 13 Basic. 1 4 9 Bern. Even if for several methods we observe comparable results for the two mesh types, for some of them the analyses are quite different. An example is given by the Basic elements with SSPRK schemes and CIP stabilization. For this method, we plot the dispersion error (49) and the stability area in Fig. 15a for the X mesh and in Fig. 15b for the T mesh. We see huge differences in P 2 and P 3 where in the former a wide region becomes unstable for δ L ≤ δ ≤ δ R and for the latter we have to decrease a lot the value of δ to obtain stable schemes.
In the case of Cubature elements with the OSS stabilization and SSPRK time integration, we have already seen in the previous section that the optimal parameters found were in a dangerous area. Repeating the stability analysis for the T mesh we see that the situation is even more complicated. In Fig. 16a we plot the analysis for the X mesh and in Fig. 16b the one for the T mesh.P 3 elements, though being stable for some parameters for the X mesh, are never stable on the T mesh. This means, that, when searching general parameters for the schemes, we have to keep in mind that different meshes leads to different results.
For completeness, we present the optimal parameters also for the T mesh in Table 16 in Appendix 2.
In general, it is important to consider more mesh types when doing this analysis. In practice, we will use the two presented above (X and T meshes). In the following, we will consider the stability region as the intersection of stability regions of both meshes.

Final Results of the Stability Analysis
Taking into consideration all the aspects seen in the previous sections, it is important to have a comprehensive result, which tells which parameters can be used in the majority of the situations. A summary of the parameters obtained for the X and T mesh is available in Appendix 2. In Table 4, instead, we present parameters obtained using the most restrictive case among different meshes and that insure a sufficiently large area of stability around them, The values denoted by * are not the optimal one, but they lay in a safer region, see Sect. 3.6. The values marked by * * cannot be used on the T mesh. "-" means that it is unstable for every parameter as explained in Sect. 3.6. These parameters can be safely used in many cases and we will validate them in the numerical sections, where, first, we validate the results of the X mesh on a linear problem on an X mesh, then we used the more general parameters in Table 4 for nonlinear problems on unstructured meshes. A special remark must be done for CubatureP 3 elements combined with the OSS and the CIP stabilizations. In Fig. 17 we see how the amplification coefficient max i ε i has always values far away from zero. For the CIP stabilization this is always true and even for theP 2 elements the stability region is very thin. As suggested in [17,38] higher order derivatives jump stabilization terms might fix this problem, but it introduces more parameters. This has not been considered here. Another remark is that the T configuration is very peculiar and, as we will see, on classical Delauney triangulations the issue seem to not affect the results. Moreover, the use of additional discontinuity capturing operators may alleviate this issue as some additional, albeit small, dissipation is explicitly introduced in smooth regions.
In Sect. 3.9, we propose to add an additional stabilization term for these unstable schemes, i.e., CubatureP 3 elements and OSS or CIP stabilization techniques. This term is based on viscous term [2,30,36,41] and allows to stabilize numerical schemes for any mesh configuration.
For the OSS stabilization we observe a similar behavior in Fig. 17. The stability that we see in that plot are only due to the the T mesh. Indeed, for the OSS stabilization on the X mesh there exists a corridor of stable values, which turn out to be unstable for the T mesh, see Fig. 18. In practice, also on unstructured grids we have not noticed instabilities when running with the parameters found with the X mesh. Hence, we suggest anyway some values of CFL and δ for these schemes, which are valid for the X mesh, noting that they might be dangerous for very simple structured meshes. The validation on unstructured meshes also for more complicated problems will be done in the next sections.
Overall, Table 4 gives some insight on the efficiency of the schemes. We remind that, in general, we prefer matrix free schemes, so this aspect must be kept in mind while evaluating the efficiency of the schemes. All the SUPG schemes, except when with DeC, and all the Basic element schemes have a mass matrix that must be inverted. Among the others we see that for first degree polynomials schemes the DeC with Bernstein polynomials and SUPG stabilization gives one of the largest CFL result, while for second degree polynomials the OSS Cubature SSPRK scheme seems the one with best performance and, for fourth order schemes, again the Bernstein DeC SUPG is one of the best.
In conclusion of this section, there are important points to highlight: -The extension of the Fourier analysis to the two-dimensional space leads to significantly different results with respect to the one-dimensional one. Both in terms of global stability of the schemes, and in terms of optimal parameters. Moreover, in opposition to [42], Bernstein elements with SUPG stabilization technique lead to stable and efficient schemes. Cubature elements, which were the most efficient in one-dimensional problems, have stability issues on the two-dimensional mesh topologies studied. -The complexity of the analysis in two-dimensional space is increased. This not only implies a larger number of degrees of freedom, but also more parameters to keep into account, including the angle of the advection term and the possible different configuration of the mesh. The visualization of the stability region of the time scheme as shown in Fig. 14 with the eigenvalues of the semi-discretization operators helps in understanding the effect of CFL and penalty coefficient on the stability of the scheme, only for methods of lines. This helps in choosing and optimizing the couple of parameters.

Remark 3
Another possibility to characterize the linear stability of numerical method is proposed by J. Miller [44]. This method is based on the study of the characteristic polynomial of the amplification matrix G. However, this method does not provide information about the phase ω, since it does not compute eigenvalues of G. For this reason, we choose the eigenanalysis.

Accounting for Discontinuity Capturing Corrections
The stabilization terms accounted for so far are linear stabilization operators. For more challenging simulations, additional non-linear stabilization techniques might be added to control the numerical solution in vicinity of strong non-linear fronts and/or discontinuities. We consider here the effect of adding an extra viscosity term, as in the entropy stabilization formulations proposed e.g. in [2,30,35,36,41]. We in particular look at the approach proposed in [30], and used for shallow water waves in [41,49] and in [9,28]. In this approach the viscosity is designed to provide a first order correction μ K = O(h) close to discontinuities, while for smooth enough solutions μ K = ch p+1 . Our idea is to embed this high order correction explicitly in the analysis of the previous section to provide a heuristic characterization of the fully discrete stability of the resulting stabilized formulation: find u h ∈ V p h that satisfies for any Viscosity term = 0. (53)

Note on the Stability of the Method
As it is done for previous stabilization terms in Sect. 2.1, we can characterize the accuracy of this method estimating the truncation error for a polynomial approximation of degree p.
Considering the smooth exact solution u e (t, x) of (53), for all functions ψ of class at least C 1 ( ) of which ψ h denotes the finite element projection, we obtain with C a constant independent of h. The estimate can be derived from standard approximation results applied to u e h − u e and to its derivatives, knowing that μ K = O(h p+1 ). Then, for a linear flux, periodic boundaries and taking μ K = μ constant along the mesh, we can test with v h = u h in (53), we get which can be integrated in time to obtain a bound on the L 2 norm of the solution.

The von Neumann Analysis
As we saw in Sect. 3.8, the T mesh configuration has stability issues. In particular, the numerical schemes using CubatureP 3 elements, SSPRK and DeC time integration methods, and the OSS and the CIP stabilization techniques are unstable. We propose to evaluate these schemes adding the viscosity term in (53). For the von Neumann analysis, we use μ K (u) = ch p+1 K in (53), with c ∈ R + , h K the cell diameter and p the degree of polynomial approximation. We show the plot of max i i to understand how the stability region behaves with respect to c using CubatureP 3 elements. In Fig. 19 the maximum amplification factor is represented for varying c, using the OSS stabilization technique and the SSPRK time integration method. We note that the same behaviour is observed with CIP and DeC. Plots are available online [43].
We can observe two main results. First, increasing the parameter c up to around 0.1 allows to expand the stability region. Second, when the viscosity coefficients reaches too high values, it is necessary to decrease the CFL (see Fig. 19c with μ = 0.05 and Fig. 19d with μ = 0.5 as an example).

Numerical Verification
We now perform numerical tests to check the validity of our theoretical findings. We initially focus on the structured grids, and in particular on the X mesh configuration, although similar verifications have been performed on the T mesh. We will use elements of degree p, with p up to 3, with time integration schemes of the corresponding order of accuracy to ensure an overall error of O( x p+1 ), under the CFL conditions discussed earlier (see also Table  15 in Sect. 2). As already stressed, numerical integration is performed with Gauss-Legendre formulae of the appropriate order to exactly integrate the variational form for Basic and Bernstein elements, while for Cubature elements we use those associated to the interpolation points.
The mesh used in the Fourier analysis is the basis of the one we will use in the numerical simulations. We will extend it periodically for the whole domain, see an example in Fig. 20a.
The exact solution is u ex (x, t) = u 0 (x − a x t, y − a y t) for all x = (x, y) ∈ and t ∈ R + . The initial conditions are displayed in Fig. 20b. We discretize the domain with the X mesh pattern, see Fig. 20a. To have approximately the same number of degrees of freedom for different degrees p, we use different mesh sizes for each order of accuracy: x 1 = {0.1, 0.05, 0.025} for P 1 , x 2 = 2 x 1 for P 2 , and x 3 = 3 x 1 for P 3 elements. In Fig. 21a, b, we study the error convergence for different schemes. In the x-axis the values of t are displayed, which we remind are proportional to x, and the error is plotted on the y-axis. These figures show a comparison between Cubature and Basic elements with OSS stabilization and SSPRK time integration. As we can see, the two schemes have correct   Table 15, and have larger computational costs because of the inversion of the mass matrix.
To show the main benefit of using the Cubature elements (diagonal mass matrix), we plot in Fig. 22 the computational time of Basic and Cubature elements for the SSPRK time scheme and all stabilization techniques.
As a first interesting result of numerical test, looking at the Fig. 22, we can clearly see that, for a fixed accuracy, Cubature elements obtain better computational times with respect to Basic elements. Moreover, as expected, the SUPG stabilization technique requires more Table 5 Convergence order for all schemes on linear advection test, using coefficients obtained in Table 15 Element & SUPG OSS CIP "-" means that the Fourier analysis showed that the scheme is unstable computational time as it requires the inversion of a mass matrix, even in the case where the CFL used is larger than the ones for OSS or CIP stabilization, see Table 15.
The order of accuracy reached by each simulations is shown in Table 5. The plots and all the errors are available at the repository [43].
Looking at the Table 5, we observe that almost all the stabilized schemes provide the expected order of accuracy. Exception to this rule are several P 2 discretization which reach an order of accuracy of ≈ 2.5, and all Bernstein B 3 polynomials with the DeC which reach an order of accuracy of 2. This result is very disappointing and it does not improve even adding more corrections, as suggested in [1,3]. Moreover, it has been independently verified that also in Fourier space the accuracy of DeC with Bernstein polynomials of degree 3 is only of order 2. This problem do not show up for steady problems, as there only the spatial discretization determines the order of accuracy. We will show it in Sect. 5.3, where we study also some steady vortexes. The authors still do not understand why the optimal order of accuracy is not reached. This opens doors to further research on this family of schemes.
Note that we do not show results for Bernstein elements with SSPRK technique because they are identical to Basic elements, but are more expensive because of the projection in the Bernstein element space and the interpolation in the quadrature points.
More comparisons on different grids (unstructured) will be done in Sect. 5.

Shallow Water Equations
We consider the non linear shallow water equations (no friction and constant topography): An analytical solution of this system is given by travelling vortexes [53]. We use here a vortex with compact support and in C 6 ( ) described by Fig. 23 Error for shallow water system (57) with respect to computational time for SSPRK method with Cubature (left) and Basic (right) elements and CIP and OSS stabilizations Table 6 Convergence order on shallow water, using coefficients obtained in Table 15 Element & OSS CIP where X c = (0.5, 0.5) is the initial vortex center, (h c , u c , v c ) = (1., 0.6, 0) is the far field state, r 0 = 0.45 is the vortex radius, h = 0.1 is the vortex amplitude, and the remaining paramters are defined as ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ω = π/r 0 angular wave frequency, distance from the vortex center.
We discretize the mesh with uniform square intervals of length x (see Fig. 20a), and as before we perform a grid convergence by respecting the constraint x 2 = 2 x 1 for P 2 elements and x 3 = 3 x 1 for P 3 elements. Because of the high cost of the SUPG technique, we only compare the OSS and the CIP stabilization techniques. As an example of results, we again show the benefit of using Cubature elements in Fig. 23. We can see that since the dimension of the discretized system is even larger than before (three times larger), the differences between Cubature and Basic elements are even more pronounced in the error-computational time plot.
In Table 6 we show the convergence orders for this shallow water problem with the CFL and δ coefficients found in Table 15.
The results obtained are similar to those of the linear advection case. We can also notice the P 2 discretization reaching the proper convergence order, i.e., 3, and Bernstein B 3 elements  Table 7 Convergence order for linear advection on unstructured mesh, using coefficients obtained in Table 4 Element & SUPG OSS CIP  Fig. 17 and Sect. 3.7). "-" means that the scheme is clearly unstable reaching an order of accuracy of ≈ 3 which is more satisfying than the results obtained for the linear advection test, but still disappointing knowing that we were expecting 4.

Simulations on Unstructured Meshes
We now perform numerical tests to check the validity of our theoretical findings using an unstructured mesh, and the most restrictive parameters in Table 4. These parameters make sure that we are stable for both T and X mesh configurations. The results have similar convergence rate to the tests on the structured meshes of the previous section.
The unstructured mesh used in this section is shown in Fig. 24, and it was created by the mesh generator gmsh. 1

Linear Advection Test
We use the same test case of Sect. 4.1. Convergence orders for all schemes are summarized in Table 7. We observe that all P 1 discretizations provide the proper convergence order. For P 2 discretization we spot a slight reduction of the order of accuracy, which lays for most of the schemes between 2 and ≈ 2.5 instead of being 3. For polynomials of degree 3, we observe an order reduction to 2 for the same schemes that lost the right order of accuracy also for X mesh in the previous section. In particular, we have that Bernstein B 3 polynomials with the DeC result in an order of accuracy of ≈ 2 instead of 4, as well as theP 3 discretization with the combination DeC and SUPG stabilization. As for the X mesh, the Basic P 3 discretization reach order of accuracy ≈ 4 for all stabilization techniques, as well as CubatureP 3 with SUPG and OSS stabilizations.
Also in this case, the results obtained withP 3 Cubature elements and OSS stabilization are stable as we can see from the convergence analysis. This might mean that just few unfortunate mesh configurations, as the T one, result in an unstable scheme and that, most of the time, the parameters found in Table 4 are reliable for this scheme. On the other hand, the combinatioñ P 3 and CIP gives an unstable scheme.
We compare error and computational time for all methods presented above in Fig. 25. Looking at P 2 and the P 3 discretizations, as expected, the mass-matrix free combination, i.e., Cubature elements with SSPRK and OSS, gives smaller computational costs than other combinations with Basic elements. Conversely, the SUPG technique increase the computational costs with respect to all other stabilizations for all schemes. That is why we will not use it for the next test. The plots and all the errors are available at the repository [43].

Remark 4 (Entropy viscosity)
As remarked in Sect. 3.9, we can improve the stability of some schemes (Cubature OSS) with extra entropy viscosity. Here, we test the convergence rate on the T mesh configuration, i.e., the one with more restrictive CFL conditions and most unstable. This test is performed using CubatureP 3 elements, SSPRK and DeC time integration methods, and the OSS and the CIP stabilization techniques. We solve again problem (56).  Table 9 Convergence order on shallow water for unstructured mesh, using coefficients obtained in Table 4 Element & OSS CIP  Fig. 17). "-" means that the scheme is clearly unstable Using formulation (53) and tuning stability coefficient δ, CFL and viscosity coefficient c found in Fig. 19, we obtain fourth order accurate schemes. These tuned coefficients, and the corresponding convergence orders are summarized in Table 8.
Many other formulations of viscosity terms exist in literature and can ensure convergent methods of order p + 1 (using P p elements) [30,36,41]. The majority use a nonlinear evaluation of the parameter μ K , based on the local entropy production.

Shallow Water Equations
In this section we test the proposed schemes on the test case of Sect. 4.2 with the unstructured mesh in Fig. 24. Convergence orders are summarized in Table 9.
Also for the shallow water equations, we have results that resemble the ones of the structured mesh. There are small differences in the order of accuracy in both directions in different schemes. Comparing also the computational time of all the schemes in Fig. 26, we can choose what we consider the best numerical method for these test cases: Cubature discretization with the OSS stabilization technique. This performance seems fully provided by the free massmatrix inversion, as the CFLs for the OSS technique (with SSPRK scheme) is approximately the same between Basic and Cubature elements (see Table 4).
The plots and all the errors are available at the repository [43].

Remark on the Steady Vortex Case
For completeness we consider now a steady vortex, similarly to what reported in [3] for the isentropic Euler equations. So, we consider again the traveling vortex proposed in Sect. 4.2 with t f = 0.1s. We compare the convergence orders between u c = 0 (steady case) and u c = 0.6 (unsteady case) in Tables 10 and 11. As we can see, in the steady case we obtain, without any additional viscous stabilization, the expected convergence order for all schemes, in particular for the DeC with Bernstein polynomial function. These results agree with the  "-" means that the scheme is clearly unstable ones in [3]. Comparing with the unsteady case, all the other schemes reach similar order of accuracy as obtained in Table 9. Running the test with additional corrections in DeC scheme, as often suggested in [1,3], does not improve the convergence order in the unsteady case (even with K = 50). These results show that a numerical error appears in the spatio-temporal integration part of the solution (27), which might be related to the fact that the high order derivatives are never penalized in our stabilizations and might produce some small oscillations.

Conclusion
This work shows also that the stability results obtained in the one dimensional analysis [42] can not be generalized for two dimensional problems on triangular meshes. In this direction, it could be interesting to perform the stability analysis on Cartesian quadrilateral meshes, to check whether in that situation the one dimensional results still hold true.
In the numerical test section, the order of accuracy found is not the expected one for all the methods, i.e., p + 1 using P p elements. For several cases, we reach only p + 1/2 or p. Among the schemes that are stable and with the right order of accuracy, the method that uses Cubature elements with OSS stabilization technique and SSPRK method of order 4 has proven to be the most accurate and less expensive. Secondly, comparing to the SUPG stabilization technique, very often used in the literature for hyperbolic system, we showed that other stabilization techniques such as CIP and OSS can provide the same accuracy and are cheaper in term of computational costs.
In this direction, it would be interesting to evaluate the stability of the CIP adding a additional penalty term on the jump of higher order derivatives as suggested in [3,13,17]. Moreover, it could be interesting to see the stability of Cubature elements using higher degree polynomials. Another interesting point to explore is the loss of accuracy obtained using the DeC with Bernstein third order polynomial basis functions for unsteady cases.
Finally, we provided a heuristic approach characterized by additional discontinuity capturing viscous operators such as those proposed in [30,36]. Even for smooth solutions, the very small additional dissipation introduced by these terms is enough to stabilize some of the symmetric mass-matrix-free approaches, otherwise linearly unstable. This allows to obtain interesting schemes for practical purposes. not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Cubature Elements, Definition and Construction
In this section we give a description of the Cubature finite elements [25,29]. In Fig. 27 we show theP 3 example comparing the Lagrangian nodes of Basic and Cubature elements.
As defined in Sect. 2.2.3, there are several requirements and optimization procedures in order to obtain the Cubature elements. These elements are very import in our study because they permit to obtain diagonal mass matrix, and so they decrease considerably the time of computation. We describe for p = 1, 2, 3 the basis functions of the Cubature elements.
Corresponding weights are w v i = 1 3 .

Cubature Elements of Degree 2
TheP 2 element contains 7 degrees of freedom: three at the vertices v 1 , v 2 and v 3 and three at the midpoint of the edges that we denote as e i j =       -At the internal points

Time Discretization Coefficients
In this appendix we introduce the time integration coefficients used in this work, to make the study fully reproducible. In Table 12 there are the RK coefficients, in Table 13 the SSPRK coefficients and in Table 14 the DeC coefficients.

Fourier Analysis
In this section we collect all the plots and results that are essential to show the results of this work, but for structural reasons were not put in the main text.

Mesh Types and Degrees of Freedom
We represent in Fig. 28 the mesh configurations used in the Fourier analysis and the degrees of freedom of the elements of degree 3. The red square represents the periodic elementary unit that contains the degrees of freedom of interest for the Fourier analysis. The symbol "/" means that the fourier analysis for the scheme results always in instability. The values denoted by * are not the optimal one, but they lay in a safer region, see Sect. 3.6

Fourier Analysis Results: Optimal Parameters
In this section, we put the optimal values of the stability analysis of Sect. 3.5 after the modification proposed in Sect. 3.6. In Table 15 we show the parameters for the X mesh and in Table 16 we show the parameters for the T mesh.

Fourier Analysis Results: Stability Area
Finally, we present a comparison of stability area between the T and the X mesh. This comparison if perform as before, for all wave angles θ . We choose as example the comparison using Basic element, SSPRK time integration method and the OSS stabilization technique in Fig. 29. The interested reader can access to results for all methods online [43].