Spatial eigensolution analysis of discontinuous Galerkin schemes with practical insights for under-resolved computations and implicit LES (cid:2)

The study focusses on the dispersion and diffusion characteristics of discontinuous spectral element methods - speciﬁcally discontinuous Galerkin (DG) - via the spatial eigensolution analysis framework built around a one-dimensional linear problem, namely the linear advection equation. Dispersion and dif- fusion characteristics are of critical importance when dealing with under-resolved computations, as they affect both the numerical stability of the simulation and the solution accuracy. The spatial eigensolution analysis carried out in this paper complements previous analyses based on the temporal approach, which are more commonly found in the literature. While the latter assumes periodic boundary conditions, the spatial approach assumes inﬂow/outﬂow type boundary conditions and is therefore better suited for the investigation of open ﬂows typical of aerodynamic problems, including transitional and fully turbulent ﬂows and aeroacoustics. The inﬂuence of spurious/reﬂected eigenmodes is assessed with regard to the presence of upwind dissipation, naturally present in DG methods. This provides insights into the accu- racy and robustness of these schemes for under-resolved computations, including under-resolved direct numerical simulation (uDNS) and implicit large-eddy simulation (iLES). The results estimated from the spatial eigensolution analysis are veriﬁed using the one-dimensional linear advection equation and successively by performing two-dimensional compressible Euler simulations that mimic (spatially develop- ing) grid turbulence.


Introduction
One of the biggest challenges that computational fluid dynamics (CFD) is facing is to expand its current capabilities to flow problems that are still not well captured by the prevailing numerical technologies adopted today [1] . In industry, for instance, the use of low-order numerical techniques in conjunction with approximated and steady-state-tailored numerical strategies, such as Reynolds-Averaged Navier-Stokes (RANS) approaches, is almost ubiquitous. While these tools are well established and numerically robust, they struggle to accurately describe a wide range of problems that are of practical interest in various branches of engineering and applied sciences. The underlying key of these flow problems, which usually involve relatively high Reynolds numbers (up to the order of millions) and turbulence, is their inherent multi-scale nature. This aspect, still today (and likely for many years), prevents the fully-R Dedicated to Professor Eleuterio Toro on the occasion of his 70 th birthday. * Corresponding author.
E-mail addresses: mengaldo@caltech.edu , gianmarco.mengaldo@gmail.com (G. Mengaldo). resolved direct numerical simulation (DNS) of this type of flows. The inability to accurately simulate such a wide range of problems crucial to several applications, in turn, inhibits significant improvements of computer-based predictive skills in fluid dynamics. Given these limitations, it is essential to explore alternative ways to enhance the capabilities of CFD and facilitate their adoption in the broader industrial community working in the field [1,2] .
From this perspective, high-fidelity computations relying on both implicit and explicit large-eddy simulation (LES) [3][4][5][6] have been used extensively in academia, both for model problemse.g. [7][8][9][10] -and for industry-relevant test cases -e.g. [11][12][13] . However, especially in industry, the rare use of LES-based approaches have been commonly confined to low-order numerical methods, such as finite volume and finite difference approaches, that are more likely to yield inaccurate results because of the unfavorable diffusion and dispersion properties of the underlying numerics. Most recently, on the other hand, high-order spectral element methods (SEMs) [14] , that have been used for many years in academia within the context of under-resolved simulations of turbulent flows, have been drawing the attention of several researchers and industry practitioners [15][16][17][18][19][20][21] . High-order methods are in fact known to offer numerical benefits over traditional loworder schemes especially when small-scale flow features need to be captured and propagated correctly over long distances [2,22,23] , such as in high-Reynolds number open flows. In addition, highorder SEMs are better suited for industrial LES due to their flexibility and superior resolution power per degree of freedom. In particular, the so-called under-resolved DNS (uDNS) approach used in conjunction with SEMs [24][25][26] has shown to be very promising. This approach uses numerical features to treat the under-resolved scales and can be cast within the better-known implicit largeeddy simulation (iLES) framework, although it has not been proven that the associated truncation errors resemble subgrid-scale models of classical LES approaches [5] . The use of SEMs for uDNS/iLES, hereafter SEM-iLES approach, can accurately resolve features that are impossible to grasp from more established CFD methodologies, such as RANS and DES, and can complement explicit LES methods.
SEMs provide mesh and polynomial refinement (hp-refinement) that allow the (possibly adaptive) bespoke tuning of the numerical properties, namely diffusion and dispersion, of the numerical method adopted that are essential when dealing with underresolved simulations of turbulent flows at high Reynolds number [26] . In addition, SEMs naturally permit an efficient and high-order handling of complex geometries, an aspect that is particularly crucial in many industry-relevant problems.
However, to obtain reliable and meaningful data from SEM-iLES (and more generally from under-resolved simulations), it is of fundamental importance to better understand the numerical characteristics of the underlying numerics regarding wave-like solution components. The present study investigates the dispersion and diffusion/dissipation properties of discontinuous Galerkin (DG) spectral element methods, see e.g. [14,[27][28][29] , through the spatial eigensolution analysis framework applied to the one-dimensional linear advection equation. The findings obtained through the eigensolution analysis are successively verified in linear advection simulations and via a two-dimensional test case that resembles the behaviour of grid turbulence at very high-Reynolds numbers.
In addition, the well-established connections between DG and flux reconstruction (FR) schemes [30][31][32] allows drawing a more comprehensive picture of discontinuous SEMs (DSEMs) in the context of uDNS/iLES methodologies. In particular, the specific FR scheme recovering DG is algorithmically identical to the DG scheme analyzed here, therefore the results of the spatial eigensolution analysis presented in this paper hold for both schemes and a similar analysis framework can be used to investigate the wider class of energy-stable FR (ESFR) schemes [33] .
The application of eigensolution analysis to spectral element methods is not new [24,[34][35][36][37][38][39][40] . Nevertheless, most of the dedicated literature covers only the temporal eigenanalysis approach, which (strictly speaking) is valid more effectively for periodic problems. When however inflow/outflow type boundary conditions are present, the spatial approach is to be preferred, this being the relevant analysis for numerical propagation of waves and more generally of signals in open flows. While the former approach focuses on the temporal evolution of spatially-coherent signals, the latter focuses on their spatial evolution and assumes temporal coherence. The two approaches complement each other and are equally important for a comprehensive understanding of any discretisation technique suited for advection-diffusion problems.
To the authors' knowledge, only [35] covers the spatial eigensolution analysis of DG schemes offering a systematic study. Some aspects of DG's spatial analysis are also discussed in [41] , but attention is there given to the distribution of numerical errors within individual elements and not to propagation characteristics in the frequency domain. Besides that, reference [42] also touches DG's spatial eigenanalysis, but focuses mainly on the super-convergence properties of well-resolved computations, hence being of lesser in-terest for more complex, practical applications. In summary, the present work revisits the results obtained for DG in [35] while discussing them in a more qualitative manner and highlighting the details relevant to under-resolved simulations. Another point addressed here is how different amounts of upwinding affect DG's numerical dispersion and diffusion properties. These effects can be observed in practical simulations depending on the Riemann solver employed. In addition, we remark that, by harnessing the wellestablished connections between DG and FR numerical discretisations, the results presented here hold for the particular FR scheme recovering DG. This paper is organized as follows. In Section 2 , we present the spatial eigensolution analysis framework for DG. Section 3 addresses the results of the spatial eigensolution analysis obtained for the linear advection equation, while Section 4 shows the results for a one-dimensional linear advection test case that corroborates the eigensolution analysis presented in Section 3 . In Section 5 , we propose a two-dimensional test case to demonstrate the insights obtained from the eigenanalysis in a more complex scenario. Section 6 gathers concluding remarks and highlights our findings.

Spatial eigensolution analysis of discontinuous SEM
In this section, we introduce the spatial eigensolution analysis framework adopted in this work to investigate the dispersion/diffusion properties of discontinuous SEMs. We successively apply this framework to the DG method and present the relevant aspects pertaining the connections between DG and high-order FR schemes that allow us to transfer the results found for DG directly to the particular FR scheme recovering DG methods.

Spatial eigensolution analysis framework
Let us consider the linear advection equation in one dimension, given by ∂u ∂t with suitable initial and boundary conditions, where a is the advection velocity, t represents time, x is the spatial coordinate and u is the independent variable. We assume wave-like solution components in the form where κ and ω are the wavenumber and the (angular) frequency of the considered component, respectively, while i = √ −1 is the imaginary unit. Inserting (2) into (1) yields the classical dispersion relation where we introduced variable = ω/a . In general, both ϖ and κ can be complex, but, for the sake of investigating the dispersion/diffusion properties of the DG scheme, we assume that one of them is real. In particular, from (3) , we define The first relation, (4a) , with κ real, is the basis of temporal analyses, while the second relation, (4b) , and ϖ real, is the reference in spatial analyses.

Spatial eigensolution analysis for DG
We start the spatial eigensolution analysis by formally discretizing (1) with the DG method, as shown in Moura et al. [24] . In the DG framework the computational domain is divided into nonoverlapping elements e , such that = e e . In each element e , the solution is represented by a sum of weighted basis functions u | e ∼ = P j=0 u j (t) φ j (ξ ) . (5) where φ j is chosen to be the orthonormal Legendre polynomial of degree j , defined in the standard domain st = [ −1 ; 1] . Equation (1) is then required to vanish locally, that is By substituting (5) into (6) , and exploiting the properties of the orthonormal Legendre polynomials as basis functions, we can retrieve the following expression where h is the mesh spacing -i.e. the length of element e . By expanding the expression of the solution u in the first integral and by applying the Gauss theorem to the second term on the righthand side of (7) , we then obtain and L e represent the right ( R ) and left ( L ) boundaries of element e , respectively, and ˜ u = ˜ is the inter-element numerical flux that depends on the values at the left ( u ) and right ( u ) sides of each given interface considered. In this work, we use a generalized upwind flux, that reads where S a = | a | /a is the sign of the advection velocity and β is an upwinding parameter that accounts, in particular, for standard upwind ( β = 1 ) and fully central ( β = 0 ) discretisations. Finally, by assuming an equispaced mesh and using (8) and (9) , the semidiscrete linear advection equation can be expressed in vector form, for each element of the computational domain, as where u = { u 0 , . . . , u P } T and the ij th entries of the above matrices Note that (10) links three consecutive elements and therefore suffices for a multi-element eigensolution analysis [35] , provided that all elements are linked in the same way.
After having derived the semi-discrete form of the DG method, we now use the spatial eigenanalysis framework briefly introduced in the previous section to investigate the dispersion/diffusion properties of the DG scheme. To that end, we work with element-byelement projections of (2) , consistent with the semi-discrete problem (10) . Direct projection of (2) within e = [ x e − h/ 2 , x e + h/ 2] leads to the vector of coefficients u = { u 0 , . . . , u P } T given by where x e is the middle point of the element and α = [ α 0 , . . . , α P ] T with α j (κh ) = st exp (iκhξ / 2) φ j (ξ ) d ξ . (13) Although coefficients α j above can be easily evaluated numerically, analytical formulas for them have been derived and are available in [24] . These values are however only relevant for the shape of wave-like numerical (polynomial) solutions within individual elements. Note also from (12) that the solution vectors at the left and right neighbouring elements, cf. (10) , can be expressed as If we now insert equation (12) into the semi-discrete advection problem (10) , with real ω and complex κ (following the spatial approach), we obtain where is a square matrix of size m = P + 1 , which is the number of (independent) element-wise degrees of freedom, P being the polynomial order.
In the temporal approach, for each real κh given, m complex values for ϖh are obtained directly through the eigenvalues of Δ. In the spatial approach, however, obtaining complex-valued κh from a given real ϖh requires the solution of the equivalent determinant problem from (14) det where z = exp (iκh ) and I is the identity matrix of required size.
As discussed in [35] , relation (16) leads to the solution of a characteristic polynomial which is quadratic in z and therefore admits (up to) two solutions. These are interpreted as two solution eigenmodes, one physical and another spurious. The latter usually comes in the form of a reflected mode generated e.g. at the interface between regions of different mesh spacing h , or near outlet regions where boundary conditions are applied without additional care. Examples from simple numerical experiments can be found in [35] .

Numerical evaluation of the eigensolutions
We now enter into some detailed technicalities required to solve (16) . In this study, the solution of the determinant problem (16) is achieved through a root-finding algorithm, namely, MAT-LAB's function newtzero . Although the latter benefits from a reasonable guess for the root's location, it is always possible to try random guesses for z until the desired solution is obtained.
Function newtzero typically returns the physical and spurious solutions for z altogether, for each given ϖh , hence it is important to track them separately. This can be done by realizing that physical and spurious/reflected solution components are expected to have opposite amplification behaviours, since spurious wave decay as they propagate backwards (i.e. contrary to the sign of a ). We note that the amplification behaviour of a (physical or spurious) solution is encapsulated in the sign of the imaginary component of κh , because as z = exp (iκh ) , one has: where Re (·) and Im (·) respectively denote the real and imaginary parts of a complex number. Employing the above reasoning is useful because it bypasses the numerical difficulties of inverting the complex exponential function in order to access the imaginary component of κh . One can be sure that a physical root has been found by checking whether | z | < 1. A random guess for a physical solution can be of course z = r 1 exp (2 π ir 2 ) , where r 1 and r 2 are real numbers within [0, 1]. These same criteria can be used for a spurious root if one considers (16) with an inverted definition for z , i.e. z = 1 / exp (iκh ) = exp (−iκh ) . This sign change guarantees that the spurious root found through this strategy will have an absolute value smaller than unity. This spurious root will obviously have to be inverted subsequently if one is to store the actual spurious root of (16) for a given ϖh .
Once the (two) sequences of complex z values are obtained for the chosen range of ϖh , dispersion-diffusion curves can be generated by plotting Re (κh ) and Im (κh ) versus ϖh . However, while the latter can be readily obtained as Im (κh ) = − ln | z| , the former invariably requires using a (multi-valued) complex logarithm function, as Re (κh ) = −i ln (z/ | z| ) . It is therefore advisable that dispersion curves are carefully adjusted so as to avoid mistakes in complex phase estimates (usually seen as discontinuities on dispersion curves). We note that admissible corrections in Re (κh ) have to be a multiple of π . To verify whether any correction is in fact justifiable, it is worth checking the derivative of κh with respect to ϖh , which can be evaluated without ambiguity. This can be done by noting that where the rightmost derivative can be approximated numerically for a physical or spurious root z as in which is a sufficiently small variation in ϖh . In particular, the real part of d ( κh ) /d ( ϖh ) will indicate whether the slopes of dispersion curves have been evaluated correctly. With the above guidelines, fully continuous dispersion-diffusion curves have been obtained for physical and spurious modes. These are discussed in Section 3 .

Connections with high-order flux-reconstruction schemes
While the broader class of DSEMs presents an heterogenous mixture of discretisation strategies, each with different numerical properties, it is important to draw a common picture for those schemes sharing strong similarities and emphasize their connections. In particular, the spatial eigensolution analysis presented for DG in the previous section can be directly extended to some of the FR schemes investigated in [30,31] . In fact, if we consider the socalled FR DG scheme, that recovers the nodal DG method, we can easily assess that, as long as the Radau polynomials of order P + 1 , defining the correction functions of the FR DG , satisfy the discrete orthogonality condition to all polynomials of order P − 1 or lesssee for instance equation (20) and (21) in [31] -the connections between the DG method and FR DG schemes hold. Hence, if we exactly project the orthonormal Legendre polynomials (used for the spatial eigenanalysis in the previous section) to a nodal Lagrange basis, we retrieve the required discrete orthogonality condition for the Radau polynomials. This allows us to directly extend the results of the previous section to the FR DG scheme. Also, the results for the two methods are identical when using consistent integration of the nonlinearities to avoid aliasing or mitigate aliasing-driven instabilities (see [43] ).
We emphasize that, for the spatial eigensolution analysis of the DG method, we could have used other basis functions, such as nodal Lagrange polynomials. The associated results have been confirmed numerically by the authors, and in particular it has been found that dispersion and diffusion curves of DG formulations based on different basis functions match exactly (note that consistent integration of all terms, including mass matrices, is assumed [43] ). In this case, the connections with FR DG schemes based on Lagrange polynomials become more obvious. Therefore, the derivation in the previous section holds for the FR DG scheme exactly and can constitute the basis for subsequent spatial eigensolution analyses of the broader FR framework, including, for instance, the ESFR schemes presented in Ref. [33] .

Spatial eigensolution analysis results
We start by illustrating the eigencurves obtained for DGbased linear advection with standard upwinding ( β = 1 ). Fig. 1 shows dispersion and diffusion curves for various polynomial orders, namely P = 0 , . . . , 8 . The numerical wavenumbers are denoted in the plots as κ * , indicating these are different from exact (pure advection) values due to DG's truncation error. Each axis is normalised by the number of degrees of freedom per element, m = P + 1 . Note that positive diffusion is here represented by Im (κ * h ) > 0 , contrary to temporal eigenanalysis [24] , since κ and ω have opposite signs in (2) . We recall that, in the spatial analysis framework, one actually evaluates κ * instead of ϖ, the latter being merely interpreted as κ * in a temporal framework. Regarding the range of h = h/m adopted in the plots, it is worth mentioning that there is no natural limit for this range here (contrary to the periodicity limit of κh = π observed in temporal eigenanalysis). However, from a practical standpoint, higher frequency values will eventually require infeasibly small time steps to be resolved in a computation. Therefore, we here adopt h = 4 as a limit, which leaves a margin beyond the corresponding temporal framework limit of π . Lastly, we note that for compressible Navier-Stokes and Euler simulations, Roe's original Riemann flux [44] is the one that more naturally represents the standing upwind discretisation, as the penalty (property jump) term, cf. (9) , in its flux formula is multiplied by the unmodified eigenvalues of the system of equations. Moreover, although no spurious eigenmodes are found when standard upwind fluxes are employed (as discussed below), this result does not hold for general meshes in two-or three-dimensions.
The element-wise constant discretisation ( P = 0 ) shown in Fig. 1 is particularly inaccurate and strongly diffusive, which is not surprising. As the order is increased, the frequency range around h = 0 where dispersion and diffusion errors are negligible also increases. The accuracy on a basis of total degrees of freedom employed (DOF ∝ m / h ) improves significantly as P is increased from lower orders, but this improvement becomes less significant at higher orders (say, above P = 4 ). These trends are similar to those observed in temporal eigenanalysis, but marked differences can be noted as well. For the sake of comparison, Fig. 2 shows eigencurves equivalent to those of Fig. 1 , but obtained via the temporal analysis framework.
Firstly, the diffusion curves in Fig. 1 rise much more slowly as the frequency increases. In the temporal analysis, this rise is much sharper, especially at higher orders. Secondly, overall diffusion levels are notably smaller in Fig. 1 . The maximum values of Im (κ * h ) scale approximately as (P + 1) 2 in the temporal framework [24] , whereas in the spatial framework the diffusion at large frequencies decreases with P . Thirdly, comparison of dispersion curves reveals that, in temporally evolving problems, individual waves can propagate either faster or slower than they should, whereas in spatially developing problems they can only travel faster. Note that a * /a = /κ for both temporal and spatial approaches, where a * is the modified advection velocity for a given frequency/wavenumber. Dashed lines indicate the exact linear advection behaviour. Group velocities observed numerically will also display different behaviour since c/a = d / d κ, with c denoting the group velocity (although this formula is only valid when diffusion is negligible, see e.g. [45] ). Temporal analysis' results indicate the possibility of both positive and negative group velocities, whereas the spatial analysis shows a monotonic rise in group velocity as the frequency increases.
Regarding the frequency range of negligible dissipative effects, the values of ϖ 1% marking a diffusion threshold of 1% damping factor per DOF crossed (named "the 1% rule" in previous studies [24,26] ) has been found to be remarkably close to the corresponding values of κ 1% defined from temporal analysis. In fact, points per wavelength estimates have been found to agree to within 1-2%, which for most practical purposes means that one can rely on the 1% rule estimates and tables presented in [24,26] . What remains to be established is whether the 1% rule remains as effective for under-resolved turbulence computations of spatially developing flows -as mentioned above, the dissipation rise with frequency observed from spatial analysis is not particularly sharp, especially at higher polynomial orders. This assessment is however left for future studies.
We now turn our attention to upwinding effects on dispersion and diffusion characteristics and discuss how the eigencurves are modified when the parameter β changes. It should be noted that a spurious (reflected) mode exists in addition to the physical mode shown in Fig. 1 whenever β = 1, as first pointed out by [35] . Fig. 3 shows the eigencurves obtained when nearly central fluxes ( β = 0 . 01 ) are employed, for various polynomial orders P = 1 , . . . , 5 , where physical eigencurves are shown in blue and spurious ones in red. We note that decreasing β even more towards zero yielded only very slight changes in the eigencurves. It has been found that in the central flux limit the dissipation eigencurves of physical and spurious modes become mirror images of one another, with vanishing diffusion at certain frequency ranges, occasionally separated by "dissipation bubbles", practically as shown in Fig. 3 . If bubbles can be tracked across various orders, a given bubble seems to grow and move away from the origin when P increases. As these move towards larger frequencies, new ones seem to originate and take their place.
The dispersion curves of spurious waves are always negative, which is related to their backwards propagation behaviour (contrary to the advection velocity sign). Spurious modes can originate e.g. at an elemental interface standing normal to the propagation direction when a variable mesh spacing is used across the interface [35] . Although some of the dispersion curves of spurious modes shown in [35] started at Re (κ * h ) = π instead of at the origin, these have not been found here. We believe that results such as those from [35] may well have been caused by a miscalculation of phase estimate, as discussed at the end of Section 2.3 . It can be argued that, on the grounds of consistency, a steady interface boundary condition ( ϖ → 0) should not generate a spatially varying solution component ( Re (κ * h ) = 0 ) into the domain, physical or otherwise.
Another point to be noted in Fig. 3 is that finite diffusion levels are obtained with central fluxes, contrary to what is found in temporal eigenanalysis. Note that the overall levels of diffusion observed as β → 0 are comparable to those obtained with β = 1 .
This is a surprising result in itself and means that central fluxes can effectively yield non-zero diffusion for under-resolved computations (of spatially developing flows), when ϖ can be large. This damping takes place however in space, not in time. Hence, once a spatially-decaying solution incoming from an inlet boundary travels through the whole domain, an asymptotic state can be reached whose norm does not decay in time. Numerical experiments demonstrating the non-vanishing dissipation levels at the central flux limit are considered in Section 4 .
The qualitative behaviour of the curves in Fig. 3 is relevant for DG-based computations of spatially developing flows when a central (or nearly so) numerical flux is used. This can be the case, for example, in a DNS where one wishes to avoid upwind dissipation and rely solely on viscous diffusion, as done e.g. in [46] . Also, some corrections of Riemann solvers aimed at low Mach number applications can effectively reduce them to nearly central fluxes around low-speed flow regions, cf. [47] . The dispersion curves of the physical modes shown in Fig. 3 become more accurate on a per DOF basis as the polynomial order increases, as expected. Regarding diffusion, the presence of dissipative bubbles might induce potentially undesirable non-smooth features in under-resolved computations which have significant energy content at high frequencies. For example, in the computation of a turbulent wake at sufficiently large Reynolds number, one of such bubbles could cause a non-physical dissipation over a narrow frequency range within the inertial zone of the energy spectrum.
A final important point concerns the dissipation of spurious modes. As β → 0, there will exist a range of frequencies for which the reflected mode will be subjected to negligible damping. This seems to happen regardless of the polynomial order and can have serious practical consequences. In aeroacoustic simulations, for example, spurious modes reflected far downstream of a region of interest could return and corrupt the solution significantly. Also, in a turbulent wake, coarsening the mesh could cause spurious reflections to interact with incoming turbulent structures and affect not only solution quality, but possibly also numerical stability, in case the coarsening is sufficiently abrupt.
We now consider the limit of very large β values, named "hyper-upwinding". This limit is relevant, for instance, to nearly incompressible DG simulations of spatially developing flows performed with the local Lax-Friedrichs (LLF) solver, also known as the Rusanov flux [48] . For compressible Navier-Stokes or Euler computations, the LLF solver essentially replaces | u | with | u | + c in the momentum equations, where u is the relevant component of the convection velocity and c is the speed of sound. This effectively amounts to an over-upwind factor β = (| u | + c) / | u | , which scales as 1 + Mach −1 and hence increases without bound when Mach → 0. Some effects of LLF-induced hyper-upwinding in turbulence computations have been discussed in [26,49] . Fig. 4 shows the eigencurves obtained with β = 100 for various polynomial orders P = 1 , . . . , 5 . Although this value of β is particularly representative of Mach ≈ 0.01, the curves shown do not change noticeably as β is further increased. Interestingly, the eigencurves' behaviour in the hyper-upwinding limit is not very different than that observed for central fluxes. One of the conclusions is that, at low Mach number, LLF-based acoustics simulations should suffer from spurious modes reflected throughout the domain as mesh spacing is varied.
The eigencurves' behaviour at intermediates values of β is now considered. Figs. 5 and 6 show, respectively for cases P = 2 and P = 4 , how dispersion and diffusion curves change as the amount of upwinding employed varies from β 1 to β 1. Away from these two limits, the curves exhibit a more monotonic behaviour, without bubbles or kinks, both for physical and spurious modes. It is also interesting to see how the diffusion curves of spurious modes tend to become the symmetrical pairs of their physical counterparts in both limits. Additionally, as the standard upwind discretisation is approached, spurious modes become increasingly damped. In fact, when β → 1, their diffusion curves tend to minus infinity (not shown) and spurious modes eventually cease to exist, cf. Fig. 1 .
By putting together all the plots considered in the present section, a consistently complete picture of the qualitative behaviour of different DG-based discretisations in terms of diffusion and dispersion characteristics can be drawn. The current analysis also provides useful guidelines for the design of effective interface numerical fluxes required by discontinuous spectral element methods and can help practitioners to have a better understanding of the role of spurious/reflected modes and on how DG methods perform in different types of under-resolved simulations. In addition, the study presented here can be directly applied to the FR DG scheme and provides the basis to retrieve similar estimates for other high-order FR schemes.

Numerical experiments on linear advection
The results presented in Section 3 are now further verified by a series of numerical experiments for the one-dimensional linear advection equation, (1) , on a domain = [0 , 1] split into 100 equally spaced elements (i.e. h = 0 . 01 ) with an advection velocity a = 1 and an homogeneous initial condition u (x, 0) = 0 . The boundary conditions adopted are The latter is implemented by the assumption that, at x = 2 , solution states outside should always match those inside. This effectively causes the corresponding interface flux to match the interior flux value, amounting to a standard upwind formula at the boundary in question. Hence, no spurious reflections are expected from the outlet. In what follows, a sequence of test cases are set with different values of the inlet frequency ω.      Fig. 7 . The first value, ϖ 1 , is prior to the diffusion bubble and the solution should not be damped as it enters the domain, since the diffusion at this frequency is nearly zero. For ϖ 2 , dissipation should be significant as the corresponding normalised frequency ϖ 2 is located at the bubble, cf. Fig. 7 . For ϖ 3 , although the frequency is larger, a small damping is expected (slightly stronger than at ϖ 1 ). Further increasing the frequency from ϖ 3 to ϖ 4 causes no significant change in dissipation, therefore similar decay rates (in space) can be expected. Lastly, for the final frequency considered, ϖ 5 , the strongest damping of all cases tested can be anticipated.
All the predictions above have been confirmed and the numerical results obtained are shown in Fig. 8 . Simulations have been conducted until t = 2 , so that no transient effects are present and asymptotic oscillatory states have been reached. The time-step used was sufficiently small to guarantee that time-integration errors are negligible. It is also interesting to see that, for severely under-resolved waves (in particular, frequencies beyond the bubble), the shape of the numerical solution can differ significantly from that prescribed at the inlet. This can be clearly seen for frequencies ϖ 3 and ϖ 4 . For a discussion about the shape of the numerical solution within individual elements, the reader is referred to [41] .
In summary, the present section has confirmed the finite dissipation levels obtained from the eigenanalysis in the limit of central flux discretisations of spatially developing problems. Although the associated decay rates manifest in space rather than in time, significant damping can be observed for sufficiently under-resolved waves. Non-monotonic diffusion features in frequency space have also been illustrated, in particular the effect of dissipation bubbles.

Numerical experiments in under-resolved vortical flows
To evaluate the impact of numerical diffusion and dispersion characteristics in the context of under-resolved simulations of more practical flows, we used the two-dimensional compressible Euler equations with unsteady inlet boundary conditions mimicking a passive generator (physical screen) of eddies propagating downstream. While the problem chosen is inviscid, Euler equations have been widely used in the literature to study the properties and behaviour of numerical methods in the limit of vanishing viscosity  -i.e. very high-Reynolds numbers -see for instance [26,49,50] . In fact, (i) at sufficiently high Reynolds, numerical diffusion eventually overcomes molecular viscosity in under-resolved turbulence computations, hence Euler simulations of free flows can be representative of Navier-Stokes equations in the limit of very high Reynolds; and (ii) we try to mimic grid turbulence, but we don't have real turbulence as simulations are two-dimensional, however we do have small-scale vorticity, which is sufficient to our purposes here.
This problem relates to the previous analysis as it is possible to see how upwind-type ( β = 1 ) Riemann solvers -e.g. Roesuppress reflections at topological mesh changes and at the outflow, since no spurious mode is present. On the other hand, Riemann solvers that deviate from standard upwinding ( β = 1) -e.g.
local Lax-Friedrichs -do produce spurious modes for non-uniform meshes and at the outflow, confirming the presence of the spurious mode that can be seen in the spatial eigensolution analysis of Section 2 .
In the following, we first introduce the governing equations, the flow configurations adopted and the set of cases investigated. We successively present the results obtained and make some considerations that show how the numerical experiments performed confirm the analysis in Section 2 .

Governing equations and flow configuration
The compressible Euler equations can be written, in compact form, as where q is the vector of conserved variables and F (q ) = [ f 1 (q ) , f 2 (q )] T is the vector of fluxes that governs the transport of q on a domain . For the two-dimensional problem considered here, these quantities assume the following form where ρ is the density, u and v are the velocity components in the streamwise and cross-flow directions, respectively, p is the pressure and E is the total energy per unit volume. In this work, we considered a perfect gas law for which the pressure is related to the total energy by the following expression where γ = 1 . 4 denotes the adiabatic coefficient.
A sketch of the flow configuration and of the (coarser) mesh used in the test cases now considered is given in Fig. 9 . The mesh shown is divided into two blocks of different resolution, with the downstream block being coarser than the upstream one. While the dimensions of the upstream block were L up, x × L up, y = 12 π × 2 π , the dimensions of the downstream block were L down, x × L down, y The simulations were carried out on two meshes, referred to as Mesh 1 (or uniform mesh) and Mesh 2 (or non-uniform mesh), where the number of elements (i.e. the mesh resolution) of the upstream block was identical for the two cases and equal to N up, x × N up, y = 72 × 12 (resulting in square-shaped elements), while the mesh resolution of the downstream block was different. Mesh 1 was used as a reference and featured the same streamwise grid spacing for both blocks. The number of elements in its downstream block was equal to N mesh 1 down,x × N mesh 1 down,y = 48 × 12. Mesh 2 had instead fewer elements in its downstream block (coarser mesh resolution), allowing for a spurious reflected mode in the presence of a sudden mesh resolution change. For Mesh 2 , the number of elements in its downstream block was equal to N mesh 2 down,x × N mesh 2 down,y = 12 × 12.
Free-slip wall boundary conditions are imposed at y = ±π . Inflow boundary conditions are set through the outer state considered by the Riemann solvers at the faces of the elements adjacent to the inlet. The following conditions are prescribed respectively for the density, the momentum components in the streamwise and cross-flow directions, and the total energy per unit volume:  with a fourth-order Runge-Kutta scheme at a sufficiently small time step to ensure negligible errors from the temporal discretisation. Two Riemann solvers are considered for the test cases, namely, Roe's original solver [44] and the local Lax-Friedrichs (LLF) flux, also known as the Rusanov solver [48] . These are arguably the most common inviscid fluxes used in DG simulations. Although these two solvers are known to display comparable performance for most well-resolved DG computations, especially at higher polynomial orders, their behaviour can be drastically different for under-resolved computations. This has been shown for temporally evolving flows in [26] and is demonstrated here for spatially developing ones. Finally, two Mach numbers are tested in the simulations, namely, 0.3 and 0.03. The latter is relied upon to illustrate the effects of hyper-upwinding induced by the Lax-Friedrichs solver, as discussed in Section 3 . By taking into account the two grids considered, eight test cases are evaluated in total. All simulations have been conducted with spectral/ hp element code Nektar ++ [51] . Implementation details regarding the boundary conditions can be found in [52] . It is important to note that an increased number of quadrature points ( q = 2 m per element, per dimension) has been employed in all test cases to ensure consistent integration of the cubic non-linearities of the compressible Euler equations, see e.g. [14] . This should suppress polynomial aliasing errors caused by insufficient integration accuracy [26,43] .

Results
The discussion of the results in this section aims to verify the applicability of the linear spatial eigensolution analysis carried out in Section 2 to nonlinear problems. Specifically, we want to investigate aspects of solution quality of the schemes and numerical fluxes considered as well as the influence of the compressibility through the different Mach numbers employed. Of particular interest was the effect of the numerical fluxes and the interaction between turbulent-like structures and the sharp change in mesh resolution at x = 12 π . This aspect can trigger spurious upstream propagating waves (i.e. spurious reflections) that can affect the accuracy of the solution and potentially undermine the robustness of the computation. To evaluate all these aspects, we assessed particularly the vorticity flow fields, ω z = ( ∂v ∂x − ∂u ∂y ) , for the two Mach numbers, meshes and numerical fluxes adopted.
In Fig. 10 , we show contours of ω z for Mesh 1 and Mesh 2 , using LLF (top two snapshots) and Roe (bottom two snapshots) fluxes and with Mach = 0.03. It is possible to appreciate how the flow field is different for the two numerical fluxes adopted. In particular, the LLF Riemann solver produces spurious reflections which seem to interact nonlinearly with the incoming eddies, causing oscillations that are perceived as small-scale noise. These reflections are generated either at the outlet boundary or at the interior interface where mesh spacing is changed. As these decay exponentially in space, cf. equation (2) with κ complex, they apparently vanish rapidly, but actually never disappear entirely. As a result, spurious reflections can interact with and affect the turbulent wake over long distances upstream of where they originate.
The presence of spurious small-scale noise in LLF-based DG computations has also been observed in under-resolved Euler turbulence simulations in a recent study [26] . This behaviour is related to the hyper-upwinding of the LLF flux that is exacerbated at small Mach numbers. We therefore expect that, by increasing the Mach number, the differences observed between LLF and Roe will decrease. This is indeed the case, as shown in Fig. 11 , where we again compare the LLF (top two snapshots) and Roe (bottom two snapshots) numerical fluxes on the two meshes used and for Mach = 0.30. At this Mach number, the over-upwind character of the LLF flux is not so strong, and reflections decay much faster as they propagate backwards.
In both Figs. 10 and 11 , it is possible to observe spurious oscillations from the outflow when using the LLF numerical flux. This is consistent with the fact that the LLF numerical flux is not solving the full Riemann problem associated to the compressible Euler and Navier-Stokes equations. Therefore, it is not producing a characteristic treatment of the outflow boundary condition that is otherwise achieved using the Roe Riemann solver (see also [52] ). In addition, Figs. 10 and 11 , show how the use of a coarser mesh region downstream create numerical oscillations upstream of the meshresolution change (denoted by the white vertical line in the plots) for the LLF numerical flux. This is particularly evident for Mach = 0.03 -second snapshot from the top in Fig. 10 . This interior reflection is however much milder than that observed near the outlet boundary of the uniform mesh, where the property variation between interfaces (considered by the Riemann solver) is much more abrupt. To verify that reflections emanating from the interface between the two mesh blocks was in fact independent from outlet reflections, we performed some additional LLF-based computations (not shown) using characteristic outflow boundary conditions, cf. [52] for implementation details. While spurious reflections from the outlet boundary practically vanished, interior reflections were still clearly visible just upstream of the interface where mesh resolution changes. Interestingly, both (interior and outlet) reflections are almost completely suppressed using the Roe flux (that represents a standard upwind flux, i.e. β = 1 ). As a last remark, we note that the coarsening of the second mesh block alleviates outlet reflections, which is not surprising, since coarser grids help to damp eddies and therefore act somewhat like sponge zones. However, coarsening the grid induces interior reflections which can again affect the upstream solution. Use of gradual coarsening is therefore advised (with Roe-like fluxes) for more complex simulations. The question of how gradual the coarsening should be is of course an open question and motivates future studies.
The results outlined in this section complete the picture regarding the behaviour of DG methods for open vortical flows and support the spatial eigensolution analysis presented in Section 2 . More specifically, the poor performance of the LLF Riemann solver (especially) at low Mach numbers and the presence of spurious reflections both at the outflow and at the junction between the coarser and finer mesh regions emphasizes how one should pay close attention to the Riemann solver implementation. In fact, in the cases of "under-" and "over-upwinding" (that is β = 1), the spurious numerical mode highlighted in the eigensolution analysis triggers undesired features that undermines the correct propagation of wavelike solutions in under-resolved problems. This can in turn lead to inaccurate solutions and prone-to-failure simulations. Note also that, the Roe Riemann solver, being representative of the standard upwind numerical flux (i.e. β = 1 ), shows a significantly better behaviour -i.e. negligible reflections at the junction between coarse and fine mesh regions and at the outflow -than the LLF flux for all the cases illustrated in this section.

Conclusion
In this paper, we have presented the formal spatial eigensolution analysis of the DG method applied to linear advection in one dimension. We have subsequently assessed these results with numerical experiments conducted using the linear advec-tion equation and the two-dimensional compressible Euler equations (in a model problem that mimicked grid turbulence). In both scenarios, the test cases confirmed the insights provided by the eigensolution analysis. Throughout the study, special attention was given to the implications of the spatial eigenanalysis for DG-based under-resolved simulations of spatially developing transitional and turbulent flows. In particular, the results discussed in this paper should help researchers and practitioners interested in underresolved DNS / implicit LES approaches based on DSEM (i.e. discontinuous spectral element methods). This study aligns with various recent works that investigate the suitability of DSEM-based uDNS/iLES, considering aspects of solution quality and numerical stability [17,18,21,25,26,49,53,54] .
In the spatial eigensolution analysis, we compared different levels of upwiding, dictated by the parameter β that scales the property jump term in numerical flux formulas. We highlighted that Roe-type fluxes represent closely a standard upwind approach for which spurious modes and reflections are expected to be negligible, at least for Cartesian meshes aligned with the flow (resembling a 1D propagation state). In addition, we pointed out how Riemann solvers that depart from standard upwinding suffer from spurious reflections for under-resolved computations of spatially developing flows whenever variable mesh spacing is adopted. In particular, we stressed how the local Lax-Friedrichs (LLF) solver exhibits a "hyperupwinding" behaviour at low Mach numbers, thereby allowing for the existence of spurious reflections that are barely damped.
We have also pointed out that the eigensolution analysis framework used here can be easily extended to high-order FR schemes. Additionally, we have highlighted that the spatial eigensolution results for the DG method are identical to those one would obtain for the particular FR scheme recovering DG [30,31] .
The results obtained with the two-dimensional model problem designed to mimic grid turbulence at infinite Reynolds numbers confirmed the insights achieved from the spatial eigenanalysis in the non-linear setting of a vortex-dominated under-resolved computation. From the eigenanalysis and test cases addressed in this study, the following guidelines have been gathered and should facilitate the successful application of DSEM-uDNS/iLES approaches: 1. Favour Roe-type (standard upwind) fluxes over more simplistic solvers. This should help to suppress spurious reflections and small-scale noise, especially at low Mach numbers.
2. Avoid sharp mesh coarsening for simulations of spatially developing flows. Although Roe-type fluxes have suppressed reflections in all cases investigated here, this benefit is likely to be somewhat reduced for irregular meshes typical of more complex flow problems. 3. Favour high polynomial orders, such as P > 4, as dispersion and diffusion curves better capture the behaviour of the governing equations on a per degree-of-freedom basis. Much higher polynomial orders are however not advised as their inherently small diffusion levels might eventually be insufficient to stabilize computations at very high Reynolds numbers.
Regarding the last point above, the use of very high polynomial orders might however be a good strategy for those willing to employ DG discretisations of improved robustness. This includes formulations that rely on consistent/over-integration (as used in the present study) and skew-symmetric (also called split-form) discretisations, see e.g. [55] .
Future studies should investigate other common Riemann solvers, such as HLLC/HLL, to verify whether these resemble Roe/LLF fluxes, as showed, for instance, in [26] and analyse the broader class of DSEM approaches, including ESFR schemes. Also, the effect of different mesh topologies and types of elements should be addressed, in order to understand how effective Roetype solvers are in suppressing reflections for complex meshes (including curved elements). Finally, a systematic study regarding which strategy of gradual mesh coarsening is more appropriate to minimize reflections and improve robustness in different scenarios would also help promoting the industrialisation of SEMs.