Dealiasing techniques for high-order spectral element methods on regular and irregular

High-order methods are becoming increasingly attractive in both academia and industry, especially in the context of computational ﬂuid dynamics. However, before they can be more widely adopted, issues such as lack of robustness in terms of numerical stability need to be addressed, particularly when treating industrial-type problems where challenging geometries and a wide range of physical scales, typically due to high Reynolds numbers, need to be taken into account. One source of instability is aliasing eﬀects which arise from the nonlinearity of the underlying problem. In this work we detail two dealiasing strategies based on the concept of consistent integration. The ﬁrst uses a localised approach, which is useful when the nonlinearities only arise in parts of the problem. The second is based on the more traditional approach of using a higher quadrature. The main goal of both dealiasing techniques is to improve the robustness of high order spectral element methods, thereby reducing aliasing-driven instabilities. We demonstrate how these two strategies can be eﬀectively applied to both continuous and discontinuous discretisations, where, in the latter, both volumetric and interface approximations must be considered. We show the key features of each dealiasing technique applied to the scalar conservation law with numerical examples and we highlight the main diﬀer-ences in terms of implementation between continuous and discontinuous spatial discretisations.


Introduction
Numerical stability is a fundamental requirement for tools, such as solvers used to model fluid dynamics problems.This is especially true when these tools are used in an engineering context, where reliability and robustness are key factors in its efficacy in the industrial pipeline.While low order finite element methods, which are traditionally popular in such tools, demonstrate a usually satisfying level of numerical stability due to their dissipation properties, high order spectral element methods typically exhibit low dissipation errors, and are therefore affected by a lack of stability which in turn affects their robustness and reliability.
Over the last decade, various attempts have been made to address this problem.Different formulations of the nonlinear terms, for instance, have different aliasing properties.In [1] it is shown that the skew-symmetric form of the convective term results in a reduced amplitude of the aliasing errors in comparison to the conservative and nonconservative forms.Examples of the skewsymmetric methodology for discontinuous Galerkin spectral element methods can be found in [2], [3].Approaches such as spectral vanishing viscosity [4] modify the underlying discretised operators in order to add artificial dissipation at high polynomial orders, thereby reducing high-frequency oscillations in the approximate solution, stabilising the method and making it suitable for high Reynolds number large-eddy simulations [5,6,7,8].Polynomial filtering [9,10,11,12] adopts a similar strategy, in which an interpolation-based filter is applied locally to each element at each time-step to prevent the buildup of unwanted oscillations in the solution field.An alternative route for stabilisation is represented by the variational multiscale schemes [13] which are based on the model used to approximate the coupling between unresolved and partially resolved scales [14,15,16].Finally, in [17] the impact of an exponential based modal filtering and over-integration on underresolved high-order turbulent flow computations is investigated and in [18] a DG method with guaranteed stability is obtained using a non-polynomial modified basis.
In this paper, we focus on the concept of over-integration [19,20,21,22], where over-integration refers to the use of more quadrature points than are necessary for a linear operator, but which might equally well be considered as consistent integration of the nonlinear operators.In the formulation of spectral element methods, on any given element one may choose the number of expansion modes used to represent the approximation of a function, and the number of quadrature points used to calculate numerical estimates of integrals which may appear in the formulation.
A popular choice is to couple a set of quadrature points with an equal number of nodal Lagrange polynomials defined at the same points, leading to a collocation method.There are many examples of this throughout the literature, both in terms of the more traditionally utilised continuous Galerkin (CG) and discontinuous Galerkin (DG) formulations, as well as newer extensions such as the flux reconstruction (FR) technique as presented by Huynh [23].In collocation methods, while most linear operators can be exactly integrated in this setting depending on the choice of quadrature, integrals of nonlinear terms typically incur numerical error.However, the computational efficiencies that can be attained through the use of a collocation formulation, especially given the presence of a diagonal mass matrix, often outweigh the numerical error that is incurred.
To illustrate this point, let us consider errors arising from the integrals of nonlinear products of polynomial expansions of order P .It is well-known that polynomials of order P can be exactly integrated up to machine precision, given some minimum number of quadrature points Q min [24].For example Gauss-Lobatto-Legendre quadrature, which is widely adopted in the context of spectral element methods, allows one to obtain exact values of integrals for order 2P − 3 polynomials given P quadrature points [25].However, when nonlinear quantities are calculated at the quadrature points, a collocation projection leads to an error known as polynomial aliasing being introduced into the obtained integrals, due to an insufficient number of quadrature points.This numerical error may be negligible for well-resolved simulations, such as flow simulations at low Reynolds numbers, but can become a critical issue for under-and marginallyresolved simulations that typically arises in large eddy simulations.In addition, in the presence of non-polynomial functions, such as the compressible Euler and Navier-Stokes equations fluxes written in conservative form, consistent integration cannot be applied with a specific rule as in the case of polynomial functions.In this case, it is still possible to consistently reduce aliasing errors but the number of quadrature points required might become too computationally expensive.
The aim of this paper is to illustrate effective and computationally efficient approaches based on consistent integration to resolve the aliasing issues arising in spectral/hp element discretisations.Specifically, we take into consideration three different spectral/hp element approximations: the continuous and discontinuous Galerkin methods as well as the flux reconstruction approach.For a scalar partial differential equation (PDE) these approximations can be written as: CG : where u is the approximate solution field, M is the mass matrix, L is the linear term and N is the nonlinear term.In DG and FR both the linear and the nonlinear term are composed by a volumetric term over each spectral element and an interface term on the trace space connecting the elements, which we denote with a subscript V and I, respectively.Henceforth, N includes aliasing issues arising from both nonlinear and "variable-coefficient" fluxes and to curved geometry, both of which might be described by nonlinear polynomial functions.Note that "variable-coefficient" fluxes and curved geometry, being functions of the spatial coordinates x, do not affect the linearity of the underlying PDE which depends on u.However, they introduce aliasing errors.As an example, the advection equation with spatially varying coefficients used in section 4 is still a linear equation, since the advection velocity does not depend on the solution, but the non-constant velocity coefficients introduce aliasing errors.
To address these aliasing issues we consider two different strategies both based on the concept of consistent integration.With the first approach we locally consistently integrate only the nonlinear terms of the PDE being solved without addressing aliasing sources arising from the mapping used in deformed/curved meshes (also referred to as geometrical-aliasing sources).Specifically, with this approach we can address only aliasing arising from the nonlinear or "variable coefficient" fluxes in our problem (also referred to as PDE aliasing).Eq. ( 1) can therefore be re-written as: where Q is the number of quadrature points used for the linear term and the mass matrix, while Q V and Q I indicate the number of quadrature points used for the consistent integration of the volumetric and the interface part of the nonlinear term, respectively.We note the local behaviour of the strategy, which allows us to selectively adjust the quadrature used for the nonlinear term.
In the second approach, we consistently integrate every term of the numerical discretisation which typically implies the over-integration of the linear terms.In this strategy we can address both PDE-aliasing driven instabilities as well as geometrical-aliasing sources arising from deformed/curved elements.Eq. (1) becomes: where the same number of quadrature points is used globally for all the terms.To summarise, the main contributions of this paper are: • Generalised analysis which encompasses DG, FR and CG discretisations, highlighting the contribution of interface and volumetric terms.
• Local dealiasing strategy which exploits sum-factorisation type approach to improve computational efficiency.
The paper is structured as follows.In section 2, we briefly summarise the different CG, DG and FR discretisations, and highlight the sources of aliasing errors which arise in each formulation.Section 3 outlines the Local and Global dealiasing strategies we adopt to reduce the aliasing errors and improve the stability of the different discretisations considered.In this section we also describe the implementation details which are required to apply such methods.In section 4, we present practical applications of the dealiasing techniques, highlighting key examples in compressible and incompressible flow simulations.Finally, in section 5, we outline the results and conclusions of the paper.

Aliasing sources in high-order spectral element methods
Aliasing errors in polynomial spectral element methods typically arise through under-integration of the terms appearing in the weak form of the equations to be discretised.When using Gauss-Lobatto-Legendre (GLL) quadrature points, for instance, the minimum number of quadrature points Q min necessary to exactly integrate any given polynomial of order P , u(ξ) ∈ P P , up to machine precision is In Galerkin methods, we are usually interested in evaluating the L 2 inner product of two polynomials (φ p , φ q ) in order to compute, for example, the mass matrix of our discretised problem, where φ p (ξ), φ q (ξ) ∈ P P are the so-called test or virtual functions, which in our studies are chosen to be the same as the solution expansion.By using Eq. ( 4) it is possible to calculate the minimum number of GLL quadrature points necessary for the quadrature to be exact as a function of the polynomial order P , as shown in Table 1.We note that to exactly integrate a linear problem it is necessary to use Q min = P + 2 GLL quadrature points.For nonlinear problems the number of quadrature points needed increases and for quadratic nonlinearities such as the convective term of the incompressible Navier-Stokes equations it becomes Q min = 3 2 (P + 1) whereas for cubic nonlinearities Q min = 2(P + 1) quadrature points are needed for the numerical integration to be exact.The above discussion holds true as long as the nonlinearities belong to the polynomial approximation space used for representing the discretised problem.If instead we have non-polynomial nonlinearities, it is not possible to perform an exact numerical integration using Gauss-type quadrature combined with the rules in Table 1.However, it is still possible to reduce the integration error to machine precision by choosing an integration quadrature rule with a sufficiently high degree.This occurs in, for instance, the compressible Euler and Navier-Stokes equations.When written in conservative form, the flux functions are rational and, therefore, non-polynomial.
In addition to the nonlinearities present in the equations (PDE aliasing sources), additional aliasing may arise from the use of deformed or curved elements.The Jacobian of the isoparametric mapping, which takes a standard element and projects it into the Cartesian space, is by definition a non-constant, non-linear polynomial for curved elements.This introduces a geometricalaliasing source which may corrupt the solution in a similar manner as the PDEaliasing and may eventually increase the overall degree of the nonlinearity.
In this section we consider three different high-order methods, namely the continuous Galerkin (CG) method, the discontinuous Galerkin (DG) method and the flux reconstruction (FR) approach.For each of these methods we highlight the sources of PDE-and geometrical-aliasing by considering the multidimensional conservation law where u = u(x, t) is the conserved variable and F(u) = [f (u), g(u), h(u)] T is the flux vector that governs the transport of u.
Note that the form of the advection term in Eq. ( 5) can affect the numerical stability.There is a rich literature, especially for the incompressible Navier-Stokes equations, where the convective term is represented in a skewsymmetric form in order to enhance stability despite losing conservation.In this regard, Malm, Schlatter et al. [26] detailed a possible relationship between consistent integration and the skew-symmetric form of the convection operator, pointing out that for any Galerkin-based method the use of consistent integration recovers imaginary eigenvalues (i.e.skew-symmetry of the convective operator) and therefore stability.In this work we do not consider other possible choices for improving stability other than consistent integration.Nevertheless, the connections between consistent integration and the skew-symmetric form of the convective term remain an interesting aspect which needs to be further explored (see also [2,3]).

CG formulation
The CG formulation is based on a continuous discretisation of the computational domain.In particular the approximation space as well as the space of the test functions is defined as where Ω = e=1 Ω e is the mesh representing the computational domain.The weak form of the CG method for Eq. ( 5) is In order to obtain a discrete form for this equation, we write u, F(u) and v in terms of an expansion of basis functions of order P , transforming Eq. ( 7) into its elemental weak form where φ pqr = φ pqr (x) are the basis functions, (p, q, r) ∈ N 3 denotes a tensor product ordering of the modes, Ω e is the elemental domain and (•, We evaluate the solution u at a set of Q quadrature points in each tensor product direction and denote it with the vector u e = B e u e (t), where B e = φ pqr (ξ m ) is the basis matrix, the subscript m refers to the m-th quadrature point and we have used the fact that in the standard CG approximation, both the trial and the test functions lie in the same expansion space.The elemental matrix form of Eq. (8) becomes: where W e = W e (J) is the weight matrix that depends on the elemental Jacobian and D e xi are the elemental differentiation matrices written with respect to the Cartesian coordinates.A more detailed explanation of the geometric definitions and the operations in an arbitrarily-shaped element can be found in Appendix A; the interested reader can also refer to [24].We can express the partial derivatives in x i in terms of partial derivatives in standard region coordinates ξ i and, using the chain rule, we can write: where Λ(•) is a diagonal matrix which contains the factors on its diagonal.Equation (10) The final step in constructing the CG discretisation of the problem in Eq. ( 5) involves assembling the global problem on Ω using the elemental contributions in Eq. (13).By denoting with u l the vector of all the expansion coefficients from each element u e and by relating u l with the global expansion coefficients u g through the assembly operation u e = A u g , we finally obtain the global matrix form for a standard CG approximation where • represents the block diagonal direct summation of all local elemental matrices Equation ( 14) can be discretised in time after the inversion of the global mass matrix M = A T (B e ) T W e B e A: where we indicate with D i the global advection matrices.In the CG formulation, if Eq. ( 5) is nonlinear, the PDE-aliasing arises from the advection terms D 1 f g , D 2 g g and D 3 h g .On the other hand, the geometrical aliasing arises from the inverse of the mass matrix as well as from the advection matrices.These depend on the elemental Jacobian determinants (hereafter referred to as Jacobians) through the weight matrix W e = W e (J).The advection matrices also depend on the geometric factors through W e as well as G i/j .In the case of an unstructured domain with deformed/curved elements, the Jacobians and the geometric factors are represented as non-constant polynomials, which may increase the aliasing issues if not properly resolved.

DG formulation
The DG formulation relies upon a discontinuous discretisation of the computational domain.In particular, the approximation space as well as the space of the test functions is defined as where with R being the space of real numbers.The weak form of the DG method applied to Eq. ( 5) is where the elemental inner products are defined as follows: and f (u δ + , u δ − ) is the interface flux which couples the elements of the DG discretisation.This interface flux can be computed in various ways, including simple upwind approaches or more sophisticated Riemann solvers depending on the equation being solved.It is important to note here that the interface flux is a combination of the information coming from two adjacent elements, namely u δ + and u δ − , which are in general discontinuous.The matrix form of Eq.( 18) is where b e is the surface integral introduced in Eq. ( 19): We can use similar relations to those in Eqs.(11,12) to rewrite Eq. ( 20) as follows T W e B e g e (u)+ where the diagonal matrix Λ e (G i ) has on its diagonal the geometric factors defined in Eq. ( 12).The final matrix form of the DG method can be written as follows T W e B e g e (u)+ where is the inverse of the elemental mass matrix.
In the DG method the source of aliasing errors are similar to those identified for the CG method.Specifically, PDE-aliasing arises from the advection terms whilst geometrical-aliasing emerges from the weight matrices which contain the elemental Jacobians as well as from the G i/j terms which include the geometric factors.
A separate discussion is instead needed for the interface fluxes f (u δ + , u δ − ).They are evaluated at the edges of the elements on a given set of points (which may differ from the quadrature points used for discretising the problem within the element domain) and they are in general non-polynomial functions since the contributions to f (u δ + , u δ − ) can come from the left (u δ + ) and from the right (u δ − ) element with respect to the edge considered.If Eq. ( 5) is nonlinear, then this term is another source of PDE-aliasing which introduces additional numerical errors.Note that the boundary term b e also includes the contribution of the geometric factors and therefore also geometrical-aliasing issues can arise.

FR formulation
In this formulation the approximation space is equal to that of the DG scheme and it is defined in Eq. ( 16).The matrix form of the FR approach is similar to that of the DG scheme, where the solution polynomials are piecewise discontinuous.The main difference arises from the use of the differential form of the equations as opposed to the variational form.In addition, the FR approach embeds various high-order spectral element methods including a nodal discontinuous Galerkin method, the spectral difference method and a class of energy-stable schemes [27].Some well-established connections have been found in the literature.In this paper we show how the connections in [28] hold true also applying both the dealiasing techniques presented.
The FR approach in matrix form can be represented as follows: where D ξi are the elemental differentiation matrices with respect to the standard space, Λ e (G i/j ) are the same quantities as in Eqs.(15,23), with G i/j being defined in Eq. ( 12), and the boundary term is with Φ being the derivative of a suitable correction function defined in the standard space, f C,e being the standard correction flux defined as the difference between the interface flux f and the corresponding volumetric flux evaluated at the edge of element e, f e , where we note that if the interface flux corresponds exactly to the flux coming from element e the difference is zero.This difference needs to be calculated in the standard space since the derivatives of the correction functions belong to the standard space.To obtain the difference between f and f e in the standard space, Δ f , we apply a transformation which depends on the Jacobians J as well as on the previously defined geometric terms G i/j : where T = T (J, G i/j ).For a better understanding of the required transformations see [29] whilst for a more detailed and comprehensive discussion concerning the FR approach, the interested reader can refer to Huynh [23] and to Vincent, Castonguay, and Jameson [27,30].
In terms of aliasing considerations, we note that PDE-aliasing sources are exactly the same as the DG approach, and in particular the advection term along with the interface flux f are directly responsible for introducing PDE-aliasing if Eq. ( 5) is nonlinear.Regarding geometrical aliasing, whilst the geometric terms seem different if compared to the DG approach, they are in fact identical and they produce the same geometrical aliasing sources.

Summary
In the previous sections we have highlighted the sources of aliasing in both continuous (CG) and discontinuous (DG and FR) formulations.These are due to PDE nonlinearities as well as to geometrical nonlinearities.In particular, whilst in CG the PDE nonlinearities depend on the volumetric flux only, in DG and FR another source of nonlinearity to take into account depends on the interface flux.In the following sections we perform some numerical experiments with the three methods above, showing the effects of these aliasing issues and explaining different procedures to reduce or eliminate them.Regarding the discontinuous approaches we analyse the role of the interface flux for PDE and geometrical nonlinearities.We also show that the aliasing instabilities of two different versions of DG schemes and two types of FR schemes are equivalent, which is in agreement with previous observations [28].

Dealiasing techniques
In this section we describe the two dealiasing approaches considered in this work, namely: • Local dealiasing: PDE-dealiasing through consistent integration of the nonlinear terms6 ; • Global dealiasing: PDE-and geometrical-dealiasing through consistent integration of all the terms of the discretisation.
The first takes into account only PDE-aliasing errors and it locally uses a consistent integration of the nonlinear terms of the PDE.The second takes into account both PDE and geometrical aliasing effects.

Local dealiasing
The first dealiasing approach presented involves the consistent integration of the nonlinear terms in Eq. ( 5).This approach can be applied to all the three formulations previously introduced, namely the CG, DG and FR discretisations, and addresses the PDE-aliasing issues.In the following, we first consider a onedimensional case before outlining the extension to two and three dimensions through a tensor product expansion.The implementation presented here has been carried out in the spectral/hp element library Nektar++ [31].This implementation, and in particular its tensor product extension in both structured and unstructured sub-domains, is computationally efficient for high polynomial orders due to its use of the sum factorisation technique.

Dealiasing in one dimension
Consider a 1D approximate polynomial solution u(ξ) of order P represented by an expansion in terms of nodal polynomial functions φ i (ξ), such as Lagrange polynomials, on a set of GLL points The dealiasing strategy consists of 3 steps.

A.) 1D interpolation
The solution u(ξ) is interpolated onto a larger set of GLL quadrature points where the number of additional points Q required depends on the degree of the nonlinearity.This step can be achieved through a matrix-vector multiplication where I Q P →Q is the interpolation matrix whose dimensions are Q × Q P .Note that we could use any set of Q points and not necessarily GLL points.For instance, the larger set of points could be represented by Gauss-Legendre points.

B.) Collocation product
We then evaluate the nonlinear term on the expanded set of points using the interpolated values of the solution u Q (ξ j ) obtained at the previous step.For simplicity we consider a quadratic nonlinearity so that the product is calculated as

C.) Galerkin projection
Finally, we perform a projection onto the original set of points.In the following we use a Galerkin projection (GP), which is equivalent to a least squares projection of the nonlinear term evaluated on the larger set of points onto the original set of points, so that In a similar manner to step A, this can be interpreted as a matrix-vector multiplication where GP Q→Q P is the Galerkin projection matrix with M[i, j] = φ Q P i φ Q P j dξ being the mass matrix and W[i, i] = φ Q i dξ being the diagonal Gauss quadrature weight matrix defined using Q points.Note that the dimensions of GP Q→Q P are Q P × Q.
To summarise, Fig. ( 1) depicts a schematic representation of the dealiasing procedure described above for one-dimensional problems.We note that this procedure is exact when the isoparametric mapping that describes the coordinates of a physical element Ω e is affine and therefore the Jacobian is constant.In addition, we note that if M is diagonal then M has diagonal components of the weights at P set of points.In this case we observe that the projection is an interpolation matrix with the rows scaled by the inverse of the P quadrature point weights whereas the columns are rescaled by the quadrature point integration weights.

Tensor product extension to 2D and 3D
For the sake of simplicity, in this section we consider a 2D tensor product element and we remark where necessary the differences in terms of computational costs between the 2D and the 3D tensor-product element cases.
Consider a 2D approximate polynomial solution u(ξ 1 , ξ 2 ) of order P represented through a nodal expansion basis, such as Lagrange polynomials, onto a set of GLL points In the 2D (equivalently 3D) tensor product case, the dealiasing technique consists of the same conceptual steps as the 1D case but the operations to perform are more numerous and consequently computationally more costly.However, the exploitation of the tensor-product properties allows an efficient implementation through the use of the sum-factorisation technique.

A.) 2D/(3D) interpolation
The first step involves interpolating the solution u(ξ 1 , ξ 2 ) onto a larger set of GLL quadrature points as follows: where the number of additional points Q required depends on the degree of the nonlinearity.The interpolation performed is divided into two (or three in 3D) sub-steps as shown in Fig. is the interpolation along ξ 2 .This implementation allows a reduction of the floating point operations required.In this case the overall number of floating point operations required is while by doing a 2D interpolation for a non-tensor product basis the floating point operations needed becomes O(Q 2 × Q 2 P ).Following a similar procedure for the 3D case gives a number of floating point operations that is O(Q×Q 2 P +Q 2 ×Q 2 P +Q 3 ×Q P ) versus O(Q 3 P ×Q 3 ) operations otherwise required for a full 3D interpolation.Further details of this type of sum factorisation can be found in [24].
Note that also in this case we could have used any larger set of points (such as Gauss-Legendre points) and not necessarily GLL points.

B.) Collocation product
The second step involves evaluating the nonlinear term on the expanded set of points using the interpolated values of the solution u Q (ξ 1 , ξ 2m ) obtained at the previous step.As with the 1D case, we evaluate the product of the interpolated values u Q (ξ 1 , ξ 2m ) using a collocation projection as follows where we again consider a quadratic nonlinearity for the sake of simplicity.

C.) Galerkin projection
The third step involves performing a projection of the nonlinear term onto the original set of points.To do so we use a Galerkin projection (GP).Similar to the first step, we split this operation into two sequential suboperations as shown in Figure 3. First (Figure 3(a)) we perform a Galerkin projection in direction ξ 2 and successively (Figure 3(b)) we perform a Galerkin projection in direction ξ 1 .
(ξ1 , ξ2 m ) In Figure 4 we show an overview of the steps above for the 2D/3D tensorproduct case.A similar tensor-product based approach can also be used for triangles in 2D as well as prismatic and tetrahedral elements in 3D as reported in Appendix C.

Global dealiasing
In order to address both PDE-and geometrical-aliasing issues, one can oversample all the terms present in the numerical formulation.For each of the CG, DG and FR formulations considered, we perform the following steps: Note that we do not make any distinctions between the 1D or the 2D/3D cases, nor between different element shapes, since here we do not explicitly exploit the tensor product advantages (although some operators may still have a tensor product form) used in the previous sections for the PDE-dealiasing through consistent integration of the nonlinear terms.This approach is similar in nature to the one presented by Kirby and Karniadakis [19].However, here we additionally perform an over-sampling of the geometric terms, the Jacobian and the geometric factors.Based on the greater number of floating operations to be performed, this dealiasing technique is more computationally expensive than the consistent integration of the nonlinear terms presented in the previous section but can address both PDE-and geometrical-aliasing issues.

Numerical results
In this section, we consider the impact of the dealiasing schemes outlined above, when applied to continuous and discontinuous discretisations, from two perspectives.First, we perform a quantitative study of the effects of these dealiasing techniques using a simple transport advection equation for which a periodic solution is known, but where the advection velocity is a nonlinear polynomial in space.We then consider various cases including deformed elements, the contribution of boundary terms in the dealiasing process, and reinforce the inequalities which lead to exact integration of discretised nonlinear functionals.We finally consider qualitative effects, giving examples of how dealiasing can enhance the stability of both compressible flows and incompressible Large Eddy Simulations.Even if, in the examples presented, the dealiasing techniques show a stabilising effect, dealiasing does not guarantee stability and in other flow studies the authors experienced instabilities even applying a consistent integration of the flux functions.In this case, for stabilising the simulation, they added an artificial viscosity of the type described in [32] and [7].
In our experiments, we chose Gauss-Lobatto-Legendre points for consistent integration.The use of a Gauss-Legendre quadrature is more efficient for reducing the aliasing issues; however, Gauss-Lobatto-Legendre quadrature is less expensive because the interpolants do not need to be evaluated at the elemental boundaries to project flux terms into the volume.

PDE aliasing
To understand the effect of the PDE-dealiasing techniques presented in the previous section, we consider the transport equation where u is the conserved variable, a x and a y are the advection velocities along the x and y directions, (x b , y b ) denote the boundaries of the numerical domain, P adv is the polynomial order of the flux and g(t) is a time dependent periodic function g(t) = cos(πt/T ) (39) on the time time interval [0, T ].The solution u evolves in time in such a way that the initial data is recovered at every period T , that is We select a period T = 1 and a final time 4T so that the final solution is equal to the initial condition in order to calculate the error.This strategy is the same as that adopted in [33] and it is widely used to evaluate properties of numerical methods when the exact solution is not available.
In the following, we first perform a comparative study between the Local and Global dealiasing approaches on the case where the advection velocities is linear, so that P adv = 1.This will highlight the connections between the two dealiasing approaches and ensures that some general properties hold true for all the results presented in the following sections.We then show how increasing P adv , and therefore the order of the advection velocities, can lead to a more important role of the interfaces in the case of the discontinuous discretisations.In all the simulations we use a tensor-product Lagrangian polynomial basis of order P = 4.

Comparative study of the different dealiasing techniques
In this section we show the results obtained using Local and Global dealiasing techniques for the test case in Eq. ( 38) when P adv = 1 on a regular grid of quadrilaterals, so that no geometrical aliasing effects are present.
For each simulation we use a forward Euler time integration scheme and the time step was chosen to be sufficiently small in order to consider the temporal error negligible.We remark that any explicit time integration scheme such as 2 nd or 4 th order Runge-Kutta schemes can be used.In space we consider each of the CG, DG and FR discretisations in turn, where the polynomial order is set to P = 4 and the initial condition is given by a collocation projection (i.e. using Q = 5 quadrature points).Throughout the following 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, N is the total number of quadrature points and w i are the quadrature weights.
In Table 2 we show the L 2 errors for the different dealiasing strategies applied to the CG, DG and FR formulations for the mesh shown in Figure 6(a).Here, 'LMM' stands for lumped mass matrix, which is the diagonal matrix obtained when Q = 5 quadrature points are used to obtain a collocation scheme, while 'EMM' refers to the exact mass matrix obtained for higher quadrature orders.This table illustrates the saturation of the L 2 error up to machine precision for Q = 6 as well as the equivalence between Local and Global dealiasing when the number of quadrature points used for the nonlinear terms are the same in the two approaches, since the geometric terms for a regular grid are constant.We also note that the equivalence of the DG and FR schemes presented in [28] holds when using a larger number of quadrature points than a collocation method.
Regarding the dealiasing of the interface fluxes, we note that we obtain a saturation of the error up to machine precision.This is because the interface fluxes for the mesh in Figure 6(a) are polynomial functions since their values come from just one side with respect to a given boundary 7 .If we instead consider the mesh in Figure 6(b), we are in a case where the interface fluxes come from both sides of a given interface as shown in Figure 5; therefore, they do not lie in the polynomial space and error saturation up to machine precision cannot be expected as presented in Table 3.This means that, in general, it is difficult to fully control the PDE-aliasing arising from the interfaces because, in this case the functions representing the interface fluxes lie outside a polynomial space.However, the error does decrease as the number of quadrature points is increased.
Table 3: L 2 errors for the different dealiasing strategies applied to the DG and FR formulations: Case B mesh as in Figure 6(b).

Role of dealiasing for higher order advection velocities
In this section, we show the results using advection velocities with nonconstant coefficients (i.e.spatially varying polynomials), so that the order P adv is chosen to be greater than one in Eq. (38).The numerical parameters are chosen to be identical to those in the previous section, and we use the grid depicted in Figure 6(a) to avoid non-saturation of the error when dealiasing the interfaces.
In Figure 7, each point represents the difference between the L 2 error calculated using (Q V , Q I ) quadrature points and the error calculated using (Q V + 1, Q I + 1) quadrature points in the LMM case 8 .We note that the machineprecision saturation of the error satisfies the inequality: where Q is the minimum number of quadrature points to exactly integrate the PDE-aliasing for a given polynomial advection velocity possessing a nonlinearity of order P adv and for a given expansion order P .In Appendix B.1, Figure B.17 shows the actual L 2 errors vs. (Q V , Q I ) for the LMM case while Tables B.6, B.7, B.8 and B.9, show the tabulated values of the L 2 errors for P adv = 1, 2, 3, 4, which also include the EMM case.From these results it is possible to see how, for the discontinuous formulations, the use of Local dealiasing on both the volumetric and interface fluxes was, for almost all of the cases considered, the best choice in terms of L 2 error among the dealiasing techniques applied, except for EMM P adv = 4 where instead Local dealiasing of only the volumetric flux provided a better result.To further quantify the effects of interface dealiasing, in Table 4 we show the difference between the L 2 errors obtained for Local dealiasing of only the volumetric flux, and Local dealiasing of both volumetric and interface fluxes, taken at the point where the error saturates according to Figure 7 for each value of P adv .We see from this data, how the difference between performing only volumetric dealiasing and volumetric plus interface dealiasing increases as P adv increases except for the last case (P adv = 4).The gap between using a lumped (LMM) and an exact (EMM) mass matrix is notable, and is due to the implicit Global dealiasing in the case of EMM which in turn reduces the effects of the interface dealiasing (through the use of a quadrature that can integrate a higher polynomial order).

Link between geometrical-and PDE-aliasing
Now that the aliasing effects of the discretisation of the PDE have been examined, in this section we show how geometrical-aliasing can be linked to PDE-aliasing, with some differences arising from the fact that geometrical nonlinearities play a slightly different role in the underlying formulation compared to PDE nonlinearities.
We consider a mesh composed of a single standard quadrilateral element Ω st = [−1, 1] 2 , where we curve either only the top edge, or both top and left edges.For each configuration, we applied different curvatures to the edges by means of Legendre polynomials of order P = 1, 2, 3 and 4. In Figure 8(a) and (8(b)) we show two examples for the mesh of order four.
In each case, the isoparametric mapping, and therefore the Jacobian determinant, is represented by an expansion of basis functions over the standard element, so that for example a P = 4 expansion contains 5 × 5 = 25 different unique modes.Figure 9 shows the magnitude of each modal coefficients of the Jacobian determinant for the 4 th order case when projected onto an orthonormal basis of order 5 (leading to 36 coefficients), in order to emphasise that the polynomial space is sufficiently large to capture the Jacobian mapping.In particular, we note that the Jacobian for both cases is effectively of fourth order, but in the case of Figure 9(a) we have nonzero modes only in the x-direction up to the fifth coefficient, whereas in the case of Figure 9(b) the nonzero modes span both x− and y− directions again up to the fifth coefficient.The sixth coefficient in both directions is zero, showing adequate support for the Jacobian determinant expansion.The test case we consider is identical to that of Eq. ( 38), aside from the use of spatially-constant advection velocities where g(t) is defined in Eq. (39).We choose a period T = 0.5 and a final time 40T .For the time-integration we used a 2 nd -order Runge-Kutta scheme while the polynomial order was P = 14 in order to have a sufficiently resolved problem.The initial condition was applied using a collocation projection.Throughout the results in this subsection, we used a Global dealiasing technique in order to target the geometrical aliasing.We first present estimates of the quadrature necessary for correctly integrating a deformed or curved mesh.For doing so, we take into account the leading order within the problem under investigation, which is the mass matrix for both CG and DG discretisations.In Figure 10 each point denotes the difference between two mass matrices, the first obtained for Q and the second for Q + 1 quadrature points and it is calculated as follows Table B.10 in Appendix B.2 quantifies the results presented in Figure 10.We note that the mass matrix does not change (approximately up to machine precision) after we have reached the minimum number of quadrature points to exactly integrate the polynomial order used for describing the geometry, given by the inequality where Q is the minimum number of quadrature points to exactly integrate the mass matrix for a given polynomial of order P order describing the geometry and for a given expansion order P .Eq. ( 45) is identical to Eq. ( 42) with the only exception being that P geom now refers to the polynomial order of the geometry deformation.The above result is consistent with the Gauss-Lobatto-Legendre quadrature rules.Figure 11 shows the differences in terms of L 2 errors as defined in Eq. ( 44) instead of differences between mass matrices.The main points to note are the following: • CG reaches a saturation of the error approximately up to eleven or twelve digits for all of the meshes.This is due to the tolerances used for inverting the mass matrix which, in this case, do not achieve machine precision.
• DG and FR show machine-precision saturation of the error for the linear mesh (i.e.J ∈ P 1 ) whilst the saturation happens at a higher level as the mesh order increases.This behaviour is due to the interface fluxes which come from both the left (internal domain) and right (boundary condition) sides of a given interface, since the normals are spatially varying across the curved edge(s).Therefore, for complex geometries, it is difficult to fully control the geometrical aliasing arising from the interfaces.
• The non-saturation of the error to machine precision for each of the spatial discretisations is not due to rational functions which do not lie in the polynomial space.In fact, the geometric terms encounter cancellations which lead to just one geometric factor which does belong to the polynomial space considered.Thus, it can be fully represented within the polynomial representation used for these test cases.
As a final comment, we note that the connections between DG and FR established in [28] also hold for a deformed/curved mesh.

Flow applications
Having considered the quantitative effects of the Local and Global dealiasing techniques, in this section we consider two examples that highlight how stability can be enhanced by using appropriate dealiasing strategies.We begin by considering a compressible Euler test case, and show that the lack of dealiasing can lead to the buildup of numerical entropy in the solution field.Secondly, we consider the LES of a NACA 0012 wingtip, and show how dealiasing can increase the robustness of the simulation.As previously stated, however, the consistent integration does not imply stability and in some simulations the high frequency modes of the nonlinear flux functions are not eliminated through dealiasing.In this case, other strategies need to be applied (eventually in conjunction with consistent integration) such as the use of artificial viscosity or modal filtering.

Compressible inviscid subsonic flow past a cylinder
In this section, we present an inviscid compressible flow past a cylinder governed by the compressible Euler equations.We set the Mach number to 0.2, use a FR-DG scheme for the spatial discretisation and a 4 th -order Runge-Kutta scheme for the time-integration.We use a circular mesh composed of 9612 quadrilateral elements and the cylinder is described by third order splines in order to achieve a high order description of the geometry.In Figure 12 we show the mesh used for simulations with the Jacobian determinant distribution.We use two different polynomial orders: P = 2 and P = 3 and the Global dealiasing technique.Figure 12 indicates that the elements of the mesh are of good quality and not highly distorted, although a third order curvature is applied to the cylinder wall.We therefore expect PDE aliasing to play a bigger role than geometrical aliasing.In Figure 13(a) we show the numerically-generated entropy for the case P = 2, Q = 3 without dealiasing techniques applied.We can clearly see a buildup of entropy from the posterior stagnation point, which is sensitive to aliasing-driven instabilities due to the strong gradients of the solution in this region.Note that ideally for an inviscid simulation, the entropy should be numerically zero.In this case the simulation diverged after a few timesteps; however, by applying the Global dealiasing techniques using Q = 4, we were able to reduce the aliasing errors.In Figure 13(b) we show a snapshot of the solution field taken at the same value of time as Figure 13(a).We note that the numerically-generated entropy at the posterior stagnation point is no longer present.The simulation for P = 3 and Q = 4 depicted in Figure 14 was also unstable.As in the previous case (i.e.P = 2, Q = 3), the use of the Global dealiasing technique stabilised the simulation.In Table 5 we report the maximum and minimum values of entropy across the computational domain for all the cases considered.As we can see, the numerically-generated entropy, which in this case can be seen as a measure of the numerical error, roughly saturates at Q = 6 for P = 2 and at Q = 8 for P = 3. Clearly the saturation is not up to machine precision because of the high-order description of the geometry and the deformed elements throughout the domain which introduce additional 'nonlinearities' into the problem.Also, the right-hand side of the compressible Euler equations is composed by rational functions.This in turn means that the right-hand side cannot be represented in a polynomial space and, therefore, a machine-precision saturation of the solution cannot be expected.However, a minimisation of the relative error can be seen as Q is increased and this may also indicate that the PDE-aliasing is potentially dominant for this problem.
Note that, for both the polynomial orders considered, we also applied the Local dealiasing approach but the results are not shown here, as they are similar to those obtained using the Global dealiasing technique.This further reinforces the role of the PDE-aliasing for this problem.

Incompressible viscous flow past a wing tip
In this section we show the results of an incompressible viscous flow past a NACA 0012 wingtip, originally studied experimentally by Chow et al. [34].
The simulation was obtained using the CG approach, with a velocity correction scheme being used to discretise the incompressible Navier-Stokes equations [35].While the experimental Reynolds number was set at Re = 4.6 × 10 6 , initial simulations were performed at a lower Reynolds number of Re = 1.2 × 10 6 .Local dealiasing, when combined with a spectral vanishing viscosity to dampen high-frequency energy buildup, yielded a stable simulation.However, in order to increase the Reynolds number and bring the simulation in line with experimental conditions, we found that without an increase in the global quadrature order (i.e.Global dealiasing technique), the simulation quickly became unstable.This highlights the role of geometric as well as PDE aliasing under these conditions for the mesh and polynomial resolution considered.
In Figure 15 we illustrate the dynamics of the flow by showing the isocontours of helicity for a P = 3 approximation, where the domain is represented using prismatic elements around the boundary and tetrahedral elements in the rest of the domain.In Figures 16(a) and 16(b), for the same simulation, we represent the aliasing errors.The top image shows up to 30% aliasing error with respect to the magnitude of the nonlinear terms while the bottom image shows regions of aliasing near the wing surface where we observe up to 600% error, highlighting the crucial nature of dealiasing in this problem.

Summary and conclusions
We have presented two dealiasing approaches based on the concept of consistent integration of the nonlinear terms or over-integration if the integration order is based on exact integration of linear operators.The first technique presented, namely the Local approach, targets the PDE-aliasing sources which arise from the nonlinearities contained in the PDE itself.This approach exploits the tensor product structure of the elements to reduce the number of floating point operations, making it computationally efficient at higher polynomial orders.The application of Local dealiasing is particularly useful when we have localised nonlinearities such as in the incompressible Navier-Stokes equations if geometric aliasing is not significant.The second technique, namely the Global approach, involves the use of higher quadrature order than is typically necessary for linear PDEs.This technique targets both PDE-and geometrical-aliasing sources which arise when deformed or curved elements are used, as often occurs in industrially-relevant problems where complex geometries are present.In contrast to the first approach, this technique is computationally more intensive, but it can address all the aliasing sources arising in any spectral element method when coupled with a sufficiently large number of quadrature points.We applied both the dealiasing techniques to different high-order spectral element methods, showing the main implementation differences between continuous (CG) and discontinuous (DG and FR) discretisations.
Based on the above analysis and supported by numerical examples we can state that: 1.The consistent integration rules presented in [19] are generally true and, for GLL points, can be complemented by the following expression where Q is the minimum number of quadrature points to exactly integrate the highest-degree of nonlinearity P order within the problem considered and for a given expansion order P exp .Note that in the above expression P order can be a combination of geometrical and PDE nonlinearities and therefore the overall degree of the nonlinearity can be higher than that dictated by the PDE itself.2. In a discontinuous discretisation, it is generally not possible to fully control (i.e. up to machine precision) all aliasing effects, since the boundary terms which dictate interface contributions introduce non-polynomial functions into the problem being solved.In addition, geometrical aliasing affects a discontinuous discretisation more than a continuous, due to the geometric factors appearing in the boundary term.3. Geometrical dealiasing, in contrast to PDE aliasing, is responsible for modifying the mass matrix and it may, therefore, change the numerical dissipation and dispersion properties of the discretisation employed as shown in [36] and [37].4. The connections between DG and FR presented in [28] hold also for the dealiasing strategies described.5.The two dealiasing strategies can be used in a complementary manner and when applied to challenging applications they have proven effective and have increased the numerical stability of the simulations although they do not guarantee stability.This result is in agreement with the findings in the recent work of Malm, Schlatter et al. [26] where a possible connection between numerical stability and consistent integration has been shown for the incompressible Navier-Stokes equation (i.e.polynomial nonlinearities).
We additionally note that interface dealiasing may play a bigger role as the order of the interface flux function increases, thus becoming comparable or more important than the volumetric dealiasing contribution.

Integration within a general-shaped region
To integrate over an element Ω e of a two-dimensional domain, we transform this region into the standard region Ω S and we have: where J 2D is the two-dimensional Jacobian defined as The partial derivative can be evaluated using the maps defined in Eqs.(A.1, A.2, A.3).If the element is straight-sided and has parallel sides, then the partial derivatives and therefore the Jacobian are constant.For deformed elements the Jacobian can be evaluated and stored at the quadrature points.This represents the Jacobian as a polynomial function and can, therefore, increase the polynomial order of the integrand.The integral is then evaluated by an appropriate quadrature rule.An analogous procedure is used for a three-dimensional element where the three-dimensional Jacobian J 3D can be evaluated from the partial derivatives given by the mappings.

Differentiation within a general-shaped region
To differentiate over an element Ω e we use the chain rule: where d is the dimension of the domain.For d = 2: The terms on the right-hand side are known from the mapping.We can now evaluate the derivatives as all the partial derivatives can be expressed in terms of differentials with respect to ξ 1 and ξ 2 .where the collapsed coordinates (η 1 , η 2 ) ∈ [−1, 1] 2 are given by the Duffy transformation The use of this collapsed coordinate system within the standard quadrilateral region [−1, 1] 2 leads to a tensor product of quadrature points within the standard triangular region, as shown in figure C. 18.Although this transformation is singular at the top vertex of the triangle, the use of Gauss-Radau quadrature in the η 2 direction, in which the top vertex is excluded, leads to a formulation without geometric singularities.Additionally, Gauss-Radau points with a choice of α = 1 and β = 0 naturally incorporate the Jacobian term which appears in integrals over the standard triangular region when weighted by a constant factor of 1 2 [24].In this setting, the sum-factorised dealiasing methods described in section 3 can be utilised for triangular elements by applying the stated tensor product logic to the grid of (η 1 , η 2 ) quadrature points, without applying the Duffy transformation to obtain the desired dealiasing effect in the simplex region.Mathematically, this approach projects the nonlinear terms down to the tensor product space spanned by the (η 1 , η 2 ) space which is typically richer than the space spanned by the triangular expansion in the (ξ 1 , ξ 2 ) space.However, since the inner product using the collapsed coordinates spans the (η 1 , η 2 ) space, this is sufficient to ensure that the nonlinear product is correctly integrated.This logic extends to three dimensions for both prismatic elements and tetrahedra when coupled with suitable extensions of the Duffy transformation to a hexahedral region.

(
B e ) T W e B e d u e dt + +(B e ) T W e D e x1 f e (u) + D e x2 g e (u) + D e x3 h e (u) = 0,

(
B e ) T W e B e d u e dt − (D e x1 B e ) T W e B e f e (u)+ −(D e x2 B e ) T W e B e g e (u)+ −(D e x3 B e ) T W e B e h e (u) + b e = 0, u P (ξi) ∈ P P ξi ξi ξj ξj f P (ξi) = u(ξi) n ∈ P P

Figure 1 :
Figure 1: Conceptual flow chart of the Local dealiasing approach through consistent integration of the nonlinear terms for 1D problems.

Figure 4 :
Figure 4: Conceptual flow chart of the Local dealiasing approach for the 2D/3D tensor-product case.
A.) Over-sample the solution polynomial onto a larger set of points, u P → u Q with Q > Q P ; B.) Evaluate the flux onto the larger set of points F(u Q P ) → F(u Q ); C.) Over-sample the geometric factors ∂xi ∂ξj from which the Jacobian determinant J and geometric factors Λ(G i ) are calculated; D.) Over-sample the differentiation matrices D ξi ; E.) (DG/FR only) Over-sample the boundary terms b e or b e ; F.) Assemble the right-hand side onto the larger set of points; G.) Project the calculated solution back to the original set of coefficients through a Galerkin projection.

Vn = ax nx 0 ≥Figure 5 :
Figure 5: Upwind flux at the interface between two elements: Aliasing arising from mixed information coming from the left and the right element.

Figure 6 :
Figure 6: Meshes for the 2D linear advection test cases: Case A in Figure 6(a) and Case B in Figure 6(a).

P adv = 4 Figure 7 :
Figure 7: Differences between the L 2 errors obtained with (Q I , Q V ) and (Q I + 1, Q V + 1) vs (Q I , Q V ) using the Local dealiasing technique for the case of LMM.
(a) Top edge curved (b) Top and left edges curved

Figure 8 :
Figure 8: Examples of 4 th order meshes employed to investigate geometrical aliasing.

Figure 9 :
Figure 9: Coefficients in an orthonormal expansion basis for the 4 th order meshes showed in Fig.(8).

Figure 10 :
Figure 10: L 2 norm of the difference between mass matrices calculated with Q and Q + 1 GLL points for the four orders considered and using both the mesh configurations shown in Figure 8. Case A corresponds to Figure 9(a) and Case B corresponds to Figure 9(b).

4 Figure 11 :
Figure 11: Differences between the L 2 errors obtained with Q and Q + 1 GLL points for the four orders considered using the mesh in Figure8(a).

Figure 12 :
Figure 12: Jacobian distribution on the mesh used for the compressible inviscid simulation of a cylinder.

Table 1 :
Number of GLL quadrature points for the GLL quadrature to be exact up to machine precision as a function of the polynomial order of the integrand.

Table 2 :
L 2 errors for the different dealiasing strategies applied to the CG, DG and FR formulations: Case A mesh as in Figure

Table 4 :
Percentage difference between the L 2 errors for Local (Q V ), Local (Q I ) and Local (Q V , Q I ) at error saturation.

Table 5 :
Maximum and minimum values across the mesh of numerically-generated entropy for both P = 2 and P = 3.