An Analysis of Solution Point Coordinates for Flux Reconstruction Schemes on Tetrahedral Elements

The flux reconstruction (FR) approach offers an efficient route to high-order accuracy on unstructured grids. In this work we study the effect of solution point placement on the stability and accuracy of FR schemes on tetrahedral grids. To accomplish this we generate a large number of solution point candidates that satisfy various criteria at polynomial orders ℘=3,4,5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\wp = 3,4,5$$\end{document}. We then proceed to assess their properties by using them to solve the non-linear Euler equations on both structured and unstructured meshes. The results demonstrate that the location of the solution points is important in terms of both the stability and accuracy. Across a range of cases it is possible to outperform the solution points of Shunn and Ham for specific problems. However, there appears to be a degree of problem-dependence with regards to the optimal point set, and hence overall it is concluded that the Shunn and Ham points offer a good compromise in terms of practical utility.


Introduction
Theoretical studies and numerical experiments suggest that unstructured high-order methods can provide solutions to otherwise intractable fluid flow problems within the vicinity of complex geometries. In 2007 Huynh proposed the flux reconstruction (FR) approach [1], a unifying framework for high-order schemes that encompasses several existing methods while simultaneously admitting an efficient implementation. Using FR it is possible to recover nodal discontinuous Galerkin (DG) methods of the type described by Hesthaven and Warburton [2], and the spectral difference schemes, original proposed by Kopriva and Kolias [3] and later popularised by Sun et al. [4]. Unlike traditional DG methods, which are based on a weak formulation, and hence require integration, the FR approach is based on the differential form of the governing system. As a consequence, implementations of FR forego having to perform numerical quadrature within each element. This not only reduces complexity, but can also lead to decreased computational cost.
One problem with the FR approach is that for a non-linear flux function aliasing driven instabilities may develop. The severity of these instabilities depends upon the degree to which the flux is under-resolved within each element. For large scale-resolving flow simulations it is often the case that features are under-resolved due to resource constraints or modelling approaches. For example when undertaking an implicit large-eddy simulation of a turbulent flow the numerics do not resolve all of the physical scales of the problem. Hence, there is a degree of underesolution. Given that there is an increasing desire amongst practitioners to undertake such scale-resolving simulations in an industrial context it is important to address these aliasing phenomena effectively.
It has been demonstrated, both theoretically [5] and empirically [6][7][8] that the degree of aliasing driven instabilities depends upon the location of the solution points inside each element. Specifically, it has been found that placing points at the abscissa of strong Gaussian quadrature rules has a positive impact on their performance.
In this paper we will assess the performance of a variety of symmetric quadrature rules when used as solution points inside of tetrahedra. In Sect. 2 we outline the requirements that a set of points must satisfy in order to be used as solution points. Empirical and theoretical motivations for wanting solution points to be at the abscissa of quadrature rules is also presented. Our candidate point sets at polynomial orders ℘ = 3, 4, 5 are presented in Sect. 3. Two numerical test cases, both based around solving the non-linear Euler equations, are introduced in Sect. 4. These test cases are then used to evaluate the performance of the candidate points with the results being presented in Sect. 5. Numerical stability across the two test cases is assessed in Sect. 6. Finally, in Sect. 7, conclusions are drawn.

Solution Point Requirements
Hesthaven and Warburton [2] identify three key criteria that solution point sets for nodal DG must satisfy. The first is that the points define a well conditioned nodal basis. This property ensures any nodal expansions based on the points will be well suited to the task of polynomial interpolation. The second is that the points be arranged symmetrically inside of each element. This eliminates any potential for the solution to be biased towards one region of an element depending on the global-to-local mapping. The third requirement is that the number of solution points must be equivalent to the rank of the polynomial basis used to represent the solution. For tetrahedra this fixes the number of points to be a tetrahedral number given by where ℘ is the polynomial order. In FR and nodal DG solution points are not just used for representing the discontinuous solution within each element; they are also used to construct a polynomial representation of the flux. However, should the flux function be non-linear then the construction of this polynomial gives rise to an implicit collocation projection from the space of the flux, which may be non-polynomial, into the space of the solution. The development of aliasing driven instabilities is closely related to how this projected polynomial compares to the L 2 optimal polynomial obtained via a full L 2 projection of the flux [5]. The differences between these two polynomials in both FR and nodal DG can be significant [2] and hence have a significant impact on the non-linear stability of these schemes.
Hence, an additional fourth requirement for solution points is that they also minimise the L 2 error arising from performing a collocation projection of non-linear functionals.

Polynomial Interpolation and Unisolvency
Consider an arbitrary scalar function u(r) on the reference tetrahedron T as shown in Fig. 1. By expanding this function in terms of a nodal basis set we can construct an interpolating polynomial where u δ (r) ≈ u(r), {r i } are a set of solution points, and { i (r)} is a nodal basis set with the property that i (r j ) = δ i j with δ i j being the Kronecker delta. Starting with an orthonormal polynomial basis set {ψ i (r)} on T such that T ψ i (r)ψ j (r) d 3 r = δ i j we can construct a nodal basis set by first computing the elements of the generalised Vandermonde matrix and then taking as required. We note here that in constructing the nodal basis set we have required the generalised Vandermonde matrix to be invertible, which is the case if det V = 0. A necessary, but not sufficient, condition for this is that all of the points must be distinct. However, it is known that in two and three dimensions only certain sets of distinct points give rise to a invertible Vandermonde matrix; a requirement which is termed unisolvency [9,10]. To illustrate this take {r i } to be a set of distinct points with a non-singular Vandermonde matrix and consider arbitrarily relabelling a pair of points. The effect of this relabelling is to interchange two columns in the Vandermonde matrix. From the properties of the determinant Fig. 2 Origins of non-unisolvency in multiple dimensions as illustrated on a triangle this will cause the sign of the determinant to flip. If this interchange is performed continuously, with the two points following different non-intersecting paths as shown in Fig. 2, it is evident from the intermediate value theorem that there must be an intermediate arrangement where the determinant is zero. Hence, while the points are all distinct they can not be used to construct a nodal basis set.

Quadrature Rules
In 2011 Castonguay et al. [6] used the FR approach to solve the Euler equations on triangular grids. As a starting point they employed the α-optimised points of Hesthaven and Warburton [2] as the set of solution points. These points are constructed with the goal of minimising the Lebesgue constant which serves as a quantification of how good a set of points are for the purposes of polynomial interpolation. However, while evaluating the performance of these points on a isentropic Euler vortex test problem the reference point set poor performance, in terms of both accuracy and stability, was observed. Noting this, and motivated by the results of Jameson et al. [5], the authors proceeded to use the abscissa of the (mildly asymmetric) quadrature rules presented in [11] as the solution points. This change gave rise to a marked improvement in both stability and accuracy. These results were subsequently confirmed by Witherden and Vincent [8] in 2014 who analysed the performance triangular quadrature rules when used as solution points. Inside of tetrahedra Williams [7] compared the three dimensional α-optimised points against the tetrahedral quadrature rules of Shunn and Ham [12]. Once again the quadrature points were found to greatly outperform the α-optimised points with regards to both stability and accuracy. A visual comparison of the third order α-optimised and Shunn-Ham points can be seen in Fig. 3 These results contribute to a strong body of empirical evidence suggesting that solution points should be taken to be the abscissa of quadrature rules. In the context of numerical integration quadrature rules are defined based on their ability to integrate polynomial functions up to a given degree inside of a domain T . Specifically, an N p point quadrature rule of strength ϕ is capable of exactly integrating any polynomial function p(x) ∈ T of maximal degree ϕ as

. A major difference between the
where {x i } are the abscissa of the rule and {ω i } the associated quadrature weights. Generating such rules however, especially those with a prescribed number of points, is non-trivial.

Candidate Point Sets
A popular approach for identifying quadrature rules inside of two-and three-dimensional domains is based around expressing Eq. 6 as a system of non-linear equations. Given a suitable set of polynomials { p i (x)} ∈ T it is possible to analytically evaluate the integrals on the left hand side of Eq. 6. This can then be approached as a least squares optimisation problem which is non-linear with respect to the abscissa and linear with respect to the weights. By seeding the non-linear least squares problem with a random initial condition it is capable of identifying a large number of distinct quadrature rules. Polyquad [13] is a software package which employs the Levenberg-Marquardt algorithm to identify fully symmetric quadrature rules of strength ϕ which a prescribed number of points inside of a variety of domains, including tetrahedra. Using Polyquad v1.0.0 a large number of symmetric quadrature rules with a tetrahedral number of points were generated. At each order Polyquad was run for 20 min on an Intel Core i7-4820K CPU. By design Polyquad ensures that for all points x i ∈ T . Although Polyquad does perform some degree of rule deduplication an additional post-processing step was added which required that no two rules had the same set of quadrature weights given a tolerance of ± 1 %. Finally, the remaining unique rules were checked for unisolvency by computing the determinant of the Vandermonde matrix. A summary of the final set of rules can be seen in Table 1.
We note here that our candidate point sets have the same quadrature strengths as those of Shunn and Ham [12].

Euler Equations
The three dimensional Euler equations govern the flow of an inviscid compressible fluid. They are both time-dependent and non-linear. Expressed in conservative form they read where Here ρ is the mass density of the fluid, v = (v x , v y , v z ) T is the fluid velocity vector, E is the total energy per unit volume, and p is the pressure. For a perfect gas the pressure and total energy are related by with γ 1.4. To solve such a system using the FR approach it is necessary to chose both a time integration scheme and an approximate Riemann solver to evaluate the normal fluxes at element interfaces.

Test Cases
For the purposes of evaluating the performance of our point sets two test cases with analytical solutions are considered. Manufactured Sinusoidal Solution Through the introduction of a source term it is possible for the non-linear Euler equations to admit a sinusoidal solution. Specifically, by taking where θ = k(x + y + z) − ωt with k = π, ω = π and a = 3, and adding this to the right hand side of Eq. 7 leads to a system that admits solutions of the form For this test case the computational domain is taken to be = [−1, 1] 3 with fully periodic boundary conditions. Structured tetrahedral meshes were generated on this domain by taking a structured hexahedral mesh with N 3 elements and splitting each hexahedron into six tetrahedra for a total element count of N E = 6N 3 . For the purposes of our experiments six meshes with N = 4, 5, 6, 8, 12, 16 were employed. A cutaway of the N = 5 mesh can be seen in Fig. 4. With these we can introduce a characteristic spacing for mesh N at polynomial order ℘ as where | | is the volume of the domain.
The meshes employed at each order are listed in Table 2. As discussed previously, here we are most interested in the performance of point sets for under-resolved simulations, which are often encountered in practice. The coarsest grids utilised for each ℘ were chosen to achieve absolute errors in the range of 10 −1 -10 0 . This is comparable to the magnitude of the solution, see Eq. 11, and thus can be considered significantly under-resolved. For each ℘ three further successively finer grids were used to give a total of four grids per ℘ in total. Isentropic Euler Vortex Following [6,7,14] we also consider the propagation of an isentropic Euler vortex in a free-stream. The initial conditions for this numerical experiment, in primitive form, are given by p(x, y, z, where f = (1 − x 2 − y 2 )/2R 2 , S is the strength of the vortex, M is the free-stream Mach number and R is the radius of the vortex. These conditions give rise to a vortex that is  [14]. The impact of this is mitigated by the observation that the exponentially-decaying vortex has a radius which is far smaller than the extent of the domain. Neglecting these effects the analytic solution of the system at a time t is simply a translation of the initial conditions. A series of four unstructured tetrahedral meshes were generated with N E = 837, 2165, 5076, 9571 elements, respectively. A cutaway of the coarsest N E = 837 mesh can be seen in Fig. 5. All four grids were employed at each of the three polynomial orders.

Error Estimation
Given that the analytic solution is available for both test cases an L 2 error can be defined as where Q δ (x, t) is a numerical quantity, Q(x, t) is the associated analytical quantity, and M i (x) is a local-to-global coordinate mapping inside of element i with J i (x) being the associated Jacobian determinant. In the third step each integral has been approximated by using a quadrature rule with abscissa {x j } and weights {w j }. So long as the rule used-which need not be symmetric or consist of a tetrahedral number of points-is of adequate strength then this will be a very accurate approximation of the true L 2 error. For the purposes of our experiments a 59-point rule of quadrature strength 9 was employed. To assess the error for the manufactured sinusoidal solution test case the quantity was taken to be the energy E(x, t) whereas for the isentropic Euler vortex test case the density ρ(x, t) was chosen. Denote the set of errors for rule i at polynomial order ℘ as {σ (℘) i j } where j runs over the four meshes employed at this polynomial order for the test case in question. Point sets where one or more of the simulations diverged are considered to have an infinite error. Taking the geometric mean of these numbers we find which gives us an overall error we can rank. The use of the geometric, as opposed to the arithmetic, mean allows us to account for the fact that magnitude of the absolute error varies greatly across the four grids. A key property of the geometric mean is that which permits us to examine the overall relative performance of our points.

Other Details
To completely define the proposed numerical experiment it is also necessary to specify the time marching algorithm, the correction functions, the approximate Riemann solver, and the choice of flux points along each edge. In this study a TVD-RK3 scheme of Gottlieb [15] is chosen for time marching. When performing the simulations it is important that the time step t is chosen such that the temporal error is negligible when compared to the spatial error. Correction functions were chosen to recover the nodal DG scheme described in Hesthaven and Warburton [2]. For computing the numerical fluxes at element interfaces a Rusanov type Riemann solver, as presented in Sun et al. [4], is employed. Finally, on the faces of the tetrahedra the flux points are taken to be the Williams Shunn points of [16].

Results
For each polynomial order all derived point sets were used to solve the test problem define in Sect. 4 with a modified version of the PyFR solver v0.2.3 [17]. Differences between the reference version of PyFR v0.2.3 and the one employed here can be found in the electronic supplementary material.

Manufactured Sinusoidal Solution
All simulations were initialised with the analytic solution at t = 0 and run up until t = 100 with a t = 5 × 10 −4 at which point L 2 errors of energy were computed in accordance with Eq. 18. At each order the top 18 rules were compared and contrasted against the Shunn-Ham When the order is increased to ℘ = 4 the results, shown in Fig. 7, are more varied. Although little improvement is observed on the finer meshes where N = 12, 8, 6, the majority of our points can be seen to outperform the Shunn-Ham points on the coarsest N = 5 mesh. The top ranked point set here has an overall error which is 90.5 % that of the reference set. The results at ℘ = 5 can be seen in Fig. 8. Improvements are observed on all four grids with the top-ranked point set having a relative error which is 64.4 % that of the Shunn-Ham points. As in the case of ℘ = 3 the top-ranked point set is found to outperform the Shunn-Ham points on all of the grids.

Isentropic Euler vortex
All simulations were initialised with the analytic solution at t = 0 and run up until t = 20 with a t = 3 × 10 −4 at which point L 2 errors of density were computed in accordance with Eq. 18. As with the manufactured sinusoidal solution test case at each order the top 18 rules were compared and contrasted against the Shunn-Ham points. Plots showing the absolute and relative errors, as a function of the characteristic mesh spacing, can be seen in Figs. 9, 10, and 11. Looking at the plots it is evident that when ℘ = 3 all of our new point sets represent an improvement over the reference set of Shunn-Ham in all but the most under-resolved grid. The top-ranked point set has an overall error which is 78.1 % of the Shunn-Ham points. When the order is increased to ℘ = 4 the results, shown in Fig. 7, are more varied. The majority of rules can be seen to under-perform the Shunn-Ham points with the best overall rule having an overall error which is 1.01 times that of the Shunn-Ham points. The results at ℘ = 5 can be seen in Fig 8. Here we note that only five of the generated point sets successfully ran to completion on all four of the grids. However, of these five point sets the majority are observed to outperform the Shunn-Ham points, especially on the coarser grids. The topranked point set is found to have an overall error which is 68.1 % that of the Shunn-Ham points.

Comparisons
A cross comparison between the top-ranked point sets for each of the two test cases can be seen in Table 3. For both ℘ = 3 and ℘ = 4 the top-ranked point set is different between the two test cases. However, at ℘ = 3 it is observed that the top-ranked Euler vortex point set represents an improvement over the Shunn-Ham points on both of the test cases. In the case of ℘ = 5 the top-ranked set is identical for both cases with a similar overall reduction in error being observed for both of the problems.
At each polynomial order the top ranked point set is provided as electronic supplementary material.

Stability
In order to assess the stability of various point sets it is interesting to tabulate the number of rules that diverged in either one, or both, of the two test cases. This can be seen in Table 4.

Conclusion
In this paper we have investigated how solution point placement contributes towards the accuracy and stability of FR schemes on tetrahedral elements. The performance of these points has been assessed with two test cases at three different polynomial orders on both structured and unstructured grids. At all three orders of accuracy the top-ranked point set is observed to differ between the two test cases. This suggests that, for a given polynomial order, there is no set of solution points that are universally optimal across both test cases.
For the manufactured sinusoidal solution test case when ℘ = 3, 5 we have identified point sets which outperforms the equivalent Shunn-Ham points on all four grids. At order ℘ = 4 a point set has been identified which outperform the Shunn-Ham points on three out of the four grids. For the isentropic Euler vortex case when ℘ = 5 the top-ranked point set outperforms the equivalent Shunn-Ham points on all four grids. At order ℘ = 3 the topranked rule outperforms the Shunn-Ham points on three out of the four grids for an overall improvement across all four grids. At order ℘ = 4 the top-ranked point set only outperforms the Shunn-Ham points on two out of the four grids with no overall improvement.
From this we conclude that while it is indeed possible to outperform the solution points of Shunn-Ham it is not possible to do so consistently. As such the Shunn-Ham points should be considered a good compromise in terms of practical utility when using FR schemes on both structured and unstructured tetrahedral grids.

Data Statement
Data relating to the results in this manuscript can be downloaded as Electronic Supplementary Material under a CC-BY-NC-ND 4.0 license.