On the Connections Between Discontinuous Galerkin and Flux Reconstruction Schemes: Extension to Curvilinear Meshes

This paper investigates the connections between many popular variants of the well-established discontinuous Galerkin method and the recently developed high-order flux reconstruction approach on irregular tensor-product grids. We explore these connections by analysing three nodal versions of tensor-product discontinuous Galerkin spectral element approximations and three types of flux reconstruction schemes for solving systems of conservation laws on irregular tensor-product meshes. We demonstrate that the existing connections established on regular grids are also valid on deformed and curved meshes for both linear and nonlinear problems, provided that the metric terms are accounted for appropriately. We also find that the aliasing issues arising from nonlinearities either due to a deformed/curved elements or due to the nonlinearity of the equations are equivalent and can be addressed using the same strategies both in the discontinuous Galerkin method and in the flux reconstruction approach. In particular, we show that the discontinuous Galerkin and the flux reconstruction approach are equivalent also when using higher-order quadrature rules that are commonly employed in the context of over- or consistent-integration-based dealiasing methods. The connections found in this work help to complete the picture regarding the relations between these two numerical approaches and show the possibility of using over- or consistent-integration in an equivalent manner for both the approaches.


Introduction
Popularity of high-order compact discretisations using spectral/hp element spatial approximations is rising rapidly and their deployment to industrial-type problems is becoming a real possibility. Spectral/hp methods can achieve an arbitrary order of accuracy and they are significantly less dissipative than the more traditional low-order methods making them a key tool in some areas of computational fluid dynamics, such as in the numerical simulation of unsteady flows. On the other hand, spectral/hp element discretisations are still affected by a lack of numerical stability which makes them not highly reliable unlike low-order methods. One of the most popular approaches to discretise hyperbolic conservation laws is based on the discontinuous Galerkin (DG) method firstly introduced by Reed and Hill [1]. Usually, the DG method is formulated in a weak form where the domain is divided into non-overlapping elements and on each of these elements a basis of polynomials and a set of quadrature points are chosen to calculate the integrals arising from the weak formulation of the problem. Some of the most efficient and ubiquitous forms of the DG method are the so-called nodal DG schemes, whereby Lagrange interpolants are combined with a set of nodal solution points on a given element [2]. We note that in a general nodal DG scheme, given an expansion on each element in terms of P Lagrange interpolants, we may perform quadrature on a separate set of Q quadrature points. If P = Q and the distribution of solution and quadrature points is identical, then we recover the so-called discontinuous Galerkin spectral element method (DG SEM ) which diagonalises the mass matrix, allowing for further computational optimisation [3,4].
In contrast to the DG method which makes use of a weak representation of the equations, a set of schemes based upon the differential formulation have been recently introduced. They offer another route which avoids the need for calculating integrals and therefore defining quadrature rules, making these schemes attractive from an implementation point of view. The first of these, namely the spectral-difference scheme, was introduced by Kopriva and Kolias [5] in 1996 and was extended to quadrilateral and triangular elements by Liu et al. [6] in 2006. Most recently, the flux reconstruction (FR) method was presented by Huynh [7].
The FR approach encapsulates various energy stable discretisations where the particular scheme recovered is dictated by a single real-valued parameter and is referred as Vincent-Castonguay-Jameson-Huynh (VCJH) scheme [8]. The FR approach allows one to recover not only differential-type schemes such as a particular SD method, but also integral-type schemes such as nodal DG schemes.
A general connection between FR and nodal DG schemes has already been examined in other works [7][8][9]. A more detailed inspection of the connections between these two schemes was carried out in [10] where it has been shown how for linear and nonlinear advection equations various equivalences exist between DG and FR schemes on regular grids. However, the extension to irregular meshes and the associated aliasing issues have not yet been explored. We consider a nodal DG scheme with Q > P and two DG SEM schemes. In the following with the subscript SEM 1 we mean that a collocation quadrature rule is used for the inner product of the advection term of the conservation law equation.
The relationships found in this paper are summarized in Table 1. In particular, we demonstrate that the equivalences found in [10] between the DG SEM scheme with lumped mass indicates that the schemes are equivalent, whereas ✗ indicates differences between the schemes * The equivalence holds true for Gauss-Lobatto-Legendre points only matrix (LMM) and the FR HU scheme introduced in [7] (also referred to as FR g2 ) as well as the equivalences between the DG SEM scheme with exact mass matrix (EMM) 2 and the FR DG scheme (also presented in [7]) hold true for irregular tensor-product meshes (i.e. deformed/curved tensor-product grids). These equivalences are complemented with some numerical experiments on regular (used as a baseline) and irregular grids for linear and nonlinear problems.
We also show numerically that, using Q > P, where Q is the number of quadrature points and P are the Lagrange interpolants inside each element, further extends the equivalences between DG and FR DG . Specifically, we found that DG (Q > P) and FR DG(Q > P) are identical when using Q > P for linear and nonlinear flux functions as well as for regular and irregular tensor-product meshes 3 . This indicates that polynomial aliasing sources for these two schemes are identical and can therefore be addressed using equivalent strategies, such as the consistent integration (through additional quadrature points) of the nonlinearities arising either from the equations themselves or from the geometry [11,12].
We finally present some results related to the computational time required by the FR and DG schemes. This paper is organised as follows. In Sect. 2 we prove theoretically the connections in the first four columns of Table 1, in Sect. 3 we further assess the theoretical work with some numerical experiments for both linear and nonlinear problems and we show numerically how DG (Q > P) and FR DG(Q > P) are identical and, finally, in Sect. 4 we draw the conclusions.

Theory
We describe the DG and FR schemes in the context of a 2D scalar conservation law. We assume that the quadrilateral elements are deformed/curved (i.e. possess spatially varying Jacobians) and we prove that the FR DG scheme and DG SEM method with an exact mass matrix evaluation are equivalent on irregular grids (i.e. spatially varying Jacobians) when the same polynomial approximation of the geometry is employed. We also demonstrate that the FR HU scheme and the DG SEM method with lumped mass matrix are equivalent on irregular quadrilateral grids.

2D Scalar Conservation Law on Irregular Quadrilateral Grids
Consider the following 2D conservation law where f = f (u) and g = g(u) are the advection fluxes in the x and y directions respectively. The domain is partitioned into N non-overlapping, quadrilateral elements n which can be deformed or curved Each quadrilateral element n can be mapped into a reference quadrilateral s = [−1, 1] 2 in the transformed space x = (ξ, η) using the generic mapping In the case of curvilinear elements we can use the standard isoparametric mapping which, for an arbitrary-shaped curved-sided two-dimensional element, is: where φ p and φ q are the same basis functions used for representing the solution. Analogous mappings can be adopted for three-dimensional curvilinear elements. For additional details the interested reader can refer to [13].

FR DG Scheme as DG SEM Method with EMM
To demonstrate that the DG method is identical to FR DG in the case of irregular grids we start multiplying Eq. (1) by a two-dimensional Lagrange polynomial (x) and integrating the equation over a local element n n ∂u δ where f δ n is the elemental polynomial flux function and u δ n is the elemental polynomial solution. We successively perform an integration by parts, substitute the boundary terms by the numerical interface fluxes and integrate by parts once more in order to derive the strong form of the DG method where n is the normal, f δ I n is the numerical interface flux and the superscript D indicates discontinuous polynomial flux functions which are directly evaluated at the solution points. We now substitute the operations on the local element n onto the reference element s in order to highlight the metric terms where b DG n is the boundary term which couples the elements of the spatial discretisation together. We can re-express Eq. (7) in a form more amenable to implementation and which keeps all the geometric terms as standard polynomials as: where we used the following relationships The expression of b DG n on the reference element s depends on the edge considered and for the bottom edge is where Using similar relations for the other three edges of s , Eq. (12) becomes We now consider the FR DG approach where we define the following transformations between the local element n and the reference element s : where the polynomial metric terms J 2D n , ∂ x ∂ξ , ∂ x ∂η , ∂ y ∂ξ , ∂ y ∂η can be evaluated, for instance, from Eq. (4). The FR approach in a reference element can be written as We now expand Eq. (13) as follows where b F R n is defined as with i being the derivatives of the correction functions which will be defined later. For the bottom edge of the reference quadrilateral element we can write where n 0 = (0, 1) and G −1 and J 2D n are defined in Eq. (9). Note that the relation in Eq. (16) is necessary to transform the flux jump from the physical space to the reference element (for additional details see [14]). As final form we obtain Using similar relations for the other edges of the reference quadrilateral element, Eq. (14) becomes Now we require the residual to be orthogonal to a set of smooth functions where we used the relation (ξ ) = i (ξ ) j (η). As shown by Huynh [7]: where 3 and 0 are the right Radau polynomials of order P + 1 on [−1, 1] which vanish at the right and top edges ξ = η = 1. This is due to the orthogonality of the right Radau polynomial of order P + 1 to all polynomials of order up to P − 1. Analogously, we can write which is identical to Eq. (12) and proves the statement.

Remark
The equivalence above holds true for any point distribution used for defining the Lagrange basis.

FR HU Scheme as DG SEM Method with LMM
Consider Eq. (12) with solution coefficients represented by Lagrange polynomials at (P + 1)×(P +1) GLL points. If we additionally choose a GLL quadrature for the DG SEM scheme, it is well known that only polynomials up to order 2P − 1 are exactly integrated (to machine precision). The associated elemental mass matrix will contain a numerical quadrature error since each integrand is a polynomial of order 2P. However, because the Lagrangian basis possesses the property that m (ξ n ) = δ m,n where δ m,n represents the Kronecker delta function, this quadrature rule gives a diagonal mass matrix, where w i, j are the weights of the GLL quadrature rule. With this quadrature rule the boundary integrals simplify and each flux jump modifies only its own boundary (edge) solution point. Equation (19) can therefore be written as If we now consider a collocation quadrature rule with solution values located at GLL points we obtain the same diagonal mass matrix as before, M = diag(w i, j J 2D n i, j ). The volumetric integrals associated to Eqs. (12) and (19) are evaluated in the same way and, therefore, in order for the two schemes to be the same we need to find a correction function which transforms the FR boundary term into the DG one. Specifically, if we use as correction functions those recovering the FR HU scheme, which can be found in [8], their derivatives vanish at all solution points except at the left/bottom edges if evaluating 3 (ξ) / 0 (η) or at the right/top edges if evaluating 1 (ξ) / 2 (η) (the zeros of are at GLL points). Thus, in the FR HU scheme the corrective flux modifies only the edge points. We can adopt the following transformation: because the derivatives of the correction functions in vector form evaluated at the (P + 1) × (P + 1) GLL points are the following: The boundary term of FR HU scheme is equal to the boundary term of the DG SEM scheme with lumped mass matrix and therefore the two schemes are equivalent.
Remark Note that this equivalence holds true only for GLL points.

Numerical Experiments
In this section we present the numerical results obtained for both linear and nonlinear problems. We used two different mesh configurations: a single-element mesh (also referred as Mesh A) whose Jacobian determinant is shown in Fig. 1 and a multi-element mesh (also referred as Mesh B) whose Jacobian determinants are depicted in Fig. 2. The curved edges in Mesh A are described through a parabolic function f 1 , which, defined within a onedimensional reference element, is where A 1 is a constant which determines the amplitude of the deformation and was set to zero for the mesh in Fig. 1a (i.e. we obtained a regular mesh), while for the meshes in Fig. 1b, c A 1 was equal to 0.2 and 0.4 respectively. Mesh B was generated using a sinusoidal function f 2 where A 2 is a constant which set the deformation amplitude and was set to zero for the mesh in Fig. 2a (i.e. we obtained a regular mesh), while for the meshes in Fig. 2b, c A 2 was equal to 0.1 and 0.2 respectively. Note that similar meshes were used in [15] for studying the accuracy and efficiency of several discontinuous high-order formulations on curved quadrilateral elements.
For both mesh configurations we tested the linear advection equation and the nonlinear compressible Euler equations. The test cases considered span all the possible combinations in Table 1.
We tested two different set of points, Gauss-Lobatto-Legendre (GLL) and Gauss-Legendre (GL) where the interpolations required at the boundaries for GL points were performed in a consistent manner for the various DG and FR schemes considered. Note that we show the results for the regular and most deformed grids on GLL points only for the sake of brevity. The conclusions we can draw from the results obtained on GL points and on the grids in the Figs. 1b and 2b are identical to those which can be drawn from the results presented here. Some additional results concerning GL points and the two grids in the Figs. 1b and 2b are reported in the "Appendix".
Throughout the numerical results we will make extensive use of relative L 2 errors defined as where u i is the numerical solution calculated at each quadrature point, u i,exact is the exact solution and N is the total number of quadrature points and w i are the weights.

Linear Problem
The first series of numerical experiments was carried out on the two-dimensional linear advection equation where τ x = τ y = 2 for Mesh A in Fig. 1 and τ x = τ y = 0.5 for Mesh B in Fig. 2. The initial conditions were applied using a collocation projection for FR HU , DG SEM , FR DG while we used a higher order projection (employing an additional quadrature point) for FR DG(Q > P) and DG (Q > P) . For both mesh configurations, we applied exact boundary conditions and an explicit 4th-order Runge-Kutta scheme for the time-integration with a final time T = 2s and a time-step sufficiently small to consider the temporal error negligible. We employed four different polynomial orders to discretise the solution. Specifically, for the meshes in Fig. 1 we used P = 10, 11, 12 and 13, whilst, for the meshes in Fig. 2, we applied P = 5, 6, 7 and 8. The main parameters used are summarized in Table 2. Table 2 Parameters used for the linear problem: Mesh A refers to the meshes represented in Fig. 1; Mesh B refers to the meshes represented in Fig. 2 a τ x , τ y T P T-I   Figure 3 represents the L 2 error vs the polynomial order obtained using GLL points for the various schemes tested where Fig. 3a refers to the regular mesh represented in Figs. 1a and 3b refers to the mesh depicted in Fig. 1c. Table 4 in the "Appendix" quantifies up to sixteen digits the results obtained for P = 10, for both GLL and GL points and all the three grids in Fig. 1.

Mesh A
The results for the other polynomial orders are not tabulated for the sake of compactness and because they provide the same information as those which can be obtained from Table 4. The following equivalences hold true for all the polynomial orders considered: • FR HU and DG SEM with lumped mass matrix on GLL points; • FR DG and DG SEM with exact mass matrix on any point distribution; • FR DG(Q > P) and DG (Q > P) on any point distribution.
The first two results confirm the demonstrations presented in the previous section. In addition, the last result indicates that these connections are valid also when using additional quadrature points than a standard collocated spectral/hp element method. Also, the equivalences are up to machine precision for all the deformation levels considered as shown in Table 4. From a polynomial aliasing perspective, where the polynomial aliasing is introduced by the curved geometry, it is interesting to note how the magnitude of the error increases as the deformation level of the mesh increases. Specifically, we can see that using either FR DG(Q > P) or DG (Q > P) provides a different result (generally a better L 2 error) for the deformed meshes. 4  Table 5 in the "Appendix" quantifies up to sixteen digits the results obtained for both GLL and GL points for P = 5 and all the three grids in Fig. 2. The connections presented for the single-element mesh configuration in the previous subsection hold true also for the multi-element mesh configuration.
From a polynomial aliasing point of view it is possible to note how the differences between the pairs FR DG , DG SEM with exact mass matrix and FR DG(Q > P) , DG (Q > P) are still present although less marked.

Nonlinear Problem
The second series of numerical experiments were undertaken using the two-dimensional compressible Euler equations which can be written as follows where q = (ρ, ρu, ρv, E) is the vector of the conserved variables and f i 1 = f i 1 (q) and f i 2 = f i 2 (q) are the vectors of the inviscid fluxes, In the above, ρ is the density, u and v are the velocity components in the x-and y-direction respectively, p is the pressure and E is the total energy for which we assumed the perfect gas law. Note that γ denotes the constant ratio of specific heats of the gas. Table 3 Parameter used for the nonlinear problem: Mesh A refers to the meshes represented in Fig. 1; Mesh B refers to the meshes represented in Fig. 2 β As test case we considered the isentropic vortex problem whose initial conditions on a two-dimensional grid can be written as follows where with β = 5, R = 1 and γ = 1.4. For both mesh configurations, we used an explicit 4th-order Runge-Kutta scheme for the time-integration with a final time T = 2s and a timestep sufficiently small to consider the temporal error negligible. We employed four different polynomial orders to discretise the solution. Specifically, for the meshes in Fig. 1 we used P = 10, 11, 12 and 13, while, for the meshes in Fig. 2, we applied P = 5, 6, 7 and 8. The main parameters used are summarized in Table 3. Figure 5 represents the L 2 error associated to the density (ρ) vs the polynomial order obtained using GLL points for the various schemes tested where Fig. 5a refers to the regular mesh represented in Figs. 1a and 5b refers to the mesh in Fig. 1c. Table 6 in the "Appendix" quantifies up to sixteen digits the results obtained for both GLL and GL points for P = 10 and all the three grids in Fig. 1. Also in the case of a nonlinear problem the same equivalences presented before maintain their validity up to machine precision as shown in Table 6. However, considerations on polynomial aliasing are different because in this case the aliasing sources arise both from the equations themselves, which are nonlinear, and from the geometry. Significant, in particular, is the gap between the pair FR DG , DG SEM and the pair FR DG(Q > P) , DG (Q > P) . Figure 6 represents the L 2 error associated to the density (ρ) versus the polynomial order obtained using GLL points for the various schemes tested where Fig. 6a refers to the regular mesh represented in Figs. 2a and 6b refers to the mesh in Fig. 2c. Table 7 in the "Appendix" quantifies up to sixteen digits the results obtained for both GLL and GL points for P = 5 and all the three grids in Fig. 2. Identical considerations can be made for both connections and polynomial aliasing issues as in the previous subsection.

Comparison of Computational Time
In this subsection we present the comparison between the various schemes considered in terms of computational time. In particular, we show the results for the linear problem on the multi-element mesh (mesh B) for both the regular case (A 2 = 0.0) and for the deformed case (A 2 = 0.2). Note also that we show the results for both GLL and GL points. The results on the single element mesh (mesh A) as well as the results for the nonlinear case are not presented since the conclusions we can draw are identical to the one we can obtain from the results presented. Figure 7 shows the CPU-time required to perform a time-step for both GLL and GL point distributions and for four different schemes on the regular mesh (A 2 = 0.0). We can note that the FR HU scheme provides a similar computational time per time-step if compared to the DG SEM -LMM scheme and that the FR DG(Q > P) scheme has a similar computational cost if compared to the DG (Q > P) scheme. Given that the DG and FR pair of schemes considered are mathematically identical one expects similar if not equal computational times per time-step. The differences observed are mainly due to the implementation of the two approaches and in particular to their different data structures. In fact, the FR approach has been implemented in the standard element space while the DG approach has been implemented in the local element space. We can also observe that the collocated schemes, FR HU and DG SEM -LMM, provides more efficient algorithms, while using the non-collocated counterparts, FR DG(Q > P) and DG (Q > P) , gives a higher computational time per time-step. This is due to the lower number of operations required per time-step when adopting a collocation projection than using a higher quadrature. Also, it can be seen how the use of GL points (Fig. 7b) provides a higher computational time per time-step than the use of GLL points (Fig. 7a). When using GL points, in fact, it is necessary to interpolate the solution and the fluxes at the boundaries of the element. This additional step increases the operation count thus the overall computational costs of a given time-step. Figure 8 shows the same results as Fig. 7, with the only difference that in this case the mesh considered is the deformed one (A 2 = 0.2). We can see that the overall cost of a given time-step is higher than the regular-mesh case. In this case in fact the operation count is larger due to the point-wise nature of the geometric factors. Also in this case it is possible to observe a higher computational costs per time-step for the case of GL points (Fig. 8b) if compared to GLL points (Fig. 8a) and the collocated schemes (FR HU and DG SEM -LMM) perform better than the non-collocated ones (FR DG(Q > P) and DG (Q > P) ).

Summary
In this section we presented several numerical experiments for both a linear and a nonlinear problem on two mesh configurations. For each mesh configuration we applied three different deformation levels, the first was a regular mesh, which was used as a baseline, while the other two were meshes with incremental deformation levels. In particular we showed how the identities shown in Table 1 and proved in Sect. 2 are further assessed (up to machine precision) across all the tests performed. An additional result shown in this section is the equivalence of FR DG(Q > P) and DG (Q > P) schemes which implies that the aliasing issues arising in FR DG and DG SEM are numerically the same and can be alleviated using a higher quadrature Q > P, i.e. consistent integration of the nonlinearities arising in the problem (although the concept of integration could be seen as out of context for an FR scheme because there are no integrals present in the formulation of this approach). Note that the use of a better point distribution with a more powerful quadrature, such as the GL points, can also alleviate aliasing issues. However, when the nonlinearity cannot be fully described by the GL quadrature, then also in this case one could use additional quadrature points for consistently integrate the nonlinearity sources. In addition, we showed the computational costs associated to the FR and DG schemes for both the GLL and GL point distributions on the linear advection problem. Although no attempt was made to optimise the underlying algorithms of the two numerical approaches, it can be seen that both provide similar computational costs as expected. The differences observed can be imputed to the different data structures between the FR and the DG approach. Specifically, the first was implemented in the standard element space (ξ , η), while the second on the local element space (x, y).

Conclusions
In this paper we established that the connections between various discontinuous Galerkin methods and high-order FR schemes are also valid on quadrilateral irregular/curved meshes. In addition, we found that the polynomial aliasing sources for the FR and DG schemes considered are identical and can be addressed using the same techniques. Specifically, we established that the FR and DG schemes taken into account are identical when a higher-order quadrature rule, usually adopted when applying over-integration-based (also referred to as consistent-integration) dealiasing strategies, is employed.
The schemes considered for the DG approach are the DG SEM with a lumped mass matrix (LMM) and a collocation projection of the solution and of the inner product, the DG SEM with an exact mass matrix (EMM) and a collocation projection of the solution and of the inner product and the DG (Q > P) with exact mass matrix (EMM) and an additional quadrature point for representing the solution and performing the inner product (i.e. using a Galerkin projection of the solution and of the inner product). For what concerns the FR schemes we took into account the FR HU scheme, the FR DG scheme and FR DG(Q > P) scheme which used an additional quadrature point to represent the solution (i.e. using a Galerkin projection of the solution).
We found that the connections between discontinuous Galerkin methods and high-order flux reconstruction schemes explored in [10] for regular grids hold true also for irregular/curvilinear meshes. In particular we mathematically proved the equivalences between FR HU and DG SEM with lumped mass matrix and the connections between FR DG and DG SEM with exact mass matrix. Both demonstrations were further assessed by numerical experiments on two different mesh configurations at different grid deformation levels and for both linear and nonlinear problems. In addition, we showed numerically how these connections are valid when applying a higher quadrature Q > P (i.e. using more quadrature points than a collocation projection: FR DG(Q > P) and DG (Q > P) ) whose use can be crucial for stabilising problems with significant polynomial aliasing issues such as under-or marginally-resolved problems (e.g. high Reynolds number simulations). This result indicates that the aliasing sources for FR DG and DG SEM are identical and can be addressed using identical strategies.
The latter result is particularly significant because it shows how the machinery for the DG method for tackling aliasing issues and improving the numerical stability of this class of schemes can be directly deployed to the FR schemes.
Finally, we showed that, in our current code-base, the computational costs are similar when considering a pair of identical FR and DG schemes. The differences observed are mainly imputable to the different implementation-related to the data structures-of the two approaches.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Tabulated Numerical Results
In this appendix we show some tabulated values of the numerical experiments carried out for both the linear and the nonlinear problem on all the grids and for both the point distributions considered.