On the behaviour of fully-discrete flux reconstruction schemes

In this study we employ von Neumann analyses to investigate the dispersion, dissipation, group velocity, and error properties of several fully-discrete flux reconstruction (FR) schemes. We consider three FR schemes paired with two explicit Runge–Kutta (ERK) schemes and two singly diagonally implicit RK (SDIRK) schemes. Key insights include the dependence of high-wavenumber numerical dissipation, relied upon for implicit large eddy simulation (ILES), on the choice of temporal scheme and time-step size. Also, the wavespeed characteristics of fully-discrete schemes and the relative dominance of temporal and spatial errors as a function of wavenumber and time-step size are investigated. Salient properties from the aforementioned theoretical analysis are then demonstrated in practice using linear advection test cases. Finally, a Burgers turbulence test case is used to demonstrate the importance of the temporal discretization when using FR schemes for ILES.


Introduction
Unsteady scale-resolving computational fluid dynamics (CFD) simulations, such as large eddy simulation (LES) and direct numerical simulation (DNS) rely on two different types of discretization when using the method of lines. A spatial discretization is required to obtain the divergence of the flux functions in space, and a temporal discretization is required to advance the solution in time. Therefore, the numerical properties of particular space-time discretizations rely on the underlying properties of both the spatial and temporal schemes that are being used, and how they interact. Recent advances have led to the development of several high-order accurate unstructured spatial discretizations, which are suitable for unsteady flows around complex geometries. These include the discontinuous Galerkin (DG) [1][2][3] and spectral difference (SD) [4,5] methods, amongst others. Recently, the flux reconstruction (FR) approach was proposed by Huynh [6] as a unifying framework for high-order unstructured numerical methods. FR does not describe a single scheme, but is in fact a family of schemes. Huynh [6] described several linearly stable schemes, including one equivalent to a collocation nodal discontinuous Galerkin method, henceforth referred to as FRDG, and another equivalent to an energy-stable SD method, henceforth referred to as FRSD. High-order methods, such as the FR, DG, and SD approaches, are particularly appealing for simulations of complex unsteady flows [7][8][9][10]. They have been found to provide more accurate solutions with fewer degrees of freedom and reduced computational cost relative to industry-standard second-order schemes [11]. The FR approach has been shown to be particularly accurate for scale-resolving simulations of complex turbulent flows, including DNS and implicit LES (ILES) [7,11,8,9,12].
Unsteady simulations using the FR, DG, and SD approaches are typically performed with high-order accurate temporal discretizations. Explicit Runge-Kutta [13,14] (ERK) schemes are, perhaps, the most popular choice for problems that are numerically non-stiff. They are able to achieve high levels of accuracy with a low cost per timestep. However, ERK schemes are only conditionally stable [15], meaning the time-step size is limited. Problems with significant numerical stiffness often necessitate the use of A-stable and L-stable implicit schemes. Popular methods include the implicit Runge-Kutta and backward differentiation formula (BDF) schemes [15]. Since there are no A-stable BDF schemes with order greater than two [15], we are only interested in implicit Runge-Kutta schemes for the current study. Of the available implicit Runge-Kutta schemes, the singly diagonally implicit Runge-Kutta (SDIRK) methods are particularly appealing [15]. They require the solution of a relatively small system of equations at each stage, when compared to more general implicit Runge-Kutta methods. They also allow for a quasi-Newton approach to be employed for non-linear problems, whereby the same Jacobian matrix can be used for several stages within a time-step, or across several time-steps, reducing the number of costly preconditioning operations [7,9]. The most popular methods for scale-resolving simulations of turbulent flows appear to be third-and fourth-order ERK schemes and third-order SDIRK schemes [16,9,17,12]. We also consider a popular five-stage fourth-order SDIRK scheme in the current study [15].
Previous studies of the FR approach have investigated the behaviour of various energy-stable FR (ESFR) schemes using semi-discrete von Neumann analysis [6,18,19]. This type of analysis provides insights into the dispersion, dissipation, numerical error, and group velocity properties of a semi-discrete scheme in the absence of temporal errors. For example, Huynh [6] investigated the behaviour of semi-discrete FRDG, FRSD, and similar schemes using von Neumann analysis. Similarly, Vincent et al. [18] and Vermeire and Vincent [19] performed von Neumann analysis for families of provably-stable semi-discrete ESFR schemes. However, none of the aforementioned studies investigated the behaviour of fully-discrete FR schemes. Yang et al. [20] showed that fully-discrete DG schemes can exhibit significantly different behaviour from their semi-discrete forms. They provided analytical expressions for the dispersion and dissipation behaviour of several schemes, with emphasis on their behaviour in the asymptotic lowwavenumber limit. The current study will follow the approach of Yang et al. [20] to investigate the fully-discrete behaviour of FRDG, FRSD, and a particularly appealing ESFR scheme identified by Vermeire and Vincent [19]. For the current study we restrict our analysis to FR schemes of polynomial degree k = 4, which is a common choice for scale-resolving simulations of turbulent flows [9,17,12,21,22]. However, the procedure described herein is readily generalized to all polynomial degrees. It is expected that results from this study will inform the use of fully-discrete FR schemes for DNS and LES of unsteady turbulent flows. The behaviour of fully discrete schemes at the highest-resolved wavenumbers is of particular interest for ILES, since ILES relies on numerical dissipation at these wavenumbers to dissipate energy from the flow, acting as a subgrid scale model [9]. Therefore, while Yang et al. [20] were primarily interested in the behaviour of schemes in the asymptotic limit, we are interested in the behaviour of FR schemes for all resolvable wavenumbers and implications on future large-scale ILES simulations. This paper is structured as follows. In Section 2 we present a formulation of the FR approach following Huynh [6]. We will also summarize the correction function formulation for polynomial degree k = 4 provided by Vincent et al. [23], which yields stable-symmetric-conservative ESFR correction functions. In Section 3 we will describe fully-discrete von Neumann analysis, which follows from the work of Yang et al. [20]. We will then present results for several popular FR schemes with ERK and SDIRK time-stepping. In Section 4 we will investigate the behaviour of these fully-discrete schemes via numerical experiments including advection of a sine wave, a marginally-resolved Gaussian bump, and ILES of viscous Burgers turbulence. Finally, we will present conclusions based on our analysis and numerical experiments, along with implications for the application of FR to large-scale ILES.

Preliminaries
The flux reconstructions (FR) approach, first introduced by Huynh [6], is used to spatially discretize a 1D conservation law of the form in a domain Ω, where u = u(x, t) is the solution with a given initial distribution u(x, 0) = u 0 , t is time, f = f (u) is the flux, and x is the spatial coordinate. We split the domain Ω into a mesh where Ω n is one of N e elements in the domain. For the sake of simplicity, we prefer to perform all operations in a computational space where each element exists on the interval I ∈ [−1, 1]. For the current study we use the linear mapping where x n is the mesh node corresponding to the left face of Ω n , Γ n is the mapping function, and ξ is the location in reference space. This mapping can be inverted according to which is also linear. This allows us to construct the transformed solution in computational space according tô and a transformed flux in computational space according tô The system of equations in computational space can then be written as Following Huynh [6], the solution within each element is represented by a degree k polynomial, which is allowed to be discontinuous at the interface between elements. This polynomial is supported by nodal basis functions generated at k + 1 solution points. Therefore, the solution within each element in computational space can be approximated aŝ whereû δ n =û δ n (ξ, t) is the polynomial representation of the solution within an element,û δ n,l =û δ n,l (t) is the value of the solution at solution point l, andφ l =φ l (ξ ) is its corresponding nodal basis function in reference space. For the one-dimensional case these basis functions are the Lagrange polynomialŝ 2.2. Constructing the flux Following Huynh [6], a polynomial representation of a discontinuous flux functionf δ D n =f δ D n (x, t) can be constructed using the same polynomial basis as the solution according tô where the superscript D denotes that this flux, like the solution, is allowed to be discontinuous at the interface between elements andf δ n,l =f δ n,l (t) is the flux evaluated at the solution points. We notice that the flux between elements must be continuous to ensure global conservation. However, since the solution is generally discontinuous between elements, so to will be the flux. Huynh [6], proposed that we generate a globally C 0 continuous flux by applying a correction denoted byf δC n to the discontinuous flux in each element such thatf δ n =f δ D n +f δC n , wheref δ n is the globally continuous flux function referred to in Eq. (7). The flux correction is computed according tô are common interface fluxes computed at the flux points between elements by an appropriate Riemann solver using the interpolated values u − L , u + L , u − R , and u + R of the solution from the neighbouring elements at each edge. The functions g L = g L (ξ ) and g R = g R (ξ ) are so-called correction functions, which are degree k + 1 polynomials with the constraints The gradient of the continuous flux function can then be found according to which is of degree k and in the same polynomial space asû n . The behaviour of an FR scheme depends on the choice of correction functions g L and g R for the left and right hand boundaries, respectively. Huynh [6] originally introduced several correction functions including one that recovered a collocation-based nodal DG scheme, one that recovered an energy-stable SD scheme, and one that recovered a socalled g 2 scheme. Vincent et al. [24] later introduced a single-parameter family of energy-stable flux reconstruction schemes, whose semi-discrete dispersion and dissipation behaviour was later studied [18].

Extended range of energy stable schemes
Vincent et al. [23] later introduced a multi-parameter range of energy stable flux reconstruction (E-ESFR) correction functions [24]. Like the single parameter schemes, all of the E-ESFR schemes are provably stable for linear advection. Schemes with k = 1 and k = 2 are single parameter families of schemes. However, for k ≥ 3 the E-ESFR schemes were found to have multiple parameters.
Here we will summarize the process used to generate the E-ESFR schemes. First we define D as the nodal differentiation operator according to whereφ j is the corresponding nodal basis function at solution point j. Similarly, the Vandermonde matrix V can be formulated as were L j (ξ i ) is a Legendre polynomial of degree j, normalized to unity at ξ = 1. Using the nodal differentiation operator and Vandermonde matrix, we can generate the modal differentiation operator Also, a modal mass matrixM can be generated according tõ We denote gradients of the correction functions given as vectors of modal coefficients for the left and right hand side, respectively, asg ξ L andg ξ R . These can be converted to nodal coefficients by and Vincent et al. [23] proved that FR correction functions are stable provided whereL andQ is a real square matrix of dimension k + 1 that satisfies They are also symmetric if whereJ and conservative if Following this procedure yields an under-determined system of equations. Solutions to this system yield linearly stable, symmetric, and conservative correction functions. For k = 4 the extended range of schemes has two parameters q 0 and q 1 [23]. Specifically, the most general form of the correction function that satisfies the given constraints is [23] where and 8 15 Also, due to the symmetry condition, By setting q 1 = 0 we can recover the original single parameter family of schemes. Subsequently, setting q 0 = q 0DG = 0 recovers the FRDG scheme and q 0 = q 0S D = 8/45 recovers the FRSD scheme. Vermeire and Vincent [19] identified a scheme with q = [q 0 , q 1 ] = [0.05, 0.01] that was particularly stable for marginally-resolved ILES of the Taylor-Green vortex test case. In this study we undertake fully-discrete analysis of these three schemes, which have been used previously for ILES [9,19,21].

Fully-discrete analysis
Fully-discrete von Neumann analysis was performed by Yang et al. [20] for Runge-Kutta DG (RKDG) schemes and Lax-Wendroff DG (LWDG) schemes with a focus on the asymptotic (well-resolved) regime. In the current study we extend their approach to FR schemes and all resolvable wavenumbers. Specifically, we consider the linear advection equation where L x is the length of the domain. We discretize the domain into elements of equal length h. Without loss of generality, we assume for this analysis that h = 1. We can express the nodal coefficients of the solution nodal basis u δ n,l as a vector u δ n for each element where n = 1 . . . N e . We also use an upwind flux at the interface between elements. A partial discretization of the FR scheme can be written as [6] where is an array of the nodal basis coefficients u δ n,l at the solution points and is the spatial discretization operator. Each sub-block L i, j generates the contributions of element j to the residual in element i. By applying an appropriate temporal discretization, the solution can be updated according to where is the full space-time discretization of the scheme and, for linear equations, is independent of the current solution. It follows that the behaviour of any scheme depends on both the spatial and temporal discretizations that are employed in the generation of M.
For any element u δ n the solution at the next time step can in general be determined from We seek solutions of the form where is α a non-zero constant vector of length k + 1, κ is a prescribed wavenumber, and ω is a numerically obtained frequency. A solution to this problem is given by the exact dispersion relation ω = κ, where ω is the frequency. By splitting the numerical frequency ω into real and imaginary parts ω r and ω i , respectively, we can write It follows that when ω i < 0 the scheme has numerical dissipation, and when ω r ̸ = κ the scheme has numerical dispersion error. Furthermore, we expect the scheme to remain stable for any time-step τ provided ω i ≤ 0. By substituting Eq. (47) into Eq. (46) and dividing by e iκ x n e −iωt we can write which can be rewritten compactly as Solutions to Eq. (50) are obtained by solving the classical eigenvalue problem for specified values of κ. Therefore, for a given κ we first compute the eigenvalues λ = e −iωτ of N . Then, the numerical frequency for each eigenvalue ω can be found as By then splitting ω into real and imaginary parts we can analyse the dispersion and dissipation properties, respectively, of the fully-discrete scheme. Each prescribed wavenumber κ will yield k + 1 admissible eigenvalue/eigenvector pairs as solutions to this problem [18]. According to Van den Abeele, each of these admissible solutions belongs to a true wavenumber in the range of −(k + 1)π ≤ κ ≤ (k + 1)π [25]. In this study we follow the approach of Vincent et al. [18] to identify true wavenumbers for each of these k + 1 admissible eigenvalues/eigenvectors. Unlike, Yang et al. [20], who were primarily interested in the low-wavenumber behaviour of DG schemes, we are interested in the behaviour of FR schemes for all permissible wavenumbers and possible implications in the context of ILES. For all schemes there exists a cut-off wavenumber κ c = (k + 1)π , which represents the maximum wavenumber that can be represented by the solution points. Similarly, there exists a cut-off for the real component of ω, which is given as ω c r = ±π/τ . This arises from the finite time sampling performed by the temporal scheme, which imposes a limit on the minimum and maximum permissible frequencies. The group velocity, which is the speed at which features of a given wavenumber propagate in space, is obtained from the slope of the dispersion curve and is given as v g = ∂ω r /∂κ [26]. We can also define the numerically obtained local and global orders of accuracy of the scheme as and respectively, where κ δ R is a specific wavenumber and E T = |ω(κ) − ω(κ)|. This approach approximates the order of accuracy of the scheme in the vicinity of κ δ R .

Explicit schemes
For a classical three-stage third-order ERK scheme the fully-discrete operator matrix is [13,14] where I is the identity matrix. Plots of dispersion, dissipation, group velocity, and numerical error are shown in Fig. 1 for the FRDG scheme with τ = 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0τ CFL , where τ CFL is the maximum stable time-step. These are chosen since large time-steps are typically preferential in minimizing the computational cost of a simulation. A total of N e = 4 elements are used to avoid imposing periodicity on the fully-discrete stencil. A reference profile is given for a semi-discrete scheme, which assumes no error from the temporal scheme. Similar plots are shown in Appendix A for the FRSD scheme and the q = [0.05, 0.01] scheme. The dispersion profiles for the all three methods shows that their behaviour at low wavenumbers is only weakly dependent, in an absolute sense, on the time-step size. However, in the vicinity of κ ≈ 4π increasing the time-step significantly reduces ω r . In contrast, near the maximum resolvable wavenumber κ ≈ κ c , increasing the time-step significantly increases ω r beyond the line κ = ω r , which describes the exact dispersion relation. This suggests The dissipation profiles for all three methods show that numerical dissipation increases in the vicinity of κ = 4π when the timestep is increased. However, the numerical dissipation initially increases at κ = κ c for small timestep sizes, and then decreases until it matches ω i = 0 at τ = τ CFL , defining the stability limit of the scheme. All three schemes initially exceed the stability limit between 4π < κ < 5π as the time-step τ is increased. Also, when τ ≈ τ CFL there is a plateau at high-wavenumbers with little numerical dissipation. This suggests that high-wavenumber features could be allowed to propagate nearly undamped when the timestep is close to the stability limit. This type of behaviour is not favourable for ILES, where numerical dissipation is being relied upon to drain energy from these under-resolved scales [9].
The group velocity for all three schemes shows nearly constant v g = 1 for well-resolved wavenumbers, regardless of the time-step size. The group velocity is uniformly decreased in the vicinity of κ = 4π with increasing time-step size. The behaviour of the group velocity near κ = κ c is less clear, and no obvious trend is apparent as a function of the timestep size for any of the schemes. However, one important conclusion is that the absolute value of the group velocity at high wavenumbers can significantly exceed the exact relation v g = 1. These results suggest that it is possible for high-wavenumber structures to propagate rapidly through the domain, several times faster than the physical wavespeed.
The numerical error E T as a function of wavenumber is found to have a strong dependence on the timestep size for all three schemes. For the FRDG scheme there is a small region near κ ≈ 10 that is relatively unaffected by the timestep size, however, the error of both the low-and high-wavenumber features are strongly influenced. It is also clear that the local order of accuracy of the FRDG scheme degrades from the expected B T = 2k + 2 to B T = 4 for well-resolved waves, which is the local order of accuracy of the RK3 scheme. Observations for the FRSD and the q = [0.05, 0.01] schemes are similar. Both achieve their expected rates of super-convergence for moderately-resolved wavenumbers of B T = 2k + 1 and B T = 2(k − 1) + 1, respectively. For well-resolved wavenumbers they are also reduced to B T = 4, matching the expected local order of accuracy of RK3.

RK4
For a classical four-stage fourth-order ERK scheme the fully-discrete operator matrix is [13,14] Plots of dispersion, dissipation, group velocity, and numerical error are shown in Fig. 2 for the FRDG scheme with N e = 5 and τ = 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0τ CFL . These are chosen since large time-steps are typically preferential in minimizing the computational cost of a simulation. A reference profile is given for a semi-discrete scheme, which assumes no error from the temporal scheme. Similar plots are shown in Appendix A for the FRSD and q = [0.05, 0.01] schemes.
The dispersion profiles for all three schemes show that ω r decreases with increasing τ for κ 4π . It then has a sharp inflection point near κ ≈ 4.5π . As τ is increased the frequency increases until it exceeds the ω c r limit of the RK4 scheme. Features with higher wavenumbers are, therefore, found to be aliased into the negative frequency range for large time-steps. Unlike the RK3 scheme, ω r = 0 at the cut-off wavenumber κ c for all time-step sizes, suggesting that high-wavenumber components will not oscillate in time.
Dissipation profiles for all three schemes show that increasing the time-step size uniformly decreases the amount of dissipation at the cut-off wavenumber κ c . This continues until it exceeds the stability condition at τ = τ CFL . In the vicinity of κ ≈ 4.5π the amount of dissipation is increased for relatively small time-step sizes. However, further increasing the time-step size causes the numerical dissipation to then uniformly decrease at all of the smallest resolved scales. These results suggest that when the time-step approaches τ CFL , the fully-discrete schemes could permit nearly undamped high-frequency waves similar to the RK3 scheme.
The group velocity behaviour of all three schemes is similar. Increasing the time-step size tends to decrease the group velocity for well-resolved wavenumbers, up to κ ≈ 4π . Above this wavenumber the behaviour is less clear. However, it is apparent that large group velocities can be expected to occur at these highest resolved wavenumbers, particularly in the vicinity of κ ≈ 4.5π . The peak group velocities for the FRDG and q = [0.05, 0.01] schemes are similar, while the peaks for the FRSD schemes are somewhat reduced in comparison. However, all of these Table 1 Butcher tableau for the SDIRK 3,3 scheme. significantly exceed the v g = 1 expected from the exact dispersion relation. These results suggest that it is possible for high-wavenumber structures to propagate rapidly through the domain when using the RK4 scheme with large time-steps.
Plots of the error E T show similar behaviour to previous results with the RK3 scheme. The error at low wavenumbers is dominated by the temporal term, which achieves the expected local order of accuracy of 5. There is a wide band around κ ≈ 10 that is relatively unaffected by the timestep size for all of the schemes. The FRDG, FRSD, and q = [0.05, 0.01] schemes achieve their expected local orders of accuracy of B T = 2k + 2 and B T = 2k + 1, and B T = 2(k − 1) + 1 in this region. The error at the highest resolved wavenumbers near κ c shows a strong dependence on the time-step size, as would be expected based on the behaviour of the dispersion and dissipation profiles described previously.

Implicit schemes
The Butcher tableau for the SDIRK 3,3 scheme is given in Table 1. Since this scheme is implicit [15], it requires the solution of several linear systems of equations per time-step. This requires several matrix inversions to generate M, unlike the explicit RK3 and RK4 schemes. As a consequence, M is in general dense for implicit schemes, suggesting that the influence of upstream boundary conditions cannot be entirely eliminated. However, the magnitude of the values in each off-diagonal block decreases with increasing distance from their corresponding diagonal block, meaning distant elements are only weakly coupled. For the time-step sizes used in the current study it was found that when N e ≥ 32 all terms in the most distant off diagonal blocks were on the order of machine precision. Hence, N e = 32 was used for all the SDIRK 3,3 analysis to approximate an effectively infinite domain. It should be noted that for larger time-step sizes a greater number of elements may be required, since the coupling with upstream elements will be greater. Plots of dispersion, dissipation, group velocity, and numerical error are shown in Fig. 3 for the FRDG scheme with non-dimensional time-step sizes τ = 1/32, 1/16, 1/8, 1/4, 1/2, and 1. It should be noted that these are significantly larger than the time-steps used for the explicit schemes, which were a fraction of their corresponding τ CFL . This is because implicit schemes are typically used with larger time-steps due to their improved linear stability properties and to offset the increased computational cost per step. Results for the FRSD and q = [0.05, 0.01] schemes are provided in Appendix A.
Plots of dispersion behaviour for all three schemes show similar trends. Increasing the time-step size tends to decrease the frequency ω r for all of the schemes at low wavenumbers. Further increasing the time-step size reduces the wavenumber until the solution exceeds the ω c r limit of the temporal scheme. This results in ω r being aliased into the negative frequency range. For smaller time-step sizes the dispersion relation does not hit the critical ω c r limit, and so continues up to κ = 5π .
Plots of the dissipation error for all three schemes show that increasing the time-step size tends to increase the amount of numerical dissipation in the well-and moderately-resolved wavenumber regions. It is expected that this will result in significant damping in these regions. In contrast, increasing the time-step size is observed to decrease the amount of numerical dissipation at the highest resolved wavenumbers.
Plots of the group velocity show that increasing the time-step size tends to decrease the magnitude of v g below the exact dispersion relation. This is most pronounced in the moderately-resolved region. It is expected that this will result in different wavespeeds for well-and moderately-resolved structures, with moderately-resolved structures propagating more slowly. Similar to the RK3 scheme, large group velocities are observed in the vicinity of κ = 5π , which suggests that high wavenumber structures could be advected at large speeds. However, unlike the RK3 scheme, there is still significant numerical dissipation in this region when using large time-steps suggesting these structures will be rapidly dissipated.  Table 2 Butcher tableau for the SDIRK 5,4 scheme.
Plots of the numerical error E T show that for small time-step sizes there is a large band around κ ≈ 10 that achieve the expected semi-discrete local orders of accuracy of B T = 2k + 2 and B T = 2k + 1, and B T = 2(k − 1) + 1 for the FRDG, FRSD, and q = [0.05, 0.01] schemes, respectively. Increasing the time-step size increases the amount of error in the low-and high-wavenumber regions, until the error for the fully-discrete scheme is dominated by the temporal error of the SDIRK 3,3 scheme with B T = 4. Large jumps in error are also apparent where the methods exceed the ω c r limit of the temporal scheme. 5,4 Similar to the SDIRK 3,3 scheme, we use N e = 32 elements for analysing the behaviour of fully-discrete schemes using SDIRK 5,4 described by the Butcher tableau provided in Table 2 [15]. Plots of dispersion, dissipation, group velocity, and numerical error are shown in Fig. 4 for the FRDG scheme with non-dimensional time-step sizes τ = 1/32, 1/16, 1/8, 1/4, 1/2, and 1. Results for the FRSD and q = [0.05, 0.01] schemes are provided in Appendix A.

SDIRK
The dispersion profiles show good agreement between the fully-discrete and semi-discrete schemes for low wavenumbers and all time-steps considered. Similar to the SDIRK 3,3 scheme, large time-steps result in the dispersion relation exceeding ω c r and getting aliased into negative frequencies. Similar to the RK4 scheme, all of the fully-discrete SDIRK 5,4 schemes are found to have ω r = 0 at κ c for the time-step sizes considered in the current study.
Plots of dissipation error show that, at low-wavenumbers, it is less sensitive to the time-step size than the SDIRK 3,3 scheme. However, at high wavenumbers the numerical dissipation still decreases with increasing time-step size. These results suggest that low-wavenumber content will be exposed to significantly less dissipation when compared to the SDIRK 3,3 scheme, while the high-wavenumber content will still have significant dissipation.
Plots of group velocity show that increasing the time-step size decreases the group velocity for all schemes in the low-and moderate-wavenumber regimes. This results in under-prediction of the group velocities for each of the three schemes. It is expected that this will result in different wavespeeds for well-and moderately-resolved structures, with moderately-resolved structures propagating more slowly. Combined with the relatively low dissipation at these wavenumbers, it can be expected that these slow-moving structures will be relatively undamped. Similar to the RK4 scheme, large group velocities are observed in the vicinity of κ = 4π , which suggests that high wavenumber structures could be advected at large speeds. However, unlike the RK4 scheme, there is still significant numerical dissipation in this region when using large time-steps.
Plots of the numerical error E T show that for small time-step sizes there is a large band around κ ≈ 10 that achieve the expected semi-discrete local orders of accuracy of B T = 2k + 2 and B T = 2k + 1, and B T = 2(k − 1) + 1 for the FRDG, FRSD, and q = [0.05, 0.01] schemes, respectively. Increasing the time-step size increases the amount of error in the low-and high-wavenumber regions, until the error for the fully-discrete scheme is dominated by the temporal error of the SDIRK 5,4 scheme with B T = 5. Large jumps in error are also apparent where the methods exceed the ω c r limit of the temporal scheme.

Sine wave
Results from von Neumann analysis suggest that while moderately-resolved waves can achieve super-convergence, it is expected that the numerical error of well-resolved waves is limited by the order of accuracy of the temporal scheme. To further study this, we consider the relative error of several fully-discrete schemes for an advecting sine wave. We consider a domain of size Ω ∈ [0, 20] after Vincent et al. [18]. A time-dependent boundary condition is used for the upstream boundary u δ (0, t) = sin(π t/2), which has an exact real frequency ω = π/2. No boundary condition is required at the downstream location, since a fully upwind flux is used, and the solution was initialized according to u δ (x, 0) = 0. By uniformly refining the mesh the expected non-dimensional wavenumber κ is determined from κ = ωh. Therefore, by refining the mesh with a constant CFL number we can explore the behaviour of a global error metric as a function of κ.
We use the same approach as Vincent et al. [18] for assessing numerical error. Since the forcing function is cyclical, it advects four units of length downstream before repeating. Therefore, we can evaluate an exact L 2 norm of the difference between the solution u δ at x ∈ [0, 4] and x ∈ [12,16] according to where E is the error metric. The integration for this error metric is performed analytically and piecewise on each element. This approach isolates errors associated with the discretization approach, as opposed to any errors associated with boundary or initial conditions [18]. Each simulation is run to a final time t = 500 to allow any initial transients to be eliminated. Tables of the numerical error for the k = 4 FRDG, FRSD, and q = [0.05, 0.01] schemes are shown in Table 3, Table 4, and Table 5, respectively. Results are included for κ = 2π through κ = π/16 for the RK3, RK4, SDIRK 3,3 , and SDIRK 5,4 temporal schemes. Bold numbers represent the maximum achieved global orders of accuracy A T . We find that all schemes achieve some level of super-accuracy (A T > k + 1) for moderately-resolved wavenumbers. For example, with RK4 and SDIRK 5,4 the observed maximum global orders of accuracy are approximately A T = 2k + 1, A T = 2k, and A T = 2(k − 1) for the FRDG, FRSD, and q = [0.05, 0.01] schemes, respectively. These match the expected orders of accuracy from semi-discrete analysis [6,18,19], and also match the local order of accuracy B T for small time-step sizes shown by the von Neumann analysis in this study. However, for well-resolved wavenumbers the observed global order decays to the global order of the time-stepping scheme. For example, with RK3 and SDIRK 3,3 the fully-discrete schemes approach a global order of A T = 3 and with RK4 and SDIRK 5,4 the fully-discrete schemes approach a global order of A T = 4. These results confirm observations from the fully-discrete von Neumann analysis. The order of accuracy of fully-discrete FR schemes can be limited by the order of accuracy of the temporal scheme, particularly in the limit of small wavenumbers. This also agrees with the observations of Yang et al. [20], who predicted this type of asymptotic behaviour with DG schemes.

Marginally-resolved Gaussian bump
For typical grid-filtered LES and ILES simulations using high-order schemes, turbulent structures exist down to the κ c resolution limit of the spatial discretization. Therefore, these types of simulations are inherently marginallyresolved and the behaviour of fully-discrete schemes across all resolvable wavenumbers is of interest, rather than just their behaviour in the well-resolved regime. To study this behaviour we consider marginally resolved advection of a Gaussian bump. We consider a domain of size Ω ∈ [−20, 20] with periodic boundaries and N e = 40 elements, each of length h = 1. The exact initial condition is a collocation projection of u δ (x, 0) = e −10x 2 with Gauss points used as the solution points. The configuration of this test case is similar to that used by Trefethen to study the dispersion behaviour of finite-difference schemes [26].

Explicit schemes
Since explicit schemes are typically run close to their stability limits, we are interested in the behaviour of the RK3 and RK4 schemes in the vicinity of τ ≈ τ CFL . Therefore, we run simulations at the maximum stable time-step τ = τ CFL and at τ = 0.5τ CFL . Each case was initially run as close as possible to, without exceeding, a final time of t = 1 using an integer number of time-steps. They were then run to a final time of t = 4000 to investigate error for lower-wavenumber features, which yields a total of 100 advection cycles through the domain.
Results are shown in Fig. 5 with RK3 and in Fig. 6 with RK4 for the FRDG, FRSD, and q = [0.05, 0.01] schemes, respectively, at t ≈ 1. It is clear that all of the solutions using both temporal schemes possess high wavenumber Table 3 Error and order of accuracy for an advecting sine wave using a k = 4 FRDG scheme with τ = 0.05h. Bold numbers denote maximum observed order of accuracy.  Table 4 Error and order of accuracy for an advecting sine wave using a k = 4 FRSD scheme with τ = 0.05h. Bold numbers denote maximum observed order of accuracy.  Table 5 Error and order of accuracy for an advecting sine wave using a k = 4 q = [0.05, 0.01] scheme with τ = 0.05h. Bold numbers denote maximum observed order of accuracy. and high group-velocity features that are advected out in front of the Gaussian bump. These structures are nearly undamped, and propagate large distances downstream. For both RK3 and RK4, the FRDG scheme appears to have the highest group velocity for these waves when τ = τ CFL , followed by the q = [0.05, 0.01] scheme, and finally followed by the FRSD scheme. This order, and the maximum observed wave speeds, agree well with the maximum observed group velocities from von Neumann analysis for each of the schemes in the limit τ = τ CFL . Therefore, we have demonstrated that the high group velocity structures predicted previously with von Neumann analysis are observed in practice. Also, in the limiting case of τ ≈ τ CFL they are nearly undamped and can advect far downstream without significant dissipation. Results are shown in Fig. 7 using RK3 and RK4 for the FRDG, FRSD, and q = [0.05, 0.01] schemes, respectively, at the final solution time t = 4000. Decreasing the time-step size was found to reduce the amount of dissipation error for all spatial discretizations using both RK3 and RK4. This agrees with observations from the von Neumann analysis of all three spatial discretizations, which showed reduced dissipation in the well-resolved region for decreasing timestep size. Another observation from von Neumann analysis is that the error in the group velocity of well-resolved waves increases with increasing time-step size. This is also observed in the numerical experiments, where non-physical oscillations ahead of and behind the main bump profile exhibit increased group velocity error with increasing time- step size. We also observe that the advection error for the FRDG scheme is more biased towards fast travelling waves, while the error of the FRSD and q = [0.05, 0.01] schemes are more biased towards slow travelling waves, which are formed in the lee of the Gaussian bump.

Implicit schemes
Implicit schemes are typically used for stiff systems of equations, where relatively large time-steps are taken to offset the increased computational cost per step relative to explicit schemes. Therefore, we are interested in the behaviour of SDIRK 3,3 and SDIRK 5,4 for relatively large time-step sizes. We run simulations with the upper-range of time-steps used for the von Neumann analysis, τ = 1/8, 1/4, 1/2, and 1. Each simulation is run to a final time t = 40, which corresponds to exactly one cycle through the domain.
Results are shown in Fig. 8 for both the SDIRK 3,3 and SDIRK 5,4 schemes. We observe that increasing the timestep with SDIRK 3,3 tends to dissipate the Gaussian bump profile for all of the spatial discretizations. It also tends to reduce its advection speed, evident from the downstream shift in the solution peak relative to the reference solution. Furthermore, there is a trailing oscillation behind the bump, which shifts further downstream with increasing τ . These results are consistent with observations from von Neumann analysis, including the reduced group velocities and increased dissipation at moderate wavenumbers when using large time-step sizes. Results with SDIRK 5,4 show that large amplitude waves are produced and trail behind the Gaussian bump when using large time-steps. While both the SDIRK 3,3 and SDIRK 5,4 schemes reduced the group velocity of moderately-resolved waves in the context of von Neumann analysis, the SDIRK 5,4 scheme did not introduce significant additional numerical dissipation. Therefore, these waves are able to persist for much longer periods of time than with the SDIRK 3,3 scheme. Another observation is that for small time-step sizes there is a significant difference between the behaviour of the FRDG, FRSD, and q = [0.05, 0.01] schemes. The FRDG tends to have a forward dispersion error for high-wavenumber structures, with large peaks being pushed in front of the Guassian bump. The FRSD scheme appears to have a more rearward dispersion error, with oscillations being reduced in front of the peak but increasing behind, relative to the FRDG scheme. This trend continues with the q = [0.05, 0.01] scheme, with the forward dispersion error continuing to decrease while the rearward dispersion error increases. These results are consistent with behaviour expected from previous observations of the dispersion and group velocity properties of all three spatial discretizations.

Viscous Burgers turbulence
Forced Burgers turbulence is a commonly used test case to study the behaviour of numerical methods for LES and ILES. For example, Adams et al. [27] used Burgers turbulence to investigate ILES using an adaptive local deconvolution method (ALDM), Yanan and Wang [28] used it to study the behaviour of LES and ILES using the FR scheme, and Moura et al. [29] used it to study under-resolved turbulence using DG and spectral h/p methods. The viscous Burgers equation is defined as where µ = 10 −4 is a constant viscosity term and s(x, t) is a time and space dependent forcing function. The common interface solution and gradient for the viscous terms are obtained using a local DG (LDG) approach for the FR scheme following Huynh [30]. The forcing function is specified following Adams et al. [27], where s(x, t) is defined in wavenumber space as where ψ is the wavenumber, A = 0.01, and −π ≤ φ ≤ π is chosen randomly for every wavenumber and at every time step. We consider all three k = 4 schemes with Ω ∈ [−π, π] using periodic boundaries and N e = 20 elements. The non-linear inviscid flux is anti-aliased at 2k + 1 Gauss points to eliminate aliasing errors at high wavenumbers.
Each simulation was initialized with u δ 0 (x, 0) = 0 and was advanced in time using the RK4 scheme. While the time step stability limit is well-defined for linear advection, for non-linear systems, such as the viscous burgers equation, it is less well-defined. Since we are interested in the behaviour of FR schemes near the stability limit, we define τ max = 1.85, 1.88, and 2.05 as the maximum time step size for which >50% of simulations remained stable for the FRDG, FRSD, and q = [0.05, 0.01] schemes, respectively. We then perform simulations at τ = 0.90, 0.92, 0.94, 0.96, 0.98, and 1.00τ max to investigate their behaviour with the RK4 scheme near their stability limit. Each simulation Fig. 9. Energy spectra for viscous Burgers turbulence using FRDG and RK4.  is advanced in time to t = 5000, allowing them to reach a fully developed state. Instantaneous energy spectrã are then time-averaged over a further 10 000 time units, whereũ(ψ) is the Fourier transform of the solution. Since the solution is dependent on a randomly parametrized source term, we also perform ensemble averaging over 100 independent simulations for each time step size. Results in the form of energy spectra are shown in Figs. 9-11 for all three schemes. A reference solution computed with a small time-step is shown demonstrating the expected behaviour based on a semi-discrete scheme. It is clear from these results, and has been shown previously [29,9], that the FR schemes can be used for ILES without the addition of an explicit subgrid scale model. Numerical dissipation, concentrated at high-wavenumbers acts as an implicit model for the un-resolved scales. However, we notice that when the time-step is increased near to τ max this dissipation mechanism can be lost and, rather than being effectively dissipated, there can be a build-up of energy at the highest resolved wavenumbers. This correlates with the dissipation characteristics shown in the von Neumann analysis studies with RK4, which demonstrated that numerical dissipation at the highest-resolved wavenumbers is diminished as the time-step approaches the theoretical stability limit τ CFL .

Conclusions
We have investigated the fully-discrete dispersion, dissipation, group velocity, and error behaviour of three different k = 4 FR schemes with RK3, RK4, SDIRK 3,3 , and SDIRK 5,4 time-stepping. In particular, we have considered the FRDG, FRSD, and q = [0.05, 0.01] schemes, which have been used previously for implicit large eddy simulation [9,22,19]. The behaviour of these fully-discrete schemes was found to be dependent on the time-stepping method being used, and the time-step size.
The fully-discrete explicit schemes with RK3 and RK4 exhibited a strong dependence on the time-step size, particularly for large wavenumbers κ 3π . The RK3 schemes showed increasing dissipation in the vicinity of κ = 4π with increasing time-step size for all spatial discretizations. However, numerical dissipation at the highest resolved wavenumbers decreased with increasing time-step for RK3, which was also observed in this region for RK4. In the context of ILES, this behaviour near the stability limit could be unfavourable, since numerical dissipation at the highest resolved wavenumbers is relied-upon to act as a simple subgrid scale dissipation mechanism. Von Neumann analysis also predicted that high-velocity waves can be generated using these ERK schemes, significantly exceeding the physical wavespeed. Lastly, numerical error at both high-and low-wavenumbers was found to be strongly dependent on the time-step size. For low-wavenumbers, the expected local semi-discrete super-convergence rates of B T = 2k + 2, 2k + 1, and 2(k − 1) + 1 for the FRDG, FRSD, and q = [0.05, 0.01] schemes, respectively, were lost and the local order of accuracy of the temporal schemes was observed instead. However, there was still a region near κ ≈ 10 in which super-convergence was maintained.
The fully-discrete implicit schemes with SDIRK 3,3 and SDIRK 5,4 also exhibited a strong dependence on the time-step size. However, since these schemes are L-Stable their behaviour for relatively large time-step sizes was investigated. With increasing time-step size SDIRK 3,3 tended to increase the amount of numerical dissipation at moderate-wavenumbers and decrease dissipation at high-wavenumbers. In addition, SDIRK 5,4 was found to primarily decrease the amount of dissipation at high-wavenumbers, leaving the dissipation at moderate-wavenumbers relatively unaffected. Both SDIRK 3,3 and SDIRK 5,4 were also found to reduce the group-velocities for moderate-and wellresolved waves as the time-step was increased. Finally, as with the explicit schemes, the numerical error at low-and high-wavenumbers was strongly dependent on the time-step size and, for very large time-step sizes, the error for all wavenumbers was found to be dominated by the temporal scheme and super-convergence properties were lost.
Numerical experiments confirmed the aforementioned results. Advection of a sine wave showed that the superconvergence properties of FR schemes, previously identified from semi-discrete von Neumann analysis [18,19], were indeed observed for moderately-resolved cases. However, this order of accuracy asymptotically approached the order of accuracy of the temporal scheme with mesh refinement at constant CFL number, corroborating the error properties observed from the current dispersion analysis and previous studies [20]. Advection of a Gaussian bump using the ERK schemes with large time-steps demonstrated the existence of high-speed nearly undamped waves, which significantly exceeded the physical wavespeed, and were also predicted by fully-discrete von Neumann analysis. Long-period simulations of the Gaussian bump test case also demonstrated the importance that choice of spatial and temporal schemes has on numerical dispersion and dissipation errors and damping of non-physical waves. Finally, ILES of Burgers turbulence demonstrated that linear von Neumann analysis of fully-discrete FR schemes can provide insight into their behaviour for more complicated non-linear simulations. Results from this study demonstrate the importance of considering the fully-discrete behaviour of FR schemes used for ILES. Future work will focus on investigating the relative behaviour of these different schemes for applied ILES of turbulent flows (see Figs. 12-19).