# Fully-discrete spatial eigenanalysis of discontinuous spectral element
methods: insights into well-resolved and under-resolved vortical flows

Niccolò Tonicello<EMAIL_ADDRESS>Rodrigo C. Moura  Guido Lodato
Gianmarco Mengaldo  Normandie Université, INSA et Université de Rouen, St.
Etienne du Rouvray (Rouen), France  Instituto Tecnológico de Aeronáutica, São
José dos Campos, Brazil  National University of Singapore, Singapore

###### Abstract

This study presents a comprehensive spatial eigenanalysis of fully-discrete
discontinuous spectral element methods, now generalizing previous spatial
eigenanalysis that did not include time integration errors. The influence of
discrete time integration is discussed in detail for different explicit Runge-
Kutta (1st to 4th order accurate) schemes combined with either Discontinuous
Galerkin (DG) or Spectral Difference (SD) methods, both here recovered from
the Flux Reconstruction (FR) scheme. Selected numerical experiments using the
improved SD method by Liang and Jameson [1, 2, 3] 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.
The effect of mesh regularity is also considered, where time integration
errors are found to be, in the case of irregular grids, less pronounced than
those of the spatial discretisation. Still, 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. In fact, for the under-resolved vortical flows considered, the
predominance of spatial errors made it practically impossible for time
integration errors to be distinctly identified. Nevertheless, for well-
resolved nonlinear simulations, the effect of time integration errors could
still be recognized. 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 behavior of
numerical errors in practical under-resolved nonlinear problems, including
under-resolved turbulence computations.

###### keywords:

High-order methods, spectral element methods, flux reconstruction,
discontinuous Galerkin, spectral difference, eigenanalysis.

††journal: Journal of Computational Physics

##  1 Introduction

Along with the increasing computational power experienced in the last few
decades, the use of Computational Fluid Dynamics (CFD) is becoming a more and
more common tool in industrial aerodynamic analysis and design [4]. The
expanding capabilities of modern computers have fuelled the development of a
large number of innovative numerical high-order schemes for fluid dynamics,
whose accuracy clearly surpasses that of more popular low order methods. In
particular, much interest has been focused on spectral element methods due to
their geometric flexibility and promising computational efficiency [5]. These
include the Discontinuous Galerkin (DG) [6, 7, 8], Spectral Difference (SD)
[9, 10, 11] and Flux Reconstruction (FR) [12, 13] schemes. The use of high-
order numerical schemes has shown promising results in the simulation of
turbulent flows, in both Direct Numerical Simulation (DNS) [14, 15, 16, 17]
and Large-Eddy Simulations (LES) [18, 19, 20, 21, 22, 23]. Due to the rapidly
developing field of high-order methods, especially in fluid dynamics, the
research community is constantly looking for more efficient and reliable tools
to quantify the peculiar numerical properties of such schemes. Indeed, for
under-resolved simulations (such as LES), the role played by the numerical
scheme can be of fundamental importance. For instance, the same computation,
performed with different methods, can lead to dramatically different outcomes
in realistic flow configurations [24, 21, 20].

Spectral analyses based on the numerical discretisation of the linear
advection equation are a common tool in the evaluation of the intrinsic
numerical properties of high-order methods for advection dominated flows. Such
techniques are particularly simple to implement and computationally
inexpensive. However simplistic, they are capable of providing useful
information regarding the scheme, such as numerical dispersion and diffusion.
The most widely used spectral analysis technique is commonly referred to as
_temporal eigenanalysis_ [25, 26, 27, 28, 29, 30], which quantifies how
specific spatial oscillations will evolve over time under periodic conditions
in space. The complementary analysis technique, called _spatial eigenanalysis_
[31, 32, 33, 34, 35], considers instead how wave-like solutions of given
frequency behave while propagating in space under inflow-outflow conditions.
Even if similar and intimately connected from a theoretical perspective, these
two approaches can provide very different insights on the general behaviour of
the numerical scheme. The former is more suitable for temporally evolving
problems, whereas the latter is more appropriate for spatially developing
problems.

Although such techniques are able to reveal useful numerical properties, they
have been often categorised as too simplistic. In fact, a direct connection
between the numerical discretisation of the linear advection equation and more
complex three-dimensional nonlinear simulations can be everything but trivial.
The high-order schemes’ research community has also produced substantial works
on generalised forms of spectral analyses. This includes analyses of linear
advection on multiple dimensions [25, 36, 37], of linear advection problems
with non-constant velocity [38, 39, 35] and even of alternative eigenanalysis
techniques, such as non-modal analysis [40] and combined-mode analysis [41].
Nonetheless, all of these approaches commonly involve a semi-discrete form of
the linear advection equation, where only spatial numerical errors are
considered. The main goal of the present work is focused on a fully-discrete
formulation of spatial eigenanalysis where also temporal errors are taken into
account. Despite the fact that some recent work has been presented regarding
fully-discrete temporal eigenanalysis [42, 43, 44], the full generalisation of
spatial eigenanalysis including a temporal discretisation remains, as of yet,
unexplored. Even if spatial errors play a critical role in the simulation of
turbulent flows, the interaction with the temporal counterparts is everything
but obvious. The temporal discretisation can significantly affect the spatial
errors, leading to much different results in terms of accuracy and stability
of the numerical simulation.

The ultimate goal of simplified spectral analyses is to predict, within a
reasonable range, how a certain scheme, for a given resolution, would behave
in practical flow configurations. For instance, it is possible to estimate
which grid resolution allows one to reliably perform Implicit LES (ILES),
thereby reducing the computational burden otherwise incurred. In the same way,
temporal discretisations would benefit from large time steps, while
introducing negligible numerical errors. In other words, a trade-off between
accuracy and computational cost needs to be found both spatially and
temporally. Finding such delicate balance using fully three-dimensional
simulations can be extremely expensive and even pointless. Ideally, efficient
and reliable simplified spectral analysis techniques aim at providing
precisely this kind of information at very low computational cost.

The paper is organised as follows. In section 2, we introduce the fully-
discrete eigenanalysis framework which serves both temporal and spatial
approaches. This framework is applied to the Flux Reconstruction method, which
can recover both the discontinuous Galerkin and Spectral difference schemes.
In section 3, we present the eigenanalysis results for the fully-discrete
spatial approach, where we focus on the interpretation of the diffusion curves
and the contributions of the different (transmitted/reflected) modes. The
observations will be specifically classified into well-resolved and under-
resolved frequency ranges, providing a direct link with DNS and LES
approaches, respectively. In section 4, we present the numerical experiments
conducted to assess the validity of the eigenanalysis. In particular, first a
discretisation of the linear advection equation will be considered, followed
by more complex inviscid flows based on the Euler equations. In line with the
theoretical discussion, numerical results will be arranged in well-resolved
and under-resolved flows. Finally, in section 5, we outline the key
conclusions of this work.

##  2 Eigenanalysis framework

In this section, we introduce the fully-discrete eigenanalysis framework,
which serves both the temporal (section 2.1) and the spatial (section 2.2)
approaches. The classical fully-discrete temporal eigenanalysis (that was
originally reported by Vermeire and Vincent [43]) is presented first, as a way
of introducing the key building blocks for the fully-discrete spatial
eigenanalysis.

###  2.1 Classical fully-discrete temporal eigenanalysis

All the relevant steps of temporal eigenanalysis will be herein introduced
using the FR scheme, which is well-known for its ability to recover different
methods with a specific choice of its correction functions [45, 46, 47]. The
present study considers two such methods in particular, namely DG and SD. In
fact, this strategy allows one to perform subsequent eigenanalysis studies for
any other method that FR is capable of recovering.

In the standard temporal analysis, the linear advection equation is
discretised looking for wave-like solutions in order to study the numerical
description of specific wavenumbers and their evolution over time. Then,
dispersion and dissipation properties of any scheme follow directly from the
corresponding eigensolutions. Let us now suppose the case of a constant
unitary advection velocity:

$\frac{\partial{u}}{\partial{t}}+\frac{\partial{u}}{\partial{x}}=0.$ (1)

This equation admits plane wave solutions of the form

$u(x,t)=e^{\iota(\theta x-\omega t)},\quad\text{with $\iota=\sqrt{-1}$},$ (2)

provided that the angular frequency $\omega=\omega(\theta)$ is such that

$\textrm{Re}(\omega)=\theta\qquad\textrm{and}\qquad\textrm{Im}(\omega)=0,$ (3)

where $\theta$ is a real-valued wavenumber chosen at the initial condition (in
the case of standard temporal eigenanalysis). Equation (3) provides what are
respectively known as dispersion and diffusion relations for the exact plane
wave solution. A numerical discretisation of the linear advection equation is
now introduced on a grid of equally spaced elements of fixed width
$h_{n}=h=1$. Under such assumptions, a FR spatial discretisation of equation
(1) within the standard element $\Omega_{n}$ can be written as

$\frac{\mathrm{d}{\hat{u}_{n}}}{\mathrm{d}{\textrm{t}}}=-2\sum_{j=0}^{N}\hat{u}_{n}\frac{\mathrm{d}{l_{j}}}{\mathrm{d}{\hat{x}}}(\hat{x}_{i})-(2\hat{f}^{I}_{L}-2\hat{u}_{L})\frac{\mathrm{d}{g_{L}}}{\mathrm{d}{\hat{x}}}(\hat{x}_{i})-(2\hat{f}^{I}_{R}-2\hat{u}_{R})\frac{\mathrm{d}{g_{R}}}{\mathrm{d}{\hat{x}}}(\hat{x}_{i}),$
(4)

where $N$ is the order of the solution polynomial,
$\hat{u}_{L}=\hat{u}_{n}(-1,t)$ and $\hat{u}_{R}=\hat{u}_{n}(+1,t)$.

To examine a sufficiently general framework, the classical set of numerical
fluxes for the linear advection equation is considered:

$\hat{f}_{L}^{I}=(1-\alpha)\hat{u}_{n-1}(+1,t)+\alpha\hat{u}_{n}(-1,t)$ (5)

and

$\hat{f}_{R}^{I}=(1-\alpha)\hat{u}_{n}(+1,t)+\alpha\hat{u}_{n+1}(-1,t).$ (6)

Note that $\alpha=0$ leads to a standard upwind scheme, whereas $\alpha=0.5$
leads to a centred scheme. In the following sections, unless stated otherwise,
standard upwind fluxes will be considered. The $\alpha$ parameter can be used
to mimic different levels of interface upwinding, which can be directly linked
to more complex numerical fluxes for the Euler or Navier-Stokes equations (via
Riemann solvers). Equation (4) can be written in matrix form as:

$\frac{\mathrm{d}{\hat{\textbf{u}}^{n}}}{\mathrm{d}{\textrm{t}}}=-2\textbf{D}\hat{\textbf{u}}^{n}-(2\hat{f}^{I}_{L}-2\textbf{l}^{T}\hat{\textbf{u}}^{n})\textbf{g}_{L}-(2\hat{f}^{I}_{R}-2\textbf{r}^{T}\hat{\textbf{u}}^{n})\textbf{g}_{R},$
(7)

where $\hat{\textbf{u}}^{n}_{i}=\hat{u}(\hat{x}_{i})$,
$\textbf{D}_{ij}=\frac{\mathrm{d}{l_{j}}}{\mathrm{d}{\hat{x}}}(\hat{x}_{i})$,
$\textbf{g}^{L/R}_{i}=\frac{\mathrm{d}{g_{L/R}}}{\mathrm{d}{\hat{x}}}(\hat{x}_{i})$
and $\textbf{r}_{i}=l_{i}(1)$, $\textbf{l}_{i}=l_{i}(-1)$.

In order to reduce the infinite dimensional problem to an element-wise
framework, Bloch wave-like type of solutions are sought, namely:

$\hat{\textbf{u}}^{n}=e^{\iota(\tilde{\theta}x_{n}/h-\tilde{\omega}t)}\hat{\textbf{v}},$
(8)

where tilde accent denotes the numerical counterparts of continuum variables
and $x_{n}$ indicates the left boundary of the element $n$ in the global
reference frame. Once $\hat{\textbf{u}}^{n}$ is properly defined, due to
periodicity, it is straightforward to obtain a closed form for all the fluxes
expressed in equation (4):

$\displaystyle\hat{f}_{L}^{I}$
$\displaystyle=(1-\alpha)e^{-\iota\tilde{\theta}}$
$\displaystyle\textbf{r}^{T}\hat{\textbf{u}}^{n}$
$\displaystyle+\alpha\textbf{l}^{T}$ $\displaystyle\hat{\textbf{u}}^{n},$ (9)
$\displaystyle\hat{f}_{R}^{I}$ $\displaystyle=(1-\alpha)$
$\displaystyle\textbf{r}^{T}\hat{\textbf{u}}^{n}$
$\displaystyle+\alpha\textbf{l}^{T}e^{\iota\tilde{\theta}}$
$\displaystyle\hat{\textbf{u}}^{n}.$

Finally, it can be seen that, due to solution periodicity, the problem has now
been reduced to an element-wise formulation, _i.e._ , the numerical fluxes now
depend only on the local solution $\hat{\textbf{u}}^{n}$. Hence, index $n$ can
be dropped and the final discretised equation can be written as

$\frac{\mathrm{d}{\hat{\textbf{u}}}}{\mathrm{d}{\textrm{t}}}=-2\bigg{(}\textbf{D}\hat{\textbf{u}}+((1-\alpha)e^{-\iota\tilde{\theta}}\textbf{r}^{T}\hat{\textbf{u}}+\alpha\textbf{l}^{T}\hat{\textbf{u}}-\textbf{l}^{T}\hat{\textbf{u}})\textbf{g}_{L}+((1-\alpha)\textbf{r}^{T}\hat{\textbf{u}}+\alpha\textbf{l}^{T}e^{\iota\tilde{\theta}}\hat{\textbf{u}}-\textbf{r}^{T}\hat{\textbf{u}})\textbf{g}_{R}\bigg{)},$

or, in a more compact way:

$\frac{\mathrm{d}{\hat{\textbf{u}}}}{\mathrm{d}{\textrm{t}}}=-2(\textbf{C}^{-}e^{-\iota\tilde{\theta}}+\textbf{C}^{0}+\textbf{C}^{+}e^{\iota\tilde{\theta}})\hat{\textbf{u}},$
(10)

where
$\textbf{C}^{0}=\textbf{D}-(\alpha-1)\textbf{g}_{L}\textbf{l}^{T}-\alpha\textbf{g}_{R}\textbf{r}^{T}$,
$\textbf{C}^{-}=(1-\alpha)\textbf{g}_{L}\textbf{r}^{T}$ and
$\textbf{C}^{+}=\alpha\textbf{g}_{R}\textbf{l}^{T}$. When FR is made to
recover specific methods (such as DG or SD), matrices C above are what will
differentiate them, but the structure of the problem does not change. In
either case, the discretised system is now reduced to the following linear
dynamical system:

$\frac{\mathrm{d}{\hat{\textbf{u}}}}{\mathrm{d}{\textrm{t}}}=\textbf{H}(\tilde{\theta})\hat{\textbf{u}}.$
(11)

The approach herein outlined represents a classical tool in numerical
analysis, particularly used for the theoretical study of novel high-order
schemes. Notice that the present analysis is often called _semi-discrete_
temporal eigenanalysis due to the absence of time discretisation. Indeed,
taking advantage of the analytical form of wave-like solutions, the dynamical
system defined by equation (11) can be reduced to the following eigenvalue
problem:

$-\iota\tilde{\omega}\hat{\textbf{u}}=\textbf{H}(\tilde{\theta})\hat{\textbf{u}}.$
(12)

Alternatively, instead of considering the analytical form of the wave-like
solution, a numerical time integration scheme can be introduced to advance in
time. The numerical solution at the time step $t+1$ can be generally
approximated by the M-stage explicit Runge-Kutta schemes as

$\hat{\textbf{u}}^{t+1}\approx\sum_{m=0}^{M}\beta_{m}\frac{(\textbf{H}(\tilde{\theta})\Delta
t)^{m}}{m!}\hat{\textbf{u}}^{t}=\textbf{T}(\tilde{\theta},\Delta
t)\hat{\textbf{u}}^{t},$ (13)

where $\Delta t$ is the time step and $\boldsymbol{\beta}$ contains the
coefficients of the specific Runge-Kutta scheme.

This formulation is denoted as _fully-discrete_ since both spatial and
temporal discretisations are simultaneously involved. Considering the present
generalised formulation, two main closely linked analyses can be performed.
The first is relatively popular and originates from the classical stability
analysis of numerical discretisation of partial differential equations. In
particular, it is well-known that in equation (13), the combined space-time
discretisation is denoted as _stable_ if the spectral radius of the iteration
matrix T is strictly smaller than unity. Namely, for the specific case of Von
Neumann analysis, the full discretisation is considered linearly stable if,
for a given time step $\Delta t$,

$\forall\,\tilde{\theta}\in\big{[}-\pi(N+1)/h;\quad\pi(N+1)/h\big{]}:\quad\lambda(\textbf{T}(\tilde{\theta},\Delta
t))<1,$ (14)

where $\lambda(\textbf{T})$ denotes the spectral radius of the matrix T.

A complementary analysis, instead, takes into account the Bloch wave-like form
of the numerical solution and generalises the concept of dispersion/diffusion
curves in the fully-discrete analyses. In fact, exploiting the Bloch-like form
of the approximate solution, it can be shown (see [48, 43]) that

$\hat{\textbf{u}}^{t+1}=e^{-\iota\tilde{\omega}\Delta t}\hat{\textbf{u}}^{t}.$
(15)

Consequently, the fully-discrete version of the dispersion/diffusion
eigenanalysis can be written as:

$e^{-\iota\tilde{\omega}\Delta
t}\hat{\textbf{u}}^{t}=\textbf{T}(\tilde{\theta},\Delta
t)\hat{\textbf{u}}^{t}.$ (16)

Notice, in fact, that for $\Delta t\to 0$ we have

T $\displaystyle\rightarrow$ $\displaystyle\textbf{I}+\textbf{H}\Delta
t+O(\Delta t^{2}),$ $\displaystyle e^{-\iota\tilde{\omega}\Delta t}$
$\displaystyle\rightarrow$ $\displaystyle 1-\iota\omega\Delta t+O(\Delta
t^{2})$

and the problem reduces to

$\displaystyle(1-\iota\tilde{\omega}\Delta t+O(\Delta
t^{2}))\hat{\textbf{u}}^{t}$ $\displaystyle=(\textbf{I}+\textbf{H}\Delta
t+O(\Delta t^{2}))\hat{\textbf{u}}^{t}$
$\displaystyle-\iota\tilde{\omega}\Delta t\hat{\textbf{u}}^{t}$
$\displaystyle=\textbf{H}\Delta t\hat{\textbf{u}}^{t}+O(\Delta t^{2})$
$\displaystyle-\iota\tilde{\omega}\hat{\textbf{u}}^{t}$
$\displaystyle=\textbf{H}\hat{\textbf{u}}^{t}$

which coincides with classical semi-discrete eigenvalue problem.

In a similar way with respect to the classical analysis, considering (16), it
is possible to compute the frequencies $\tilde{\omega}$ for any given value of
$\tilde{\theta}$ and $\Delta t$. Therefore, the choice of the time integration
scheme and, in particular, of a stable time step $\Delta t$, will affect the
dispersion and diffusion curves. Such curves can be interpreted in the same
manner as in the semi-discrete analysis and significant deviations from the
semi-discrete dispersion/diffusion curves can crucially affect the numerical
properties of the spatial discretisation.

We finally note that, had we chosen to work in the equations above with
unspecified values of the advection velocity $a$ and the mesh spacing $h$, the
relevant time-step factor appearing in equation (13) would have been
$\overline{\Delta t}=a\Delta t/h$—namely, the Courant–Friedrichs–Lewy (CFL)
number—instead of simply $\Delta t$. Hence, we remark that, although we chose
to work with $\Delta t$ for the sake of simplicity, one can actually have in
mind the CFL number whenever $\Delta t$ values are mentioned in the remainder
of the paper.

###  2.2 Fully-discrete spatial eigenanalysis

In the classical temporal eigenanalysis, for any real-valued wavenumber
$\tilde{\theta}$, the numerical angular frequency $\tilde{\omega}$ can be
obtained through the eigenvalue problem in equation (16). The frequency
$\tilde{\omega}$ will describe the numerical evolution of the specified
wavenumber $\tilde{\theta}$ over time. It is however possible to invert the
problem by imposing a time frequency $\tilde{\omega}$ and subsequently seek
for wavenumbers $\tilde{\theta}$ that solve the inverted problem. Such
approach is commonly referred to as _spatial eigenanalysis_ and it represents
a very useful tool for non-periodic problems since it aims at studying the
dynamics of prescribed time frequencies over space.

Within the framework of semi-discrete eigenanalysis, the equation defining the
relationship between $\tilde{\theta}$ and $\tilde{\omega}$ is the same, namely
equation (12), and can be written in an alternative form as

$\det(\textbf{H}(\tilde{\theta})+\iota\tilde{\omega}\textbf{I})=0.$ (17)

In this case, $\tilde{\omega}$ is a real-valued prescribed scalar and the
nonlinear equation (17) is solved in the unknown $\tilde{\theta}$. This
nonlinear equation can be particularly complex to solve and, except for some
very particular cases, cannot be analytically inverted. Considering the fully-
discrete space-time discretisation, an equivalent equation which defines the
relation between $\tilde{\omega}$ and $\tilde{\theta}$ can be derived from
equation (16), namely

$\det(\textbf{T}(\tilde{\theta},\Delta t)-e^{-\iota\tilde{\omega}\Delta
t}\textbf{I})=0.$ (18)

Note of course that, in the temporal eigenanalysis, the addition of a temporal
discretisation does not significantly complicate the overall problem. In fact,
the matrix T can be easily computed once $\tilde{\theta}$ and $\Delta t$ are
prescribed. Subsequently, the eigenvalue problem can be solved in
$e^{-\iota\tilde{\omega}\Delta t}$, which can be finally inverted to obtain
the exact value of $\tilde{\omega}$. However, in the spatial eigenanalysis,
the specific matrix T influences the general structure of the nonlinear
equation to be solved in $\tilde{\theta}$. In particular, if the simple
explicit Euler integration scheme is considered, then equation (18) becomes:

$\det(\textbf{I}+\Delta t\textbf{H}(\tilde{\theta},\Delta
t)-e^{-\iota\tilde{\omega}\Delta t}\textbf{I})=0.$ (19)

If $z$ is defined as $z=e^{\iota\tilde{\theta}}$, the final nonlinear equation
in $z$ can be written as

$\det((1-e^{-\iota\tilde{\omega}\Delta t})\textbf{I}-2\Delta
t(\textbf{C}^{-}z^{-1}+\textbf{C}^{0}+\textbf{C}^{+}z))=0.$ (20)

Again, notice that for $\Delta t\to 0$, the diagonal term
$1-e^{-\iota\tilde{\omega}\Delta t}$ tends to $\iota\tilde{\omega}\Delta t$
and the fully-discrete formulation goes back to the classical semi-discrete
form:

$\det(\iota\tilde{\omega}\textbf{I}-2(\textbf{C}^{-}z^{-1}+\textbf{C}^{0}+\textbf{C}^{+}z))=0,$
(21)

which is exactly equation (17). On the other hand, for a RK22 scheme, the
formulation becomes much more complex as higher powers of $z$ appear in the
final equation:

$\displaystyle\det\bigg{(}$ $\displaystyle(1-e^{-\iota\tilde{\omega}\Delta
t})\textbf{I}+$ (22) $\displaystyle+\Delta
t^{2}\big{(}\textbf{C}^{-}\textbf{C}^{-}\big{)}z^{-2}+$
$\displaystyle+\big{(}\Delta
t^{2}(\textbf{C}^{-}\textbf{C}^{0}+\textbf{C}^{0}\textbf{C}^{-})-2\Delta
t\textbf{C}^{-}\big{)}z^{-1}+$ $\displaystyle+\big{(}\Delta
t^{2}(\textbf{C}^{-}\textbf{C}^{0}+\textbf{C}^{0}\textbf{C}^{0}+\textbf{C}^{+}\textbf{C}^{-})-2\Delta
t\textbf{C}^{0}\big{)}$ $\displaystyle+\big{(}\Delta
t^{2}(\textbf{C}^{+}\textbf{C}^{0}+\textbf{C}^{0}\textbf{C}^{+})-2\Delta
t\textbf{C}^{+}\big{)}z+$ $\displaystyle+\Delta
t^{2}\big{(}\textbf{C}^{+}\textbf{C}^{+}\big{)}z^{2}\bigg{)}=0.$

This way, the order of the nonlinear equation increases significantly and
additional roots for $z$ are expected to be found (especially for time
discretisation schemes of higher order). The role of multiple wave modes in
spatial eigenanalysis has been discussed in more detail in [34]. As properly
explained in the next section, only one solution will be denoted as the
physical one, whereas the others are considered spurious.

##  3 Theoretical Results

In this section, we show the results associated with the fully-discrete
spatial eigenanalysis just introduced. The theoretical findings of the fully-
discrete temporal eigenanalysis previously reported by Vermeire and Vincent
[43] were correctly reproduced within the present implementation too. These
have been used for validation purposes only and, for the sake of brevity, are
not shown here.

It is important to highlight that, from equation (2), the numerical wave-like
solutions depend on the numerical equivalent of $\theta$, i.e.
$\tilde{\theta}$ (whose imaginary part makes the diffusion curve plots below)
according to

$u\propto\exp[-\mbox{Im}(\tilde{\theta})x]\exp[\iota\mbox{Re}(\tilde{\theta})x-\omega
t]\mbox{ ,}$ (23)

where reflected waves should be considered with $\Delta x<0$, thus being
stable for $\mbox{Im}(\tilde{\theta})<0$.

In the following, we consider high-order FR-based DG/SD methods in space,
coupled with explicit Runge-Kutta schemes of increasing order of accuracy (and
increasing number of stages): RK11, RK22, RK33 and RK45. Hence, we start with
the $1^{\mathrm{st}}$-order accurate, one-stage Runge-Kutta scheme (more
commonly known as the Explicit Euler method) until finally addressing the
well-known $4^{\mathrm{th}}$-order accurate, five-stages Runge-Kutta scheme by
Carpenter and Kennedy [49]. This allows for a detailed exposition of the
fully-discrete eigenanalysis results, whose complexity grows with the number
of stages of the temporal scheme, as shown below.

In fig. 1, we show the dissipation curve for the $5^{\mathrm{th}}$-order FR-SD
standard upwind method when coupled with RK11. This curve is shown for
increasing values of the time step $\Delta t$, namely $\Delta
t=0.01,0.02,\dots,0.1$, with dissipation values becoming more negative as
$\Delta t$ increases. Note that, given a $\Delta t$ value, only one
eigenmode/curve exists in this case, corresponding to a single transmitted
physical wave. This is because, from eq. 20, the characteristic polynomial in
the unknown $z$ has the same order of the semi-discrete analysis (compare with
eq. 21). We recall that, in the semi-discrete analysis [32], only the physical
transmitted mode exists in case of standard upwinding ($\alpha=0$ in eq. 10),
whereas one reflected spurious wave exists in addition if $\alpha\neq 0$. As
for the curves in fig. 1, negative dissipation values indicate in fact anti-
dissipative behaviour (convective instability), since we are dealing with a
physical transmitted wave, cf. eq. 23. Hence, fig. 1 essentially re-affirms
that RK11 (explicit Euler) is only marginally stable for pure advection [50].
Interestingly, for moderately large $\Delta t$, waves regain stability beyond
a certain frequency $\omega$, probably due to increased dissipation from the
spatial discretisation at large frequencies. These features will be tested in
the next section, as they highlight how the fully-discrete eigenanalysis is
useful in revealing non-intuitive characteristics stemming from the coupling
between spatial and temporal discretisations.

Figure 1: Dissipation curve for the $5^{\mathrm{th}}$-order FR-SD standard
upwind scheme coupled with RK11 (explicit Euler) denoting a single transmitted
physical mode, shown for $\Delta t=0.01,0.02,\dots,0.1$ (with dissipation
values becoming more negative as $\Delta t$ increases). The plot on the right
shows a closer look near the origin. These re-affirm that RK11 is convectively
unstable for pure advection. Figure 2: Dissipation curves for the
$2^{\mathrm{nd}}$-order FR-SD standard upwind scheme coupled with RK22 for
increasing values of $\tau=\Delta t/\Delta t_{\max}$. The plots on the right
show the same left curves, but for a larger frequency range. The semi-discrete
($\tau=0$) result also appears when $\tau>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 $\omega^{c}=\pi/\Delta t$.

Now we consider the classical $2^{\mathrm{nd}}$-order explicit Runge-Kutta
(RK22) scheme, this time coupled for simplicity with the
$2^{\mathrm{nd}}$-order FR-SD standard upwind method. The resulting
discretisation is stable for $\Delta t<\Delta t_{\max}$, where $\Delta
t_{\max}$ is the critical time-step size at which instabilities first appear.
In fig. 2, we show the corresponding dissipation curves for increasing values
of $\tau=\Delta t/\Delta t_{\max}$. When $\tau=0$, the results match those of
the semi-discrete analysis, where only the transmitted physical mode exists in
case of standard upwinding. For $\tau>0$, a spurious transmitted mode appears
in addition to the transmitted physical one, not to be confused with the thin
continuous curve which simply shows the semi-discrete ($\tau=0$) physical
curve as reference (note that the physical curve is shown in red, whereas the
spurious one is shown as a thick black curve). The existence of exactly two
modes is guaranteed by the characteristic polynomial in eq. 22 being only
quadratic in $z$, since matrix $\textbf{C}^{+}$ becomes zero in case of
standard upwinding ($\alpha=0$ in eq. 10). This feature can be traced
algebraically to the number of stages of the time discretisation adopted. In
general, the number of modes (in case of standard upwinding) matches the
number of stages of the chosen explicit Runge-Kutta scheme. The physical mode
can be identified as the one _anchored_ to the origin of the dissipation
plots, thus featuring vanishing dissipation levels for well-resolved waves (as
expected in simulations of pure advection).

Returning to fig. 2, we note that, as $\tau$ increases, the physical mode
becomes less dissipative overall (at least for $\omega h(N+1)^{-1}>0.5$, say)
until finally touching the horizontal axis, whose crossing (for $\tau>1$)
means anti-dissipation (convective instability). The spurious mode, despite
being _infinitely dissipated_ for $\tau\rightarrow 0$, also has its
dissipation reduced as $\tau$ increases; in fact, for moderately large $\tau$,
its dissipation levels are comparable to those of the physical mode. Small
dissipation is typically undesirable for spurious modes, since it allows them
to survive longer as artificial oscillations in a simulation, potentially
affecting not only solution quality, but also numerical stability [33].
Curiously, as $\tau\rightarrow 1$, the spurious mode’s dissipation gradually
approaches that of the physical mode, almost matching it at $\tau=1$. This
behaviour closely resembles that of reflected spurious modes which exist for
either very low or very high upwinding levels (i.e. either centred fluxes or
_hyper-upwinding_), cf. [32], except that reflected modes have dissipation
values of opposite sign when compared to those of their corresponding physical
modes. In any case, for the present scheme, both physical and spurious modes
touch the horizontal axis of zero dissipation when $\tau=1$, becoming
therefore unstable for $\Delta t>\Delta t_{\max}$.

Still in fig. 2, note that the spurious mode’s dissipation curve (e.g. when
considered until $\omega^{c}=\pi/\Delta t$, such that $\omega^{c}h(N+1)^{-1}$
is marked as a vertical dashed line) resembles a mirror image of the physical
mode’s curve. What is more generally true, though, is that the spurious curve
is simply a shifted replication of the physical curve, while both are periodic
functions of period $2\omega^{c}$. Interestingly, the periodicity of
eigencurves in spatial analysis has not been noted in the literature so far.
This is because $\omega^{c}\rightarrow\infty$ when $\tau\rightarrow 0$, cf.
fig. 2, consistent with the fact that the spurious mode becomes _infinitely
dissipated_ for $\tau\rightarrow 0$, as already mentioned. Therefore,
periodicity is only visible in fully-discrete spatial analysis. This is in
contrast with temporal analysis, where one sees periodicity already in the
semi-discrete case, as well as the shifted replication property amongst
multiple eigencurves [29, 51].

Although the cut-off value $\omega^{c}=\pi/\Delta t$ can be rightfully
recognized as the Nyquist frequency imposed by the discrete time-stepping even
outside of the fully-discrete spatial framework [43], it is only in the fully-
discrete spatial analysis that $\omega^{c}$ exhibits visual correlation with
the periodicity of the eigencurves, as already noticed in fig. 2. In fact, the
normalized version of $\omega^{c}$, namely $\omega^{c}h(N+1)^{-1}$, is the
counterpart of the wavenumber limit $\theta h(N+1)^{-1}=\pi$ typically used in
temporal analysis. Hence, it seems appropriate to plot spatial analysis curves
only until $\omega^{c}$. Nevertheless, as it will become clear in subsequent
plots, it so happens that the actual limit observed for explicit Runge-Kutta
schemes of $R$ stages is $R\,\omega^{c}/2$. Interestingly, the conclusion is
that multiple stages seem to affect the temporal resolution as if multiple
single-stage time steps of size $\delta t=2\Delta t/R$ had been employed
instead. This seems akin to how multiple element-wise polynomial modes
increase the spatial resolution of a single element, as if it had _sub-
elements_.

A first obvious consequence of adopting $R\,\omega^{c}/2$ as the right-most
value of $\omega$ when plotting eigencurves is that this bound tends to
infinity for $\tau\rightarrow 0$. Indeed, in the semi-discrete spatial
analysis, every value of $\omega$ is allowed, in principle. Secondly, such
bound gets closer and closer to zero for larger time steps, as only lower and
lower frequencies can be solved by the temporal discretisation when the time
step increases. These effects are highlighted in fig. 3, where only the
physical mode is represented for different values of $\tau$.

Figure 3: Dissipation curves of the physical transmitted mode at increasing
values of $\tau$ (from $\tau=0$ in black until $\tau=1$ in yellow). Results
obtained for the $2^{\mathrm{nd}}$-order FR-SD standard upwind discretisation
coupled with RK22.

We now consider the $5^{\mathrm{th}}$-order FR-SD standard upwind scheme
coupled with either RK33 (fig. 4) or RK45 (fig. 5). Dissipation curves for
$\tau=0.5$ and $\tau=1$ are respectively shown on the left and right plots
within these figures. As previously mentioned, the number of roots in the
respective characteristic polynomials matches the number of stages in the
adopted time scheme. This is why RK33 (fig. 4) yields three transmitted modes,
whereas RK45 (fig. 5) yields five. The correspondence between the number of
eigenmodes and stages in spatial analysis resembles that between eigenmodes
and element-wise polynomial modes in temporal analysis (in line with the
previous comment regarding the interpretation of stages as sub-divisions of
the actual time step $\Delta t$). In figs. 4 and 5, the black and red vertical
dashed lines indicate, respectively, $\omega^{c}$ and $R\,\omega^{c}/2$, the
latter matching half the period of the physical mode, as expected. Note that
the physical mode is found to be symmetric with respect to the red dashed line
(actual Nyquist frequency). Moreover, all eigenmodes are observed to be
shifted replications of the physical mode, similar to temporal analysis. This
is sometimes useful for the identification of the physical mode.

The physical curves of figs. 4 and 5 have been isolated and are shown in fig.
6 for different values of $\tau$. For comparison, the equivalent plots for FR-
DG are given in fig. 7. As one can see, the general trends observed for FR-SD
are also seen for FR-DG. Some differences, however, can be discerned. For
example, since the CFL limit is more restrictive for FR-DG than for FR-SD
(both under standard upwinding), the dissipative curves of FR-DG extend to
higher frequencies (recall that $\omega^{c}=\pi/\Delta t$). Moreover, as
already pointed out in previous works [31, 32, 33, 35], FR-DG is spatially
less dissipative than FR-SD, overall. A more complete comparison between FR-SD
and FR-DG is shown in the Appendix (cf. figs. 28 and 29).

Figure 4: Dissipation curves for the $5^{\mathrm{th}}$-order standard upwind
FR-SD coupled with RK33 at $\tau=0.5$ (left) and $\tau=1$ (right). The
physical mode is shown as a continuous red curve. The vertical black dashed
line marks $\omega^{c}$, whereas the vertical red dashed line marks
$R\,\omega^{c}/2$ (indicating half a period).

Figure 5: Dissipation curves for the $5^{\mathrm{th}}$-order standard upwind
FR-SD coupled with RK45 at $\tau=0.5$ (left) and $\tau=1$ (right). The
physical mode is shown as a continuous red curve. The vertical black dashed
line marks $\omega^{c}$, whereas the vertical red dashed line marks
$R\,\omega^{c}/2$ (indicating half a period).

Figure 6: Dissipation curves of the physical transmitted mode at increasing
values of $\tau$ (from $\tau=0$ in black until $\tau=1$ in yellow). Results
obtained for the $5^{\mathrm{th}}$-order FR-SD standard upwind discretisation
coupled with either RK33 (left) or RK45 (right).

Figure 7: Dissipation curves of the physical transmitted mode at increasing
values of $\tau$ (from $\tau=0$ in black until $\tau=1$ in yellow). Results
obtained for the $5^{\mathrm{th}}$-order FR-DG standard upwind discretisation
coupled with either RK33 (left) or RK45 (right).

To quantify the mutual influence between time and spatial discretisations, the
approximation error can be evaluated separately for both temporal and spatial
approaches, since the exact solution is known. In particular, for the temporal
approach, the error $E_{T}$ will read

$E_{T}=|\tilde{\omega}(\theta)-\omega(\theta)|,$ (24)

whereas, in the spatial analysis the error $E_{S}$ can be evaluated as

$E_{S}=|\tilde{\theta}(\omega)-\theta(\omega)|.$ (25)

Notice, of course, that the estimates above take into account both dispersive
and diffusive errors. Both errors have been evaluated for the
$5^{\mathrm{th}}$-order standard upwind FR-SD method coupled with the RK33
scheme. Results are shown in fig. 8. Similar to what has been found in [43],
for small values of $\theta$, the accuracy of the scheme is reduced by the
time integration method. In particular, the expected super-convergence of the
FR-SD scheme (of order $2N+1=9$) reduces to a $4^{\mathrm{th}}$-order
convergence. For a certain intermediate range of wavenumbers, instead, the
theoretical accuracy of the scheme is preserved. The smaller the prescribed
time step, the longer such range will penetrate into the small wavenumbers
region. The same applies in the spatial analysis: for small frequencies, the
scheme converges with order $4$ and for intermediate frequencies with order
$9$. The differences between the two plots are restricted to the high
frequencies/wavenumbers region. In fact, since the imaginary parts of both
$\tilde{\omega}$ and $\tilde{\theta}$ are negligible for small
frequencies/wavenumbers, the temporal and spatial approaches are equivalent,
whereby $\tilde{\omega}$ and $\tilde{\theta}$ simply switch axes.

Figure 8: Approximation errors $E_{T}$ and $E_{S}$ obtained, respectively,
from the temporal (left) and spatial (right) frameworks for a
$5^{\mathrm{th}}$-order FR-SD method coupled with the RK33 scheme.

As already discussed in previous works on spatial eigenanalysis [31, 32],
whenever the numerical interface flux departs from the standard upwind
condition, twice as many roots will exist in the relevant characteristic
polynomial. The new roots are typically found to represent wave-like modes
that propagate in the opposite direction (reflected modes) and are
consequently considered spurious. As already shown for upwind fluxes, for any
given value of $\omega$, the number of roots depends on the number of stages
of the time integration scheme. In the same way, for non-upwind schemes, the
total number of roots is expected to be twice the number of stages. In the
Appendix (cf. figs. 30 and 31), plots showing all such eigencurves are given.
Here, to avoid unnecessary confusion, only the physical modes are shown. The
case considered is that of nearly centered flux ($\alpha=0.49$), whose
dissipation curves are shown in fig. 9, again for the $5^{\mathrm{th}}$-order
FR-SD method coupled with either RK33 or RK45, as previously shown for the
standard upwind case in fig. 6.

Figure 9: Dissipation curves of the physical transmitted mode at increasing
values of $\tau$ (from $\tau=0$ in black until $\tau=1$ in yellow). Results
obtained for the $5^{\mathrm{th}}$-order FR-SD with nearly centered flux
($\alpha=0.49$) coupled with either RK33 (left) or RK45 (right).

Similarly to the standard upwind case, for increasing values of $\tau$,
numerical dissipation decreases in the high frequency region. For sufficiently
high values of $\tau$, the curve for $\mathrm{Im}(\tilde{\theta})$ is expected
to _touch_ the horizontal axis with the potential consequent onset of anti-
dissipative behaviours, leading to convective instability. This is clearly
visible in fig. 9 for the RK33 scheme (left). Less so for the RK45 scheme, for
which a similar behaviour is expected for high frequencies which are well
beyond the plotted range. For almost centered fluxes, in the spatial analysis
framework, the presence of _dissipative bubbles_ has been already observed in
previous works [32, 35]. In other words, significant dissipation appears
concentrated in certain frequency ranges, whereas other ranges can be almost
completely unaffected. In fig. 9, it can be noticed that, for increasing
values of $\tau$, the separation of scales is less strong and the dissipative
curves are characterised by a more monotonic behaviour. Non-monotonic
dissipation profiles in frequency space might potentially lead to undesirable
non-smooth features in simulations with broadband energy content. It is
expected that having neighbouring frequencies so differently affected can
deteriorate the distribution of dissipation among scales and impact accuracy
and stability of the numerical simulation.

Consequently, only well-resolved wavenumbers/frequencies behave quite
similarly in temporal/spatial approaches, whereas large differences are
observed in the high wavenumbers/frequencies region. For practical turbulent
simulations using high-order discontinuous spectral element methods, depending
on the local level of resolution, different ranges of scales and frequencies
are of particular interest. For highly resolved flows, such as in DNS, the
correlation with spectral analysis will be mainly restricted to the well-
resolved regions (i.e. small wavenumbers/frequencies). Instead, reducing the
resolution of the scheme and passing from DNS to LES of turbulent flows,
higher wavenumber and frequency regions in the plots will provide more useful
information.

##  4 Numerical experiments

In order to validate the theoretical framework introduced in previous
sections, a series of numerical experiments is presented to quantify the
interplay between temporal and spatial errors for practical flows. First, we
consider the discretisation of the one-dimensional linear advection using the
SD scheme. We then take into account more complex numerical simulations of the
Euler equations in order to evaluate the influence of the time step on both
well-resolved and under-resolved flows. In all the Euler computations
presented herein, if not stated differently, the original Roe flux has been
employed as a representative of the standard upwind flux. All the simulations
presented in the following subsections have been performed using the improved
SD method by Liang and Jameson [1, 2, 3].

###  4.1 One-dimensional linear advection equation

The first numerical test consists in the discretisation of the one-dimensional
linear advection equation with an oscillating inlet boundary condition. This
simulation allows to monitor the general behaviour of dissipative bubbles when
the prescribed time step approaches the CFL limit. Tests are conducted using
almost centered numerical fluxes and RK33 or RK45 time integration schemes for
different injected frequencies. Accordingly, fig. 10 shows the same curves
presented in fig. 9 in which the tested frequencies are indicated in order to
appreciate the expected level of dissipation of each test. The relevant
equivalent frequencies are $\omega_{1}=1$, $\omega_{2}=2.75$ and
$\omega_{3}=3.5$ (only for the RK33 test), which are highlighted in fig. 10
using dashed vertical lines. Close to $\omega_{2}$ and $\omega_{3}$, the
expected levels of numerical dissipation are highlighted using blue dots.

The domain $\Omega=[0,1]$ is split into $100$ equally spaced elements
($h=0.01$). The frequency of the inlet signal can then be directly linked to
the equivalent frequency $\omega h/(N+1)$. The simulation matches the
theoretical analysis: a $5^{\mathrm{th}}$ order SD discretisation, coupled
with a RK33/RK45 time integration schemes has been employed. To see the
influence of the temporal errors, increasing values of $\tau$ have been
considered.

(a) RK33

(b) RK45

Figure 10: Physical dissipative curves varying $\tau$ (FR-SD) using almost
centered numerical fluxes. Dashed, vertical lines denote the frequencies
$\omega_{1}$, $\omega_{2}$ and $\omega_{3}$ (only for RK33). For $\omega_{2}$
and $\omega_{3}$ the dots indicate the levels of dissipation expected for the
specific values of $\tau$ used in the simulations.

First, in fig. 11, the velocity signal $u(x)$ is plotted along the $x$
direction for $\omega=\omega_{1}$ and for different values of $\tau$ using a
RK33 time marching scheme. According to the theoretical analysis presented in
the first part of the paper, for increasing values of $\tau$, larger values of
dissipation are expected (see fig. 10). Such conclusion is confirmed by the
numerical simulations shown in fig. 11. The amplitude of the signal $u(x)$
tends, in fact, to decrease for larger time steps.

(a) $\tau=0.1$

(b) $\tau=0.5$

(c) $\tau=0.9$

Figure 11: Solution signal $u(x)$ at prescribed frequency $\omega_{1}$ for
different values of $\tau$ (RK33).

The same experiment has been repeated using a RK45 scheme and the results are
shown in fig. 12. The same conclusions apply: increasing the time step
magnitude leads to an increase of numerical dissipation. Notice, furthermore,
that for $\tau=0.9$ the signal is more dissipated using the RK33 scheme rather
than the RK45 one. In fact, from fig. 10, it can be noticed that the
dissipative curves, in proximity of $\omega_{1}$, tend to increase a bit more
using the RK33 scheme rather than the RK45 one.

(a) $\tau=0.1$

(b) $\tau=0.5$

(c) $\tau=0.9$

Figure 12: Solution signal $u(x)$ at prescribed frequency $\omega_{1}$ for
different values of $\tau$ (RK45).

The same comparison has been performed for the frequency $\omega_{2}$, which
is located in between dissipative bubbles. More or less the same behaviour is
expected, where larger values of $\tau$ lead to stronger numerical
dissipation. The solution $u(x)$ is shown for the RK33 scheme in fig. 13 and
for the RK45 scheme in fig. 14. The only interesting difference is that, in
agreement with the theoretical analysis, for $\tau=0.9$ the RK33 scheme is
less dissipative than the RK45. For $\tau\approx 0$ and $\tau=0.5$, instead,
the levels of numerical dissipation are quite similar.

(a) $\tau=0.1$

(b) $\tau=0.5$

(c) $\tau=0.9$

Figure 13: Solution signal $u(x)$ at prescribed frequency $\omega_{2}$ for
different values of $\tau$ (RK33).

(a) $\tau=0.1$

(b) $\tau=0.5$

(c) $\tau=0.9$

Figure 14: Solution signal $u(x)$ at prescribed frequency $\omega_{2}$ for
different values of $\tau$ (RK45).

In the RK33 case, the additional frequency $\omega_{3}$ has been considered.
At this frequency, in fact, a different behaviour with respect to the low
frequencies is expected: for increasing values of $\tau$, numerical
dissipation is expected to decrease instead of increase like for all the other
frequencies. Such behaviour is confirmed in fig. 15 where for $\tau$ very
close to the CFL limit the numerical dissipation is significantly lower. In
similarity with fully-discrete temporal eigenanalysis, the highest frequencies
are the first to experience instability. Then, in a general sense, the present
test case suggests that for a reasonably large, low frequency region, the
influence of the time integration scheme is essentially dissipative and larger
time steps cause stronger dissipation. In the high-frequency region, instead,
the behaviour is the opposite: in the semi-discrete analysis such frequencies
are usually dumped very fast, whereas, approaching the CFL limit, numerical
dissipation decrease significantly, until, eventually, it changes sign leading
the simulation to be unstable. Consequently, in similarity with respect the
fully-discrete temporal eigenanalysis, the largest frequencies are responsible
for the CFL instability.

(a) $\tau=0.1$

(b) $\tau=0.5$

(c) $\tau=0.9$

Figure 15: Solution signal $u(x)$ at prescribed frequency $\omega_{3}$ for
different values of $\tau$ (RK33).

Finally it is interesting to study the behaviour of a scheme that, according
to the fully-discrete spatial eigenanalysis herein presented, is expected to
be unstable for any choice of time-step. The same computational grid and order
of approximation have been chosen. In similarity with previous tests, a series
of frequencies have been simulated. In fig. 16 the dissipative curves for a
$5^{\mathrm{th}}$ order FR-SD discretisation coupled with an explicit Euler
time-integration is shown. In particular, in order to quantify the instability
of such space/time discretisation, the following equivalent frequencies have
been considered: $\omega_{0}=0.1$ and $\omega_{1}=1.3$. The low frequency test
is meant to investigate the instability of the prescribed discretisation for
any choice of computational time step. We expect that for any choice of time-
step the prescribed sinusoidal signal will grow while propagating in the
domain. For very low frequencies the increase will likely be very small, but
nonetheless still present.

Figure 16: Physical dissipative curves varying $\tau$ (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.

For the second frequency, instead, a different behaviour is expected to take
place: for sufficiently small time steps the numerical scheme should dissipate
the inlet oscillations, but, for increasing time steps, numerical dissipation
should decrease, until it becomes anti-dissipation, leading the numerical
solution to grow while propagating in the domain. The sudden change of
behaviour, from dissipative to anti-dissipative is predicted to take place at
$\Delta t\approx 0.01$ from the theoretical analysis. The intersection between
the dissipative curve for $\Delta t\approx 0.01$ and the $x$-axis is
highlighted in fig. 16. Consequently, after proper rescaling by the grid size
of the present simulation ($h=0.01$), we expect to observe such change of
trend for time steps around $\Delta t=1.0\times 10^{-4}$. Three time-step
values are therefore tested near $\Delta t=1.0\times 10^{-4}$, as shown in
fig. 17.

Figure 17: Numerical solution of $\omega_{0}=0.1$ using explicit Euler time
integration scheme. Top, $\Delta t=9.825\times 10^{-5}$; Center $\Delta
t=9.925\times 10^{-5}$; Bottom $\Delta t=1.0925\times 10^{-4}$.

The numerical solution behaves as expected, where a sudden change from
dissipative to anti-dissipative character is observed for increasing time
steps. For $\Delta t=9.925\times 10^{-5}$, the prescribed signal at the inlet
is barely affected by numerical dissipation. For smaller time steps, numerical
dissipation increases and the sinusoidal signal is expected to decrease in
amplitude as numerically confirmed in the top plot of fig. 17. Finally, for
time steps slightly larger, instead, negative values of the imaginary part of
$\tilde{\theta}$ arise and the solution is expected to grow over space, as
confirmed in the bottom plot of fig. 17. Numerical results are then in very
good agreement with the theoretical analysis herein presented which is able to
properly predict the effect of numerical dissipation in a fully-discrete
framework where both spatial and temporal errors contribute on the overall
accuracy and regularity of the numerical solution.

Regarding the simulation of the lowest frequency, in order to allow a
sufficient large number of periods within the unitary domain, the number of
elements has been increased to $1000$. In fig. 18 the results of the simple
advection of the sinusoidal signal with lowest frequency is shown for the
explicit Euler scheme with a considerably smaller time step, namely $\Delta
t=1.0\times 10^{-6}=5\times 10^{-3}h/(N+1)$. On the top, the solution is
probed at $x=0.5$, on the bottom, instead, the numerical solution at the end
of the simulation is shown. As predicted by the theoretical 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 line $u=1$.

Figure 18: Numerical solution of $\omega_{1}=1.3$ using explicit Euler time
integration scheme with $\Delta t=1.0\times 10^{-6}$. Top, probed solution at
$x=0.5$; Bottom, numerical solution at the end of the simulation. The red
horizontal line indicates the iso-line $u=1$. Figure 19: Closer look of the
solution of fig. 18 around $u=1$ (red line).

The same simulation has been performed using an even smaller time step $\Delta
t=1.0\times 10^{-7}=5\times 10^{-4}h/(N+1)$. The probed and final solution of
the numerical simulation is shown in fig. 20. Even if barely perceptible, the
solution tends to grow while advected along $x$. The same detailed view of
fig. 19 is reproduced in fig. 21. In this zoomed visualisation, it is possible
to observe more clearly the increase in amplitude of the sinusoidal signal
while propagating into the domain, whereas no growth is observed monitoring
the probed solution.

Figure 20: Numerical solution of $\omega_{1}=1.2$ using explicit Euler time
integration scheme with $\Delta t=1.0\times 10^{-7}$. Top, probed solution at
$x=0.5$; Bottom, numerical solution at the end of the simulation. The red
horizontal line indicates the iso-line $u=1$. Figure 21: Closer look of the
solution of fig. 20 around $u=1$ (red line).

The above results confirm that the explicit Euler scheme is marginally stable
for pure advection, a well-known deficiency of this method [50]. Still, for
under-resolved problems (i.e. large frequencies), the fully-discrete scheme
revealed to be stable, probably due to the increased upwind dissipation acting
against the anti-dissipation of the explicit Euler method. In either case,
once the solution establishes itself in space, it does not grow over time — in
line with what is expected within the spatial analysis framework: either
stability or convective instability.

###  4.2 Accurately resolved flows

In the theoretical analysis, presented in the first part of the paper, the
discussion was organized making distinction between high and low frequencies,
considering more or less extended upper bounds for the frequency $\omega$. In
the low frequency regions, the reduced order of accuracy for time steps
approaching the CFL limit was observed. The convergence results already
reported in the temporal eigenanalysis by Vermeire et al. [43] were correctly
reproduced within the present fully-discrete spatial eigenanalysis framework
(see fig. 8). Such feature highlights the similarity between temporal and
spatial eigenanalysis for well-resolved wavenumbers/frequencies respectively.
For high frequencies, instead, the fully-discrete spatial eigenanalysis has
shown significant variations from the semi-discrete analysis. In the present
subsection we will perform a series of numerical experiments on accurately
resolved nonlinear inviscid flows in order to validate the theoretical
analysis for well-resolved frequencies. A second subsection will instead be
dedicated to under-resolved flows which are commonly associated to high-
frequencies which are only partially resolved by the numerical scheme. In
particular, as a representative example of the spatial eigenanalysis, a
variation of the classic isentropic vortex simulation [52] have been
considered. Such formulation is based on an inflow-type condition and is
characterised by spatial grid inhomogeneity due to a sudden change of mesh
resolution in the middle of the domain. Clearly, such configuration is
particularly representative of spatial eigenanalysis, where an oscillatory
flow field is imposed at the inlet boundary and let evolve along the
streamwise direction. Finally, in terms of turbulent flows, the study of
highly accurate fields, where all the scales and frequencies are sufficiently
well-resolved by the numerical scheme, is directly linked to the Direct
Numerical Simulation (DNS) approach.

####  4.2.1 Inlet-outlet advection of isentropic Euler vortex

A similar numerical set-up with respect to the classical isentropic vortex
test case has been considered in order to avoid the use periodic boundary
conditions, which are classically associated to temporal eigenanalysis. A
formulation considering inlet/outlet boundaries is, instead, a much more
representative example of spatial eigenanalysis, where a temporal frequency is
imposed to the system at the inlet and the evolution of the numerical solution
is studied while propagating through the space.

At the initial time, two vortices are centered at $x=1/4L$ and $x=3/4L$, where
$L=40$ denotes the full horizontal length of the domain. Each vortex is
characterised by the smooth solution:

$\displaystyle\rho$
$\displaystyle=\bigg{(}1-\frac{\beta^{2}(\gamma-1)}{16\gamma\pi^{2}}e^{2(1-r^{2})}\bigg{)}^{\frac{1}{\gamma-1}},$
(26) $\displaystyle u$
$\displaystyle=u_{0}-\frac{\beta(y-y_{c})}{2\pi}e^{1-r^{2}},$ $\displaystyle
v$ $\displaystyle=v_{0}+\frac{\beta(x-x_{c})}{2\pi}e^{1-r^{2}},$
$\displaystyle p$ $\displaystyle=\rho^{\gamma},$

where $r=\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}}$, with $x_{c}$ and $y_{c}$
indicating the coordinates of the vortex’s center, and $\gamma=1.4$ indicates
the specific heat ratio. The parameter $\beta$ defines the strength of the
vortex and, for this particular study, a value equal to $5$ has been chosen,
providing a relatively strong vortex. Subsequently, the solution is advected
by the constant velocity field $(u_{0},v_{0})=(5,0)$ while a constant stream
of new isentropic vortices is imposed on the left boundary. After a full
period has passed, two new vortices have entered the domain from the inlet
boundary, whereas the two initial ones have left the domain through the outlet
boundary.

The smooth analytical solution of the isentropic vortex test case provides an
ideal set-up to evaluate the convergence of numerical schemes for the
discretisation of Euler equations. The numerical solution of Euler equations
has been considered on a series of increasingly refined grids. Namely, the
number of elements of the different meshes reads $N_{\mathrm{el}}=N_{x}\times
N_{y}=2N_{y}\times N_{y}$ with $N_{y}=10,20,40,80,160$. In order to highlight
numerical errors, both regular and deformed grids have been considered. An
example of deformed mesh is shown in fig. 22. Furthermore, in order to study
the evolution of the initial solution while propagating in space, a sudden
change on streamwise grid resolution has been placed in the middle of the
domain (see fig. 22). This way, 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 errors associated to the spatial discretisation,
relatively high orders of approximation have been considered. In particular
$6^{\mathrm{th}}$ and $8^{\mathrm{th}}$ order accurate Spectral Difference
simulations have been performed, coupled with a $3^{\mathrm{rd}}$ order Runge-
Kutta scheme as time integration scheme. In order to evaluate the influence of
temporal numerical errors, the time step has been chosen such that
$\tau=0.1,0.5,0.9$.

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

The goal of the present test case consists in showing the essential
equivalence between spatial and temporal approaches in the well-resolved
regions (the former for well-resolved wavenumbers, the latter for well-
resolved frequencies). When we prescribe the stream of vortices through the
inlet condition, given the large number of degrees of freedom (DOF) employed
in the tests (high polynomial order), we are in fact working in the low-
frequency, well-resolved region of spatial analysis. The $L^{2}$ error of the
density field is easily evaluated and a convergence study can be performed
similarly to the classical test case (whose results are omitted for the sake
of brevity).

The convergence plots for the $6^{\mathrm{th}}$ and $8^{\mathrm{th}}$ order
simulations in 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. In
particular, the expected convergence rate arising from the spatial
discretisation is eventually reduced to the order of the time integration
scheme. When meshes of stronger deformation are considered, the same behaviour
is observed, but at smaller grid sizes. This is believed to happen because of
the augmented spatial numerical errors, which can hide the effect of large
time steps.

(a) N=5

(b) N=7

Figure 23: Density’s convergence when using RK33. Black line, $\tau=0.1$;
purple line $\tau=0.5$; orange line $\tau=0.9$. Solid line, structured grid;
dashed line, deformed grid.

In order to illustrate the influence of the time step on the evolution of the
vortex, the absolute value of the thermodynamic entropy
$s=\log(\rho^{-\gamma}p)$ has been employed as an indicator of numerical
errors. Of course, since an inviscid, shock-free, Euler flow is being
considered, the entropy should be theoretically conserved everywhere in the
domain. Any deviation from a constant entropy state is then entirely
associated to numerical artifacts. The entropy fields at the final time using
$\tau=0.9$ and $\tau=0.99$ are shown in fig. 24, for the $6^{\mathrm{th}}$
order discretisation coupled with RK33. It is interesting to notice that in
most of the domain the entropy fields of the two simulations are quite
similar. The only relevant differences are represented by the highly irregular
spurious oscillations close to the line $y=0$. These are stronger in the
second part of the domain, as they are connected with spurious transmitted
waves (only barely dissipated as $\tau\rightarrow 1$). Smaller spurious waves
are reflected into the first part of the domain too, but are more quickly
dissipated as they propagate backwards. In agreement with the work by Vermeire
et al. [43], a non-negligible influence of the time step choice is observed
only very close to the CFL limit.

Figure 24: Absolute value of the thermodynamic entropy
$s=\ln(\rho^{-\gamma}p)$ throughout the domain for the $6^{\mathrm{th}}$ order
simulation coupled with RK33 using $\tau=0.9$ (top) and $\tau=0.99$ (bottom).
Vertical white line denotes the location of the mesh coarsening.

###  4.3 Under-resolved flows

As mentioned in the previous section, a first comparative study between
spatial and temporal eigenanalysis can be performed considering well-resolved
flow fields. In such conditions, both theories should return the expected
decrease of spatial accuracy for large time steps close to the CFL limit. In
the previous section such behaviour has been numerically confirmed. In this
second test, instead, high frequencies and wavenumbers (i.e. under-resolved
flow fields) will be considered. In contrast to the previous section, where
well-resolved fields have been naturally associated to DNS, now the behaviour
of the numerical scheme represents Large Eddy Simulations (LES, and in
particular Implicit LES), where only the large-scale dynamics is adequately
resolved by the number of DOF employed. According to the theoretical analysis,
for cases employing standard upwind fluxes, large time step values should
considerably decrease the overall level of numerical dissipation, affecting
the simulation both in terms of stability and accuracy. To perform these tests
while focusing on spatial analysis errors, a two-dimensional duct flow has
been considered.

####  4.3.1 Two-dimensional grid turbulence

As a typical sample of non-periodic test case, a simple duct flow with
oscillating inlet velocity has been considered. The domain consists in a
$20\pi\times 2\pi$ rectangle. Although two-dimensional and inviscid, the
present test case is sufficiently insightful to quantify the relevant role of
the numerical scheme for under-resolved flows, as already shown in many
different works [32, 33, 34, 35].

The inlet boundary condition reads

$\displaystyle\rho$ $\displaystyle=\rho_{\infty},$
$\displaystyle\boldsymbol{\mathrm{u}}$
$\displaystyle=u_{\infty}\begin{pmatrix}1+A\sin(Ky)\sin(\Omega t)\\\ 0\\\
\end{pmatrix},$ $\displaystyle p$ $\displaystyle=p_{\infty},$

where $\rho_{\infty}=1$, $u_{\infty}=1$ and
$p_{\infty}=(\rho_{\infty}u_{\infty}^{2})/(\gamma\mathrm{Ma}^{2})$ in order to
get a desired value of inflow Mach number. The parameters defining the inflow
perturbations have been set as $A=1/2$, $K=5$, $\Omega=1$. The prescribed Mach
number is set equal to $0.03$. For this particular choice, the flow is then
essentially incompressible. Non-reflective Roe flux-based boundary conditions
are employed at the outlet, whereas top and bottom boundaries are modeled as
free-slip walls. The results have been obtained with a $6^{\rm{th}}$-order SD
scheme coupled with RK33.

In order to trigger spurious numerical oscillations, a sudden change in
streamwise mesh spacing is imposed at $x=12\pi$, as shown in fig. 25. The
number of elements along the cross-flow direction is $N_{\rm el}^{y}=12$.
Along the streamwise direction, this number is $N^{x_{1}}_{\rm el}=72$ in the
first block and $N^{x_{2}}_{\rm el}=12$ in the second block.

Figure 25: Example of computational grid used in the duct flow simulation.

A typical flow field is illustrated in fig. 26, where the two-dimensional
nature of the flow leads to a well-known reverse cascade process, with small
vortices eventually merging into larger ones. Regardless, the mechanism by
which kinetic energy is exchanged across scales is significantly influenced by
the underlying numerical errors, whereby the present under-resolved test case
is insightful with regards to implicit LES.

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

Although Roe-type solvers have been found to successfully suppress spurious
oscillations coming from the outflow boundary and from the line of mesh size
change [32, 33, 35], these results typically rely on simulations that are
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 theory) or not.
The test case has been considered for $\tau=0.5$, $0.9$ and $0.99$.

Similar to the previous test case (isentropic vortex), the time step influence
is very small from a _macroscopic_ point of view. All cases computed with
$\tau<0.95$ (approximately) behaved quite similarly. Only in very close
proximity of the CFL limit, spurious oscillations arose from the numerical
solution. As in the previous case, such behaviour is qualitatively represented
in fig. 27 using the thermodynamic entropy field. Again, spurious transmitted
waves are seen when $\tau$ 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 for this test problem).

Figure 27: Entropy field. The white vertical line denotes the change of mesh
size interface.

Only for values of $\tau$ close to unity, spurious entropy oscillations take
place. Due to the structured nature of the mesh, spurious oscillations arise
diffusively in the whole domain. In fact, a large number of cells reaches
critical values of $\Delta t$ simultaneously, producing spurious entropy
oscillations. The same conclusions obtained from isentropic vortex case seem
to apply here: the spatial under-resolution hides/overshadows the influence of
the prescribed time step. The influence of the time integration scheme is then
noticeable only very close to instability, whereas almost any time step
fulfilling the CFL condition does not impact the overall accuracy of the
numerical solution (already limited by the spatial discretisation).

##  5 Conclusions

A fully-discrete generalisation of the spatial eigenanalysis has been
presented and applied to the Spectral Difference (SD) and Discontinuous
Galerkin (DG) schemes. A detailed discussion has been presented on spatial
dissipative curves and their dependence on the time step size for 1st- to 4th-
order explicit Runge-Kutta time integration schemes. The behaviour for both
well-resolved and under-resolved frequencies has been considered. In the
fully-discrete spatial eigenanalysis, it has been shown that in the low
frequency range, in similarity with the fully-discrete temporal eigenanalysis,
the theoretical spatial convergence rate degrades to the local order of
accuracy of the time integration scheme. In the high wavenumbers/frequencies
regions, instead, temporal and spatial eigenanalyses have significantly
different behaviours. Both are influenced by the time step size and in both
cases numerical dissipation generally decreases for larger time steps
employing upwind fluxes. In the case of almost centered fluxes, the
dissipative bubbles, typical of spatial eigenanalysis of high-order methods,
tend to be smoother for increasing time step size. Regions that experience
negligible numerical dissipation in the semi-discrete analysis are
characterised by stronger diffusion in the fully-discrete generalisation.

The SD discretisation of the one-dimensional linear advection equation has
been used to validate the theoretical findings. Numerical experiments have
shown good agreement with the theory where the time step size has shown to
enhance numerical dissipation for certain frequencies and suppressed it for
others. Nonlinear simulations using Euler equations have been performed to
quantify the interplay between spatial and temporal errors for more complex
dynamics. In agreement with the previous discussion, the Euler simulations
have been grouped in well-resolved and under-resolved flows. An inlet-outlet
formulation of the classical isentropic vortex test case have been used to
validate fully-discrete temporal and spatial eigenanalyses for well-resolved
flows. As predicted by the theory, the expected convergence rate, for
sufficiently large time steps, is significantly polluted by temporal errors.

In the case of under-resolved flows, instead, the influence of the time step
size has shown to be negligible. Spurious oscillations become large enough to
visibly pollute the solution only for very large time steps approaching the
CFL limit. In under-resolved flows, in fact, the spatial errors dominate their
temporal counterparts, overshadowing the influence of the time integration
scheme. In conclusion, the theoretical framework introduced in the first part
of the paper is not always representative of more complex nonlinear dynamics.
For well-resolved flows, the slower convergence has been observed in both the
theoretical framework and in the Euler flow simulations. Consequently, the
choice of the time step size plays a more relevant role in highly accurate
simulations such as DNS of turbulent flows. Instead, in the presence of
significant spatial under-resolution (like in LES), temporal errors are more
likely overshadowed by spatial errors.

## Acknowledgements

The use of the SD solver originally developed by Antony Jameson’s group at
Stanford University is gratefully acknowledged. This work was granted access
to the high-performance computing resources of CRIANN. The PhD scholarship of
the first author is founded by the University of Rouen Normandie. RCM
acknowledges support from São Paulo Research Foundation (FAPESP) via grant
2020/10910-8. GL acknowledges support from the _Agence Nationale de la
Recherche_ (ANR) under grant number ANR-18-CE05-0030. GM acknowledges support
from National University of Singapore start-up grant WBS R-265-000-A36-133.
The authors also thank David Kopriva, Tristan Montoya and David Zingg for
fruitful discussions.

## Appendix

A series of plots for the dissipation of the FR-DG scheme are given in Fig.
28, whereas plots for the FR-SD scheme are given in Fig. 29. The SD curves
look similar to those of DG, but are moderately shifted up (stronger
dissipation, note the vertical axis) and to the left (toward lower
frequencies).

In Fig. 30, the evolution of the multiple modes of the (nearly) centred FR-SD
with RK33 and RK45 are shown for $N=2$. Similar plots for $N=4$ are shown in
Fig. 31.

Figure 28: Dissipation of the FR-DG scheme for RK33 (left) and RK45 (right) at
polynomial orders $N=1,2,3,6$ (top to bottom). Figure 29: Dissipation of the
FR-SD scheme for RK33 (left) and RK45 (right) at polynomial orders $N=1,2,3,6$
(top to bottom). Figure 30: Numerical dissipation for a $3^{\rm{rd}}$ 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). Figure 31:
Numerical dissipation for a $5^{\rm{th}}$ 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).

## References

  * [1] C. Liang, S. Premasuthan, A. Jameson, High-order accurate simulation of low-mach laminar flow past two side-by-side cylinders using spectral difference method, Computers & structures 87 (11-12) (2009) 812–827.
  * [2] C. Liang, A. Jameson, Z. J. Wang, Spectral difference method for compressible flow on unstructured grids with mixed elements, Journal of Computational Physics 228 (8) (2009) 2847–2858.
  * [3] A. Jameson, A proof of the stability of the spectral difference method for all orders of accuracy, J. Sci. Comput. 45 (1) (2010) 348–358.
  * [4] J. Slotnick, A. Khodadoust, J. Alonso, D. Darmofal, W. Gropp, E. Lurie, D. Mavriplis, CFD vision 2030 study: A path to revolutionary computational aerosciences, Tech. rep., NASA Langley Research Center, Hampton, Virginia 23681-2199, NASA Contractor Report 218178 (Mar. 2014).
  * [5] P. E. Vincent, A. Jameson, Facilitating the adoption of unstructured high-order methods amongst a wider community of fluid dynamicists, Mathematical Modelling of Natural Phenomena 6 (3) (2011) 97–140.
  * [6] J. S. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods: algorithms, analysis, and applications, Springer Science & Business Media, 2007\.
  * [7] B. Cockburn, C. Shu, The local discontinuous Galerkin finite element method for convection-diffusion systems, SIAM J. Numer. Anal. 35 (1998) 2440–2463.
  * [8] B. Cockburn, C. Shu, The Runge-Kutta discontinuous Galerkin finite element method for conservation laws V: Multidimensional systems, J. Comput. Phys. 141 (1998) 199–224.
  * [9] D. A. Kopriva, J. H. Kolias, A conservative staggered-grid chebyshev multidomain method for compressible flows, Journal of computational physics 125 (1) (1996) 244–261.
  * [10] Y. Liu, M. Vinokur, Z. Wang, Spectral difference method for unstructured grids I: basic formulation, J. Comput. Phys. 216 (2) (2006) 780–801.
  * [11] Z. Wang, Y. Liu, G. May, A. Jameson, Spectral difference method for unstructured grids II: Extension to the Euler equations, J. Sci. Comput. 32 (1) (2007) 45–71.
  * [12] H. T. Huynh, A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods, in: 18th AIAA Computational Fluid Dynamics Conference, 2007, p. 4079.
  * [13] P. E. Vincent, P. Castonguay, A. Jameson, A new class of high-order energy stable flux reconstruction schemes, Journal of Scientific Computing 47 (1) (2011) 50–72.
  * [14] J.-B. Chapelier, M. De La Llave Plata, F. Renac, E. Lamballais, Evaluation of a high-order discontinuous galerkin method for the dns of turbulent flows, Computers & Fluids 95 (2014) 210–226.
  * [15] M. de la Llave Plata, V. Couaillier, M.-C. Le Pape, On the use of a high-order discontinuous galerkin method for dns and les of wall-bounded turbulence, Computers & Fluids 176 (2018) 320–337.
  * [16] F. Bassi, L. Botti, A. Colombo, A. Crivellini, A. Ghidoni, F. Massa, On the development of an implicit high-order discontinuous galerkin method for dns and implicit les of turbulent flows, European Journal of Mechanics-B/Fluids 55 (2016) 367–379.
  * [17] B. Krank, M. Kronbichler, W. A. Wall, Direct numerical simulation of flow over periodic hills up to $\mathrm{Re}_{H}=10,595$, Flow, turbulence and combustion 101 (2) (2018) 521–551.
  * [18] G. Lodato, N. Tonicello, B. Pinto, Large-eddy simulation of bypass transition on a zero-pressure-gradient flat plate using the spectral-element dynamic model, Flow, Turbulence and Combustion (2021) 1–30.
  * [19] N. Tonicello, G. Lodato, L. Vervisch, Analysis of high-order explicit les dynamic modeling applied to airfoil flows, Flow, Turbulence and Combustion (2021) 1–28.
  * [20] P. Fernandez, N.-C. Nguyen, J. Peraire, On the ability of discontinuous Galerkin methods to simulate under-resolved turbulent flows, arXiv preprint arXiv:1810.09435.
  * [21] R. C. Moura, G. Mengaldo, J. Peiró, S. J. Sherwin, On the eddy-resolving capability of high-order discontinuous Galerkin approaches to implicit LES/under-resolved DNS of Euler turbulence, Journal of Computational Physics 330 (2017) 615–623.
  * [22] M. de la Llave Plata, E. Lamballais, F. Naddei, On the performance of a high-order multiscale dg approach to les at increasing reynolds number, Computers & Fluids 194 (2019) 104306.
  * [23] G. Mengaldo, D. Moxey, M. Turner, R. C. Moura, A. Jassim, M. Taylor, J. Peiro, S. Sherwin, Industry-relevant implicit large-eddy simulation of a high-performance road car via spectral/hp element methods, SIAM Review 63 (4) (2021) 723–755.
  * [24] C. Cox, W. Trojak, T. Dzanic, F. Witherden, A. Jameson, Accuracy, stability, and performance comparison between the spectral difference and flux reconstruction schemes, Computers & Fluids 221 (2021) 104922.
  * [25] S. K. Lele, Compact finite difference schemes with spectral-like resolution, Journal of computational physics 103 (1) (1992) 16–42.
  * [26] C. Bogey, C. Bailly, A family of low dispersive and low dissipative explicit schemes for flow and noise computations, Journal of Computational physics 194 (1) (2004) 194–214.
  * [27] K. Van den Abeele, C. Lacor, Z. J. Wang, On the stability and accuracy of the spectral difference method, Journal of Scientific Computing 37 (2) (2008) 162–188.
  * [28] P. E. Vincent, P. Castonguay, A. Jameson, Insights from von neumann analysis of high-order flux reconstruction schemes, Journal of Computational Physics 230 (22) (2011) 8134–8154.
  * [29] R. C. Moura, S. J. Sherwin, J. Peiró, Linear dispersion–diffusion analysis and its application to under-resolved turbulence simulations using discontinuous Galerkin spectral/hp methods, Journal of Computational Physics 298 (2015) 695–710.
  * [30] J. Vanharen, G. Puigt, X. Vasseur, J.-F. Boussuge, P. Sagaut, Revisiting the spectral analysis for high-order spectral discontinuous methods, Journal of Computational Physics 337 (2017) 379–402.
  * [31] F. Q. Hu, H. L. Atkins, Eigensolution analysis of the discontinuous Galerkin method with nonuniform grids: I. one space dimension, Journal of Computational Physics 182 (2) (2002) 516–545.
  * [32] G. Mengaldo, R. Moura, B. Giralda, J. Peiró, S. Sherwin, Spatial eigensolution analysis of discontinuous Galerkin schemes with practical insights for under-resolved computations and implicit LES, Computers & Fluids 169 (2018) 349–364.
  * [33] G. Mengaldo, D. De Grazia, R. C. Moura, S. J. Sherwin, Spatial eigensolution analysis of energy-stable flux reconstruction schemes and influence of the numerical flux on accuracy and robustness, Journal of Computational Physics 358 (2018) 1–20.
  * [34] R. C. Moura, M. Aman, J. Peiró, S. J. Sherwin, Spatial eigenanalysis of spectral/hp continuous Galerkin schemes and their stabilisation via DG-mimicking spectral vanishing viscosity for high Reynolds number flows, Journal of Computational Physics 406 (2020) 109112.
  * [35] N. Tonicello, G. Lodato, L. Vervisch, A comparative study from spectral analyses of high-order methods with non-constant advection velocities, Journal of Scientific Computing 87 (3) (2021) 1–38.
  * [36] W. Trojak, R. Watson, A. Scillitoe, P. G. Tucker, Effect of mesh quality on flux reconstruction in multi-dimensions, Journal of Scientific Computing 82 (3) (2020) 1–36.
  * [37] G. Sáez-Mischlich, J. Sierra-Ausín, J. Gressier, The spectral difference raviart-thomas method for two and three-dimensional elements and its connection with the flux reconstruction formulation, arXiv preprint arXiv:2105.08632.
  * [38] J. Manzanero, G. Rubio, E. Ferrer, E. Valero, Dispersion-dissipation analysis for advection problems with nonconstant coefficients: Applications to discontinuous Galerkin formulations, SIAM Journal on Scientific Computing 40 (2) (2018) A747–A768.
  * [39] J. Manzanero, G. Rubio, E. Ferrer, E. Valero, D. A. Kopriva, Insights on aliasing driven instabilities for advection equations with application to gauss-lobatto discontinuous Galerkin methods (2017).
  * [40] P. Fernandez, R. C. Moura, G. Mengaldo, J. Peraire, Non-modal analysis of spectral element methods: Towards accurate and robust large-eddy simulations, Computer Methods in Applied Mechanics and Engineering 346 (2019) 43–62.  doi:10.1016/j.cma.2018.11.027.
  * [41] M. Alhawwary, Z. Wang, A combined-mode Fourier analysis of DG methods for linear parabolic problems, SIAM Journal on Scientific Computing 42 (6) (2020) A3825–A3858.
  * [42] H. Yang, F. Li, J. Qiu, Dispersion and dissipation errors of two fully discrete discontinuous galerkin methods, Journal of Scientific Computing 55 (3) (2013) 552–574.
  * [43] B. Vermeire, P. Vincent, On the behaviour of fully-discrete flux reconstruction schemes, Computer Methods in Applied Mechanics and Engineering 315 (2017) 1053–1079.
  * [44] M. Alhawwary, Z. J. Wang, Fourier analysis and evaluation of dg, fd and compact difference methods for conservation laws, Journal of Computational Physics 373 (2018) 835–862.
  * [45] D. De Grazia, G. Mengaldo, D. Moxey, P. Vincent, S. Sherwin, Connections between the discontinuous galerkin method and high-order flux reconstruction schemes, International journal for numerical methods in fluids 75 (12) (2014) 860–877.
  * [46] G. Mengaldo, Discontinuous spectral/hp element methods: development, analysis and applications to compressible flows, Ph.D. thesis, Imperial College London (2015).
  * [47] G. Mengaldo, D. De Grazia, P. E. Vincent, S. J. Sherwin, On the connections between discontinuous galerkin and flux reconstruction schemes: extension to curvilinear meshes, Journal of Scientific Computing 67 (3) (2016) 1272–1292.
  * [48] K. Asthana, A. Jameson, High-order flux reconstruction schemes with minimal dispersion and dissipation, Journal of Scientific Computing 62 (3) (2015) 913–944.
  * [49] M. H. Carpenter, C. A. Kennedy, Fourth-order 2N-storage Runge-Kutta schemes, Tech. rep., NASA Langley Research Center (1994).
  * [50] H. Lomax, T. H. Pulliam, D. W. Zingg, Fundamentals of Computational Fluid Dynamics, Springer Science & Business Media, 2013.
  * [51] R. C. Moura, S. J. Sherwin, J. Peiró, Eigensolution analysis of spectral/hp continuous Galerkin approximations to advection–diffusion problems: Insights into spectral vanishing viscosity, Journal of Computational Physics 307 (2016) 401–422.
  * [52] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in: Advanced numerical approximation of nonlinear hyperbolic equations, Springer, 1998, pp. 325–432.