An extended range of stable-symmetric-conservative Flux Reconstruction correction functions

The Flux Reconstruction (FR) approach offers an efficient route to achieving high-order accuracy on unstructured grids. Addi-tionally, FR offers a flexible framework for defining a range of numerical schemes in terms of so-called FR correction functions. Recently, a one-parameter family of FR correction functions were identified that lead to stable schemes for 1D linear advection problems. In this study we develop a procedure for identifying an extended range of stable, symmetric, and conservative FR correction functions. The procedure is applied to identify ranges of such correction functions for various orders of accuracy. Numerical experiments are undertaken, and the results found to be in agreement with the theoretical findings.


Introduction
High-order methods for computational aerodynamics on unstructured grids offer the promise of increased accuracy at reduced cost, within the vicinity of complex engineering geometries. As such they have garnered continued interest over the past decades. However, to-date, their 'real-world' adoption in both industry and academia remains limited [1]. In 2007 Huynh proposed the Flux Reconstruction (FR) approach to high-order methods [2]. Based on a differential form of the governing system, it is hoped FR (also referred to as Lifting Collocation Penalty [3] or Correction Procedure via Reconstruction [4]) will facilitate adoption of high-order methods amongst a wider community of fluid dynamicists.
Various properties of FR schemes, including their dispersion and dissipation characteristics [5,6], their associated Courant-Friedrichs-Lewy (CFL) limit [2,5], and their fundamental stability [7], are all determined in full or in part by the form of their associated FR correction functions. These correction functions act to lift inter-element flux jumps from the boundary into the interior of each element. Building on the work of Huynh [2] and Jameson [8], Vincent, Castonguay and Jameson recently identified a one-parameter family of correction functions that lead to stable FR schemes for 1D linear advection problems [7]. Identification of these correction functions, henceforth referred to as Vincent-Castonguay-Jameson-Huynh (VCJH) correction functions, provided significant insight into stability properties various FR schemes. However, further work is required in order to determine a full specification of the necessary and sufficient conditions that should be imposed on correction functions in order to guarantee stability.
In this study we develop a procedure for identifying an extended range of stable, symmetric, and conservative FR correction functions. The procedure is applied to identify ranges of such correction functions for various orders or accuracy. In all cases the original one-parameter VCJH correction functions are found to be a sub-set of the extended ranges. Numerical experiments are undertaken in order to verify the theoretical findings.

Overview
FR schemes are similar to nodal DG schemes, which are arguably the most popular type of unstructured high-order method (at least in the field of computational aerodynamics). Like nodal DG schemes, FR schemes utilise a high-order (nodal) polynomial basis to approximate the solution within each element of the computational domain, and like nodal DG schemes, FR schemes do not explicitly enforce inter-element solution continuity. However, unlike nodal DG schemes, FR methods are based solely on the governing system in a differential form. A description of the FR approach in 1D is presented below. For further information see the original paper of Huynh [2].

Preliminaries
Consider solving the following 1D scalar conservation law within an arbitrary domain Ω, where x is a spatial coordinate, t is time, u = u(x, t) is a conserved scalar quantity and f = f (u) is the flux of u in the x direction. Additionally, consider partitioning Ω into N distinct elements, each denoted Ω n = {x|x n < x < x n+1 }, such that The FR approach requires u is approximated in each Ω n by a function u δ n = u δ n (x, t), which is a polynomial of degree k within Ω n , and identically zero elsewhere. Additionally, the FR approach requires f is approximated in each Ω n by a function f δ n = f δ n (x, t), which is a polynomial of degree k + 1 within Ω n , and identically zero elsewhere. Consequently, when employing the FR approach, a total approximate solution u δ = u δ (x, t) and a total approximate flux f δ = f δ (x, t) can be defined within Ω as where no level of inter-element continuity in u δ is explicitly enforced. However, f δ is required to be C 0 continuous at element interfaces. Note the requirement that each f δ n is one degree higher than each u δ n , which consequently ensures the divergence of f δ n is of the same degree as u δ n within Ω n .

Implementation
From an implementation perspective, it is advantageous to transform each Ω n to a standard element Ω S = {x|−1 ≤ x ≤ 1} via the mappinĝ which has the inverse Having performed such a transformation, the evolution of u δ n within any individual Ω n (and thus the evolution of u δ within Ω) can be determined by solving the following transformed equation within the standard element Ω S ∂û δ ∂t + is a polynomial of degree k, is a polynomial of degree k + 1, and J n = (x n+1 − x n )/2. The FR approach to solving Eq. (2.6) within the standard element Ω S can be described in five stages. The first stage involves representingû δ in terms of a nodal basis as follows: where l i are Lagrange polynomials defined as 10) x i (i = 0 to k) are k + 1 distinct solution points within Ω S , andû δ i =û δ i (t) (i = 0 to k) are values ofû δ at the solution pointsx i .
The second stage of the FR approach involves constructing a degree k polynomialf δ D =f δ D (x, t), defined as the approximate transformed discontinuous flux within Ω S . Specifically,f δ D is obtained via a collocation projection at the k + 1 solution points, and can hence be expressed aŝ are simply values of the transformed flux at each solution pointx i (evaluated directly from the approximate solution). The fluxf δ D is termed discontinuous since it is calculated directly from the approximate solution, which is in general discontinuous between elements.
The third stage of the FR approach involves evaluating the approximate solution at either end of the standard element Ω S (i.e. atx = ±1). These values, in conjunction with analogous information from adjoining elements, are then used to calculate numerical interface fluxes. The exact methodology for calculating such numerical interface fluxes will depend on the nature of the equations being solved. For example, when solving the Euler equations one may use a Roe type approximate Riemann solver [9], or any other two-point flux formula that provides for an upwind bias. In what follows the numerical interface fluxes associated with the left and right hand ends of Ω S (and transformed appropriately for use in Ω S ) will be denotedf δ I L andf δ I R respectively. The penultimate stage of the FR approach involves constructing the degree k + 1 polynomialf δ , by adding a correction fluxf δC =f δC (x, t) of degree k + 1 tof δ D , such that their sum equals the transformed numerical interface flux atx = ±1, yet in some sense followsf δ D within the interior of Ω S . In order to definef δC such that it satisfies the above requirements, consider first defining degree k +1 correction functions g L = g L (x) and g R = g R (x) to approximate zero (in some sense) within Ω S , as well as satisfying g L (−1) = 1, g L (1) = 0, (2.12) g R (−1) = 0, g R (1) = 1, (2.13) and g L (x) = g R (−x). (2.14) A suitable expression forf δC can now be written in terms of g L and g R aŝ Using this expression, the degree k + 1 approximate transformed total fluxf δ within Ω S can be constructed from the discontinuous and correction fluxes as follows: The final stage of the FR approach involves evaluating the divergence off δ at each solution pointx i using the expression These values can then be used to advanceû δ in time via a suitable temporal discretisation of the following semidiscrete expression

Comments
The nature of a particular FR scheme depends on three factors, namely the location of the solution pointsx i , the methodology for calculating the interface fluxesf δ I L andf δ I R , and the form of the correction functions g L and g R . Huynh [2] showed previously that a collocation based nodal DG scheme is recovered in 1D if the correction functions g L and g R are the right and left Radau polynomials respectively. Also, Huynh [2] showed that SD type methods can be recovered (at least for a linear flux function) if the correction functions g L and g R are set to zero at a set of k points within Ω S (located symmetrically about the origin).
Several additional forms of g L and g R were also suggested by Huynh [2], leading to the development of new schemes with various stability and accuracy properties. Building on this work, and the study of Jameson [8], Vincent, Castonguay and Jameson recently identified a one-parameter family of VCJH correction functions that lead to stable FR schemes for 1D linear advection problems [7].

Preliminaries
If the f (u) = au where a is a constant scalar (i.e. if the flux is linear), then Eq. (2.18) can be written as Eq. (3.1) can be written in matrix form as follows: and On defining a Vandermonde matrix V as where L j (x) is a Legendre polynomial of degree j (normalised to unity atx = 1), one can multiplying through Eq. (3.2) from the left by V −1 to obtain and thus which can be written as are vectors of modal Legendre expansion coefficient for the solution, left correction function derivative, and right correction function derivative respectively, and is the modal Legendre differentiation matrix.

Stability
Theorem 1. For all k, 1D FR correction functions are stable for a linear flux if whereM is the modal Legendre mass matrix defined as Q is a real square matrix of dimension k + 1 that satisfies andL andR are defined as Proof. On multiplying Eq. (3.9) from the left byũ δT (M +Q) one obtainŝ Eq. (3.15) implies thatM +Q is symmetric, hence (3.19) can be written as Eq. (3.16) implies thatQD is anti-symmetric and hencẽ and Eq. (3.17) implies that Hence, Eq. (3.20) can be written as which using the fact that can be written as and hence On rewriting Eq. (3.26) in terms of physical space quantities from the nth element one obtains and hence d dt where u δ n is a vector of the physical solution at the solution points inside the nth element, and f δ I n and f δ I n+1 are physical numerical interface fluxes evaluated at x n and x n+1 respectively. If the numerical flux at each internal interface x n (1 ≤ n ≤ N − 1) is defined to have the form where 0 ≤ κ ≤ 1 (with κ = 0 recovering a fully upwind scheme, and κ = 1 recovering a central scheme), and if for simplicity the domain Ω is assumed to be periodic such that then summing Eq. (3.28) over all elements leads to hence the broken norm ∥u δ ∥ will remain bounded, and hence all norms of the solution will remain bounded via equivalence of norms in a finite dimensional space. whereJ

Symmetry
on multiplying from the right byL one obtains and hencẽ where, as defined,gx L [0] andgx R [0] are the zero mode coefficients in a Legendre expansion of the left and right correction function derivatives respectively.
Proof. Using the orthogonality of Legendre polynomials, Eq. (3.42) implies Hence the schemes will be conservative.

Summary
For all k, correction functions defined by Eqs. (3.12) and (3.13) will be stable, symmetric and conservative provided the conditions defined by Eqs.

Recovery of Vincent-Castonguay-Jameson-Huynh schemes
If q 1 = 0, then Eq. (4.3) collapses to and Eq. (4.6) collapses to These define VCJH correction functions for k = 3, parameterised by q 0 , with q 0 = 0 recovering a DG scheme, q 0 = 3/14 recovering the energy-stable SD scheme described by Jameson [8], and q 0 = 8/21 recovering the g 2 described by Huynh [2]. We note that this description of a VCJH correction function, in terms of a modal Legendre expansion of its derivative, is similar to that presented previously by Huynh [10].

Numerical experiments
Numerical experiments were undertaken to demonstrate that, when k = 3, correction functions defined in terms of q 0 and q 1 via Eqs.  The shaded grey area bounded by a black dashed line highlights the theoretically stable region in q 0 − q 1 space. Solid circles denote schemes that were found to be numerically stable. Hollow circles denote schemes that were found to be numerically unstable.
were used as solution points within each element. Periodic boundary conditions were applied at the ends of Ω, and the following Gaussian profile was prescribed within Ω at t = 0 u(x, 0) = e −20x 2 . (4.10) Time integration was performed using an explicit low-storage five-stage fourth-order Runge-Kutta method [11]. A scheme was deemed to be numerically unstable if the solution at any solution point attained a value of 1000 or greater before t = 300. Otherwise the scheme was deemed to be numerically stable. A plot illustrating which of the schemes were found to be numerically unstable, and which were found to be numerically stable, is shown in Fig. 1. Results of the numerical experiments are in agreement with the theoretical results of Section 4.1, since all schemes within the theoretically stable region of q 0 −q 1 space, defined by Eq. (4.6), were found to be numerically stable. Additionally, it can be seen that all schemes outside of the theoretically stable region were found to be numerically unstable.
where the specific forms of f (q 1 , q 2 ), g(q 0 , q 1 , q 2 ) and h(q 1 , q 2 ) are omitted for brevity. In order for Eq. (5.3) to satisfy the conservation conditions defined by Eq. (3.42) it is required that q 2 f (q 1 , q 2 ) = 0, which requires q 2 = 0 so that the denominator in (5.4) is non-zero. Substituting Eq. (5.2) into Eqs. (3.12) and (3.13) with q 2 = 0 leads to  Fig. 6. Also, for reference, the differential form of the norm defined by Eq. (3.32) when k = 4 is given in Appendix B, B.2.

Recovery of Vincent-Castonguay-Jameson-Huynh schemes
If q 1 = 0, then Eq. (5.5) collapses to gx L [4] = − 9 9 q 0 + 2 , (5.9) Fig. 2. Plot comparing theoretical and numerical results when k = 4. The shaded grey area bounded by a black dashed line highlights the theoretically stable region in q 0 − q 1 space. Solid circles denote schemes that were found to be numerically stable. Hollow circles denote schemes that were found to be numerically unstable. and Eq. (5.8) collapses to These define VCJH correction functions for k = 4, parameterised by q 0 , with q 0 = 0 recovering a DG scheme, q 0 = 8/45 recovering the energy-stable SD scheme described by Jameson [8], and q 0 = 5/18 recovering the g 2 described by Huynh [2]. We note that this description of a VCJH correction function, in terms of a modal Legendre expansion of its derivative, is similar to that presented previously by Huynh [10].

Numerical experiments
Numerical experiments were undertaken to demonstrate that, when k = 4, correction functions defined in terms of q 0 and q 1 via Eqs. (5.5) and (5.6) result in stable FR schemes if q 0 and q 1 satisfy the constraints defined by Eq. (5.8).
Specifically, an equispaced sampling of 451 different schemes within a region of q 0 − q 1 space bounded by −1 ≤ q 0 ≤ 4 and −1 ≤ q 1 ≤ 1 were used to solve Eq. (2.1) with the flux function defined by Eq. (4.9). For each of the 451 numerical experiments the setup was identical to that described in Section 4.3. A scheme was deemed to be numerically unstable if the solution at any solution point attained a value of 1000 or greater before t = 300. Otherwise the scheme was deemed to be numerically stable. A plot illustrating which of the schemes were found to be numerically unstable, and which were found to be numerically stable, is shown in Fig. 2. Results of the numerical experiments are in agreement with the theoretical results of Section 5.1, since all schemes within the theoretically stable region of q 0 − q 1 space, defined by Eq. (5.8), were found to be numerically stable. Additionally, it can be seen that all schemes outside of the theoretically stable region were found to be numerically unstable.

Derivation
For reference, when k = 5 Substituting Eq. (6.2) into Eqs. (3.12) and (3.13) leads to and (c) q 2 = 0.2. Fig. 3. Plots comparing theoretical and numerical results when k = 5. The shaded grey areas bounded by black dashed lines highlight the theoretically stable regions of q 0 − q 1 space with fixed q 2 = −0.2 (a), q 2 = 0.0 (b), and q 2 = 0.2 (c). Solid circles denote schemes that were found to be numerically stable. Hollow circles denote schemes that were found to be numerically unstable.
In summary, when k = 5, correction functions defined in terms of q 0 , q 1 , and q 2 via Eqs. (6.3) and (6.4) will result in stable FR schemes if q 0 , q 1 , and q 2 satisfy the constraints defined by Eq. (6.6). Examples of such correction functions are shown in Appendix A, Figs. 7-9. Also, for reference, the differential form of the norm defined by Eq. (c) q 2 = 0.2. Fig. 4. Plots comparing theoretical and numerical results when k = 6. The shaded grey areas bounded by black dashed lines highlight the theoretically stable regions of q 0 − q 1 space with fixed q 2 = −0.2 (a), q 2 = 0.0 (b), and q 2 = 0.2 (c). Solid circles denote schemes that were found to be numerically stable. Hollow circles denote schemes that were found to be numerically unstable.

Recovery of Vincent-Castonguay-Jameson-Huynh schemes
If q 1 = q 2 = 0, then Eq. (6.3) collapses to and Eq. (6.6) collapses to These define VCJH correction functions for k = 5, parameterised by q 0 , with q 0 = 0 recovering a DG scheme, q 0 = 5/33 recovering the energy-stable SD scheme described by Jameson [8], and q 0 = 12/55 recovering the g 2 described by Huynh [2]. We note that this description of a VCJH correction function, in terms of a modal Legendre expansion of its derivative, is similar to that presented previously by Huynh [10].

Numerical experiments
Numerical experiments were undertaken to demonstrate that, when k = 5, correction functions defined in terms of q 0 , q 1 and q 2 via Eqs. (6.3) and (6.4) result in stable FR schemes if q 0 , q 1 and q 2 satisfy the constraints defined by Eq. (6.6).
Specifically, an equispaced sampling of 1353 different schemes within a region of q 0 − q 1 − q 2 space bounded by −1 ≤ q 0 ≤ 4, −1 ≤ q 1 ≤ 1, and −0.2 ≤ q 2 ≤ 0.2 were used to solve Eq. (2.1) with the flux function defined by Eq. (4.9). For each of the 1353 numerical experiments the setup was identical to that described in Section 4.3. A scheme was deemed to be numerically unstable if the solution at any solution point attained a value of 1000 or greater before t = 300. Otherwise the scheme was deemed to be numerically stable. Plots illustrating which of the schemes were found to be numerically unstable, and which were found to be numerically stable, are shown in Fig. 3. Results of the numerical experiments are in agreement with the theoretical results of Section 6.1, since all schemes within the theoretically stable region of q 0 −q 1 −q 2 space, defined by Eq. (6.6), were found to be numerically stable. Additionally, it can be seen that all schemes outside of the theoretically stable region were found to be numerically unstable.

Derivation
For reference, when k = 6  By inspection, the most general form ofQ that simultaneously satisfies the stability conditions defined by Eqs. (3.15) and (3.16), and the symmetry condition defined by Eq. (3.34), is Substituting Eq. (7.2) into Eqs. (3.12) and (3.13) leads to where the specific forms of f (q 1 , q 2 , q 3 ), g(q 0 , q 1 , q 2 , q 3 ) and h(q 1 , q 2 , q 3 ) are omitted for brevity. In order for Eq. (7.3) to satisfy the conservation conditions defined by Eq. (3.42) it is required that q 3 f (q 1 , q 2 , q 3 ) = 0, which requires q 3 = 0 so that the denominator in (7.4) is non-zero. Substituting Eq. (7.2) into Eqs.  These define VCJH correction functions for k = 6, parameterised by q 0 , with q 0 = 0 recovering a DG scheme, q 0 = 12/91 recovering the energy-stable SD scheme described by Jameson [8], and q 0 = 7/39 recovering the g 2 described by Huynh [2]. We note that this description of a VCJH correction function, in terms of a modal Legendre expansion of its derivative, is similar to that presented previously by Huynh [10].

Numerical experiments
Numerical experiments were undertaken to demonstrate that, when k = 6, correction functions defined in terms of q 0 , q 1 and q 2 via Eqs. (7.5) and (7.6) result in stable FR schemes if q 0 , q 1 and q 2 satisfy the constraints defined by Eq. (7.8).
Specifically, an equispaced sampling of 1353 different schemes within a region of q 0 − q 1 − q 2 space bounded by −1 ≤ q 0 ≤ 4, −1 ≤ q 1 ≤ 1, and −0.2 ≤ q 2 ≤ 0.2 were used to solve Eq. (2.1) with the flux function defined by Eq. (4.9). For each of the 1353 numerical experiments the setup was identical to that described in Section 4.3. A scheme was deemed to be numerically unstable if the solution at any solution point attained a value of 1000 or greater before t = 300. Otherwise the scheme was deemed to be numerically stable. Plots illustrating which of the schemes were found to be numerically unstable, and which were found to be numerically stable, are shown in Fig. 4. Results of the numerical experiments are in agreement with the theoretical results of Section 7.1, since all schemes within the theoretically stable region of q 0 −q 1 −q 2 space, defined by Eq. (7.8), were found to be numerically stable. Additionally, it can be seen that all schemes outside of the theoretically stable region were found to be numerically unstable.

Conclusions
Building on the work of Huynh [2] and Jameson [8], Vincent, Castonguay and Jameson recently identified a one-parameter family of VCJH correction functions that lead to stable FR schemes for 1D linear advection problems [7]. In this study we developed a procedure for identifying an extended range of stable, symmetric, and conservative FR correction functions. The procedure was applied to identify ranges of such correction functions for various orders of accuracy. In all cases the original one-parameter VCJH correction functions were found to be a subset of the extended ranges. Numerical experiments were undertaken, and the results found to be in agreement with the theoretical findings. Future studies should extend the approach presented here to simplex elements in 2D and 3D.         When k = 6,M +Q can be decomposed as where ϵ 1 = ϵ 2 = ϵ 3 = q 2 22050 , ϵ 4 = ϵ 5 = q 1 − 9q 2 1786050 , ϵ 6 = q 0 − 11q 1 + 54q 2 216112050 . (B.11) Hence, ∥u δ ∥ can be written in differential form as x n (u δ n ) 2 +