Computers and Fluids 266 (2023) 106060

Contents lists available at ScienceDirect

Computers and Fluids

journal homepage: www.elsevier.com/locate/compfluid

Fully-discrete spatial eigenanalysis of discontinuous spectral element

methods: Insights into well-resolved and under-resolved vortical flows
Niccolò Tonicello a ,∗, Rodrigo C. Moura b , Guido Lodato a , Gianmarco Mengaldo c
a
Normandie Université, INSA et Université de Rouen, St. Etienne du Rouvray (Rouen), France
b
Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil
c
National University of Singapore, Singapore

ARTICLE INFO ABSTRACT

Keywords: This study presents a comprehensive spatial eigenanalysis of fully-discrete discontinuous spectral element
High-order methods methods, now generalising previous spatial eigenanalysis that did not include time integration errors. The
Spectral element methods influence of discrete time integration is discussed in detail for different explicit Runge–Kutta (1st to 4th order
Flux reconstruction
accurate) schemes combined with either Discontinuous Galerkin (DG) or Spectral Difference (SD) methods,
Discontinuous galerkin
both here recovered from the Flux Reconstruction (FR) scheme. Selected numerical experiments using the
Spectral difference
Eigenanalysis
improved SD method by Liang et al. (2009) [53,54] and Jameson (2010) [55] are performed to quantify
the influence of time integration errors on actual simulations. These involve test cases of varied complexity,
from one-dimensional linear advection equation studies to well-resolved and under-resolved inviscid vortical
flows. When simulations are well-resolved, the overall order of accuracy of the (fully-discrete) method of
choice is limited to that of the time integration scheme. Moreover, it is shown that, while both well-resolved
and under-resolved simulations of linear problems correlate well with the eigenanalysis prediction of time
integration errors, the correlation can be much worse for under-resolved nonlinear problems as observed via
numerical experiments. In fact, in the numerical simulation of under-resolved vortical flows, the predominance
of spatial errors made it practically impossible for time integration errors to be distinctly identified. As a
result, the eigenanalysis predictions are expected to hold (even if partially) in direct numerical simulations
of turbulence. This highlights that the interaction between space and time discretisation errors is more
complex than otherwise anticipated, contributing to the current understanding about when eigenanalysis can
effectively predict the behaviour of numerical errors in practical under-resolved nonlinear problems, including
under-resolved turbulence computations.

1. Introduction to the rapidly developing field of high-order methods, especially in

fluid dynamics, the research community is constantly looking for more
Along with the increasing computational power experienced in the efficient and reliable tools to quantify the peculiar numerical properties
last few decades, the use of Computational Fluid Dynamics (CFD) of such schemes. Indeed, for under-resolved simulations (such as LES),
is becoming a more and more common tool in industrial aerody-
the role played by the numerical scheme can be of fundamental impor-
namic analysis and design [1]. The expanding capabilities of modern
tance. For instance, the same computation, performed with different
computers have fuelled the development of a large number of innova-
tive numerical high-order schemes for fluid dynamics, whose accuracy methods, can lead to dramatically different outcomes in realistic flow
clearly surpasses that of more popular low order methods. In particular, configurations [17,18,21].
much interest has been focused on spectral element methods due to Spectral analyses based on the numerical discretisation of the lin-
their geometric flexibility and promising computational efficiency [2]. ear advection equation are a common tool in the evaluation of the
These include the Discontinuous Galerkin (DG) [3–5], Spectral Differ- intrinsic numerical properties of high-order methods for advection
ence (SD) [6–8] and Flux Reconstruction (FR) [9,10] schemes. The dominated flows. Such techniques are particularly simple to imple-
use of high-order numerical schemes has shown promising results in ment and computationally inexpensive. However simplistic, they are
the simulation of turbulent flows, in both Direct Numerical Simula-
capable of providing useful information regarding the scheme, such as
tion (DNS) [11–14] and Large-Eddy Simulations (LES) [15–20]. Due

∗ Corresponding author.
E-mail address: ntonicel@sissa.it (N. Tonicello).

https://doi.org/10.1016/j.compfluid.2023.106060
Received 30 August 2022; Received in revised form 23 June 2023; Accepted 14 September 2023
Available online 23 September 2023
0045-7930/© 2023 Elsevier Ltd. All rights reserved.
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

numerical dispersion and diffusion. The most widely used spectral anal- complex inviscid flows based on the Euler equations. In line with
ysis technique is commonly referred to as temporal eigenanalysis [22– the theoretical discussion, numerical results will be arranged in well-
27], which quantifies how specific spatial oscillations will evolve over resolved and under-resolved flows. Finally, in Section 5, we outline the
time under periodic conditions in space. The complementary analysis key conclusions of this work.
technique, called spatial eigenanalysis [28–32], considers instead how
wave-like solutions of given frequency behave while propagating in 2. Eigenanalysis framework
space under inflow-outflow conditions. Spatial analysis can be used to
predict problems without spatial periodicity at the boundaries. On the In this section, we introduce the fully-discrete eigenanalysis frame-
other hand, periodicity in time at a fixed point in space is required. work, which serves both the temporal (Section 2.1) and the spatial
This notwithstanding, the spatial eigenanalysis can be considered as (Section 2.2) approaches. The classical fully-discrete temporal eigen-
a suitable strategy to tackle problems with more general boundary analysis (that was originally reported by Vermeire and Vincent [40]) is
conditions not enforcing periodicity in space. presented first, as a way of introducing the key building blocks for the
Even if similar and intimately connected from a theoretical per- fully-discrete spatial eigenanalysis.
spective, these two approaches can provide very different insights on
the general behaviour of the numerical scheme. The former is more 2.1. Classical fully-discrete temporal eigenanalysis
suitable for temporally evolving problems, whereas the latter is more
appropriate for spatially developing problems. All the relevant steps of temporal eigenanalysis will be herein
Although such techniques are able to reveal useful numerical prop- introduced using the FR scheme, which is well-known for its ability
erties, they have been often categorised as too simplistic. In fact, a to recover different methods with a specific choice of its correction
direct connection between the numerical discretisation of the linear functions [43–45]. The present study considers two such methods in
advection equation and more complex three-dimensional nonlinear particular, namely DG and SD. In fact, this strategy allows one to
simulations can be everything but trivial. The high-order schemes’ re- perform subsequent eigenanalysis studies for any other method that FR
search community has also produced substantial works on generalised
is capable of recovering.
forms of spectral analyses. This includes analyses of linear advection
In the standard temporal analysis, the linear advection equation
on multiple dimensions [22,33,34], of linear advection problems with
is discretised looking for wave-like solutions in order to study the
non-constant velocity [32,35,36] and even of alternative eigenanaly-
numerical description of specific wavenumbers and their evolution over
sis techniques, such as non-modal analysis [37] and combined-mode
time. Then, dispersion and dissipation properties of any scheme follow
analysis [38]. Nonetheless, all of these approaches commonly involve
directly from the corresponding eigensolutions. Let us now suppose the
a semi-discrete form of the linear advection equation, where only
case of a constant unitary advection velocity:
spatial numerical errors are considered. The main goal of the present
work is focused on a fully-discrete formulation of spatial eigenanalysis 𝜕𝑢 𝜕𝑢
+ = 0. (1)
where also temporal errors are taken into account. Despite the fact that 𝜕𝑡 𝜕𝑥
some recent work has been presented regarding fully-discrete temporal This equation admits plane wave solutions of the form
eigenanalysis [39–42], the full generalisation of spatial eigenanalysis
𝑢(𝑥, 𝑡) = 𝑒𝜄(𝜃𝑥−𝜔𝑡) , with 𝜄2 = −1, (2)
including a temporal discretisation remains, as of yet, unexplored.
Even if spatial errors play a critical role in the simulation of turbulent provided that the angular frequency 𝜔 = 𝜔(𝜃) is such that
flows, the interaction with the temporal counterparts is everything
but obvious. The temporal discretisation can significantly affect the Re(𝜔) = 𝜃 and Im(𝜔) = 0, (3)
spatial errors, leading to much different results in terms of accuracy where 𝜃 is a real-valued wavenumber chosen at the initial condition
and stability of the numerical simulation. (in the case of standard temporal eigenanalysis). Eq. (3) provides what
The ultimate goal of simplified spectral analyses is to predict, within are respectively known as dispersion and diffusion relations for the
a reasonable range, how a certain scheme, for a given resolution, would exact plane wave solution. A numerical discretisation of the linear
behave in practical flow configurations. For instance, it is possible to advection equation is now introduced on a grid of equally spaced
estimate which grid resolution allows one to reliably perform Implicit elements of fixed width ℎ𝑛 = ℎ = 1. Under such assumptions, a FR
LES (ILES), thereby reducing the computational burden otherwise in- spatial discretisation of Eq. (1) within the standard element 𝛺𝑛 can be
curred. In the same way, temporal discretisations would benefit from written as
large time steps, while introducing negligible numerical errors. In other
d̂
𝑢𝑛 ∑𝑁
d𝑙𝑗 d𝑔 d𝑔
words, a trade-off between accuracy and computational cost needs to = −2 ̂
𝑢𝑛 (̂ 𝑢𝐿 ) 𝐿 (̂
𝑥𝑖 ) − (2𝑓̂𝐿𝐼 − 2̂ 𝑢𝑅 ) 𝑅 (̂
𝑥𝑖 ) − (2𝑓̂𝑅𝐼 − 2̂ 𝑥𝑖 ), (4)
be found both spatially and temporally. Finding such delicate balance dt 𝑗=0
d̂
𝑥 d̂
𝑥 d̂
𝑥
using fully three-dimensional simulations can be extremely expensive
where 𝑁 is the order of the solution polynomial, ̂ 𝑢𝑛 (−1, 𝑡) and
𝑢𝐿 = ̂
and even pointless. Ideally, efficient and reliable simplified spectral
𝑢𝑅 = ̂
̂ 𝑢𝑛 (+1, 𝑡), 𝑙𝑗 is the Lagrange polynomial defined on the point 𝑥𝑗 ,
analysis techniques aim at providing precisely this kind of information
and 𝑔𝐿∕𝑅 are called correction functions (see [9,10,32] for additional
at very low computational cost.
details on the Flux Reconstruction scheme).
The paper is organised as follows. In Section 2, we introduce the
To examine a sufficiently general framework, the classical set of
fully-discrete eigenanalysis framework which serves both temporal and
numerical fluxes for the linear advection equation is considered:
spatial approaches. This framework is applied to the Flux Reconstruc-
tion method, which can recover both the discontinuous Galerkin and 𝑓̂𝐿𝐼 = (1 − 𝛼)̂
𝑢𝑛−1 (+1, 𝑡) + 𝛼̂
𝑢𝑛 (−1, 𝑡) (5)
Spectral difference schemes. In Section 3, we present the eigenanalysis
results for the fully-discrete spatial approach, where we focus on the and
interpretation of the diffusion curves and the contributions of the differ- 𝑓̂𝑅𝐼 = (1 − 𝛼)̂
𝑢𝑛 (+1, 𝑡) + 𝛼̂
𝑢𝑛+1 (−1, 𝑡). (6)
ent (transmitted/reflected) modes. The observations will be specifically
classified into well-resolved and under-resolved frequency ranges, pro- Note that 𝛼 = 0 leads to a standard upwind scheme, whereas 𝛼 = 0.5
viding a direct link with DNS and LES approaches, respectively. In leads to a centred scheme. In the following sections, unless stated
Section 4, we present the numerical experiments conducted to assess otherwise, standard upwind fluxes will be considered. The 𝛼 parameter
the validity of the eigenanalysis. In particular, first a discretisation of can be used to mimic different levels of interface upwinding, which can
the linear advection equation will be considered, followed by more be directly linked to more complex numerical fluxes for the Euler or

2
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Navier–Stokes equations (via Riemann solvers). Eq. (4) can be written radius of the iteration matrix 𝐓 is strictly smaller than unity. Namely,
in matrix form as: for the specific case of Von Neumann analysis, the full discretisation is
𝐮𝑛
d̂ considered linearly stable if, for a given time step 𝛥𝑡,
𝐮𝑛 − (2𝑓̂𝐿𝐼 − 2𝐥𝑇 ̂
= −2𝐃̂ 𝐮𝑛 )𝐠𝐿 − (2𝑓̂𝑅𝐼 − 2𝐫 𝑇 ̂
𝐮𝑛 )𝐠𝑅 , (7) [ ]
dt
∀ 𝜃̃ ∈ −𝜋(𝑁 + 1)∕ℎ; 𝜋(𝑁 + 1)∕ℎ ∶ 𝜆(𝐓(𝜃, ̃ 𝛥𝑡)) < 1, (14)
d𝑙 𝐿∕𝑅 𝐿∕𝑅 d𝑔
where ̂ 𝐮𝑛𝑖 = ̂ 𝑥𝑖 ), 𝐃𝑖𝑗 = d̂𝑥𝑗 (̂
𝑢(̂ 𝑥𝑖 ), 𝐠𝑖 = d̂ 𝑥𝑖 ) and 𝐫𝑖 = 𝑙𝑖 (1),
(̂
𝑥 where 𝜆(𝐓) denotes the spectral radius of the matrix 𝐓.
𝐥𝑖 = 𝑙𝑖 (−1).
A complementary analysis, instead, takes into account the Bloch
In order to reduce the infinite dimensional problem to an element-
wave-like form of the numerical solution and generalises the concept
wise framework, Bloch wave-like type of solutions are sought, namely:
of dispersion/diffusion curves in the fully-discrete analyses. In fact,
exploiting the Bloch-like form of the approximate solution, it can be
̃
𝐮𝑛 = 𝑒𝜄(𝜃𝑥𝑛 ∕ℎ−𝜔𝑡)
̂ ̃ ̂
𝐯, (8) shown (see [40,46]) that

where tilde accent denotes the numerical counterparts of continuum 𝐮𝑡+𝛥𝑡 = 𝑒−𝜄𝜔𝛥𝑡
̂ ̃ ̂𝑡
𝐮. (15)
variables and 𝑥𝑛 indicates the left boundary of the element 𝑛 in the
Consequently, the fully-discrete version of the dispersion/diffusion
global reference frame. Once ̂ 𝐮𝑛 is properly defined, due to periodicity,
eigenanalysis can be written as:
it is straightforward to obtain a closed form for all the fluxes expressed
in Eq. (4): 𝑒−𝜄𝜔𝛥𝑡
̃ ̂𝑡
𝐮𝑡 .
̃ 𝛥𝑡)̂
𝐮 = 𝐓(𝜃, (16)
̃
𝑓̂𝐿𝐼 = (1 − 𝛼)𝑒−𝜄𝜃 𝑇 ̂𝑛
𝐫 𝐮 + 𝛼𝐥 𝑇 ̂𝑛
𝐮 , Notice, in fact, that for 𝛥𝑡 → 0 we have
̃
(9)
𝑓̂𝑅𝐼 = (1 − 𝛼) 𝐫𝑇 ̂
𝐮𝑛 + 𝛼𝐥𝑇 𝑒𝜄𝜃 𝐮𝑛 .
̂ 𝐓→ 𝐈 + 𝐇𝛥𝑡 + 𝑂(𝛥𝑡2 ),
Finally, it can be seen that, due to solution periodicity, the problem has 𝑒−𝜄𝜔𝛥𝑡
̃
→ 1 − 𝜄𝜔𝛥𝑡 + 𝑂(𝛥𝑡2 )
now been reduced to an element-wise formulation, i.e., the numerical
and the problem defined by Eq. (16) reduces to
fluxes now depend only on the local solution ̂ 𝐮𝑛 . Hence, index 𝑛 can be
dropped and the final discretised equation can be written as ̃ + 𝑂(𝛥𝑡2 ))̂
(1 − 𝜄𝜔𝛥𝑡 𝐮𝑡 = (𝐈 + 𝐇𝛥𝑡 + 𝑂(𝛥𝑡2 ))̂
𝐮𝑡
(
d̂𝐮 ̃ 𝐮𝑡 + 𝑂(𝛥𝑡2 )
̃ 𝐮𝑡 = 𝐇𝛥𝑡̂
−𝜄𝜔𝛥𝑡̂
= −2 𝐃̂ 𝐮 + ((1 − 𝛼)𝑒−𝜄𝜃 𝐫 𝑇 ̂
𝐮 + 𝛼𝐥𝑇 ̂
𝐮 − 𝐥𝑇 ̂
𝐮)𝐠𝐿
dt ̃ 𝐮𝑡 = 𝐇̂
−𝜄𝜔̂ 𝐮𝑡 , as 𝛥𝑡 → 0,
)
̃
+ ((1 − 𝛼)𝐫 𝑇 ̂
𝐮 + 𝛼𝐥𝑇 𝑒𝜄𝜃 ̂
𝐮 − 𝐫𝑇 ̂
𝐮)𝐠𝑅 , which coincides with classical semi-discrete eigenvalue problem.
In a similar way with respect to the classical analysis, considering
or, in a more compact way: (16), it is possible to compute the frequencies 𝜔̃ for any given value
d̂
𝐮 ̃ ̃ of 𝜃̃ and 𝛥𝑡. Therefore, the choice of the time integration scheme and,
= −2(𝐂− 𝑒−𝜄𝜃 + 𝐂0 + 𝐂+ 𝑒𝜄𝜃 )̂
𝐮, (10) in particular, of a stable time step 𝛥𝑡, will affect the dispersion and
dt
diffusion curves. Such curves can be interpreted in the same manner as
where 𝐂0 = 𝐃 − (𝛼 − 1)𝐠𝐿 𝐥𝑇 − 𝛼𝐠𝑅 𝐫 𝑇 , 𝐂− = (1 − 𝛼)𝐠𝐿 𝐫 𝑇 and 𝐂+ = 𝛼𝐠𝑅 𝐥𝑇 .
in the semi-discrete analysis and significant deviations from the semi-
When FR is made to recover specific methods (such as DG or SD),
discrete dispersion/diffusion curves can crucially affect the numerical
matrices 𝐂 above are what will differentiate them, but the structure
properties of the spatial discretisation.
of the problem does not change. In either case, the discretised system
We finally note that, had we chosen to work in the equations
is now reduced to the following linear dynamical system:
above with unspecified values of the advection velocity 𝑎 and the mesh
d̂𝐮 ̃ 𝐮 with 𝐇 = −2(𝐂− 𝑒−𝜄𝜃̃ + 𝐂0 + 𝐂+ 𝑒𝜄𝜃̃ ). spacing ℎ, the relevant time-step factor appearing in Eq. (13) would
= 𝐇(𝜃)̂ (11)
dt have been 𝛥𝑡 = 𝑎𝛥𝑡∕ℎ—namely, the Courant–Friedrichs–Lewy (CFL)
The approach herein outlined represents a classical tool in numerical number—instead of simply 𝛥𝑡. Hence, we remark that, although we
analysis, particularly used for the theoretical study of novel high-order chose to work with 𝛥𝑡 for the sake of simplicity, one can actually have
schemes. Notice that the present analysis is often called semi-discrete in mind the CFL number whenever 𝛥𝑡 values are mentioned in the
temporal eigenanalysis due to the absence of time discretisation. In- remainder of the paper.
deed, taking advantage of the analytical form of wave-like solutions,
the dynamical system defined by Eq. (11) can be reduced to the 2.2. Fully-discrete spatial eigenanalysis
following eigenvalue problem:
̃ 𝐮. In the classical temporal eigenanalysis, for any real-valued wave-
−𝜄𝜔̂
̃ 𝐮 = 𝐇(𝜃)̂ (12) ̃ the numerical angular frequency 𝜔̃ can be obtained through
number 𝜃,
Alternatively, instead of considering the analytical form of the wave- the eigenvalue problem in Eq. (16). The frequency 𝜔̃ will describe the
like solution, a numerical time integration scheme can be introduced numerical evolution of the specified wavenumber 𝜃̃ over time. It is
to advance in time. The numerical solution at the time step 𝑡 + 1 can be however possible to invert the problem by imposing a time frequency
generally approximated by the M-stage explicit Runge–Kutta schemes 𝜔̃ and subsequently seek for wavenumbers 𝜃̃ that solve the inverted
as problem. Such approach is commonly referred to as spatial eigenanalysis
and it represents a very useful tool for non-periodic problems since
∑
𝑀
̃
(𝐇(𝜃)𝛥𝑡)𝑚
𝐮𝑡+𝛥𝑡 =
̂ 𝛽𝑚 𝐮𝑡 = 𝐓(𝜃,
̂ 𝐮𝑡 ,
̃ 𝛥𝑡)̂ (13) it aims at studying the dynamics of prescribed time frequencies over
𝑚=0
𝑚! space.
where 𝛥𝑡 is the time step and 𝜷 contains the coefficients of the specific Within the framework of semi-discrete eigenanalysis, the equation
Runge–Kutta scheme. defining the relationship between 𝜃̃ and 𝜔̃ is the same, namely Eq. (12),
This formulation is denoted as fully-discrete since both spatial and and can be written in an alternative form as
temporal discretisations are simultaneously involved. Considering the ̃ + 𝜄𝜔𝐈)
det(𝐇(𝜃) ̃ = 0. (17)
present generalised formulation, two main closely linked analyses can
be performed. The first is relatively popular and originates from the In this case, 𝜔̃ is a real-valued prescribed scalar and the nonlinear
classical stability analysis of numerical discretisation of partial differ- equation (17) is solved in the unknown 𝜃.̃ This nonlinear equation can
ential equations. In particular, it is well-known that in Eq. (13), the be particularly complex to solve and, except for some very particular
combined space–time discretisation is denoted as stable if the spectral cases, cannot be analytically inverted. Considering the fully-discrete

3
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 1. Dissipation curve for the 5th-order FR-SD standard upwind scheme coupled with RK11 (explicit Euler) denoting a single transmitted physical mode, shown for
𝛥𝑡 = 0.01, 0.02, … , 0.1 (with dissipation values becoming more negative as 𝛥𝑡 increases). The plot on the right shows a closer look near the origin. These re-affirm that RK11
is convectively unstable for pure advection.

space–time discretisation, an equivalent equation which defines the 3. Theoretical results

relation between 𝜔̃ and 𝜃̃ can be derived from Eq. (16), namely
In this section, we show the results associated with the fully-discrete
̃ 𝛥𝑡) − 𝑒−𝜄𝜔𝛥𝑡
det(𝐓(𝜃, ̃
𝐈) = 0. (18) spatial eigenanalysis just introduced. The theoretical findings of the
fully-discrete temporal eigenanalysis previously reported by Vermeire
Note of course that, in the temporal eigenanalysis, the addition of a
and Vincent [40] were correctly reproduced within the present imple-
temporal discretisation does not significantly complicate the overall
mentation too. These have been used for validation purposes only and,
problem. In fact, the matrix 𝐓 can be easily computed once 𝜃̃ and 𝛥𝑡
for the sake of brevity, are not shown here.
are prescribed. Subsequently, the eigenvalue problem can be solved in It is important to highlight that, from Eq. (2), the numerical wave-
𝑒−𝜄𝜔𝛥𝑡
̃ , which can be finally inverted to obtain the exact value of 𝜔. ̃ like solutions depend on the numerical equivalent of 𝜃, i.e. 𝜃̃ (whose
However, in the spatial eigenanalysis, the specific matrix 𝐓 influences imaginary part makes the diffusion curve plots below) according to
the general structure of the nonlinear equation to be solved in 𝜃. ̃ In
̃ exp[𝜄Re(𝜃)𝑥
𝑢 ∝ exp[−Im(𝜃)𝑥] ̃ − 𝜔𝑡] , (23)
particular, if the simple explicit Euler integration scheme is considered,
then Eq. (18) becomes: where reflected waves should be considered with 𝑥 < 0, thus being
stable for Im(𝜃) ̃ < 0.
̃ 𝛥𝑡) − 𝑒−𝜄𝜔𝛥𝑡
det(𝐈 + 𝛥𝑡𝐇(𝜃, ̃
𝐈) = 0. (19)
In the following, we consider high-order FR-based DG/SD methods
̃ in space, coupled with explicit Runge–Kutta schemes of increasing
If 𝑧 is defined as 𝑧 = 𝑒𝜄𝜃 , the final nonlinear equation in 𝑧 can be written
as order of accuracy (and increasing number of stages): RK11, RK22,
RK33 [47–49] and RK45 [50]. Hence, we start with the 1st-order
det((1 − 𝑒−𝜄𝜔𝛥𝑡
̃
)𝐈 − 2𝛥𝑡(𝐂− 𝑧−1 + 𝐂0 + 𝐂+ 𝑧)) = 0. (20) accurate, one-stage Runge–Kutta scheme (more commonly known as
the Explicit Euler method) until finally addressing the well-known
Again, notice that for 𝛥𝑡 → 0, the diagonal term 1 − 𝑒−𝜄𝜔𝛥𝑡̃ tends 4th-order accurate, five-stages Runge–Kutta scheme by Carpenter and
to 𝜄𝜔𝛥𝑡
̃ and the fully-discrete formulation goes back to the classical Kennedy [50]. This allows for a detailed exposition of the fully-discrete
semi-discrete form: eigenanalysis results, whose complexity grows with the number of
stages of the temporal scheme, as shown below.
̃ − 2(𝐂− 𝑧−1 + 𝐂0 + 𝐂+ 𝑧)) = 0,
det(𝜄𝜔𝐈 (21)
In Fig. 1, we show the dissipation curve for the 5th-order FR-
which is exactly Eq. (17). On the other hand, for a RK22 scheme SD standard upwind method when coupled with RK11. This curve
is shown for increasing values of the time step 𝛥𝑡, namely 𝛥𝑡 =
[47], the formulation becomes much more complex as higher powers
0.01, 0.02, … , 0.1, with dissipation values becoming more negative as
of 𝑧 appear in the final equation:
( 𝛥𝑡 increases. Note that, given a 𝛥𝑡 value, only one eigenmode/curve
exists in this case, corresponding to a single transmitted physical wave.
det (1 − 𝑒−𝜄𝜔𝛥𝑡
̃
)𝐈
This is because, from Eq. (20), the characteristic polynomial in the
2
( )
− − −2
+ 𝛥𝑡 𝐂 𝐂 𝑧 unknown 𝑧 has the same order of the semi-discrete analysis (compare
( ) with Eq. (21)). We recall that, in the semi-discrete analysis [29], only
+ 𝛥𝑡2 (𝐂− 𝐂0 + 𝐂0 𝐂− ) − 2𝛥𝑡𝐂− 𝑧−1
( 2 − 0 ) (22) the physical transmitted mode exists in case of standard upwinding
+ 𝛥𝑡 (𝐂 𝐂 + 𝐂0 𝐂0 + 𝐂+ 𝐂− ) − 2𝛥𝑡𝐂0 (𝛼 = 0 in Eq. (10)), whereas one reflected spurious wave exists in
( )
+ 𝛥𝑡2 (𝐂+ 𝐂0 + 𝐂0 𝐂+ ) − 2𝛥𝑡𝐂+ 𝑧 addition if 𝛼 ≠ 0. As for the curves in Fig. 1, negative dissipation values
)
( ) indicate in fact anti-dissipative behaviour (convective instability), since
+ 𝛥𝑡2 𝐂+ 𝐂+ 𝑧2 = 0. we are dealing with a physical transmitted wave, cf. Eq. (23). Hence,
Fig. 1 essentially re-affirms that RK11 (explicit Euler) is only marginally
This way, the order of the nonlinear equation increases significantly stable for pure advection [51]. Interestingly, for moderately large 𝛥𝑡,
and additional roots for 𝑧 are expected to be found (especially for waves regain stability beyond a certain frequency 𝜔, probably due to
time discretisation schemes of higher order). The role of multiple increased dissipation from the spatial discretisation at large frequen-
wave modes in spatial eigenanalysis has been discussed in more detail cies. These features will be tested in the next section, as they highlight
in [31]. As properly explained in the next section, only one solution how the fully-discrete eigenanalysis is useful in revealing non-intuitive
will be denoted as the physical one, whereas the others are considered characteristics stemming from the coupling between spatial and tem-
spurious. poral discretisations. Now we consider the classical 2nd-order explicit

4
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 2. Dissipation curves for the 2nd-order FR-SD standard upwind scheme coupled with RK22 for increasing values of 𝜏 = 𝛥𝑡∕𝛥𝑡max . The plots on the right show the same left
curves, but for a larger frequency range. The semi-discrete (𝜏 = 0) result also appears when 𝜏 > 0 for reference as a thin continuous curve. The physical mode is shown in red,
whereas the spurious transmitted mode is shown as a thick black curve. A vertical dashed line indicates the cut-off frequency 𝜔𝑐 = 𝜋∕𝛥𝑡.

Runge–Kutta (RK22) scheme, this time coupled for simplicity with the modes is guaranteed by the characteristic polynomial in Eq. (22) being
2nd-order FR-SD standard upwind method. The resulting discretisation only quadratic in 𝑧, since matrix 𝐂+ becomes zero in case of standard
is stable for 𝛥𝑡 < 𝛥𝑡max , where 𝛥𝑡max is the critical time-step size at upwinding (𝛼 = 0 in Eq. (10)). This feature can be traced algebraically
which instabilities first appear. In Fig. 2, we show the corresponding to the number of stages of the time discretisation adopted. In general,
dissipation curves for increasing values of 𝜏 = 𝛥𝑡∕𝛥𝑡max . When 𝜏 = 0, the number of modes (in case of standard upwinding) matches the
the results match those of the semi-discrete analysis, where only the number of stages of the chosen explicit Runge–Kutta scheme. The
transmitted physical mode exists in case of standard upwinding. For 𝜏 > physical mode can be identified as the one anchored to the origin of
0, a spurious transmitted mode appears in addition to the transmitted the dissipation plots, thus featuring vanishing dissipation levels for
physical one, not to be confused with the thin continuous curve which well-resolved waves (as expected in simulations of pure advection).
simply shows the semi-discrete (𝜏 = 0) physical curve as reference Returning to Fig. 2, we note that, as 𝜏 increases, the physical mode
(note that the physical curve is shown in red, whereas the spurious becomes less dissipative overall (at least for 𝜔ℎ(𝑁 +1)−1 > 0.5, say) until
one is shown as a thick black curve). The existence of exactly two finally touching the horizontal axis, whose crossing (for 𝜏 > 1) means

5
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 3. Dissipation curves of the physical transmitted mode at increasing values of 𝜏 (from 𝜏 = 0 in black until 𝜏 = 1 in yellow). Results obtained for the 2nd-order FR-SD standard
upwind discretisation coupled with RK22.

anti-dissipation (convective instability). The spurious mode, despite to how multiple element-wise polynomial modes increase the spatial
being infinitely dissipated for 𝜏 → 0, also has its dissipation reduced resolution of a single element, as if it had sub-elements.
as 𝜏 increases; in fact, for moderately large 𝜏, its dissipation levels A first obvious consequence of adopting 𝑅 𝜔𝑐 ∕2 as the right-most
are comparable to those of the physical mode. Small dissipation is value of 𝜔 when plotting eigencurves is that this bound tends to infinity
typically undesirable for spurious modes, since it allows them to survive for 𝜏 → 0. Indeed, in the semi-discrete spatial analysis, every value of
longer as artificial oscillations in a simulation, potentially affecting 𝜔 is allowed, in principle. Secondly, such bound gets closer and closer
not only solution quality, but also numerical stability [30]. Curiously, to zero for larger time steps, as only lower and lower frequencies can
as 𝜏 → 1, the spurious mode’s dissipation gradually approaches that be solved by the temporal discretisation when the time step increases.
of the physical mode, almost matching it at 𝜏 = 1. This behaviour These effects are highlighted in Fig. 3, where only the physical mode
closely resembles that of reflected spurious modes which exist for either is represented for different values of 𝜏.
very low or very high upwinding levels (i.e. either centred fluxes or We now consider the 5th-order FR-SD standard upwind scheme
hyper-upwinding ), cf. [29], except that reflected modes have dissipation coupled with either RK33 (Fig. 4) or RK45 (Fig. 5). Dissipation curves
values of opposite sign when compared to those of their corresponding for 𝜏 = 0.5 and 𝜏 = 1 are respectively shown on the left and right
physical modes. In any case, for the present scheme, both physical plots within these figures. As previously mentioned, the number of
and spurious modes touch the horizontal axis of zero dissipation when roots in the respective characteristic polynomials matches the number
𝜏 = 1, becoming therefore unstable for 𝛥𝑡 > 𝛥𝑡max . of stages in the adopted time scheme. This is why RK33 (Fig. 4)
Still in Fig. 2, note that the spurious mode’s dissipation curve (e.g. yields three transmitted modes, whereas RK45 (Fig. 5) yields five.
when considered until 𝜔𝑐 = 𝜋∕𝛥𝑡, such that 𝜔𝑐 ℎ(𝑁 + 1)−1 is marked The correspondence between the number of eigenmodes and stages
as a vertical dashed line) resembles a mirror image of the physical in spatial analysis resembles that between eigenmodes and element-
mode’s curve. What is more generally true, though, is that the spurious wise polynomial modes in temporal analysis (in line with the previous
curve is simply a shifted replication of the physical curve, while both comment regarding the interpretation of stages as sub-divisions of the
are periodic functions of period 2𝜔𝑐 . Interestingly, the periodicity of actual time step 𝛥𝑡). In Figs. 4 and 5, the black and red vertical dashed
eigencurves in spatial analysis has not been noted in the literature lines indicate, respectively, 𝜔𝑐 and 𝑅 𝜔𝑐 ∕2, the latter matching half
so far. This is because 𝜔𝑐 → ∞ when 𝜏 → 0, cf. Fig. 2, consistent the period of the physical mode, as expected. Note that the physical
with the fact that the spurious mode becomes infinitely dissipated for mode is found to be symmetric with respect to the red dashed line
𝜏 → 0, as already mentioned. Therefore, periodicity is only visible in (actual Nyquist frequency). Moreover, all eigenmodes are observed
fully-discrete spatial analysis. This is in contrast with temporal analysis, to be shifted replications of the physical mode, similar to temporal
where one sees periodicity already in the semi-discrete case, as well as analysis. This is sometimes useful for the identification of the physical
the shifted replication property amongst multiple eigencurves [26,52]. mode.
Although the cut-off value 𝜔𝑐 = 𝜋∕𝛥𝑡 can be rightfully recognised The physical curves of Figs. 4 and 5 have been isolated and are
as the Nyquist frequency imposed by the discrete time-stepping even shown in Fig. 6 for different values of 𝜏. For comparison, the equivalent
outside of the fully-discrete spatial framework [40], it is only in the plots for FR-DG are given in Fig. 7. As one can see, the general trends
fully-discrete spatial analysis that 𝜔𝑐 exhibits visual correlation with the observed for FR-SD are also seen for FR-DG. Some differences, however,
periodicity of the eigencurves, as already noticed in Fig. 2. In fact, the can be discerned. For example, since the CFL limit is more restrictive for
normalised version of 𝜔𝑐 , namely 𝜔𝑐 ℎ(𝑁 +1)−1 , is the counterpart of the FR-DG than for FR-SD (both under standard upwinding), the dissipative
wavenumber limit 𝜃ℎ(𝑁 + 1)−1 = 𝜋 typically used in temporal analysis. curves of FR-DG extend to higher frequencies (recall that 𝜔𝑐 = 𝜋∕𝛥𝑡).
Hence, it seems appropriate to plot spatial analysis curves only until 𝜔𝑐 . Moreover, as already pointed out in previous works [28–30,32], FR-
Nevertheless, as it will become clear in subsequent plots, it so happens DG is spatially less dissipative than FR-SD, overall. A more complete
that the actual limit observed for explicit Runge–Kutta schemes of 𝑅 comparison between FR-SD and FR-DG is shown in the Appendix (cf.
stages is 𝑅 𝜔𝑐 ∕2. Interestingly, the conclusion is that multiple stages Figs. 29 and 30).
seem to affect the temporal resolution as if multiple single-stage time To quantify the mutual influence between time and spatial discreti-
steps of size 𝛿𝑡 = 2𝛥𝑡∕𝑅 had been employed instead. This seems akin sations, the approximation error can be evaluated separately for both

6
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 4. Dissipation curves for the 5th-order standard upwind FR-SD coupled with RK33 at 𝜏 = 0.5 (left) and 𝜏 = 1 (right). The physical mode is shown as a continuous red curve.
The vertical black dashed line marks 𝜔𝑐 , whereas the vertical red dashed line marks 𝑅 𝜔𝑐 ∕2 (indicating half a period).

Fig. 5. Dissipation curves for the 5th-order standard upwind FR-SD coupled with RK45 at 𝜏 = 0.5 (left) and 𝜏 = 1 (right). The physical mode is shown as a continuous red curve.
The vertical black dashed line marks 𝜔𝑐 , whereas the vertical red dashed line marks 𝑅 𝜔𝑐 ∕2 (indicating half a period).

Fig. 6. Dissipation curves of the physical transmitted mode at increasing values of 𝜏 (from 𝜏 = 0 in black until 𝜏 = 1 in yellow). Results obtained for the 5th-order FR-SD standard
upwind discretisation coupled with either RK33 (left) or RK45 (right).

temporal and spatial approaches, since the exact solution is known. In 4th-order convergence. For a certain intermediate range of wavenum-
particular, for the temporal approach, the error 𝐸𝑇 will read bers, instead, the theoretical accuracy of the scheme is preserved. The
smaller the prescribed time step, the longer such range will penetrate
𝐸𝑇 = |𝜔(𝜃)
̃ − 𝜔(𝜃)|, (24)
into the small wavenumbers region. The same applies in the spatial
whereas, in the spatial analysis the error 𝐸𝑆 can be evaluated as analysis: for small frequencies, the scheme converges with order 4 and
for intermediate frequencies with order 9. The differences between the
𝐸𝑆 = |𝜃(𝜔)
̃ − 𝜃(𝜔)|. (25) two plots are restricted to the high frequencies/wavenumbers region.
Notice, of course, that the estimates above take into account both In fact, since the imaginary parts of both 𝜔̃ and 𝜃̃ are negligible for
dispersive and diffusive errors. Both errors have been evaluated for small frequencies/wavenumbers, the temporal and spatial approaches
the 5th-order standard upwind FR-SD method coupled with the RK33 are equivalent, whereby 𝜔̃ and 𝜃̃ simply switch axes. As already dis-
scheme. Results are shown in Fig. 8. Similar to what has been found cussed in previous works on spatial eigenanalysis [28,29], whenever
in [40], for small values of 𝜃, the accuracy of the scheme is reduced the numerical interface flux departs from the standard upwind con-
by the time integration method. In particular, the expected super- dition, twice as many roots will exist in the relevant characteristic
convergence of the FR-SD scheme (of order 2𝑁 + 1 = 9) reduces to a polynomial. The new roots are typically found to represent wave-like

7
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 7. Dissipation curves of the physical transmitted mode at increasing values of 𝜏 (from 𝜏 = 0 in black until 𝜏 = 1 in yellow). Results obtained for the 5th-order FR-DG standard
upwind discretisation coupled with either RK33 (left) or RK45 (right).

Fig. 8. Approximation errors 𝐸𝑇 and 𝐸𝑆 obtained, respectively, from the temporal (left) and spatial (right) frameworks for a 5th-order FR-SD method coupled with the RK33
scheme at increasing values of 𝜏 (from 𝜏 = 0 in black until 𝜏 = 1 in yellow).

Fig. 9. Dissipation curves of the physical transmitted mode at increasing values of 𝜏 (from 𝜏 = 0 in black until 𝜏 = 1 in yellow). Results obtained for the 5th-order FR-SD with
nearly centered flux (𝛼 = 0.49) coupled with either RK33 (left) or RK45 (right).

modes that propagate in the opposite direction (reflected modes) and touch the horizontal axis with the potential consequent onset of anti-
are consequently considered spurious. As already shown for upwind dissipative behaviours, leading to convective instability. This is clearly
fluxes, for any given value of 𝜔, the number of roots depends on the visible in Fig. 9 for the RK33 scheme (left). Less so for the RK45 scheme,
number of stages of the time integration scheme. In the same way, for which a similar behaviour is expected for high frequencies which
for non-upwind schemes, the total number of roots is expected to be are well beyond the plotted range. For almost centered fluxes, in the
twice the number of stages. In the Appendix (cf. Figs. 31 and 32), plots spatial analysis framework, the presence of dissipative bubbles has been
showing all such eigencurves are given. Here, to avoid unnecessary already observed in previous works [29,32]. In other words, significant
confusion, only the physical modes are shown. The case considered is dissipation appears concentrated in certain frequency ranges, whereas
that of nearly centered flux (𝛼 = 0.49), whose dissipation curves are other ranges can be almost completely unaffected. In Fig. 9, it can
shown in Fig. 9, again for the 5th-order FR-SD method coupled with be noticed that, for increasing values of 𝜏, the separation of scales
either RK33 or RK45, as previously shown for the standard upwind case is less strong and the dissipative curves are characterised by a more
in Fig. 6. Similarly to the standard upwind case, for increasing values monotonic behaviour. Non-monotonic dissipation profiles in frequency
of 𝜏, numerical dissipation decreases in the high frequency region. space might potentially lead to undesirable non-smooth features in
For sufficiently high values of 𝜏, the curve for Im(𝜃) ̃ is expected to simulations with broadband energy content. It is expected that having

8
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

neighbouring frequencies so differently affected can deteriorate the dis- comparison has been performed for the frequency 𝜔2 , which is located
tribution of dissipation among scales and impact accuracy and stability in between dissipative bubbles. More or less the same behaviour is ex-
of the numerical simulation. pected, where larger values of 𝜏 lead to stronger numerical dissipation.
Consequently, only well-resolved wavenumbers/frequencies behave The solution 𝑢(𝑥) is shown for the RK33 scheme in Fig. 13 and for
quite similarly in temporal/spatial approaches, whereas large differ- the RK45 scheme in Fig. 14. The only interesting difference is that, in
ences are observed in the high wavenumbers/frequencies region. For agreement with the theoretical analysis, for 𝜏 = 0.9 the RK33 scheme
practical turbulent simulations using high-order discontinuous spectral is less dissipative than the RK45. For 𝜏 ≈ 0 and 𝜏 = 0.5, instead, the
element methods, depending on the local level of resolution, differ- levels of numerical dissipation are quite similar. In the RK33 case, the
ent ranges of scales and frequencies are of particular interest. For additional frequency 𝜔3 has been considered. At this frequency, in fact,
highly resolved flows, such as in DNS, the correlation with spectral a different behaviour with respect to the low frequencies is expected:
analysis will be mainly restricted to the well-resolved regions (i.e. for increasing values of 𝜏, numerical dissipation is expected to decrease
small wavenumbers/frequencies). Instead, reducing the resolution of instead of increase like for all the other frequencies. Such behaviour
the scheme and passing from DNS to LES of turbulent flows, higher is confirmed in Fig. 15 where for 𝜏 very close to the CFL limit the
wavenumber and frequency regions in the plots will provide more numerical dissipation is significantly lower. In similarity with fully-
useful information. discrete temporal eigenanalysis, the highest frequencies are the first to
experience instability. Then, in a general sense, the present test case
4. Numerical experiments suggests that for a reasonably large, low frequency region, the influence
of the time integration scheme is essentially dissipative and larger time
In order to validate the theoretical framework introduced in previ- steps cause stronger dissipation. In the high-frequency region, instead,
ous sections, a series of numerical experiments is presented to quantify the behaviour is the opposite: in the semi-discrete analysis such fre-
the interplay between temporal and spatial errors for practical flows. quencies are usually dumped very fast, whereas, approaching the CFL
First, we consider the discretisation of the one-dimensional linear ad- limit, numerical dissipation decrease significantly, until, eventually, it
vection using the SD scheme. We then take into account more complex changes sign leading the simulation to be unstable. Consequently, in
numerical simulations of the Euler equations in order to evaluate the similarity with respect the fully-discrete temporal eigenanalysis, the
influence of the time step on both well-resolved and under-resolved largest frequencies are responsible for the CFL instability. Finally it
flows. In all the Euler computations presented herein, if not stated
is interesting to study the behaviour of a scheme that, according to
differently, the original Roe flux has been employed as a representative
the fully-discrete spatial eigenanalysis herein presented, is expected
of the standard upwind flux. All the simulations presented in the
to be unstable for any choice of time-step. The same computational
following subsections have been performed using the improved SD
grid and order of approximation have been chosen. In similarity with
method by Liang and Jameson [53–55].
previous tests, a series of frequencies have been simulated. In Fig. 16
the dissipative curves for a 5th order FR-SD discretisation coupled with
4.1. One-dimensional linear advection equation
an explicit Euler time integration is shown. In particular, in order to
quantify the instability of such space/time discretisation, the following
The first numerical test consists in the discretisation of the one-
equivalent frequencies have been considered: 𝜔0 = 0.1 and 𝜔1 = 1.3.
dimensional linear advection equation with an oscillating inlet bound-
The low frequency test is meant to investigate the instability of the
ary condition. This simulation allows to monitor the general behaviour
prescribed discretisation for any choice of computational time step. We
of dissipative bubbles when the prescribed time step approaches the
expect that for any choice of time-step the prescribed sinusoidal signal
CFL limit. Tests are conducted using almost centered numerical fluxes
will grow while propagating in the domain. For very low frequencies
(𝛼 = 0.49, upwinding coefficient) and RK33 or RK45 time integration
the increase will likely be very small, but nonetheless still present. For
schemes for different injected frequencies. Accordingly, Fig. 10 shows
the second frequency, instead, a different behaviour is expected to take
the same curves presented in Fig. 9 in which the tested frequencies
are indicated in order to appreciate the expected level of dissipation of place: for sufficiently small time steps the numerical scheme should
each test. The relevant equivalent frequencies are 𝜔1 = 1, 𝜔2 = 2.75 dissipate the inlet oscillations, but, for increasing time steps, numerical
and 𝜔3 = 3.5 (only for the RK33 test), which are highlighted in Fig. 10 dissipation should decrease, until it becomes anti-dissipation, leading
using dashed vertical lines. Close to 𝜔2 and 𝜔3 , the expected levels of the numerical solution to grow while propagating in the domain. The
numerical dissipation are highlighted using blue dots. sudden change of behaviour, from dissipative to anti-dissipative is
The domain 𝛺 = [0, 1] is split into 100 equally spaced elements predicted to take place at 𝛥𝑡 ≈ 0.01 from the theoretical analysis. The
(ℎ = 0.01). The frequency of the inlet signal can then be directly intersection between the dissipative curve for 𝛥𝑡 ≈ 0.01 and the 𝑥-axis is
linked to the equivalent frequency 𝜔ℎ∕(𝑁 + 1). The simulation matches highlighted in Fig. 16. Consequently, after proper rescaling by the grid
the theoretical analysis: a 5th order SD discretisation, coupled with size of the present simulation (ℎ = 0.01), we expect to observe such
a RK33/RK45 time integration schemes has been employed. To see change of trend for time steps around 𝛥𝑡 = 1.0 × 10−4 . Three time-step
the influence of the temporal errors, increasing values of 𝜏 have been values are therefore tested near 𝛥𝑡 = 1.0×10−4 , as shown in Fig. 17. The
considered. First, in Fig. 11, the velocity signal 𝑢(𝑥) is plotted along numerical solution behaves as expected, where a sudden change from
the 𝑥 direction for 𝜔 = 𝜔1 and for different values of 𝜏 using a RK33 dissipative to anti-dissipative character is observed for increasing time
time marching scheme. According to the theoretical analysis presented steps. For 𝛥𝑡 = 9.925 × 10−5 , the prescribed signal at the inlet is barely
in the first part of the paper, for increasing values of 𝜏, larger values affected by numerical dissipation. For smaller time steps, numerical
of dissipation are expected (see Fig. 10). Such conclusion is confirmed dissipation increases and the sinusoidal signal is expected to decrease in
by the numerical simulations shown in Fig. 11. The amplitude of the amplitude as numerically confirmed in the top plot of Fig. 17. Finally,
signal 𝑢(𝑥) tends, in fact, to decrease for larger time steps. The same for time steps slightly larger, instead, negative values of the imaginary
experiment has been repeated using a RK45 scheme and the results part of 𝜃̃ arise and the solution is expected to grow over space, as
are shown in Fig. 12. The same conclusions apply: increasing the time confirmed in the bottom plot of Fig. 17. Numerical results are then
step magnitude leads to an increase of numerical dissipation. Notice, in very good agreement with the theoretical analysis herein presented
furthermore, that for 𝜏 = 0.9 the signal is more dissipated using the which is able to properly predict the effect of numerical dissipation
RK33 scheme rather than the RK45 one. In fact, from Fig. 10, it can be in a fully-discrete framework where both spatial and temporal errors
noticed that the dissipative curves, in proximity of 𝜔1 , tend to increase contribute on the overall accuracy and regularity of the numerical
a bit more using the RK33 scheme rather than the RK45 one. The same solution.

9
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 10. Physical dissipative curves varying 𝜏 (FR-SD) using almost centered numerical fluxes. Dashed, vertical lines denote the frequencies 𝜔1 , 𝜔2 and 𝜔3 (only for RK33). For
𝜔2 and 𝜔3 the dots indicate the levels of dissipation expected for the specific values of 𝜏 used in the simulations.

Fig. 11. Solution signal 𝑢(𝑥) at prescribed frequency 𝜔1 for different values of 𝜏 (RK33).

Fig. 12. Solution signal 𝑢(𝑥) at prescribed frequency 𝜔1 for different values of 𝜏 (RK45).

Fig. 13. Solution signal 𝑢(𝑥) at prescribed frequency 𝜔2 for different values of 𝜏 (RK33).

10
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 14. Solution signal 𝑢(𝑥) at prescribed frequency 𝜔2 for different values of 𝜏 (RK45).

Fig. 15. Solution signal 𝑢(𝑥) at prescribed frequency 𝜔3 for different values of 𝜏 (RK33).

Fig. 16. Physical dissipative curves varying 𝜏 (FR-SD) using Euler time integration. On the left, a closer look near the origin. Vertical lines indicate the frequencies that have been
numerically simulated.

Regarding the simulation of the lowest frequency, in order to allow marginally stable for pure advection, a well-known deficiency of this
a sufficient large number of periods within the unitary domain, the method [51]. Still, for under-resolved problems (i.e. large frequencies),
number of elements has been increased to 1000. In Fig. 18 the results the fully-discrete scheme revealed to be stable, probably due to the
of the simple advection of the sinusoidal signal with lowest frequency increased upwind dissipation acting against the anti-dissipation of the
is shown for the explicit Euler scheme with a considerably smaller time explicit Euler method. In either case, once the solution establishes itself
step, namely 𝛥𝑡 = 1.0×10−6 = 5×10−3 ℎ∕(𝑁 +1). On the top, the solution in space, it does not grow over time — in line with what is expected
is probed at 𝑥 = 0.5, on the bottom, instead, the numerical solution within the spatial analysis framework: either stability or convective
at the end of the simulation is shown. As predicted by the theoretical instability.
analysis, the solution tends to grow over space. Such behaviour is even
more evident in Fig. 19, where the solution is zoomed close to the 4.2. Accurately resolved flows
line 𝑢 = 1. The same simulation has been performed using an even
smaller time step 𝛥𝑡 = 1.0 × 10−7 = 5 × 10−4 ℎ∕(𝑁 + 1). The probed In the theoretical analysis, presented in the first part of the paper,
and final solution of the numerical simulation is shown in Fig. 20. the discussion was organised making distinction between high and
Even if barely perceptible, the solution tends to grow while advected low frequencies, considering more or less extended upper bounds for
along 𝑥. The same detailed view of Fig. 19 is reproduced in Fig. 21. the frequency 𝜔. In the low frequency regions, the reduced order
In this zoomed visualisation, it is possible to observe more clearly the of accuracy for time steps approaching the CFL limit was observed.
increase in amplitude of the sinusoidal signal while propagating into The convergence results already reported in the temporal eigenanal-
the domain, whereas no growth is observed monitoring the probed ysis by Vermeire et al. [40] were correctly reproduced within the
solution. The above results confirm that the explicit Euler scheme is present fully-discrete spatial eigenanalysis framework (see Fig. 8). Such

11
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 17. Numerical solution of 𝜔0 = 0.1 using explicit Euler time integration scheme. Top, 𝛥𝑡 = 9.825 × 10−5 ; Center 𝛥𝑡 = 9.925 × 10−5 ; Bottom 𝛥𝑡 = 1.0925 × 10−4 .

Fig. 18. Numerical solution of 𝜔1 = 1.3 using explicit Euler time integration scheme with 𝛥𝑡 = 1.0 × 10−6 . Top, probed solution at 𝑥 = 0.5; Bottom, numerical solution at the end
of the simulation. The red horizontal line indicates the iso-line 𝑢 = 1.

Fig. 19. Closer look of the solution of Fig. 18 around 𝑢 = 1 (red line).

feature highlights the similarity between temporal and spatial eigen- theoretical analysis for well-resolved frequencies. A second subsection
analysis for well-resolved wavenumbers/frequencies respectively. For will instead be dedicated to under-resolved flows which are commonly
high frequencies, instead, the fully-discrete spatial eigenanalysis has associated to high-frequencies which are only partially resolved by
shown significant variations from the semi-discrete analysis. In the the numerical scheme. In particular, as a representative example of
present subsection we will perform a series of numerical experiments the spatial eigenanalysis, a variation of the classic isentropic vortex
on accurately resolved nonlinear inviscid flows in order to validate the simulation [56] has been considered. Such formulation is based on

12
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 20. Numerical solution of 𝜔1 = 1.2 using explicit Euler time integration scheme with 𝛥𝑡 = 1.0 × 10−7 . Top, probed solution at 𝑥 = 0.5; Bottom, numerical solution at the end
of the simulation. The red horizontal line indicates the iso-line 𝑢 = 1.

Fig. 21. Closer look of the solution of Fig. 20 around 𝑢 = 1 (red line).

an inflow-type condition and is characterised by spatial grid inhomo- At the initial time, two vortices are centered at 𝑥 = 1∕4𝐿 and
geneity due to a sudden change of mesh resolution in the middle of 𝑥 = 3∕4𝐿, where 𝐿 = 40 denotes the full horizontal length of the
the domain. Clearly, such configuration is particularly representative domain. Each vortex is characterised by the smooth solution:
of spatial eigenanalysis, where an oscillatory flow field is imposed ( ) 1
at the inlet boundary and let evolve along the streamwise direction. 𝛽 2 (𝛾 − 1) 2(1−𝑟2 ) 𝛾−1
𝜌= 1− 𝑒 ,
Finally, in terms of turbulent flows, the study of highly accurate fields, 16𝛾𝜋 2
where all the scales and frequencies are sufficiently well-resolved by 𝛽(𝑦 − 𝑦𝑐 ) 1−𝑟2
𝑢 = 𝑢0 − 𝑒 , (26)
the numerical scheme, is directly linked to the Direct Numerical Simu- 2𝜋
lation (DNS) approach. It is worthwhile mentioning that for the Euler 𝛽(𝑥 − 𝑥𝑐 ) 1−𝑟2
𝑣 = 𝑣0 + 𝑒 ,
equations the characteristic speed of the system might change in the 2𝜋
domain and, differently with respect to the linear advection case, the 𝑝 = 𝜌𝛾 ,
grid size is not uniform as well. Consequently, the maximum allowed √
where 𝑟 = (𝑥 − 𝑥𝑐 )2 + (𝑦 − 𝑦𝑐 )2 , with 𝑥𝑐 and 𝑦𝑐 indicating the co-
time step depends on the local flow characteristics and grid size. In ordinates of the vortex’s center, and 𝛾 = 1.4 indicates the specific
order to estimate this quantity, in the test cases herein considered, we heat ratio. The parameter 𝛽 defines the strength of the vortex and, for
simply picked increasing constant time steps until reaching the point this particular study, a value equal to 5 has been chosen, providing a
of instability. Once the maximum time step is so defined, we evaluated relatively strong vortex. Subsequently, the solution is advected by the
fractions of it in order to study the influence of time integration. constant velocity field (𝑢0 , 𝑣0 ) = (5, 0) while a constant stream of new
Finally, in order to accentuate time integration errors,√ we evaluated an isentropic vortices is imposed on the left boundary. After a full period
estimate of the maximum allowed time step as (‖𝐮‖+ 𝛾𝑝∕𝜌)−1 ℎ∕(𝑁 +1) has passed, two new vortices have entered the domain from the inlet
and verified for each test case that a sufficiently extensive region of the boundary, whereas the two initial ones have left the domain through
domain was characterised by large values of this quantity. the outlet boundary.
The smooth analytical solution of the isentropic vortex test case pro-
4.2.1. Inlet-outlet advection of isentropic Euler vortex vides an ideal set-up to evaluate the convergence of numerical schemes
A similar numerical set-up with respect to the classical isentropic for the discretisation of Euler equations. The numerical solution of
vortex test case has been considered in order to avoid the use peri- Euler equations has been considered on a series of increasingly refined
odic boundary conditions, which are classically associated to temporal grids. Namely, the number of elements of the different meshes reads
eigenanalysis. A formulation considering inlet/outlet boundaries is, 𝑁el = 𝑁𝑥 × 𝑁𝑦 = 2𝑁𝑦 × 𝑁𝑦 with 𝑁𝑦 = 10, 20, 40, 80, 160. In order
instead, a much more representative example of spatial eigenanalysis, to highlight numerical errors, both regular and deformed grids have
where a temporal frequency is imposed to the system at the inlet and been considered. An example of deformed mesh is shown in Fig. 22.
the evolution of the numerical solution is studied while propagating Furthermore, in order to study the evolution of the initial solution while
through the space. propagating in space, a sudden change on streamwise grid resolution

13
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 22. Example of deformed computational grid used in the isentropic vortex simulation.

has been placed in the middle of the domain (see Fig. 22). This way, 4.3. Under-resolved flows
the three-fold influence of inlet, outlet and mesh coarsening can be
probed in a single test. Finally, in order to keep small the numerical As mentioned in the previous section, a first comparative study
errors associated to the spatial discretisation, relatively high orders of between spatial and temporal eigenanalysis can be performed consider-
approximation have been considered. In particular 6th and 8th order ing well-resolved flow fields. In such conditions, both theories should
accurate Spectral Difference simulations have been performed, coupled return the expected decrease of spatial accuracy for large time steps
with a 3rd order Runge–Kutta scheme as time integration scheme. In close to the CFL limit. In the previous section such behaviour has been
order to evaluate the influence of temporal numerical errors, the time
numerically confirmed. In this second test, instead, high frequencies
step has been chosen such that 𝜏 = 0.1, 0.5, 0.9. The goal of the present
and wavenumbers (i.e. under-resolved flow fields) will be considered.
test case consists in showing the essential equivalence between spatial
In contrast to the previous section, where well-resolved fields have
and temporal approaches in the well-resolved regions (the former for
been naturally associated to DNS, now the behaviour of the numerical
well-resolved wavenumbers, the latter for well-resolved frequencies).
scheme represents Large Eddy Simulations (LES, and in particular Im-
When we prescribe the stream of vortices through the inlet condition,
plicit LES), where only the large-scale dynamics is adequately resolved
given the large number of degrees of freedom (DOF) employed in
by the number of DOF employed. According to the theoretical analysis,
the tests (high polynomial order), we are in fact working in the low-
for cases employing standard upwind fluxes, large time step values
frequency, well-resolved region of spatial analysis. The 𝐿2 error of
should considerably decrease the overall level of numerical dissipation,
the density field is easily evaluated and a convergence study can be
affecting the simulation both in terms of stability and accuracy. To
performed similarly to the classical test case (whose results are omitted
perform these analyses, the classical Taylor–Green vortex and a two-
for the sake of brevity).
dimensional duct flow have been herein considered as representative
The convergence plots for the 6th and 8th order simulations in
test cases.
conjunction with RK33 are shown in Fig. 23. These plots very closely
resemble those found for the classical (periodic) case. We note that
time steps close to the CFL limit slow down the error convergence. 4.3.1. Taylor–Green vortex test case
In particular, the expected convergence rate arising from the spatial The Taylor–Green Vortex (TGV) constitutes a well-established test
discretisation is eventually reduced to the order of the time integration case to study vortex dynamics, turbulent transition, turbulent decay and
scheme. When meshes of stronger deformation are considered, the same energy dissipation processes in a three-dimensional setting [57].
behaviour is observed, but at smaller grid sizes. This is believed to The problem consist in a cubic domain [−𝐿𝜋, 𝐿𝜋]3 with periodic
happen because of the augmented spatial numerical errors, which can boundary conditions applied to all faces. The smooth initial condition
hide the effect of large time steps. In order to illustrate the influence reads
of the time step on the evolution of the vortex, the absolute value
⎧ 𝜌 = 𝜌0 ,
of the thermodynamic entropy 𝑠 = log(𝜌−𝛾 𝑝) has been employed as ⎪ ( ) ( ) ( )
an indicator of numerical errors. Of course, since an inviscid, shock- ⎪ 𝑥 𝑦 𝑧
⎪𝑢1 = 𝑈0 sin 𝐿 cos 𝐿 cos 𝐿 ,
free, Euler flow is being considered, the entropy should be theoretically ⎪ ( ) ( ) ( )
conserved everywhere in the domain. Any deviation from a constant ⎪ 𝑥 𝑦 𝑧
⎨𝑢2 = −𝑈0 cos 𝐿 sin 𝐿 cos 𝐿 , (27)
entropy state is then entirely associated to numerical artifacts. The ⎪
entropy fields at the final time using 𝜏 = 0.9 and 𝜏 = 0.99 are shown ⎪𝑢 = 0,
in Fig. 24, for the 6th order discretisation coupled with RK33. It is ⎪ 3 [ ( ) ( )][ ( ) ]
⎪ 𝜌0 𝑈 0 2𝑥 2𝑦 2𝑧
interesting to notice that in most of the domain the entropy fields of ⎪ 𝑝 = 𝑃0 + cos + cos cos + 2 .
⎩ 16 𝐿 𝐿 𝐿
the two simulations are quite similar. The only relevant differences are
represented by the highly irregular spurious oscillations close to the Unity has been assigned to both 𝑈0 and 𝜌0 (reference velocity and
line 𝑦 = 0. These are stronger in the second part of the domain, as they density, respectively) and the initial value of the pressure 𝑃0 has been
are connected with spurious transmitted waves (only barely dissipated chosen such that the corresponding initial Mach number is equal to 0.1.
as 𝜏 → 1). Smaller spurious waves are reflected into the first part of The characteristic time scale of the problem is defined as 𝑇 = 𝐿∕𝑈0 . For
the domain too, but are more quickly dissipated as they propagate this value of the Mach number, the flow is practically incompressible.
backwards. In agreement with the work by Vermeire et al. [40], a non- The computational domain is subdivided in 323 uniform cubic elements
negligible influence of the time step choice is observed only very close and discretised with a 6th order SD scheme. A 3rd order accurate
to the CFL limit. Runge–Kutta method has been used to advance the solution in time.

14
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 23. Density’s convergence when using RK33. Black line, 𝜏 = 0.1; purple line 𝜏 = 0.5; orange line 𝜏 = 0.9. Solid line, structured grid; dashed line, deformed grid.

Fig. 24. Absolute value of the thermodynamic entropy 𝑠 = ln(𝜌−𝛾 𝑝) throughout the domain for the 6th order simulation coupled with RK33 using 𝜏 = 0.9 (top) and 𝜏 = 0.99
(bottom). Vertical white line denotes the location of the mesh coarsening.

For inviscid flows, the kinetic energy balance (under the periodic The time evolution of numerical dissipation is shown in Fig. 25.
boundary condition applied here) can be written as Before 𝑡∕𝑇 ≈ 4 the flow field is well resolved and the numerical
( ) ( ) dissipation takes very small values. As the flow field starts to transit
𝑑 1 𝜕𝑢𝑖 to turbulence, more and more unresolved scales naturally develop. The
− 𝜌𝑢𝑖 𝑢𝑖 𝑑𝑉 = − 𝑝 𝑑𝑉 + 𝜀num . (28)
𝑑𝑡 2 ∫ ∫ 𝜕𝑥𝑖 transfer of kinetic energy to high wavenumbers is then dissipated by the
⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ⏟⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏟ numerical scheme. Numerical dissipation, in fact, tends to peak around
Variation of kinetic energy Pressure-dilatation
𝑡∕𝑇 = 10, when the flow is fully turbulent.
The total amount of numerical dissipation can be estimated from In Fig. 25, numerical dissipation is plotted for three values of 𝜏
(namely, 𝜏 = 0.5, 𝜏 = 0.9 and 𝜏 = 0.95). However, in terms of 𝜀num ,
the previous balance subtracting the pressure-dilatation work from
the three simulations are almost indistinguishable between each other.
temporal variation of the kinetic energy. Moreover, due to the low Stability and accuracy of the numerical simulation are affected by the
Mach number, pressure-dilatation work is negligible and numerical time step size very abruptly approaching the CFL limit. No indication
dissipation is the main mechanism in the kinetic energy dynamics. of instability is noticeable before this sudden threshold is reached. In

15
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 25. Numerical dissipation 𝜀num for different values of 𝜏. Solid line, 𝜏 = 0.5; dashed line, 𝜏 = 0.9; dash-dotted line 𝜏 = 0.95.

other words, the gradual decrease of numerical dissipation observed in coming from the outflow boundary and from the line of mesh size
the theoretical framework does not manifest in the present numerical change [29,30,32], these results typically rely on simulations that are
simulation. well-resolved in time. The main goal of the present test is to check
whether such behaviour is affected by time step size (as predicted by
4.3.2. Two-dimensional grid turbulence theory) or not. The test case has been considered for 𝜏 = 0.5, 0.9 and
As a typical sample of non-periodic test case, a simple duct flow with 0.99.
oscillating inlet velocity has been considered. The domain consists in a Similar to the previous test case (isentropic vortex), the time step
20𝜋 × 2𝜋 rectangle. Although two-dimensional and inviscid, the present influence is very small from a macroscopic point of view. All cases
test case is sufficiently insightful to quantify the relevant role of the computed with 𝜏 < 0.95 (approximately) behaved quite similarly. Only
numerical scheme for under-resolved flows, as already shown in many in very close proximity of the CFL limit, spurious oscillations arose
different works [29–32]. from the numerical solution. As in the previous case, such behaviour is
The inlet boundary condition reads qualitatively represented in Fig. 28 using the thermodynamic entropy
field. Again, spurious transmitted waves are seen when 𝜏 is sufficiently
𝜌 = 𝜌∞ , close to unity. In addition, the line of mesh size change seems to trigger
( ) more spurious reflected waves than the outflow boundary (at least
1 + 𝐴 sin(𝐾𝑦) sin(𝛺𝑡)
𝐮 = 𝑢∞ , for this test problem). Only for values of 𝜏 close to unity, spurious
0
entropy oscillations take place. Due to the structured nature of the
𝑝 = 𝑝∞ , mesh, spurious oscillations arise diffusively in the whole domain. In
where 𝜌∞ = 1, 𝑢∞ = 1 and 𝑝∞ = (𝜌∞ 𝑢2∞ )∕(𝛾Ma2 ) in order to get a fact, a large number of cells reaches critical values of 𝛥𝑡 simultaneously,
desired value of inflow Mach number. The parameters defining the producing spurious entropy oscillations. The same conclusions obtained
from isentropic vortex case seem to apply here: the spatial under-
inflow perturbations have been set as 𝐴 = 1∕2, 𝐾 = 5, 𝛺 = 1. The
resolution hides/overshadows the influence of the prescribed time step.
prescribed Mach number is set equal to 0.03. For this particular choice,
The influence of the time integration scheme is then noticeable only
the flow is then essentially incompressible. Non-reflective Roe flux-
very close to instability, whereas almost any time step fulfilling the
based boundary conditions are employed at the outlet, whereas top
CFL condition does not impact the overall accuracy of the numerical
and bottom boundaries are modelled as free-slip walls. The results have
solution (already limited by the spatial discretisation). The fact itself
been obtained with a 6th-order SD scheme coupled with RK33.
that the presented theory does not fully predict nonlinear dynam-
In order to trigger spurious numerical oscillations, a sudden change ics brings useful information for practical CFD applications. Whereas
in streamwise mesh spacing is imposed at 𝑥 = 12𝜋, as shown in Fig. 26. classical semi-discrete analyses, although simplified, usually provide
The number of elements along the cross-flow direction is 𝑁el𝑦 = 12. insightful general knowledge on the diffusive/dispersive properties
𝑥
Along the streamwise direction, this number is 𝑁el1 = 72 in the first of spatial discretisations, much different conclusions can be drawn
𝑥
block and 𝑁el2 = 12 in the second block. A typical flow field is when time integration errors are involved. Wave propagation, in fact,
illustrated in Fig. 27, where the two-dimensional nature of the flow behaves very differently in linear and nonlinear dynamics, where in
leads to a well-known reverse cascade process, with small vortices even- the latter the full spectrum of temporal and spatial frequencies are
tually merging into larger ones. Regardless, the mechanism by which strongly coupled with each other. Future research will focus on further
kinetic energy is exchanged across scales is significantly influenced by generalisations of spectral eigenanalyses, for example by introducing
the underlying numerical errors, whereby the present under-resolved quasi-linear/nonlinear operators in order to fully capture the strong
test case is insightful with regards to implicit LES. Although Roe-type coupling between spatial and temporal errors in the simulation of
solvers have been found to successfully suppress spurious oscillations under-resolved turbulent flows.

16
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 26. Example of computational grid used in the duct flow simulation.

Fig. 27. Example of duct flow field (velocity magnitude). The white vertical line denotes the change of mesh size interface.

Fig. 28. Entropy field. The white vertical line denotes the change of mesh size interface.

5. Conclusions experiments have shown good agreement with the theory where the
time step size has shown to enhance numerical dissipation for certain
A fully-discrete generalisation of the spatial eigenanalysis has been frequencies and suppress it for others. Nonlinear simulations using
presented and applied to the Spectral Difference (SD) and Discontin- Euler equations have been performed to quantify the interplay between
uous Galerkin (DG) schemes. Whereas previous research focused on spatial and temporal errors for more complex dynamics. In agreement
fully-discrete generalisations of standard temporal eigenanalysis for
with the previous discussion, the Euler simulations have been grouped
similar schemes, in this work we presented the first fully-discrete ver-
in well-resolved and under-resolved flows. An inlet-outlet formulation
sion of the spatial eigenalaysis for high-order methods, complementing,
of the classical isentropic vortex test case have been used to vali-
in a way, the work of Vermeire and Vincent [40] on temporal Von
Neumann analysis. A detailed discussion has been presented on spatial date fully-discrete temporal and spatial eigenanalyses for well-resolved
dissipative curves and their dependence on the time step size for 1st- flows. As predicted by the theory, the expected convergence rate, for
to 4th-order explicit Runge–Kutta time integration schemes. The be- sufficiently large time steps, is significantly polluted by temporal errors.
haviour for both well-resolved and under-resolved frequencies has been
considered. In the fully-discrete spatial eigenanalysis, it has been shown In the case of under-resolved flows, instead, the influence of the
that in the low frequency range, in similarity with the fully-discrete time step size has shown to be negligible. Spurious oscillations become
temporal eigenanalysis, the theoretical spatial convergence rate de- large enough to visibly pollute the solution only for very large time
grades to the local order of accuracy of the time integration scheme. steps approaching the CFL limit. In under-resolved flows, in fact, the
In the high wavenumbers/frequencies regions, instead, temporal and spatial errors dominate their temporal counterparts, overshadowing the
spatial eigenanalyses have significantly different behaviours. Both are
influence of the time integration scheme. In conclusion, the theoretical
influenced by the time step size and in both cases numerical dissipation
framework introduced in the first part of the paper is not always
generally decreases for larger time steps employing upwind fluxes. In
representative of more complex nonlinear dynamics. For well-resolved
the case of almost centered fluxes, the dissipative bubbles, typical of
spatial eigenanalysis of high-order methods, tend to be smoother for flows, the slower convergence has been observed in both the theoretical
increasing time step size. Regions that experience negligible numerical framework and in the Euler flow simulations. Consequently, the choice
dissipation in the semi-discrete analysis are characterised by stronger of the time step size plays a more relevant role in highly accurate
diffusion in the fully-discrete generalisation. simulations such as DNS of turbulent flows. Instead, in the presence
The SD discretisation of the one-dimensional linear advection equa- of significant spatial under-resolution (like in LES), temporal errors are
tion has been used to validate the theoretical findings. Numerical more likely overshadowed by spatial errors.

17
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 29. Dissipation of the FR-DG scheme for RK33 (left) and RK45 (right) at polynomial orders 𝑁 = 1, 2, 3, 6 (top to bottom).

18
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 30. Dissipation of the FR-SD scheme for RK33 (left) and RK45 (right) at polynomial orders 𝑁 = 1, 2, 3, 6 (top to bottom).

19
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 31. Numerical dissipation for a 3rd order FR-SD discretisation using centred numerical fluxes coupled with RK33 (left) and RK45 (right) time integration schemes. The thinner
curves closer to the dashed line represent the semi-discrete result (no temporal error).

20
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

Fig. 32. Numerical dissipation for a 5th order FR-SD discretisation using centred numerical fluxes coupled with RK33 (left) and RK45 (right) time integration schemes. The thinner
curves closer to the dashed line represent the semi-discrete result (no temporal error).

CRediT authorship contribution statement Data availability

Niccolò Tonicello: Conceptualization, Methodology, Software, In- No data was used for the research described in the article.
vestigation. Rodrigo C. Moura: Formal analysis, Writing – review &
editing. Guido Lodato: Software, Writing – review & editing. Gian- Acknowledgements
marco Mengaldo: Writing – review & editing.
The use of the SD solver originally developed by Antony Jameson’s
Declaration of competing interest group at Stanford University is gratefully acknowledged. This work was
granted access to the high-performance computing resources of CRI-
The authors declare that they have no known competing finan- ANN. RCM acknowledges support from São Paulo Research Foundation
cial interests or personal relationships that could have appeared to (FAPESP), Brazil via grant 2020/10910-8. GL acknowledges support
influence the work reported in this paper. from the Agence Nationale de la Recherche (ANR), France under grant

21
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

number ANR-18-CE05-0030. GM acknowledges support from National [23] Bogey C, Bailly C. A family of low dispersive and low dissipative explicit schemes
University of Singapore, Singapore start-up grant WBS R-265-000-A36- for flow and noise computations. J Comput Phys 2004;194(1):194–214.
[24] Van den Abeele K, Lacor C, Wang ZJ. On the stability and accuracy of the
133. The authors also thank David Kopriva, Tristan Montoya and David
spectral difference method. J Sci Comput 2008;37(2):162–88.
Zingg for fruitful discussions. [25] Vincent PE, Castonguay P, Jameson A. Insights from von Neumann analysis of
high-order flux reconstruction schemes. J Comput Phys 2011;230(22):8134–54.
Appendix [26] Moura RC, Sherwin SJ, Peiró J. Linear dispersion–diffusion analysis and its ap-
plication to under-resolved turbulence simulations using discontinuous Galerkin
spectral/hp methods. J Comput Phys 2015;298:695–710.
A series of plots for the dissipation of the FR-DG scheme are given [27] Vanharen J, Puigt G, Vasseur X, Boussuge J-F, Sagaut P. Revisiting the spec-
in Fig. 29, whereas plots for the FR-SD scheme are given in Fig. 30. tral analysis for high-order spectral discontinuous methods. J Comput Phys
The SD curves look similar to those of DG, but are moderately shifted 2017;337:379–402.
[28] Hu FQ, Atkins HL. Eigensolution analysis of the discontinuous Galerkin
up (stronger dissipation, note the vertical axis) and to the left (towards
method with nonuniform grids: I. one space dimension. J Comput Phys
lower frequencies). 2002;182(2):516–45.
In Fig. 31, the evolution of the multiple modes of the (nearly) [29] Mengaldo G, Moura R, Giralda B, Peiró J, Sherwin S. Spatial eigensolution anal-
centred FR-SD with RK33 and RK45 are shown for 𝑁 = 2. Similar plots ysis of discontinuous Galerkin schemes with practical insights for under-resolved
for 𝑁 = 4 are shown in Fig. 32. computations and implicit LES. Comput & Fluids 2018;169:349–64.
[30] Mengaldo G, De Grazia D, Moura RC, Sherwin SJ. Spatial eigensolution analysis
of energy-stable flux reconstruction schemes and influence of the numerical flux
References on accuracy and robustness. J Comput Phys 2018;358:1–20.
[31] Moura RC, Aman M, Peiró J, Sherwin SJ. Spatial eigenanalysis of spec-
[1] Slotnick J, Khodadoust A, Alonso J, Darmofal D, Gropp W, Lurie E, et al. CFD tral/hp continuous Galerkin schemes and their stabilisation via DG-mimicking
vision 2030 study: A path to revolutionary computational aerosciences. Tech. spectral vanishing viscosity for high Reynolds number flows. J Comput Phys
rep., Hampton, Virginia 23681-2199: NASA Langley Research Center; 2014, 2020;406:109112.
NASA Contractor Report 218178. [32] Tonicello N, Lodato G, Vervisch L. A comparative study from spectral analyses
[2] Vincent PE, Jameson A. Facilitating the adoption of unstructured high-order of high-order methods with non-constant advection velocities. J Sci Comput
methods amongst a wider community of fluid dynamicists. Math Model Nat 2021;87(3):1–38.
Phenom 2011;6(3):97–140. [33] Trojak W, Watson R, Scillitoe A, Tucker PG. Effect of mesh quality on flux
[3] Hesthaven JS, Warburton T. Nodal discontinuous Galerkin methods: Algorithms, reconstruction in multi-dimensions. J Sci Comput 2020;82(3):1–36.
analysis, and applications. Springer Science & Business Media; 2007. [34] Sáez-Mischlich G, Sierra-Ausín J, Gressier J. The spectral difference raviart-
[4] Cockburn B, Shu C. The local discontinuous Galerkin finite element method for Thomas method for two and three-dimensional elements and its connection with
convection-diffusion systems. SIAM J Numer Anal 1998;35:2440–63. the flux reconstruction formulation. 2021, arXiv preprint arXiv:2105.08632.
[5] Cockburn B, Shu C. The Runge-Kutta discontinuous Galerkin finite element [35] Manzanero J, Rubio G, Ferrer E, Valero E. Dispersion-dissipation analysis for
method for conservation laws V: Multidimensional systems. J Comput Phys advection problems with nonconstant coefficients: Applications to discontinuous
1998;141:199–224. Galerkin formulations. SIAM J Sci Comput 2018;40(2):A747–68.
[6] Kopriva DA, Kolias JH. A conservative staggered-grid Chebyshev multidomain [36] Manzanero J, Rubio G, Ferrer E, Valero E, Kopriva DA. Insights on aliasing
method for compressible flows. J Comput Phys 1996;125(1):244–61. driven instabilities for advection equations with application to Gauss-Lobatto
[7] Liu Y, Vinokur M, Wang Z. Spectral difference method for unstructured grids I: discontinuous Galerkin methods. 2017.
[37] Fernandez P, Moura RC, Mengaldo G, Peraire J. Non-modal analysis of spectral
basic formulation. J Comput Phys 2006;216(2):780–801.
element methods: Towards accurate and robust large-eddy simulations. Comput
[8] Wang Z, Liu Y, May G, Jameson A. Spectral difference method for unstructured
Methods Appl Mech Engrg 2019;346:43–62. http://dx.doi.org/10.1016/j.cma.
grids II: Extension to the Euler equations. J Sci Comput 2007;32(1):45–71.
2018.11.027.
[9] Huynh HT. A flux reconstruction approach to high-order schemes including
[38] Alhawwary M, Wang Z. A combined-mode Fourier analysis of DG methods for
discontinuous Galerkin methods. In: 18th AIAA computational fluid dynamics
linear parabolic problems. SIAM J Sci Comput 2020;42(6):A3825–58.
conference. 2007, p. 4079.
[39] Yang H, Li F, Qiu J. Dispersion and dissipation errors of two fully discrete
[10] Vincent PE, Castonguay P, Jameson A. A new class of high-order energy stable
discontinuous Galerkin methods. J Sci Comput 2013;55(3):552–74.
flux reconstruction schemes. J Sci Comput 2011;47(1):50–72.
[40] Vermeire B, Vincent P. On the behaviour of fully-discrete flux reconstruction
[11] Chapelier J-B, M D, Renac F, Lamballais E. Evaluation of a high-order dis-
schemes. Comput Methods Appl Mech Engrg 2017;315:1053–79.
continuous Galerkin method for the DNS of turbulent flows. Comput & Fluids
[41] Vanharen J, Puigt G, Vasseur X, Boussuge J-F, Sagaut P. Revisiting the spec-
2014;95:210–26.
tral analysis for high-order spectral discontinuous methods. J Comput Phys
[12] de la Llave Plata M, Couaillier V, Le Pape M-C. On the use of a high-order
2017;337:379–402.
discontinuous Galerkin method for DNS and LES of wall-bounded turbulence.
[42] Alhawwary M, Wang ZJ. Fourier analysis and evaluation of DG, FD and compact
Comput & Fluids 2018;176:320–37.
difference methods for conservation laws. J Comput Phys 2018;373:835–62.
[13] Bassi F, Botti L, Colombo A, Crivellini A, Ghidoni A, Massa F. On the devel-
[43] De Grazia D, Mengaldo G, Moxey D, Vincent P, Sherwin S. Connections between
opment of an implicit high-order discontinuous Galerkin method for DNS and
the discontinuous Galerkin method and high-order flux reconstruction schemes.
implicit LES of turbulent flows. Eur J Mech B Fluids 2016;55:367–79.
Internat J Numer Methods Fluids 2014;75(12):860–77.
[14] Krank B, Kronbichler M, Wall WA. Direct numerical simulation of flow over [44] Mengaldo G. Discontinuous spectral/hp element methods: Development, analysis
periodic hills up to Re𝐻 = 10, 595. Flow Turbul Combust 2018;101(2):521–51. and applications to compressible flows [Ph.D. thesis], Imperial College London;
[15] Lodato G, Tonicello N, Pinto B. Large-eddy simulation of bypass transition on a 2015.
zero-pressure-gradient flat plate using the spectral-element dynamic model. Flow [45] Mengaldo G, De Grazia D, Vincent PE, Sherwin SJ. On the connections between
Turbul Combust 2021;1–30. discontinuous Galerkin and flux reconstruction schemes: extension to curvilinear
[16] Tonicello N, Lodato G, Vervisch L. Analysis of high-order explicit LES dynamic meshes. J Sci Comput 2016;67(3):1272–92.
modeling applied to airfoil flows. Flow Turbul Combust 2021;1–28. [46] Asthana K, Jameson A. High-order flux reconstruction schemes with minimal
[17] Fernandez P, Nguyen N-C, Peraire J. On the ability of discontinuous Galerkin dispersion and dissipation. J Sci Comput 2015;62(3):913–44.
methods to simulate under-resolved turbulent flows. 2018, arXiv preprint arXiv: [47] Gottlieb S, Shu C-W, Tadmor E. Strong stability-preserving high-order time
1810.09435. discretization methods. SIAM Rev 2001;43(1):89–112.
[18] Moura RC, Mengaldo G, Peiró J, Sherwin SJ. On the eddy-resolving capability [48] Butcher JC. The numerical analysis of ordinary differential equations:
of high-order discontinuous Galerkin approaches to implicit LES/under-resolved Runge-Kutta and general linear methods. Wiley-Interscience; 1987.
DNS of Euler turbulence. J Comput Phys 2017;330:615–23. [49] Hairer E, Wanner G, Solving O. II: Stiff and differential-algebraic problems. Berlin
[19] de la Llave Plata M, Lamballais E, Naddei F. On the performance of a high-order [etc.]: Springer; 1991.
multiscale DG approach to LES at increasing Reynolds number. Comput & Fluids [50] Carpenter MH, Kennedy CA. Fourth-order 2N-storage Runge-Kutta schemes. Tech.
2019;194:104306. rep., NASA Langley Research Center; 1994.
[20] Mengaldo G, Moxey D, Turner M, Moura RC, Jassim A, Taylor M, et al. [51] Lomax H, Pulliam TH, Zingg DW. Fundamentals of computational fluid dynamics.
Industry-relevant implicit large-eddy simulation of a high-performance road car Springer Science & Business Media; 2013.
via spectral/hp element methods. SIAM Rev 2021;63(4):723–55. [52] Moura RC, Sherwin SJ, Peiró J. Eigensolution analysis of spectral/hp continuous
[21] Cox C, Trojak W, Dzanic T, Witherden F, Jameson A. Accuracy, stability, and Galerkin approximations to advection–diffusion problems: Insights into spectral
performance comparison between the spectral difference and flux reconstruction vanishing viscosity. J Comput Phys 2016;307:401–22.
schemes. Comput & Fluids 2021;221:104922. [53] Liang C, Premasuthan S, Jameson A. High-order accurate simulation of low-mach
[22] Lele SK. Compact finite difference schemes with spectral-like resolution. J laminar flow past two side-by-side cylinders using spectral difference method.
Comput Phys 1992;103(1):16–42. Comput Struct 2009;87(11–12):812–27.

22
N. Tonicello et al. Computers and Fluids 266 (2023) 106060

[54] Liang C, Jameson A, Wang ZJ. Spectral difference method for compress- [56] Shu C-W. Essentially non-oscillatory and weighted essentially non-oscillatory
ible flow on unstructured grids with mixed elements. J Comput Phys schemes for hyperbolic conservation laws. In: Advanced numerical approximation
2009;228(8):2847–58. of nonlinear hyperbolic equations. Springer; 1998, p. 325–432.
[55] Jameson A. A proof of the stability of the spectral difference method for all [57] Taylor GI, Green AE. Mechanism of the production of small eddies from large
orders of accuracy. J Sci Comput 2010;45(1):348–58. ones. Proc R Soc Lond Ser A Math Phys Eng Sci 1937;158(895):499–521.

23