A flux-differencing formulation with Gauss nodes

In this work, we propose an extension of telescopic derivative operators for the DGSEM with Gauss nodes, and we prove that this formulation is equivalent to its usual matrix counterpart. Among other possible applications, this allows extending the stabilization methods already developed for Gauss-Lobatto nodes to Gauss nodes, also ensuring properties such as entropy stability while retaining their improved accuracy.


Introduction
Among the different spatial discretization frameworks, the Discontinuous Galerkin Spectral Element Method (DGSEM) [1,2] approximates a function as a set of discontinuous, high-order polynomials defined in a tessellation of the spatial domain. The DGSEM is a collocation method that represents the solution as a combination of Lagrange polynomials that shares the nodes with a certain quadrature rule, usually Gauss or Gauss-Lobatto.
Fisher et al. showed in [3] that diagonal-norm summation-by-parts (SBP) derivative operators can be rewritten in telescopic form. In this case, the action of the operator resembles a finite volume scheme, ensuring local conservation in the sense of Lax-Wendroff. In the context of the DGSEM, the derivative operator with Gauss-Lobatto nodes falls in this category, whereas the use of Gauss nodes leads to a generalized SBP operator [4,5,6,7,8].
In this work we show that the generalized formulation of Chan [8] also admits a flux-differencing form and thus, from a more practical standpoint, entropy-stable schemes can be generated by introducing certain dissipative numerical fluxes at the interfaces of the sub-elements [9,10,11].

High-order DGSEM on Gauss nodes
Beginning with a one-dimensional grid where ξ i and ω i are the nodes and weights of the Gauss quadrature as in fig. 1,ξ Fisher et al. showed in [3] that SBP derivative operators, D, can be rewritten in telescopic form. In this case, the action of the operator resembles a finite volume scheme, As stated in the introduction, the derivative operator of the DGSEM with Gauss-Lobatto nodes falls in this category, while we can use the developments of Chan [7,8] if Gauss nodes are used. Following Chan's notation, the DGSEM discretization of a generic one-dimensional conservation law, u t + f x = 0, is The mass matrix, M = diag(Jω i ), contains the Jacobians and weights of the quadrature rule in each mesh element. The projection of the interior values to the faces is represented by V f , and B = diag(−1, 1). The matrix Q is a generalized SBP operator representing the integral of the derivative of the basis functions, where l i (x) is the Lagrange polynomial associated to the node ξ i . The numerical fluxes are indicated with a star (f L and f R ), and sub-indices with a tilde mean that the magnitude has been computed on so-called entropy projected variables, Finally, F S is the matrix with all the combinations of entropy-conservative fluxes (including the entropyprojected values at the interfaces), F S ij = f S (u i , u j ), and A•B is the Hadamard product of matrices A and B. More details can be found in [7].
Imposing the equality of eq. (2) with a finite volume scheme, where ∆f represents the flux-differencing formulation of eq. (1), we find an expression for the different subcell fluxes,f , With the additional constraintsf 0 = f L andf N +1 = f R , we can compute the value off i+1 fromf i and eq. (4), obtaining an expression for the interface fluxes that resembles the one developed by Rueda-Ramírez et al. [12] for the semi-discretization,f We remark that eq. (2) is recovered whenf 0 andf i+1 are introduced in eq. (3) by construction. The last two terms of eq. (5) represent a new addition with respect to eq. (3.9) of [9], and they couple the volume and surface integrals. These new operations result in a higher computational cost. However, the use of Gauss nodes also entails a higher accuracy that can compensate it as shown in [12].
There is, however, a problem with eq. (5). For a set of N + 1 Gauss nodes there are N + 2 complementary staggered fluxes, but we have N + 3 equations. The flux at the right face,f N +1 , can be computed fromf i+1 (i = N ), but it also must be equal to f R . We can overcome this issue by proving that both equations are equivalent and our system is not over-constrained.
Now we set i = N and apply different simplifications to prove that eq. (6) is equivalent tof N +1 = f R . In this case, the third and fourth terms of the right-hand side can be further simplified by considering that, for the Lagrange interpolating polynomials, N i=0 l i (x) = 1. Substituting this allows us to "exchange" f L in the first term by f R , Written in this form we can now apply the symmetry property of the entropy-conservative numerical flux, Finally, we remark that S ij = −S ji and thus, the product of the last term is the Hadamard product of the skew-symmetric matrix S with the symmetric matrix f S ij . This is, in fact, another skew-symmetric matrix and the last term is also zero,f N +1 = f R .
Note that this proof is also applicable to the standard Gauss DGSEM scheme since the expression for the staggered fluxes is similar, simply lacking the terms containing entropy-projected quantities.

Subcell limiting
The existence of an underlying flux-differencing formula for the Gauss-DGSEM enables the use of the subcell limiting strategies presented by Rueda-Ramírez et al. [11] to improve the robustness of the method, e.g., in the presence of shocks. In particular, we propose a hybrid scheme obtained as a convex combination of the high-order Gauss-DGSEM with a first-order FV method, with m i the entries of the diagonal mass matrix. The individual fluxes are given aŝ wheref FV i is a robust first-order approximation of the flux and α i is a so-called blending coefficient, which is selected such that the resulting scheme exhibits some desired properties, e.g., positivity, non-oscillatory behavior, etc.
When using the DGSEM on Gauss-Lobatto nodes, it is possible to combine the high-order and low-order methods at the element level, i.e., with α i a piece-wise constant function that is different at each element.
This limiting technique, known as element-wise blending, does not require a flux-differencing formula and has been shown to retain conservation for any choice of α [10,13]. The proofs in [10,13]

Results
This section includes some numerical results that confirm the theoretical developments of sections 2 and 3.
In both cases we solve the Euler equations using the local Lax-Friedrichs numerical flux. In section 4.1, the non-dissipative term uses the two-point flux from Chandrashekar [14], whereas we employ a simple average in section 4.2.

Convergence
We test the numerical accuracy of the subcell approach by comparing it against the formulation obtained by Chan for Gauss nodes [8]. Considering the solution to the one-dimensional Euler equations,

Sedov blast
To illustrate the shock-capturing capacity of the hybrid DGSEM/FV method, we simulate a Sedov blast problem describing the evolution of a blast wave expanding from an initial concentration of density and pressure. For the initial condition, we assume a gas in rest, v 1 (t = 0) = v 2 (t = 0) = 0, with a Gaussian distribution of density and pressure, where we choose σ ρ = 0.25 and σ p = 0.15. Furthermore, the ambient density is set to ρ 0 = 1 and the ambient pressure to p 0 = 10 −1 . We complement the simulation domain, Ω = [−1, 1] 2 , with periodic boundary conditions, and tessellate it using K = 64 2 quadrilateral elements. We represent the solution with polynomials of degree N = 3 and run the simulation until t = 1.
To avoid non-physical oscillations in the vicinity of shocks, we use a subcell-wise limiting strategy (8) to impose a non-oscillatory behavior on the density of each node i, (9) where N (i) is the so-called low-order stencil of node i, a set containing all the neighboring nodes to i, u ij is the so-called bar state, λ max ij is an estimate of the maximum wave speed between nodes i and j, andˆ n ij is the normal vector (normalized metric term) at the interface between nodes i and j. We refer the reader to [11] for details on how to compute α with an algebraic flux correction scheme, such that (9) is guaranteed.
Using the procedure described in [11], we obtain a provisional blending coefficientα i at each node of our domain. The subcell-wise limiting strategy (8) is then applied taking the maximum blending coefficient on both sides of every interface, α i = max(α i ,α i+1 ). The same procedure is applied at element interfaces when using Gauss nodes. As explained above, the surface terms of DG and FV are equal when using Gauss-Lobatto nodes. Hence, it is not necessary to blend the fluxes at the interfaces. Figure 3 illustrates the density and blending coefficient in the domain obtained with the hybrid DGSEM/FV method using Gauss and Gauss-Lobatto nodes at time t = 0.6. Both schemes capture the expanding shock correctly while using the high-order DG method in most of the domain.

Gauss
LGL Gauss LGL We compute the mean blending coefficient at a particular time as where e ∈ [1, K] denotes the element index, K is the number of elements of the domain, i, j ∈ [0, N ] are the node indices, N is the polynomial degree,α e ij (t) is the provisional (nodal) blending coefficient of node ij of element e at time t, and V is the total area of the domain. Figure 4 shows the evolution of the mean blending coefficient and number of time steps taken for the simulation of the blast wave with CFL=0.9. For this particular example and the non-oscillatory condition on density (9), the ES Gauss scheme requires more limiting than the ES LGL scheme for the same number of degrees of freedom, polynomial degree, and CFL number. However, the LGL scheme takes more time steps to finish the simulation as it requires shorter time-step sizes.
The maximum allowable time-step size that is required for condition (9) in 2D reads [16] ∆t ≤ min .
The larger quadrature weights of the Gauss-DGSEM at the boundary nodes of the elements result in a less restrictive CFL condition. As a result, the DGSEM simulation with Gauss nodes is not only more accurate, but also computationally more efficient than the LGL simulation for the same number of degrees of freedom.

Conclusions
We have presented in this work a novel flux-differencing expression of the entropy-conservative formulation of J. Chan [7,8] in terms of staggered fluxes. This approach can be applied explicitly to the DGSEM framework with Gauss nodes. We have also proved that this telescopic formulation of the derivative operators exists with and without entropy-projected variables.
When applied to a simple test case, we have shown that the convergence properties of the flux-differencing method match those of the entropy-conservative approach used as the baseline. This allowed us to apply some subcell limiting strategies already developed for the DGSEM with Gauss-Lobatto nodes, although more work is needed in this aspect.

Appendix A. Extension to higher dimensions and curvilinear grids
The flux-differencing formula (3) extends to multiple space dimensions on curvilinear grids. We can write the Gauss-DGSEM in two dimensions as where we introduce a new notation:f (i−1,i)j denotes the high-order telescoping flux between nodes (i − 1, j) and (i, j), andf i(j−1,j) denotes the high-order telescoping flux between nodes (i, j − 1) and (i, j).
The telescoping fluxes are uniquely defined between two neighboring nodes. For instance, in the first coordinate direction they read, where the surface numerical fluxes depend on the inner entropy-projected solution, the outer solution, and the normal vector at the boundary of the element, The volume numerical fluxes are entropy-conservative two point fluxes with metric dealiasing (see, e.g., [12]), f 1S (i,k)j = f S (u ij , u kj ) · J a 1 (i,k)j , and we use the so-called contravariant metric vector, a 1 := ∇ξ, which relates the physical-frame coordinates (x, y) with the reference-frame coordinates (ξ, η), and the average operator,