∎

11institutetext: Will Trojak (corresponding author)22institutetext: Department
of Engineering, University of Cambridge, Cambridge, UK, CB2 1PZ
22email<EMAIL_ADDRESS>
33institutetext: Rob Watson 44institutetext: School of Mechanical and
Aerospace Engineering, Queen’s University Belfast, Belfast, UK, BT9 5AH
55institutetext: Ashley Scillitoe 66institutetext: Alan Turing Institute,
Kings Cross, London, UK, NW1 2DB
77institutetext: Paul G. Tucker 88institutetext: Department of Engineering,
University of Cambridge, Cambridge, UK, CB2 1PZ

# Effect of Mesh Quality on Flux Reconstruction in Multi-Dimensions

Will Trojak  Rob Watson  Ashley Scillitoe  Paul G. Tucker

(Received: 18 September 2018)

###### Abstract

Theoretical methods are developed to understand the effect of non-uniform
grids on Flux Reconstruction (FR) in multi-dimensions. The analysis reveals
that the same effect of expanding and contracting grids is seen in two
dimensions as in one dimension. Namely, that expansions cause instability and
contractions cause excess dissipation. Subsequent numerical experiments on the
Taylor-Green Vortex with jittered elements show the effect of localised
regions of expansion and contraction, with an initial increase in the kinetic
energy observed on non-uniform meshes. Some comparison is made between second-
order FR and second-order finite volume (FV). FR is found to be more resilient
to mesh deformation, however, FV is found to be more resolved when operated at
second order on the same mesh. In both cases, it is recommended that a kinetic
energy preserving/conservation formulation should be used as this can greatly
increase resilience to mesh deformation.

###### MSC:

65M60 65T99 76F65 76M10

##  1 Introduction

Since the inception of Spectral Volume methods by Wang Wang2002 , the
trajectory of high order methods has trended towards the Flux Reconstruction
method of Huynh Huynh2007  and Vincent et al. Vincent2010 . This approach
draws on the work of Finite Elements, see Brenner and Ridgway-Scott
Brenner2008 , enabling the high performance of Flux Reconstruction on
heterogeneous computing — as can be seen in the highly efficient use of vast
computing resources by Vincent et al. Vincent2017 . However, the move towards
high order was not born out of a need for more efficient use of modern HPC
environments. For example, Brandvik and Pullan Brandvik2011  showed that high
throughput could be obtained using second-order Finite Volume (FV) methods.
Instead, the main motivating factor has been the increased uptake by industry
of turbulence resolving methods, such as Large Eddy Simulation (LES), as this
allows for far better exploration of flow physics and moves towards the aim of
computational wind tunnels. The main feature of LES is the modelling of the
very smallest scales of motion, which removes the need for Direct Numerical
Simulation (DNS) levels of resolution. However, Chow and Moin Chow2003  and
Ghosal Ghosal1996  showed that, for LES, the need to keep the truncation error
small to enable the sensible use of sub-grid scale models meant that the grid
requirements were demanding. A move to higher order would mean that the
scaling of the truncation error was to a higher power of grid spacing — thus
lowering the grid requirements and decoupling the scaling of aliasing error
and truncation error. Hence, for wall-resolved LES, calculations are often
impractical unless the improved mesh resolution requirements of high order
methods are considered.

The analytical understanding of Flux Reconstruction has been explored to a
large extent in the work of Vincent et al. Vincent2011 , Jameson et al.
Jameson2012 , and Castonguay et al. Castonguay2014 , where the stability of
linear advection, advection-diffusion, and non-linear problems were presented.
The key findings were the energy stability of FR on linear problems, and the
condition for energy stability on non-linear problems. In addition, by
investigating the dispersion and dissipation characteristics of FR, the
existence of superconvergence after temporal integration and the corresponding
CFL limits were found. This work was limited to one dimension — although still
applicable, the investigation of the exact behaviour of FR in higher
dimensions has been limited, such as that of Williams and Jameson Williams2014
and Sheshadri and Jameson Sheshadri2016 . This work focused primarily on the
proof of the Sobolev type energy stability in 2D in a similar manner to that
of Hesthaven and Warburton Hesthaven2007 , alongside some numerical studies
performed for validation.

The advantage of FR — that leads to high performance on heterogeneous and
massively parallel architectures — is its unstructured and sub-domain nature.
Unstructured grids also allow far more complex geometries to be considered,
but the resulting meshes experience deformation, expansion, and contraction of
the elements. We wish to characterise the performance of FR under these
conditions, and so far the effect of linear mesh deformation on FR has been
considered in one dimension by Trojak et al. Trojak2017a . Therefore, we make
use of the seminal work of Lele Lele1992 , in which the dispersion and
dissipation of finite difference methods were considered in both one and two
dimensions. We wish to repeat this process for FR, but extend it to also
consider deformed grids.

In this paper, we present an extension to the one-dimensional analytical work
of Vincent et al. Vincent2011  and Trojak et al. Trojak2017a . This extension
will be shown for a two-dimensional case on quadrilaterals with rectilinear
mesh stretching, but could also be performed on higher dimensional hypercubes.
From the basis of this more general von Neumann analysis, the behaviour of FR
on linearly mapped meshes can be explored. The investigation has been
restricted to linear transformation as these are of key importance for complex
industrial simulations due to there fundamental nature. For example, they
occur in meshes where mesh generators have simply tessellated elements to fill
the domain. Therefore, understanding their character is key, however, we
should point to some recent work that has numerically investigated curved
meshes Mengaldo2016 ; Yu2014 .

The aim of this work is to understand the effect of moving to higher
dimensionality on key metrics governing scheme performance, such as CFL limit,
dispersion, and dissipation. Finally, the Taylor-Green vortex will be used to
understand how deformed meshes affect full Navier-Stokes calculations, with
reference calculations performed by an industrial second order finite volume
method.

##  2 Flux Reconstruction

Flux Reconstruction Castonguay2012 ; Huynh2007  (FR) applied to the linear
advection equation will form the basis of the initial investigation to be
carried out, and for the reader’s convenience, an overview of the scheme is
presented here. However, for a more detailed understanding, the reader should
consult Castonguay Castonguay2012  or Huynh Huynh2007 . This 1D scheme can be
readily converted to two and three dimensions for quadrilaterals and
hexahedrals, respectively. First, let us consider the one-dimensional
advection equation:

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

The FR method is related to the Discontinuous Galerkin (DG) method Reed1973
and makes use of the same subdivision of the domain into discontinuous sub-
domains:

$\mathbf{\Omega}=\bigcup_{n=1}^{N}{\mathbf{\Omega}_{n}}$ (2)

Within the standardised sub-domain, $\hat{\mathbf{\Omega}}\in\mathbb{R}^{d}$,
computational spatial variables are defined. When $d=1$,
$\hat{\mathbf{\Omega}}=[-1,1]$, using $\xi$ to denote the value taken. This
computational space is discretised with $(p+1)^{d}$ solution points, and
$2d(p+1)^{d-1}$ flux points, placed at the edges of the sub-domain. The
solution and flux point locations are determined using a tensor grid of a 1D
quadrature. Fig. 1(a) shows a 1D example of this. To transform from
$\mathbf{\Omega}_{n}\rightarrow\hat{\mathbf{\Omega}}$, a Jacobian $J_{n}$ is
defined such that:

$\hat{u}^{\delta}=\hat{u}^{\delta}(\xi;t)=J_{n}u^{\delta}(x;t)$ (3)

With this domain set up, we now proceed with defining the steps to construct a
continuous solution from the discontinuous segments. The first stage is to
define a local solution polynomial in $\hat{\mathbf{\Omega}}$ using Lagrange
interpolation.

$\displaystyle l_{k}(\xi)$ $\displaystyle=\prod_{i=0,i\neq
k}^{p}{\frac{\xi-\xi_{i}}{\xi_{k}-\xi_{i}}}$ (4)
$\displaystyle\hat{u}^{\delta}(\xi)$
$\displaystyle=\sum_{i=0}^{p}{\hat{u}^{\delta}_{i}l_{i}(\xi)}$ (5)

Repeating the interpolation for the discontinuous flux in
$\hat{\mathbf{\Omega}}$:

$\hat{f}^{\delta D}=\hat{f}^{\delta D}(\xi,t)=\sum_{i=0}^{p}{\hat{f}^{\delta
D}_{i}l_{i}(\xi)}$ (6)

Here we define $\hat{f}^{\delta D}$ as the transformed discontinuous flux
polynomial. Now using the Jacobian and the solution polynomials, the primitive
values can be calculated in the physical domain $\mathbf{\Omega}_{n}$:

$\overline{u}^{\delta}(x)=\frac{\hat{u}^{\delta}(\xi)}{J_{n}}=\sum_{i=0}^{p}{u^{\delta}_{i}l_{i}(\xi)}$
(7)

(a) Flux and solution point layout for $p=3$ in $\hat{\mathbf{\Omega}}$, with
corresponding left and right Huynh correction functions.

(b) Diagram showing the procedure of correcting one, central, interface.

Figure 1: Point layout in $\hat{\mathbf{\Omega}}$ for $p=3$ and cell interface
topology.

The primitive polynomial can then be interpolated to the interface and define
as $\hat{u}^{\delta}_{l}=\hat{u}^{\delta}(-1)$ and
$\hat{u}^{\delta}_{r}=\hat{u}^{\delta}(1)$. The values at the interface, $I$,
then allow for a common interface flux, $f^{\delta I}_{I}$, to be calculated
in the physical domain. The is shown graphically in Fig. 1(b). For a general
case, this is done using a Riemann solver on the primitives at the interface,
such as: Roe Roe1981 ; flux vector splitting vanLeer1982 ; or HLL Harten1983 .
In order to get a spatially continuous solution over $\mathbf{\Omega}$, the
common interface flux has to be incorporated into the solution. For FR this is
done by using a correction function to propagate the corrected flux into
$\mathbf{\Omega}_{n}$. The definition of the correction function was shown to
be important in the determination of the characteristics of FR, Vincent et al.
Vincent2011 . Primarily the correction function is a polynomial which, in one
dimension, has the boundary conditions:

$\displaystyle h_{L}(-1)=h_{R}(1)=1$ (8) $\displaystyle h_{R}(-1)=h_{L}(1)=0$
(9)

Beyond this, several sets of stable correction functions have been defined,
firstly unified by Vincent et al. Vincent2010  and later expanded in
Vincent2015 ; Trojak2018 ; Trojak2018b ; Trojak2018f . However, in this paper,
we will focus on the correction function defined by Huynh Huynh2007 , in
particular, the Huynh $g_{2}$ correction function, which is shown in Fig.
1(a). We will also consider the correction that can recover Nodal
Discontinuous Galerkin (NDG) Hesthaven2008  as it provides a good point of
comparison, due to DG relevance and maturity. It should be noted that the NDG
correction function will only recover NDG in FR for homogeneous linear flux
functions, due to the different mechanisms of aliasing. At times we will also
explore some of the effects of correction function in two dimensions and for
this, we will restrict ourselves to the Original Stable FR (OSFR) family of
Vincent2010 , for brevity. This is a one-parameter family of the correction
functions for which we will call the control parameter $\iota$.

The correction to the flux function is then calculated using the difference
between the discontinuous and common interface values and a correction
function. The correction is then defined as:

$\hat{f}^{\delta C}=(\hat{f}^{\delta I}_{L}-\hat{f}^{\delta
D}_{L})h_{L}(\xi)+(\hat{f}^{\delta I}_{R}-\hat{f}^{\delta D}_{R})h_{R}(\xi)$
(10)

and, hence the corrected continuous gradient of the flux is then:

$\displaystyle\frac{\partial\hat{f}^{\delta}}{\partial{\xi}}$
$\displaystyle=\frac{\mathrm{d}\hat{f}^{\delta
D}}{\mathrm{d}\xi}+\frac{\mathrm{d}\hat{f}^{\delta C}}{\mathrm{d}\xi}$ (11)
$\displaystyle=\sum_{j=0}^{p}\hat{f}^{\delta
D}_{j}\frac{\mathrm{d}l_{j}(\xi)}{\mathrm{d}\xi}+(\hat{f}^{\delta
I}_{L}-\hat{f}^{\delta
D}_{L})\frac{\mathrm{d}h_{L}(\xi)}{\mathrm{d}\xi}+(\hat{f}^{\delta
I}_{R}-\hat{f}^{\delta D}_{R})\frac{\mathrm{d}h_{R}(\xi)}{\mathrm{d}\xi}$ (12)

Finally, the solution is advanced in time following Eq. (13) — this can be
done via a sensible choice of temporal integration.

$\frac{\partial\hat{u}^{\delta}}{\partial{t}}=-\frac{\partial\hat{f}^{\delta}}{\partial{\xi}}$
(13)

The method detailed here was shown for simplicity in one dimension, but in
subsequent sections, the method for extending this to a high dimension will be
detailed. We will briefly state that to increase the dimensionality of the
method a tensor product is used, which is this same method as is used in
analysis and in the formal implementation of a solver for hypercube elements.

##  3 Two-Dimensional Von Neumann Analysis

The procedure for investigating the dispersion and dissipation properties of
finite element methods has been laid out in some detail by Huynh Huynh2007 ,
Hesthaven and Warburton Hesthaven2007 , and Vincent et al. Vincent2011 . It is
broadly classified as a von Neumann analysis. The procedure was, however, only
performed in 1D, with critical insight into the analytical performance of FR
when applied to more realistic problems overlooked. Extension of the analysis
to higher dimension domains was performed by Lele Lele1992  for various finite
difference schemes. This did, however, avoid the increased complexity of
finite element von Neumann analysis. To begin our extension we introduce the
2D linear advection equation:

$\displaystyle\frac{\partial u}{\partial{t}}+\nabla\cdot\mathbf{F}$
$\displaystyle=0$ (14) $\displaystyle\mathbf{F}=\begin{bmatrix}f\\\
g\end{bmatrix}$ $\displaystyle=u\mathbf{a}=\begin{bmatrix}au\\\
bu\end{bmatrix}$ (15)

Flux reconstruction then uses the superposition of the discontinuous and
corrected flux divergence, meaning Eq. (15) can be rewritten as:

$\frac{\partial\mathbf{u}_{i,j}}{\partial{t}}=-\nabla\cdot\mathbf{F}^{\delta
D}_{i,j}-\nabla\cdot\mathbf{F}^{\delta C}_{i,j}$ (16)

Taking the following definition of the Jacobian, the computational-physical
domain transformation can be defined:

$\displaystyle\mathbf{G}$ $\displaystyle=\begin{bmatrix}\frac{\partial
x}{\partial{\xi}}\>\frac{\partial y}{\partial{\xi}}\\\ \frac{\partial
x}{\partial{\eta}}\>\frac{\partial
y}{\partial{\eta}}\end{bmatrix}=\begin{bmatrix}G_{1}\>G_{2}\\\
G_{3}\>G_{4}\end{bmatrix}~{}\quad\mathrm{and}\quad J=|\mathbf{G}|$ (17)
$\displaystyle u$
$\displaystyle=J^{-1}\hat{u},\quad\mathbf{F}=J^{-1}\mathbf{G}\hat{\mathbf{F}},\quad\nabla\cdot\mathbf{F}=J^{-1}\hat{\nabla}\cdot\hat{\mathbf{F}}$
(18)

where we use $\hat{\nabla}$ to mean
$[\frac{\partial}{\partial{\xi}},\frac{\partial}{\partial{\eta}}]^{T}$ in 2D.
From the work of Huynh Huynh2007 , Castonguay Castonguay2012 , and Sheshadri
et al. Sheshadri2016a , Eq.(15) is written in two dimensions as:

$\hat{\nabla}\cdot\hat{\mathbf{F}}^{\delta
D}=\sum^{p}_{i=0}\sum^{p}_{j=0}\hat{f}^{\delta D}_{i,j}\frac{\partial
l_{i}(\xi)}{\partial{\xi}}l_{j}(\eta)+\sum^{p}_{i=0}\sum^{p}_{j=0}\hat{g}^{\delta
D}_{i,j}\frac{\partial l_{j}(\eta)}{\partial{\eta}}l_{i}(\xi)\\\ $ (19)
$\begin{split}\hat{\nabla}\cdot\hat{\mathbf{F}}^{\delta
C}=\sum^{p}_{i=0}\Bigg{(}&(\hat{f}^{\delta I}_{L,i}-\hat{f}^{\delta
D}_{L,i})\frac{\mathrm{d}h_{L,i}}{\mathrm{d}\xi}+(\hat{f}^{\delta
I}_{R,i}-\hat{f}^{\delta D}_{R,i})\frac{\mathrm{d}h_{R,i}}{\mathrm{d}\xi}+\\\
&(\hat{g}^{\delta I}_{B,i}-\hat{g}^{\delta
D}_{B,i})\frac{\mathrm{d}h_{B,i}}{\mathrm{d}\eta}+(\hat{g}^{\delta
I}_{T,i}-\hat{g}^{\delta
D}_{T,i})\frac{\mathrm{d}h_{T,i}}{\mathrm{d}\eta}\Bigg{)}\end{split}$

where we use $L$, $R$, $B$, and $T$ subscripts to mean left, right, bottom,
and top respectively. We will then impose that grid transformations are purely
rectilinear, _i.e._ $G_{2}=G_{3}=0$. This is to reduce the number of dependent
variables will still allowing an important form of grid deformation to be
investigated. We may now use Eq. (19) and convert it into a matrix form:

$\displaystyle\hat{\nabla}\cdot\hat{\mathbf{F}}^{\delta D}$
$\displaystyle=\mathbf{D}_{\xi}\hat{\mathbf{f}}^{\delta}_{i,j}+\mathbf{D}_{\eta}\hat{\mathbf{g}}^{\delta}_{i,j}$
(20) $\displaystyle\nabla\cdot\mathbf{F}^{\delta D}$
$\displaystyle=G^{-1}_{1,i,j}\mathbf{D}_{\xi}\mathbf{f}^{\delta}_{i,j}+G^{-1}_{4,i,j}\mathbf{D}_{\eta}\mathbf{g}^{\delta}_{i,j}$
(21)

To apply the correction function, we need to calculate the interface values
around the element. For the case of generalised central/upwinding with
upwinding ratio $\alpha$, the common interface fluxes may be written as:

$\displaystyle G^{-1}_{4,i,j}\hat{f}^{\delta I}_{L}$
$\displaystyle=a\big{(}\alpha G^{-1}_{4,i-1,j}\hat{u}^{\delta}_{i-1,j,R}$
$\displaystyle+(1-\alpha)G^{-1}_{4,i,j}\hat{u}^{\delta}_{i,j,L}\big{)}$ (22)
$\displaystyle G^{-1}_{4,i,j}\hat{f}^{\delta I}_{R}$
$\displaystyle=a\big{(}\alpha G^{-1}_{4,i,j}\hat{u}^{\delta}_{i,j,R}$
$\displaystyle+(1-\alpha)G^{-1}_{4,i+1,j}\hat{u}^{\delta}_{i+1,j,L}\big{)}$
(23) $\displaystyle G^{-1}_{1,i,j}\hat{g}^{\delta I}_{B}$
$\displaystyle=b\big{(}\alpha G^{-1}_{1,i,j-1}\hat{u}^{\delta}_{i,j-1,T}$
$\displaystyle+(1-\alpha)G^{-1}_{1,i,j}\hat{u}^{\delta}_{i,j,B}\big{)}$ (24)
$\displaystyle G^{-1}_{1,i,j}\hat{g}^{\delta I}_{T}$
$\displaystyle=b\big{(}\alpha G^{-1}_{1,i,j}\hat{u}^{\delta}_{i,j,T}$
$\displaystyle+(1-\alpha)G^{-1}_{1,i,j+1}\hat{u}^{\delta}_{i,j+1,B}\big{)}$
(25)

where $\alpha=1$ gives rise to upwinding and $\alpha=0.5$ produces central
difference. Hence the divergence correction can be written as:

$\begin{split}\hat{\nabla}\cdot\hat{\mathbf{F}}^{\delta
C}_{i,j}=a\alpha\bigg{(}G^{-1}_{4,i-1,j}\mathbf{h_{L}}\mathbf{l_{R}}^{T}\hat{\mathbf{u}}^{\delta}_{i-1,j}-G^{-1}_{4,i,j}\mathbf{h_{L}}\mathbf{l_{L}}^{T}\hat{\mathbf{u}}^{\delta}_{i,j}\bigg{)}&+\\\
a(1-\alpha)\bigg{(}G^{-1}_{4,i+1,j}\mathbf{h_{R}}\mathbf{l_{L}}^{T}\hat{\mathbf{u}}^{\delta}_{i+1,j}-G^{-1}_{4,i,j}\mathbf{h_{R}}\mathbf{l_{R}}^{T}\hat{\mathbf{u}}^{\delta}_{i,j}\bigg{)}&+\\\
b\alpha\bigg{(}G^{-1}_{1,i,j-1}\mathbf{h_{B}}\mathbf{l_{T}}^{T}\hat{\mathbf{u}}^{\delta}_{i,j-1}-G^{-1}_{1,i,j}\mathbf{h_{B}}\mathbf{l_{B}}^{T}\hat{\mathbf{u}}^{\delta}_{i,j}\bigg{)}&+\\\
b(1-\alpha)\bigg{(}G^{-1}_{1,i,j+1}\mathbf{h_{T}}\mathbf{l_{B}}^{T}\hat{\mathbf{u}}^{\delta}_{i,j+1}-G^{-1}_{1,i,j}\mathbf{h_{T}}\mathbf{l_{T}}^{T}\hat{\mathbf{u}}^{\delta}_{i,j}\bigg{)}&\end{split}$
(26)

where $\mathbf{h_{L}}$ is the _gradient_ of the left correction function at
the solution points and again $\mathbf{l_{L}}$ are the values of the
polynomial basis at the left interface and so on for $R$, $T$, and $B$.
Therefore, by grouping terms by their cell indexing and transforming each term
into the physical domain:

$\begin{split}\frac{\partial\mathbf{u}_{i,j}}{\partial{t}}=&-a\Big{(}G_{1,i-1,j}^{-1}\mathbf{C}_{L}\mathbf{u}^{\delta}_{i-1,j}+G_{1,i,j}^{-1}\mathbf{C}_{0\xi}\mathbf{u}^{\delta}_{i,j}+G_{1,i+1,j}^{-1}\mathbf{C}_{R}\mathbf{u}^{\delta}_{i+1,j}\Big{)}\\\
&-b\Big{(}G_{4,i,j-1}^{-1}\mathbf{C}_{B}\mathbf{u}^{\delta}_{i,j-1}+G_{4,i,j}^{-1}\mathbf{C}_{0\eta}\mathbf{u}^{\delta}_{i,j}+G_{4,i,j+1}^{-1}\mathbf{C}_{T}\mathbf{u}^{\delta}_{i,j+1}\Big{)}\end{split}$
(27)

where

$\displaystyle\mathbf{C}_{L}$
$\displaystyle=\alpha\mathbf{h_{L}}\mathbf{l_{R}}^{T}\quad\mathbf{C}_{R}$
$\displaystyle=(1-\alpha)\mathbf{h_{R}}\mathbf{l_{L}}^{T}\quad\mathbf{C}_{0\xi}$
$\displaystyle=\mathbf{D}_{\xi}-\alpha\mathbf{h_{L}}\mathbf{l_{L}}^{T}-(1-\alpha)\mathbf{h_{R}}\mathbf{l_{R}}^{T}$
(28) $\displaystyle\mathbf{C}_{B}$
$\displaystyle=\alpha\mathbf{h_{B}}\mathbf{l_{T}}^{T}\quad\mathbf{C}_{T}$
$\displaystyle=(1-\alpha)\mathbf{h_{T}}\mathbf{l_{B}}^{T}\quad\mathbf{C}_{0\eta}$
$\displaystyle=\mathbf{D}_{\eta}-\alpha\mathbf{h_{B}}\mathbf{l_{B}}^{T}-(1-\alpha)\mathbf{h_{T}}\mathbf{l_{T}}^{T}$
(29)

Finally, we are here interested in the frequency response of the system, and,
importantly to engineers and technicians, how the cell’s orientation relative
to an oncoming wave affects performance. Therefore, we impose a trial solution
of the form:

$u(x,y;t)=\exp{(ik(x\cos\theta+y\sin\theta-ct))}$ (30)

and by substitution into Eq.(15), the advection velocity, $\mathbf{a}$, can be
found, which is shown schematically in Fig. 2(a).

$\mathbf{a}=\begin{bmatrix}a\\\ b\end{bmatrix}=\begin{bmatrix}\cos\theta\\\
\sin\theta\end{bmatrix}$ (31)

(a) Schematic showing inclined plane wave passing through a cell with a
geometrically transformed rectilinear stencil.

(b) Schematic showing inclined plane wave passing through a cell with a
geometrically transformed rectilinear stencil.

Figure 2: Linear advection schematic for two and three dimensions.

The plane wave can then be projected into the computational domain and
discretised as:

$\mathbf{u}_{i,j}=\mathbf{v}\exp{\bigg{(}ik\Big{(}\big{(}0.5(\xi+1)\delta_{i}+x_{i}\big{)}\cos\theta+\big{(}0.5(\eta+1)\delta_{j}+y_{i}\big{)}\sin\theta-
ct\Big{)}\bigg{)}}$ (32)

where, for brevity, $\delta_{i}=x_{i}-x_{i-1}$ and $\delta_{j}=y_{j}-y_{j-1}$
are defined. Inserting Eq.(32) into Eq.(27), an Eigenvalue problem can be
obtained as:

$\begin{split}-ikc(k)\mathbf{v}=&-\cos\theta\Big{(}G_{1,i-1,j}^{-1}\mathbf{C}_{L}\exp{\big{(}-ik\delta_{i-1}\cos\theta\big{)}}+G_{1,i,j}^{-1}\mathbf{C}_{0\xi}+G_{1,i+1,j}^{-1}\mathbf{C}_{R}\exp{\big{(}ik\delta_{i}\cos\theta\big{)}}\Big{)}\mathbf{v}\\\
&-\sin\theta\Big{(}G_{4,i,j-1}^{-1}\mathbf{C}_{B}\exp{\big{(}-ik\delta_{j-1}\sin\theta\big{)}}+G_{4,i,j}^{-1}\mathbf{C}_{0\eta}+G_{4,i,j+1}^{-1}\mathbf{C}_{T}\exp{\big{(}ik\delta_{j}\sin\theta\big{)}}\Big{)}\mathbf{v}\end{split}$
(33)

where $\Re{k(c(k))}=\Re{(\omega)}$ and $\Im{(kc(k))}=\Im{(\omega)}$ are the
dispersion and dissipation, respectively, and $\omega$ is the modified angular
frequency response of the system. BY studying the trial solution of Eq.(30) it
can be understood that if $\Im{(\omega)}>0$ then the amplitude of the wave
will increase and vice versa. Furthermore, if $\Re{(\omega)}\neq k$ then a
wave will move at a different speed compared to the other waves inside the
packet, causing the quality of the interpolation to be affected as the
solution is advanced in time as the components that make it up move at
different speeds. An important point is the difference between phase velocity
$\omega/k$ and group velocity $\mathrm{d}\omega/\mathrm{d}k$. Phase velocity
is the speed of a wave in a packet of waves. Group velocity is the speed of
the packet. Therefore, changes to $\mathrm{d}\omega/\mathrm{d}k$, i.e can be
thought of more as a change to the physics due to the numerical method.

Equation (33) can alternatively be cast in the form of an update equation. If
initially Eq.(27) is combined with Eq.(32), then a new matrix,
$\mathbf{Q}_{i,j}$, can be defined:

$\frac{\partial\mathbf{u}_{i,j}}{\partial{t}}=\mathbf{Q}_{i,j}\mathbf{u}_{i,j}$
(34)
$\begin{split}\mathbf{Q}_{i,j}=&-J_{i,j}\cos\theta\Big{(}G_{4,i-1,j}^{-1}\mathbf{C}_{L}\exp{\big{(}-ik\delta_{i-1}\cos\theta\big{)}}+G_{4,i,j}^{-1}\mathbf{C}_{0\xi}+G_{4,i+1,j}^{-1}\mathbf{C}_{R}\exp{\big{(}ik\delta_{i}\cos\theta\big{)}}\Big{)}\\\
&-J_{i,j}\sin\theta\Big{(}G_{1,i,j-1}^{-1}\mathbf{C}_{B}\exp{\big{(}-ik\delta_{j-1}\sin\theta\big{)}}+G_{1,i,j}^{-1}\mathbf{C}_{0\eta}+G_{1,i,j+1}^{-1}\mathbf{C}_{T}\exp{\big{(}ik\delta_{j}\sin\theta\big{)}}\Big{)}\end{split}$
(35)

This definition of the semi-discrete FR operator, $\mathbf{Q}$, can then be
used to form what is called the update equation by imposing some temporal
discretisation. As such we may write:

$\displaystyle\mathbf{u}^{n+1}_{i,j}$
$\displaystyle=\mathbf{R}(\mathbf{Q}_{i,j})\mathbf{u}^{n}_{i,j}$ (36)
$\displaystyle\mathbf{R}_{33}$
$\displaystyle=\mathbf{I}+\frac{\tau\mathbf{Q}_{i,j}}{1!}+\frac{(\tau\mathbf{Q}_{i,j})^{2}}{2!}+\frac{(\tau\mathbf{Q}_{i,j})^{3}}{3!}$
(37)

where the superscript denotes the time level, and our update matrix is
$\mathbf{R}$. Shown here is also an example definition for $\mathbf{R}$ for a
3-step $3^{\mathrm{rd}}$-order Runge-Kutta time integration scheme. Finally,
in keeping with von Neumann’s theorems Isaacson1994 ; Lax1956  and Banach’s
fixed point theorem Kress1998 , the spectral radius of $\mathbf{R}$ has to be
less than or equal to 1 for stability. $\rho(\mathbf{R})\leqslant
1\>\forall\>k\in\mathbb{R}$.

In recent works by Vermerie et al. Vermeire2017a  and Trojak et al.
Trojak2018c , the Fourier analysis was extended by fully discretising the
equation. This is performed by taking Eq.(36) and again applying Eq.(30). This
results in:

$\exp{(-ik(c-1)\tau)}\mathbf{v}=\lambda\mathbf{v}=\exp{(ik\tau)}\mathbf{R}(k,\tau)\mathbf{v}$
(38)

where the time step from $n$ to $n+1$ is $\tau$. Hence, rearranging for the
modified wave speed:

$c=\frac{i\log{(\lambda)}}{k\tau}+1$ (39)

where $\lambda$ are the eigenvalues of $\exp{(ik\tau)}\mathbf{R}$. The
advantage of this further analysis is that it gives the dispersion and
dissipation relations of the full scheme as would be experienced when applied
as implicit LES.

Lastly, the linear FR operator matrix, $\mathbf{Q}_{i,j}$, can be diagonalised
as:

$\mathbf{Q}_{i,j}=ik\mathbf{W}_{i,j}\mathbf{\Lambda}_{i,j}\mathbf{W}^{-1}_{i,j}$
(40)

where $\mathbf{W}$ is a matrix of eigenvectors and $\mathbf{\Lambda}$ is a
diagonal matrix of normalised eigenvalues. From this, the weights of modes
used to reconstruct the solution can be found as:

$\mathbf{u}_{i,j}=\mathbf{W}_{i,j}\boldsymbol{\beta}_{i,j}$ (41)

where $\boldsymbol{\beta}$ is an array of mode weights used to project
$\mathbf{u}_{i,j}$ into the functional space of FR. The posedness of the
projection can then be measured for wavenumbers using the number of the matrix
of modes defined as:

$\kappa(\mathbf{W}_{i,j})=\frac{\sigma_{\mathrm{max}}(\mathbf{W}_{i,j})}{\sigma_{\mathrm{min}}(\mathbf{W}_{i,j})}$
(42)

where $\sigma(\mathbf{W}_{i,j})$ is a singular value of $\mathbf{W}_{i,j}$,
with the matrix becoming singular as $\kappa\rightarrow\infty$.

The results of this section can then be extended to n-dimensions and, in
particular, the three–dimensional case will be investigated to show the
continuation of trends with higher dimensionality. The analysis can broadly be
repeated and is excluded for brevity, but importantly the prescribed solution
is taken as:

$u=\exp{\big{(}ik(x\cos\phi\cos\theta+y\cos\phi\sin\theta+z\sin\phi-
ct)\big{)}}$ (43)

where the angles are as shown in Fig. 2(b), and hence the 3D convective
velocities for linear advection are:

$\mathbf{a}=\begin{bmatrix}\cos\phi\cos\theta\\\ \cos\phi\sin\theta\\\
\sin\phi\end{bmatrix}$ (44)

##  4 Analytical Findings

The analytical methods presented in section 3 allow us to investigate many
properties of FR, however from Eq.(35-37) it can be seen that the functional
space of $\mathbf{Q}$ is 8 dimensional, leading to the functional space of
$\rho(\mathbf{R})$ being 9 dimensional
($\tau,\gamma_{x},\gamma_{y},\Delta_{x},\Delta_{y},k,\theta,\iota,p$).
Therefore we need to restrict our investigation to some key results relating
to grid deformation. Firstly, understanding the dispersion and dissipation
($\Re{(\omega)}$ & $\Im{(\omega)}$) in 2D for both uniform and stretched grids
will be important. Secondly, we wish to understand how higher dimensionality
and grid deformations affect the temporal stability of FR through evaluation
of the CFL limits Courant1967 . Here, the dispersion and dissipation relations
will be useful in explaining the trends seen and will aid in linking this work
to that of Trojak et al. Trojak2017a . The definition of the CFL number in
higher dimensions will be taken as:

$\mathrm{CFL}_{d}=\tau\sum^{d}_{i=1}\frac{a_{i}}{\Delta_{i}}$ (45)

where $d$ is the dimensionality, $\tau$ is the time step and, $a_{i}$ and
$\Delta_{i}$ are the characteristic velocity and grid spacing in the
$i^{\mathrm{th}}$ dimension, respectively. The CFL limit is then the maximum
value of CFL at which the scheme is stable in a von Neumann sense. Finally, we
wish to understand if correction functions can be used to alleviate any
effects of deformation by understanding how the scheme properties vary with
correction function. Within this study, the link to the posed nature of the
linear system will be explored with regard to how this relates to the other
properties.

###  4.1 Review of 1D Grid Expansion

Before commencing with the Fourier/von Neumann analysis in 2D, we will give a
brief review of the behaviour exhibited in one dimension, for which more
detail can be found in Trojak et al. Trojak2017a . Figure 3 shows the results
for Huynh $g_{2}$ correction functions for $p=4$ on a few geometrically
expanding grids.

(a) Dispersion

(b) Dissipation

Figure 3: Upwinded 1D FR, $p=4$, with Huynh $g_{2}$ correction functions for
several geometric expansion ratios.

Starting by considering the dispersion, Fig. 3(a) shows that expanding grids
cause dispersion overshoot, as is seen by the dispersion relation going above
the line $\hat{\omega}=1$. While contractions cause dispersion to undershoot.
More interesting is the impact of the grid on the dissipation as is shown in
Fig. 3(b). Here expanding grids are shown to have a region where
$\Im{(\hat{\omega})}>0$ at low wavenumbers, while contracting grids in the
same wavenumber range have $\Im{(\hat{\omega})}<0$. This implies that a wave
propagating across a contracting grid is more quickly dissipated, whereas on
an expanding grid the amplitude of the wave is increased. Hence, expanding
grids are strictly unstable. This result was first presented in Trojak2017a
and was confirmed with numerical experiments.

This characteristic is obviously of importance as it could lead to the scheme
being non-conservative. Hence, the importance of studying this phenomena in
higher dimensions.

###  4.2 Effect of Grid Expansion on Dispersion and Dissipation

For the higher dimensional case, we begin by considering the dispersion and
dissipation on a uniform grid in two dimensions. We are concerned here with
the primary mode — as FR has multiple modes, this is the one that physically
represents the wave. Although, as was found by Asthana Asthana2017 , this may
not be how the energy distributes itself. We identify the physical mode as the
mode whose dispersion relation that goes through zero and dissipation relation
similar to those seen in Fig. 3(b).

The dispersion and dissipation relations are then shown in Figs. 4 and 5. It
is clear that for all orders FR becomes more dispersive and dissipative at
$\theta=45^{\circ}$. This is more easily seen for the dissipation relation,
but in the case of dispersion is displayed by the indent at $\hat{k}\approx
3\pi/4$. Furthermore, there doesn’t seem to be any widening in the range of
angles over which FR becomes more dispersive and dissipative as the order is
reduced. By comparison with the results of Lele Lele1992 , where a similar
test is performed for standard and compact difference schemes, FR shows a
comparatively smaller change in performance as the angle is varied. It is
thought that this due to the method of polynomial fitting used by FR, namely
that this implementation of FR used a tensor grid of monomials _i.e._ , the
number of solution points is $(p+1)^{d}$ and hence the monomials in the
interpolation go from $(\xi^{0}\eta^{0},\xi^{1}\eta^{0}\dots\xi^{p}\eta^{p})$.
By contrast, finite differences do not include the cross product terms, which
will become increasingly dominant as the angle is increased.

(a) $p=2$

(b) $p=3$

(c) $p=4$

(d) $p=5$

Figure 4: Primary mode dispersion for 2D upwinded FR, with Huynh $g_{2}$
corrections, at various orders. Normalised wavenumber as radial distance
(markers at $\pi/4$ intervals), and element angle of incidence as azimuthal
distance.

(a) $p=2$

(b) $p=3$

(c) $p=4$

(d) $p=5$

Figure 5: Primary mode dissipation for 2D upwinded FR, with Huynh $g_{2}$
corrections, at various orders. Normalised wavenumber as radial distance
(markers at $\pi/4$ intervals), and element angle of incidence as azimuthal
distance.

Moving on, we then consider the impact of non-uniform grids on the character
of the dispersion and dissipation. There are two cases that have been
identified from previous work as being of interest. Firstly, when the grid is
expanding, does this give the same positive dissipation in higher dimensions
as seen in the 1D case, Fig. 3(b)? Secondly, if positive dissipation is seen
in the 2D case, will an orthogonal contraction help to stabilise the grid? For
example, if $\gamma_{x}=1.1$ then will setting $\gamma_{y}=0.9$ help to reduce
the positive dissipation?

(a) Dispersion, $\gamma_{x}=1.1$, $\gamma_{y}=1$

(b) Dissipation, $\gamma_{x}=1.1$, $\gamma_{y}=1$

(c) Dispersion, $\gamma_{x}=1.1$, $\gamma_{y}=0.9$

(d) Dissipation, $\gamma_{x}=1.1$, $\gamma_{y}=0.9$

Figure 6: Two dimensional upwinded FR, $p=3$ with Huynh $g_{2}$ corrections,
for different grid expansion factors. Normalised wavenumber as radial distance
(markers at $\pi/4$ intervals), and element angle of incidence as azimuthal
distance. The solid black line on the dissipation plots is the contour of zero
dissipation.

The first of these questions is explored in Figs. 6(a) & 6(b). It is observed
in Fig. 6(b) that expanding grids do cause positive dissipation in higher
dimensions. This is more clearly displayed for some specific angles in Fig.
7(d), where the dissipation is seen to be slightly positive at low
wavenumbers. In one dimension, this behaviour was previously explained that as
a wave moves through elements of different size the group velocity
($\mathrm{d}\omega/\mathrm{d}k$) will change. For an expanding grid, this
leads to low wavenumber energy collecting in elements. The same mechanism
looks to be responsible in higher dimensions.

(a) Dispersion: $\gamma_{x}=\gamma_{y}=1$

(b) Dissipation: $\gamma_{x}=\gamma_{y}=1$

(c) Dispersion: $\gamma_{x}=1.1$, $\gamma_{y}=1$

(d) Dissipation: $\gamma_{x}=1.1$, $\gamma_{y}=1$

Figure 7: Two dimensional upwinded FR, $p=3$, with Huynh $g_{2}$ corrections
at selected incident angles.

The next question of whether a contraction orthogonal to the expansion will
help to stabilise the situation is considered in Figs. 6(c) and 6(d). From
Fig. 6(d), it can be seen that the answer is yes. However, there is still a
region of positive dissipation for $\theta<45^{\circ}$ and as the incidence
angle of the wave approaches $\theta=0$ the stabilisation brought about by the
contraction decays. The $\theta=0$ case is identical to the case of
$\gamma_{y}=1$.

For completeness, we have included plots of the dispersion and dissipation for
NDG in Fig. 8. By comparing the un-stretched and stretched results of Fig.
8(a)-8(b) to those of Fig. 8(c)-8(d) it is clear that the effect of stretching
is similar to that seen with the Huynh $g_{2}$ corrections in Fig. 6. However,
the notable between NDG and Huynh’s $g_{2}$ corrections is that NDG exhibits
more dissipation at higher wavenumbers, as well as dispersion overshoot. This
is difficult to see precisely in Fig. 8, but is indicated by the darker band
at $\hat{k}\approx 3\pi/4$. Furthermore, the location of the zero dissipation
contour is different between the two correction functions, which for NDG is at
a higher wavenumber. This is caused by the unstretched NDG having lower
dissipation at lower wavenumbers. As is seen by comparing Fig. 8(b) and Fig.
6(b).

(a) Dispersion, $\gamma_{x}=\gamma_{y}=1.0$

(b) Dissipation, $\gamma_{x}=\gamma_{y}=1.0$

(c) Dispersion, $\gamma_{x}=1.1$, $\gamma_{y}=0.9$

(d) Dissipation, $\gamma_{x}=1.1$, $\gamma_{y}=0.9$

Figure 8: Two dimensional upwinded FR, $p=3$ with NDG corrections, for
different grid expansion factors. Normalised wavenumber as radial distance
(markers at $\pi/4$ intervals), and element angle of incidence as azimuthal
distance. The solid black line on the dissipation plots is the contour of zero
dissipation.

Lastly, for this investigation into the dispersion and dissipation
characteristics of FR, we wish to make a note on the Nyquist frequency of the
elements. As was also the case for 1D stretched elements, the Nyquist
frequency has a dependency on the expansion ratio. This is found from the
harmonic mean of the 1D Nyquist frequencies, then normalised by the element
size at that angle. Hence, the normalised wavenumber is then:

$\hat{k}=k/k_{\mathrm{nq}}=k\max{\big{\\{}\cos{(\theta)},\sin{(\theta)}\big{\\}}}\bigg{(}\frac{1}{p+1}\bigg{)}\sqrt{\bigg{(}\frac{\cos{(\theta)}}{\gamma_{x}}\bigg{)}^{2}+\bigg{(}\frac{\sin{(\theta)}}{\gamma_{y}}\bigg{)}^{2}}$
(46)

###  4.3 Effect of Grid Aspect Ratio on CFL Limit

Beyond the resolution of the scheme, is the question of setting up a case and
running it on some machine. For this, knowledge of the temporal stability is
key and we will begin by looking at the effect of the relative size of an
element in $x$ and $y$ on the CFL limits of FR. In this case, the grid is not
expanding or contracting, merely the ratio of $\Delta_{x}$ to $\Delta_{y}$ is
varied.

As can be seen from Fig. 9 there is a clear impact on the CFL limit of FR when
elements are rectangular, with no change in the CFL limit when the waves are
aligned with the grid. What is evident is that the angle of incidence where
the CFL limit is smallest, for a given size ratio, is when a wave is incident
at an angle of $\tan^{-1}{(\Delta_{y}/\Delta_{x})}$. This angle corresponds to
the maximum length across the element and hence, when the wave is decomposed
into $x$ and $y$ components, the wavenumbers are lowest — _i.e._ ,
$\inf_{\theta\in\mathbb{R}}({\sup{\\{\Delta_{y}\cos{(\theta)},\Delta_{x}\sin{(\theta)}\\}})}$
when $\theta=\tan^{-1}{(\Delta_{y}/\Delta_{x})}$. As the wavenumbers will be
at their lowest necessary to form the wave, from the 1D dissipation of FR, the
dissipation will also be at its lowest. Therefore, there is less dissipation
in the spatial scheme available to counteract the negative dissipation of the
temporal scheme, hence reducing temporal stability at
$\theta=\tan^{-1}{(\Delta_{y}/\Delta_{x})}$.

(a) $p=3$

(b) $p=4$

Figure 9: Effect of varying: the relative size an element in $x$ and $y$; the
angle of incidence; and CFL in 2D for upwinded FR. The method of temporal
integration used here was RK44. ( ) $\Delta_{y}/\Delta_{x}=\tan{\theta}$

###  4.4 Effect of Grid Expansion on CFL Limit

Now we introduce to the grid an expansion or contraction in $x$ and $y$, with
varying incident angles. The results of these distortions are shown in Fig.
10. The minimum CFL limit is seen to be at $\theta=45^{\circ}$, with temporal
performance peaking as expected at $\theta=0^{\circ},90^{\circ}$. This result
agrees with that of Fig. 9. However, this also shows that the angle of minimum
CFL is only dependent on the local element shape, and in the case investigated
here the central element is always square. Furthermore, when the CFL limit in
the quasi-1D case ($\theta=0^{\circ},90^{\circ}$) is compared to the results
of Trojak et al. Trojak2017a , the CFL limit is found to be lower than the 1D
case. This may be due in part to the increased modes of the system and their
coupling leading to a less stable system. This decrease is corroborated by
numerical tests, and the analytical reasoning will follow shortly.

A second point to note is that for non-grid-aligned waves the expansion or
contraction in both components affects stability. With a contraction
orthogonal to expansion again helping to stabilise the scheme. This point is
subtle, because if the decomposition of the wave into $x$ and $y$ were
linearly independent then it would be expected that the lowest CFL limit would
dominate. However, this result demonstrates that there is coupling between the
$x$ and $y$ components — which could be used advantageously, as was discussed
earlier.

(a) $\gamma_{y}=0.8$

(b) $\gamma_{y}=0.9$

(c) $\gamma_{y}=1.0$

(d) $\gamma_{y}=1.1$

Figure 10: CFL limit for 2D linear advection with FR ($p=4$) using Huynh
correction function, showing variation with $\theta$ and $\gamma_{x}$ for some
set values of $\gamma_{y}$. Time integration is RK44.

Throughout the analytical tests in which waves were injected at incidence on a
square central element, the Nyquist wavenumber was found to be
$k_{nq,\theta}=k_{nq,0}/\cos(\theta)$ for $0^{\circ}\leqslant\theta\leqslant
45^{\circ}$. This result can be understood in one of two ways. Firstly that a
wave at an angle can draw on more points in the normal direction to form a fit
of higher wavenumbers. Or that the wave can be thought of as being a coupled
decomposition of the wave into the $x$ and $y$ directions, and although it has
been said that these are not independent, this does mean that higher
wavenumbers can be supported.

To understand further why expanding meshes are less stable and contracting
meshes more stable it can be illuminating to consider the 1D linear advection
equation, which for upwinded FR can be written as:

$\frac{\partial\mathbf{u}_{j}}{\partial
t}=-J^{-1}_{j}\mathbf{C}_{0}\mathbf{u}_{j}-J^{-1}_{j-1}\mathbf{C}_{-1}\mathbf{u}_{j-1}$
(47)

where $\mathbf{C}_{-1}$ and $\mathbf{C}_{0}$ are defined similarly to
$\mathbf{C}_{L}$ and $\mathbf{C}_{0,\xi}$ in Eq.(28). Then using Euler’s
method to temporally integrate this equation, we find:

$\mathbf{u}^{n+1}_{j}=\mathbf{u}^{n}_{j}-\tau\big{(}J^{-1}_{j}\mathbf{C}_{0}\mathbf{u}^{n}_{j}+J^{-1}_{j-1}\mathbf{C}_{-1}\mathbf{u}^{n}_{j-1}\big{)}=\big{(}\mathbf{I}-\tau
J^{-1}_{j}\mathbf{C}_{0}\big{)}\mathbf{u}^{n}_{j}-\tau
J^{-1}_{j-1}\mathbf{C}_{-1}\mathbf{u}^{n}_{j-1}$ (48)

where the superscript $n$ denotes the time step. If
$u^{n-m}_{j-m-1}(\forall\>m\in\mathbb{N}$) is then recursively substituted,
the final form is then:

$\mathbf{u}^{n+1}_{j}=\frac{1}{2}\sum^{\infty}_{m=0}(\mathbf{I}-\tau(2\gamma)^{m}\mathbf{C_{0}})(\underbrace{-2\gamma\tau\mathbf{C}_{-1}}_{\mathbf{T}})^{m}\mathbf{u}^{n-m}_{j-m}$
(49)

where we assume the solution is on a geometrically expanding grid in order to
substitute for the Jacobian — hence being only valid for linearly transformed
elements. If we consider that $\mathbf{C}_{0}$ and $\mathbf{C}_{-1}$ are
linear operators, then rather than prescribing a solution, the dynamics of
linear operators can be used Bayart2009 . So, if the mesh extends infinitely
downwind, then it is sufficient to say that Eq.(49) is stable when it is a
hypercylic orbit. Hence, the stability criterion is that
$-2\gamma\tau\mathbf{C}_{-1}=\mathbf{T}$ is a matrixable linear hypercylic
operator. The definition of which is that $\sup{\|\mathbf{T}^{n}\|}\leqslant
1,\forall\>n\in\mathbb{N}$, which in turn implies that $\rho(\mathbf{T})=1$.
What this aims to show is that the stability criterion is dependent on the
product of $\tau$ and $\gamma$, as well as on $\mathbf{C}_{0}$. Therefore, as
$\gamma$ increases the maximum stable $\tau$ decreases, for constant
$\mathbf{C}_{-1}$. This also explains that although the scheme may be formally
unstable, correct setting of $\tau$ for a given $\gamma$ can lead to a
$\mathbf{T}$ that is still hypercyclic and give a bounded solution. However
because of the $(-1)^{m}$ this does mean that if
$\mathbf{U}_{n-1}\subset\mathbb{R}$ is the set of solutions at some time step
$n-1$, then the solution $\mathbf{u}^{n}_{j}\nsubseteq\mathbf{U}_{n-1}$ _i.e_
, for a sinusoidal solution the computed value may exceed the prescribed
magnitude.

###  4.5 Effect of Grid Expansion of Fully Discrete Dispersion and Dissipation

Following on from the exploration of grid expansion on temporal stability
limits, we will present the fully discretised Fourier analysis in two
dimensions. In this investigation, we look at the dispersion and dissipation
of Huynh $g_{2}$ correction functions for several expansion ratios. Throughout
this investigation as the angle is swept through from $0-90^{\circ}$ the time
step will be held constant as this is reflective of practical applications.

(a) Dispersion: $\gamma_{x}=\gamma_{y}=1$

(b) Dissipation: $\gamma_{x}=\gamma_{y}=1$

(c) Dispersion: $\gamma_{x}=1.1$, $\gamma_{y}=1$

(d) Dissipation: $\gamma_{x}=1.1$, $\gamma_{y}=1$

(e) Dispersion: $\gamma_{x}=1.1$, $\gamma_{y}=0.9$

(f) Dissipation: $\gamma_{x}=1.1$, $\gamma_{y}=0.9$

Figure 11: Dispersion and dissipation of upwinded FR, $p=3$, with Huynh
$g_{2}$ corrections and explicit RK44 temporal integration, $\tau=0.18$. The
radial distance is the normalised wavenumber (including the effect of angle),
and the azimuthal distance is the angle of incidence to the element. The solid
black line on the dissipation plots is the contour of zero dissipation.

Figure 11 displays the dispersion and dissipation relation for Huynh $g_{2}$
correction functions when fully discretised with RK44 explicit temporal
integration. It should be noted that there are some anomalous artefacts in the
data due to the complexity of sorting and selecting the eigenvalues. Primarily
this shows that when fully discretised, the instability on expanding grids is
still present. This is then coupled to the same behaviour that was observed by
Vermerie et al. Vermeire2017a  and Trojak et al. Trojak2018c . Namely that the
gradient of the dispersion can be much larger and hence the magnitude group
velocity can be very large. Also observed, which is seen here, is that as the
explicit time step is increased the dissipation is reduced, with the largest
difference seen at higher wavenumbers. As a result the wavenumber at which
$\Im{(\hat{\omega})}=0$ (the solid black lines) in Fig 11(d) & 11(f) is not
appreciably changed compared to the semi-discrete results. However, in Fig 11
there are regions of significantly increased and decreased dissipation at high
wavenumbers. The general impact that these results show is that, when fully
discretised, FR becomes more heterogeneous.

###  4.6 Effect of Correction Functions with Grid Expansion

As was mentioned in section 1, a series of correction functions with peak
temporal stability and spatial accuracy was proposed by Vincent et al.
Vincent2011 , defined by a correction parameter, $\iota_{+}$. These correction
functions exhibited the superconvergence expected of Nodal DG (Cockburn et al.
Cockburn1999 ), however $\iota_{+}\neq\iota_{\mathrm{DG}}$ as they account for
the variation in the dispersion caused by the discrete temporal integration.
We wish to investigate if the advantage of this family of correction functions
is maintained in 2D or if an analogue can be found. Therefore, the correction
function parameter is varied for different angles and grid expansion rate, the
results of which are presented in Fig. 12.

(a) $p=3$, $\gamma_{x}=1.0$ & $\gamma_{y}=1.0$

(b) $p=4$, $\gamma_{x}=0.8$ & $\gamma_{y}=1.0$

(c) $p=4$, $\gamma_{x}=1.0$ & $\gamma_{y}=1.0$

(d) $p=4$, $\gamma_{x}=1.3$ & $\gamma_{y}=1.0$

Figure 12: CFL limit for 2D linear advection, at several orders. Varying
correction function parameter, $\iota$, and angle $\theta$. Time integration
is RK44. The value of $\iota_{+}$ is shown as a dashed black line.

Using Fig. 12(c) as an example, the peak CFL at $\iota_{+}$ can be seen
clearly in the case of $\theta=0^{\circ}$, with its peak value reduced in
comparison to the 1D case. (This was discussed in section 4.4. However, the
clear peak at $\theta=0^{\circ},90^{\circ}$ does not significantly persist as
the wave angle increases to $\theta=45^{\circ}$, with the peak becoming
substantially flattened. Therefore, the balancing effect that the modification
of correction function has on the dispersion of the scheme seems to have a
limited scope. Also, it seems there is no other correction function able to
achieve the same effect at the intermediate range of angles. Furthermore, the
persistence of the peak CFL limit is not seen as the expansion rate and order
is varied, with $\iota_{+}$, in fact, suffering the most appreciable decay in
performance as the angle is varied, when compared to the other correction
functions.

To investigate further the impact of varying correction functions on higher
dimensional problems, we will consider the projection of the solution into the
functional space of FR for linear advection. The method for understanding this
is by studying the posedness of the linear operators. The process of
projection was outlined in section 3. Presented here are the several results
showing how the posedness varies with angle, correction function, and order.

(a) Variation of mode condition number with angle and wavenumber for upwinded
2D FR, $p=2$, with Huynh $g_{2}$ correction functions.

(b) Variation of mode condition number with angle and wavenumber for upwinded
2D FR, $p=2$, with Nodal DG correction functions.

(c) Variation of mode condition number against order, $p$, for upwinded 2D FR,
$\theta=45^{\circ}$, with Huynh $g_{2}$ correction functions.

(d) Variation of mode condition number against order, $p$, for upwinded 1D FR,
with Huynh $g_{2}$ correction functions.

Figure 13: Condition number of various linear FR configurations

Several insights into the different multidimensional behaviour of FR can be
gained by studying Fig. 13. By comparison of Fig. 13(a) & 13(b) it can be seen
that Huynh’s $g_{2}$ correction function causes the projection to be more ill-
posed on average compared to Nodal DG methods. Hence, the superconvergent DG
recovering scheme of Vincent et al. Vincent2011  has decreased temporal
stability. This is because the ill-posedness indicates the sensitivity of the
reconstruction to change, a more sensitive reconstruction means that error can
result in energy being transferred to more dissipative modes. This point is
important as it means to show that ill-conditioning is not directly the
mechanism for loss. But the movement of energy to other more dissipative modes
is, and the condition number is symptomatic of that.

Furthermore, in both cases as, across all wavenumbers, the condition number
can be seen to increase with incident angle. When compared with Fig. 10 and
Fig. 12(c) it can be seen that for a given correction function the CFL limit
reaches its minimum at $\theta=45^{\circ}$, therefore a stark increase in
condition number can also cause a decrease in temporal stability, potentially
due to too much transport between modes.

A result presented by Trojak et al. Trojak2017a  was that for FR the best
points per wavelength (PPW) performance was seen for $p=4$ with the PPW
increasing for orders higher than this. A result that is exhibited by Fig.
13(c) is that the condition number for $p=5$ schemes is higher than for lower
orders and may have passed a point where increased order is out weighted by
inaccuracy in ill-conditioning. This may explain the optimal result seen and
is touching on the fundamental problem characterised by Runge’s phenomena,
that high-order may introduce accuracy through high-order but may also
introduce inaccuracy reflected by a high condition. A result that cannot be
clearly seen in either Fig. 13(a) or 13(b), however, was exhibited by FR was
for $\theta=0$ the condition number was significantly higher than that found
for 1D FR owing to the naturally poor conditioning of a 2D system acting as a
quasi-1D one. This is linked to the lower CFL number experienced in this case,
as was shown in Fig. 10. To show that the results found in this section extend
to higher dimensions, the von Neumann analysis was repeated for 3D ’hexi-
linear’ grids. The primary result of interest is the increased condition
number of the functional projection and this can be seen in Fig.14. The
message is that, as expected, the ill-posedness of the reconstruction
increases with the order. While small increases in the condition number can
give increased temporal stability, larger increases in condition number tend
to act to destabilise the coupled spatial-temporal scheme.

Figure 14: Variation of condition number with wave angles and wavenumber for
upwinded 3D FR, $p=3$, $\gamma_{x}=\gamma_{y}=\gamma_{z}=1$, using Huynh
$g_{2}$ correction functions.

##  5 Non-linear Navier-Stokes Equations with Randomised Grids

It is common within the CFD community to use the canonical Taylor-Green Vortex
(TGV) Taylor1937  test case to assess the numerics of a solver applied to the
Navier-Stokes equations with turbulence — and to that end, there is a plethora
of DNS data available for comparison Brachet1983 ; DeBonis2013 . However, this
case is quite contrived and ultimately will favour spectral or structured
methods due to the Cartesian and periodic domain, whilst also being
unrepresentative of engineering flows that are often wall bounded and/or have
complex geometries. Hence, we propose linearly deforming the elements of the
mesh by jittering the corner nodes to be more representative of real mesh
conditions. Importantly, these deformations will introduce cross
multiplication into the Jacobian, as well as local regions of expansion and
contraction.

The initial conditions of the TGV being used here are those of DeBonis
DeBonis2013 , where the character of the flow is controlled by the non-
dimensional parameters defined as:

$R_{e}=\frac{\rho_{0}U_{0}L}{\mu},\quad\quad P_{r}=0.71=\frac{\mu\gamma
R}{\kappa(\gamma-1)},\quad\quad M_{a}=0.08=\frac{U_{0}}{\sqrt{\gamma RT_{0}}}$
(50)

where we will use the standard set of free-variables for the velocity,
density, pressure, and gas characteristics:

$U_{0}=1,\quad\quad\rho_{0}=1,\quad\quad p_{0}=100,\quad\quad
R=1,\quad\quad\gamma=1.4,\quad\quad L=1$ (51)

Here, due to the solver implementation, we use a specific gas constant of
unity and hence, to achieve the required Reynolds and Prandtl numbers, the
dynamic viscosity and thermal conductivity can be set appropriately.

As has been stated, we will take the uniform periodic mesh on the domain
$\mathbf{\Omega}\in[-\pi,\pi]^{3}$, and jitter the corner nodes of the
elements that are interior to the domain. The amount of jitter is calculated
using a time seeded random number shifted to be centred about zero and scaled
by a global factor between zero and unity. The scaling factor is such that
zero gives a uniform mesh and unity could lead to edges of zero length. After
jittering, the solution points are then linearly positioned within the element
using the thin plate spline radial basis function together with the mapping
from the uniform to jittered corner nodes. This gives a linear mapping of
uniform solution points to solution points within the jittered elements.
Finally, a quality metric is needed to describe, in a single number, the
relative quality of the meshes produced. We opted for a volume ratio shape
factor, slightly redefined as:

$q_{h}=\frac{6\sqrt{\pi}V_{h}}{S_{h}^{3/2}}$ (52)

where $S_{h}$ is the surface area of the hexahedral element and $V_{h}$ is the
volume of the hexahedral elements. The quality metric, $q_{h}$, is then
defined as the ratio of the volume of the element to the volume of a sphere
with the same surface area, with $q_{h}=\sqrt{\pi/6}$ for a perfect cube. To
put this parameter into context, some example meshes are shown in Fig. 15.

(a) $q_{h}=\sqrt{\pi/6}\approx 0.7236$

(b) $q_{h}=0.7201$

(c) $q_{h}=0.7016$

Figure 15: Example slices through a 3D hexahedral mesh to illustrate the mesh
quality metric.

The statistics that will be studied here are the decay of the kinetic energy
and the enstrophy decay rate, which are defined respectively as:

$\displaystyle-\frac{\mathrm{d}E_{k}}{\mathrm{d}t}$
$\displaystyle=-\frac{1}{2\rho_{0}|\mathbf{\Omega}|}\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathbf{\Omega}}\rho(u^{2}+v^{2}+w^{2})\mathrm{d}\mathbf{x}$
(53) $\displaystyle\epsilon$
$\displaystyle=\frac{\mu}{\rho_{0}^{2}|\mathbf{\Omega}|}\int_{\mathbf{\Omega}}\rho(\boldsymbol{\omega}\cdot\boldsymbol{\omega})\mathrm{d}\mathbf{x}$
(54)

where $\boldsymbol{\omega}=\nabla\times[u,v,w]^{T}$ is vorticity and
$|\mathbf{\Omega}|$ is the domain volume.

(a) Selected turbulent kinetic energy dissipation.

(b) Variation of turbulent kinetic energy dissipation with jitter. Dashed
contour at zero dissipation.

Figure 16: Effect of jitter on turbulent kinetic energy dissipation of the TGV
($R_{e}=1600$) for FR, $p=2$, with Huynh $g_{2}$ correction functions on a
$40^{3}$ DoF mesh. Explicit time step size is $\Delta t=1\times 10^{-3}$.

(a) Selected turbulent kinetic energy dissipation.

(b) Variation of turbulent kinetic energy dissipation with jitter. Dashed
contour at zero dissipation.

Figure 17: Effect of jitter on turbulent kinetic energy dissipation of the TGV
($R_{e}=1600$) for FR, $p=4$, with Huynh $g_{2}$ correction functions on a
$120^{3}$ DoF mesh. Explicit time step is $\Delta t=1\times 10^{-3}$.

Fig. 16 & 17 shows the first of these results. First looking at Fig. 16(a) &
17(a) which shows two specific dissipation curves for an uniform and jittered
mesh. At the beginning of the simulation, there is a clear time at which the
global energy increases.

Extending these runs to cover multiple grid qualities, Fig. 16(b) & 17(b), it
is observed that as the grid quality decreases a region where turbulent
kinetic energy increases soon emerges. As time progress, energy dissipation is
again seen and the point of peak dissipation arrives early, moving from
$t\approx 8.5$ to $t\approx 7.5$. The same behaviour is seen for both $p=2$
and $p=4$. From comparison of $p=2$ and $p=4$, it seems that $p=4$ is slightly
more robust to grid deformation, as $p=4$ was able to run at $q_{h}\approx
0.7$. Whereas for $p=2$, $q_{h}$ could not be reduced much below $0.717$ for
$120^{3}$ DoF without completely diverging.

The explanation of this is believed to be that initially the regions in the
mesh that are locally expanding cause an increase in the energy due to the
positive dissipation. This was discussed in Section 4.2, that for the linear
advection equation, dissipation is positive at low wavenumbers for expanding
grids and negative for contracting grids at the same wavenumber. It is thought
that as the simulation progresses, the energy cascade of large scales to small
scales then means that more of the solution lies in the more dissipative
higher wavenumber region for both expanding and contracting grids. Therefore,
the net dissipation at a later time is higher than the uniform case and hence
the peak dissipation is earlier. This is consistent with a lower $R_{e}$ and
hence higher global dissipation. We can conclude that the stability of this
case is brought about by the physics of the Navier-Stokes and cascade of
energy from low to high wavenumbers, which sidesteps the problem of positive
dissipation of low wavenumbers on expanding meshes. This result is interesting
as it is in slight contradiction to the result of Trojak et al. Trojak2017a ,
where FR was found to be more robust to grid stretching than second order FV
for Euler’s equations. However, in that investigation, the Isentropic
Convecting Vortex was considered where there was a large convective velocity.
This may expose slightly different properties, although a full explanation is
not known.

(a) $p=2$

(b) $p=4$

Figure 18: Comparison of TGV enstrophy for $120^{3}$ degree of freedom grid
with similar $q_{h}$.

Studying the effect of jittered grids on enstrophy, shown in Fig. 18, it is
clear that as the grid is stretched the enstrophy increases. This is
indicative of an increase in the vorticity, with the rise occurring within
$t=0-1$. This is consistent with the assertion that the energy is being added
in the small scales, as at this time there are only small scales present.
Furthermore, if the energy were only added to the bulk flow, i.e. $k=0$, it
would not be seen in the enstrophy. After the initial increase, the enstrophy
returns to following the trend of uniform case. However, in the case of $p=2$,
Fig. 18(a), a larger initial increase is seen followed by a wider peak. The
wider peak is similar in character to that of the uniform case and is due to
the grid being mildly under-resolved in the $p=2$ case relative to the DNS.

To provide some reference as to how FR performs relative to an established
method we will use an edge-based Finite Volume (FV) method for comparison. The
FV method is a standard central second order method with L2 Roe smoothing
Osswald2016  for stabilisation, which has been validated previously
Scillitoe2016 . The particular FR scheme used in this comparison is $p=1$,
giving second order, the same as the FV scheme. However, this puts FR at a
significant disadvantage as its numeric characteristics at low order are
particularly poor. For example, consider the dispersion and dissipation
relations in Fig. 19, which, by comparison to the result of Lele Lele1992 ,
show that FR has noticeably lower resolving abilities when compared against a
second order FD scheme.

(a) Dispersion

(b) Dissipation

Figure 19: Dispersion and Dissipation relation for 1D upwinded FR, $p=1$, with
DG correction function.

With this in mind, we present the results of tests on various jittered grids
with a total of $170^{3}$ degrees of freedom in Fig. 20. For the uniform case,
the enstrophy clearly shows that FR is underresolved compared to FV, which is
also shown by a slightly increased rate of dissipation earlier — indicating
that the implicit filter is too narrow. If we now consider the effect of
jittering, several things may be concluded. Firstly, we were unable to run FR
with $j_{f}=0.5$ as the simulation quickly became unstable. Secondly, for
$-dE_{k}/dt$ it seems that the peak value is less sensitive with FR than with
FV, with central FV seeing some very large amplitude oscillation in
$-dE_{k}/dt$. This is likely to be rooted in the central differencing at the
interfaces, as if we change to a kinetic energy preserving formulation
Jameson2008 ; WatsonTucker2015 , as is displayed in Fig. 21, these
oscillations are removed and the sensitivity to jitter is reduced. The
enstrophy (Fig. 21(b)) seems to indicate that a large amount of what seemed to
be resolved energy may have in fact been dispersion induced fluctuations.
However, in both cases FV was able to run with grids up to $j_{f}=0.9$ and
$q_{h}=0.6382$ — not shown — it appears that in these cases the added
stability of the smoothing has greatly helped FV, especially in the central
difference case where running without smoothing caused the case to fail even
at low levels of jitter.

(a) Kinetic energy dissipation

(b) Enstrophy based dissipation

Figure 20: Comparison of FR, $p=1$ with DG correction functions with a second
order central FV scheme with L2 Roe smoothing both with $170^{3}$ degrees of
freedom and $\Delta t\approx 5\times 10^{-4}$. A reference DNS solution is
provided from Brachet1983 .

Before moving on, it must be noted that in FR we once again see a dip in
$dE_{k}/dt$, which is also present to a lesser extent in both versions of FV
tested. This is again linked to the instability caused by locally expanding
grids, but in the case of FV the dip is smaller and is aided by the use of
smoothing. The conclusion for this comparison is that FV is somewhat resilient
to degradation to the mesh quality, with the resilience coming from smoothing
for central differenced FV. This allowed highly warped meshes to be run, but
at the expense of accuracy with excess dissipation affecting the solution. KEP
was found to be far more resilient and could even run without smoothing. FR,
when run at a low order, was unsuited to this problem, but did see less
degradation compared to central FV and there is the potential for $p=1$ FR to
equally benefit on poorer quality meshes from smoothing via a different
Riemann solver, such as that of Roe Roe1981 .

The result shown in this section clearly highlight a key issue of FR, that it
works on the strong conservative form of the equation. Consequently, it relies
on the accurate interpolation and accurate calculation of gradients to make
sure that the scheme does not stray far from conservation. However, this
implementation does not strictly enforce conservation and hence when
inaccuracy in the gradient or interpolation is introduced there is no remedial
action taken. Some recent work that aims to fix this issue, see Abe et al.
Abe2015 ; Abe2016 ; Abe2018 , where the common interface values are used to
enforce conservation. Following from the results of the present work, these
extensions to FR should be considered essential.

(a) Kinetic energy dissipation

(b) Enstrophy based dissipation

Figure 21: Comparison of FR, $p=1$ with DG correction functions with a second
order KEP FV scheme with L2 Roe smoothing both with $170^{3}$ degrees of
freedom and $\Delta t\approx 5\times 10^{-4}$. A reference DNS solution is
provide from Brachet1983 .

##  6 Conclusions

Through this work, we have presented a theoretical extension of the FR von
Neumann analysis to higher dimensions. This allowed as to understand the
character of the dispersion and dissipation relations of FR as the incident
angle of a wave was varied. Differences were noted between the behaviours of
FR and finite differencing methods, primarily that FR saw a lower variation in
character with the angle of incidence. The effect of higher dimensionality on
the CFL limit was also found, with higher dimensionality causing a reduction
in the CFL limit.

Investigations were then performed on deformed meshes and theoretically, the
same behaviour was seen in two dimensions as in one. Specifically, that
expanding meshes cause instability and contracting meshes cause excess
dissipation, however, when coupled together the effects can act to cancel each
other out. Numerical experiments were subsequently performed using the Taylor-
Green vortex case, but with the element corner node positions jittered. Tests
showed $5^{\mathrm{th}}$ order to be more resilient to poor quality meshes
than $3^{\mathrm{rd}}$ order, but in both cases, the effect of localised
regions of expansion are thought to be responsible for an initial increase in
the kinetic energy of the solution. The appearance of smaller turbulent scales
within the TGV solution as time progressed then counteracted this effect, as
high wavenumbers on locally contracting regions experience excess dissipation.
Lastly, a comparison was made between FR and a second order FV method. It was
found that FR was more resilient to mesh deformation that FV methods, however,
FR is far from optimal when running at second order. In both cases, it
recommended that kinetic energy presenting/conservation methods should be used
as they will greatly increase resilience to mesh quality.

## Acknowledgements

The support of the Engineering and Physical Sciences Research Council of the
United Kingdom is gratefully acknowledged under the award reference 1750012.
This work was performed using resources provided by the Cambridge Service for
Data Driven Discovery (CSD3) operated by the University of Cambridge Research
Computing Service (http://www.csd3.cam.ac.uk/), provided by Dell EMC and Intel
using Tier-2 funding from the Engineering and Physical Sciences Research
Council (capital grant EP/P020259/1), and DiRAC funding from the Science and
Technology Facilities Council (www.dirac.ac.uk).

## References

  * (1) Abe, Y., Haga, T., Nonomura, T., Fujii, K.: On the freestream preservation of high-order conservative flux-reconstruction schemes.  Journal of Computational Physics 281, 28–54 (2015).  DOI 10.1016/j.jcp.2014.10.011.  URL http://dx.doi.org/10.1016/j.jcp.2014.10.011
  * (2) Abe, Y., Haga, T., Nonomura, T., Fujii, K.: Conservative high-order flux-reconstruction schemes on moving and deforming grids.  Computers and Fluids 139, 2–16 (2016).  DOI 10.1016/j.compfluid.2016.03.024.  URL http://dx.doi.org/10.1016/j.compfluid.2016.03.024
  * (3) Abe, Y., Morinaka, I., Haga, T., Nonomura, T., Shibata, H., Miyaji, K.: Stable, non-dissipative, and conservative flux-reconstruction schemes in split forms.  Journal of Computational Physics 353, 193–227 (2018).  DOI 10.1016/j.jcp.2017.10.007.  URL https://doi.org/10.1016/j.jcp.2017.10.007
  * (4) Asthana, K., Watkins, J., Jameson, A.: On Consistency and Rate of Convergence of Flux Reconstruction For Time-Dependent Problems.  Journal of Computational Physics 334, 367–391 (2017).  DOI 10.1016/j.jcp.2017.01.008.  URL http://linkinghub.elsevier.com/retrieve/pii/S0021999117300189
  * (5) Bayart, F., Matheron, E.: Dynamics of Linear Operators, first edn.  Cambridge University Press, New York (2009)
  * (6) Brachet, M.E., Orszag, S.A., Nickel, B.G., Morf, R.H., Frisch, U.: Small-Scale Structure of the Taylor-Green Vortex.  Journal of Fluid Mechanics 130, 411–452 (1983)
  * (7) Brandvik, T., Pullan, G.: An Accelerated 3D Navier–Stokes Solver for Flows in Turbomachines.  Journal of Turbomachinery 133(April 2011), 1–9 (2009).  DOI 10.1115/1.4001192
  * (8) Brenner, S.C., Ridgway-Scott, L.: The Mathematical Theory of Finite Element Methods, third edn.  Springer-Verlag, New York (2008)
  * (9) Castonguay, P.: High-order Energy Stable Flux Reconstruction Schemes For Fluid Flow Simulations on Unstructured Grids.  Ph.D. thesis, Stanford University (2012)
  * (10) Castonguay, P., Williams, D.M., Vincent, P.E., Jameson, A.: Energy Stable Flux Reconstruction Schemes for Advection-Diffusion Problems.  Computer Methods in Applied Mechanics and Engineering 267(1), 400–417 (2013).  DOI 10.1007/s10915-013-9780-2.  URL http://dx.doi.org/10.1016/j.cma.2013.08.012
  * (11) Chow, F.K., Moin, P.: A Further Study of Numerical Errors in Large-Eddy Simulations.  Journal of Computational Physics 184(2), 366–380 (2003).  DOI 10.1016/S0021-9991(02)00020-7
  * (12) Cockburn, B., Luskin, M., Shu, C.W., Suli, E.: Post-Processing of Galerkin Methods for Hyperbolic Problems.  In: Discontinuous Galerkin Methods (Newport, RI, 1999), pp. 291–300. Springer-Verlag, Newport (1999)
  * (13) Courant, R., Friedrichs, K., Lewy, H.: On the Partial Differnce Equations of Mathematical Physics.  IBM Journal 11(2), 215–234 (1967).  DOI 10.1147/rd.112.0215
  * (14) DeBonis, J.: Solutions of the Taylor-Green Vortex Problem Using High-Resolution Explicit Finite Difference Methods.  51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition (February 2013) (2013).  DOI 10.2514/6.2013-382.  URL http://arc.aiaa.org/doi/10.2514/6.2013-382
  * (15) Ghosal, S.: An Analysis of Numerical Errors in Large-Eddy Simulations of Turbulence.  Journal of Computational Physics 125(1), 187–206 (1996).  DOI 10.1006/jcph.1996.0088.  URL http://www.sciencedirect.com/science/article/pii/S0021999196900881
  * (16) Harten, A., Lax, P.D., van Leer, B., Vanleer, B.: On upstream differencing and Godunov-type schemes for hyperbolic conservation laws.  SIAM Review 25(1), 35–61 (1983).  DOI 10.1137/1025002
  * (17) Hesthaven, J.S., Kirby, R.M.: Filtering in Legendre spectral methods.  Mathematics of Computation 77(263), 1425–1452 (2008).  DOI 10.1090/S0025-5718-08-02110-8.  URL http://www.ams.org/jourcgi/jour-getitem?pii=S0025-5718-08-02110-8
  * (18) Hesthaven, J.S., Warburton, T.: Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, _Texts in Applied Mathematics_ , vol. 54, first edn.  Springer New York, New York, NY (2008).  DOI 10.1007/978-0-387-72067-8.  URL http://link.springer.com/10.1007/978-0-387-72067-8
  * (19) Huynh, H.T.: A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods.  In: 18th AIAA Computational Fluid Dynamics Conference, vol. 2007-4079, pp. 1–42 (2007).  DOI 10.2514/6.2007-4079.  URL http://arc.aiaa.org/doi/pdf/10.2514/6.2007-4079
  * (20) Isaacson, E., Keller, H.B.: Analysis of Numerical Methods, second edn.  John Wilery & Sons Ltd., New York (1994).  DOI 10.2307/2003280.  URL http://books.google.pl/books?id=y77n2ySMJHUC
  * (21) Jameson, A.: Formulation of kinetic energy preserving conservative schemes for gas dynamics and direct numerical simulation of one-dimensional viscous compressible flow in a shock tube using entropy and kinetic energy preserving schemes.  Journal of Scientific Computing 34(2), 188–208 (2008).  DOI 10.1007/s10915-007-9172-6
  * (22) Jameson, A., Vincent, P.E., Castonguay, P.: On The Non-Linear Stability of Flux Reconstruction Schemes.  Journal of Scientific Computing 50(2), 434–445 (2012).  DOI 10.1007/s10915-011-9490-6
  * (23) Kress, R.: Numerical Analysis, _Graduate Texts in Mathematics_ , vol. 181, first edn.  Springer New York, New York, NY (1998).  DOI 10.1007/978-1-4612-0599-9.  URL http://link.springer.com/10.1007/978-1-4612-0599-9
  * (24) Lax, P., Richtmyer, R.: Survey of the Stability of Linear Finite Difference Equations.  Communications on Pure and Applied Mathematics IX(1), 267–293 (1956).  DOI 10.1002/cpa.3160090206
  * (25) van Leer, B.: Flux-Vector Splitting for the Euler Equation.  Tech. Rep. February (1982).  DOI 10.1007/3-540-11948-5
  * (26) Lele, S.K.: Compact Finite Difference Schemes With Spectral-Like Resolution.  Journal of Computational Physics 103(1), 16–42 (1992).  DOI 10.1016/0021-9991(92)90324-R
  * (27) Mengaldo, G., De Grazia, D., Vincent, P.E., Sherwin, S.J.: On the Connections Between Discontinuous Galerkin and Flux Reconstruction Schemes: Extension to Curvilinear Meshes.  Journal of Scientific Computing 67(3), 1272–1292 (2016).  DOI 10.1007/s10915-015-0119-z
  * (28) Oßwald, K., Siegmund, A., Birken, P., Hannemann, V., Meister, A.: L2 Roe: a low dissipation version of Roe’s approximate Riemann solver for low Mach numbers.  International Journal for Numerical Methods in Fluids 81(2), 71–86 (2016).  DOI 10.1002/fld.4175.  URL http://doi.wiley.com/10.1002/fld.4175
  * (29) Reed, W.H., Hill, T.R.: Triangular Mesh Methods for the Neutron Transport Equation.  Los Alamos Report LA-UR-73-479 (836), 10 (1973)
  * (30) Roe, P.L.: Approximate Riemann Solvers , Parameter Vectors, and Difference Schemes.  Journal of Computational Physics 43, 357–372 (1981)
  * (31) Scillitoe, A.D., Tucker, P.G., Adami, P.: Numerical Investigation of Three-Dimensional Separation in an Axial Flow Compressor: The Influence of Free-Stream Turbulence Intensity and Endwall Boundary Layer State.  In: Volume 2D: Turbomachinery, p. V02DT44A025. ASME, Seoul, South Korea (2016).  DOI 10.1115/GT2016-57241.  URL http://proceedings.asmedigitalcollection.asme.org/proceeding.aspx?doi=10.1115/GT2016-57241
  * (32) Sheshadri, A.: An Analysis of Stability of the Flux Reconstruction Formulation With Applications to Shock Capturing.  Phd thesis, Stanford University (2016)
  * (33) Sheshadri, A., Jameson, A.: On the Stability of the Flux Reconstruction Schemes on Quadrilateral Elements for the Linear Advection Equation.  Journal of Scientific Computing 67(2), 769–790 (2016).  DOI 10.1007/s10915-015-0102-8
  * (34) Taylor, G.I., Green, A.E.: Mechanism of the Production of Small Eddies from Large Ones.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 158(895), 499–521 (1937).  DOI 10.1098/rspa.1937.0036.  URL http://rspa.royalsocietypublishing.org/cgi/doi/10.1098/rspa.1937.0036
  * (35) Trojak, W.: Generalised Lebesgue Stable Flux Reconstruction.  Arxiv 1805(12481v2), 1–15 (2018).  URL http://arxiv.org/abs/1805.12481
  * (36) Trojak, W.: Generalised Sobolev Stable Flux Reconstruction.  Arxiv 1804(04714), 1–19 (2018).  URL http://arxiv.org/abs/1804.04714
  * (37) Trojak, W., Watson, R., Tucker, P.G.: High Order Flux Reconstruction on Stretched and Warped Meshes.  In: AIAA SciTech, 55th AIAA Aerospace Sciences Meeting, pp. 1–12. Grapevine Texas (2017)
  * (38) Trojak, W., Watson, R., Tucker, P.G.: Temporal Stabilisation of Flux Reconstruction on Linear Problems.  In: 2018 Fluid Dynamics Conference, AIAA AVIATION Forum. American Institute of Aeronautics and Astronautics, Atlanta (2018).  DOI doi:10.2514/6.2018-4263.  URL https://doi.org/10.2514/6.2018-4263
  * (39) Trojak, W., Witherden, F.D.: A New Family of Weighted One-Parameter Flux Reconstruction Schemes pp. 1–37 (2018)
  * (40) Vermeire, B.C., Vincent, P.E.: On the Behaviour of Fully-Discrete Flux Reconstruction Schemes.  Journal of Computational Physics 315(1), 1053–1079 (2017).  DOI 10.1016/j.jcp.2016.09.034.  URL http://linkinghub.elsevier.com/retrieve/pii/S0021999116304557
  * (41) Vincent, P.E., Castonguay, P., Jameson, A.: A New Class of High-Order Energy Stable Flux Reconstruction Schemes.  Journal of Scientific Computing 47(1), 50–72 (2010).  DOI 10.1007/s10915-010-9420-z.  URL http://link.springer.com/10.1007/s10915-010-9420-z
  * (42) Vincent, P.E., Castonguay, P., Jameson, A.: Insights From von Neumann Analysis of High-Order Flux Reconstruction Schemes.  Journal of Computational Physics 230(22), 8134–8154 (2011).  DOI 10.1016/j.jcp.2011.07.013.  URL http://dx.doi.org/10.1016/j.jcp.2011.07.013
  * (43) Vincent, P.E., Farrington, A.M., Witherden, F.D., Jameson, A.: An Extended Range of Stable-Symmetric-Conservative Flux Reconstruction Correction Functions.  Computer Methods in Applied Mechanics and Engineering 296, 248–272 (2015).  DOI 10.1016/j.cma.2015.07.023.  URL http://dx.doi.org/10.1016/j.cma.2015.07.023
  * (44) Vincent, P.E., Witherden, F.D., Vermeire, B.C., Park, J.S., Iyer, A.: Towards Green Aviation with Python at Petascale.  In: International Conference for High Performance Computing, Networking, Storage and Analysis, SC, November, pp. 1–11 (2017).  DOI 10.1109/SC.2016.1
  * (45) Wang, Z.: Spectral (Finite) Volume Method for Conservation Laws on Unstructured Grids. Basic Formulation: Basic Formulation.  Journal of Computational Physics 178(1), 210–251 (2002).  DOI 10.1006/jcph.2002.7041.  URL http://www.sciencedirect.com/science/article/pii/S0021999102970415{%}5Cnhttp://pdn.sciencedirect.com/science?{_}ob=MiamiImageURL{&}{_}cid=272570{&}{_}user=929460{&}{_}pii=S0021999102970415{&}{_}check=y{&}{_}origin=article{&}{_}zone=toolbar{&}{_}coverDate=01-May-2002{&}view=c{&}originContent
  * (46) Watson, R., Tucker, P.G.: Optimising Flux Reconstruction Schemes for Large Eddy Simulation.  In: AMMCS-CAIMS, pp. 1–7 (2015)
  * (47) Williams, D.M., Jameson, A.: Energy stable flux reconstruction schemes for advection-diffusion problems on tetrahedra.  Journal of Scientific Computing 59(3), 721–759 (2014).  DOI 10.1007/s10915-013-9780-2.  URL http://dx.doi.org/10.1016/j.jcp.2013.05.007
  * (48) Yu, M., Wang, Z.J., Liu, Y.: On the Accuracy and Efficiency of Discontinuous Galerkin , Spectral Difference and Correction Procedure via Reconstruction Methods.  Journal of Computational Physics 259, 70–95 (2014).  DOI 10.1016/j.jcp.2013.11.023.  URL http://dx.doi.org/10.1016/j.jcp.2013.11.023

##  Appendix A Nomenclature

Roman
---
$a$ | convective velocity in x
$b$ | convective velocity in y
$c(k)$ | wavespeed at wavenumber $k$
$\mathbf{C}_{0\xi}$ | centre cell FR matrix in $\xi$
$\mathbf{C}_{0\eta}$ | centre cell FR matrix in $\eta$
$\mathbf{C}_{L}$ | left cell FR matrix
$\mathbf{C}_{R}$ | right cell FR matrix
$\mathbf{C}_{B}$ | bottom cell FR matrix
$\mathbf{C}_{T}$ | top cell FR matrix
$\mathbf{D}_{\xi}$ |  $\xi$ first derivative matrix
$\mathbf{D}_{\eta}$ |  $\eta$ first derivative matrix
$h_{L}\>\mathrm{\&}\>h_{R}$ | left and right correction functions
$h_{B}\>\mathrm{\&}\>h_{T}$ | bottom and top correction functions
$k_{nq}$ | solution point Nyquist wavenumber, $(p+1)/\delta_{j}$
$\hat{k}$ |  $k_{nq}$ normalised wavenumber, $[0,\pi]$
$l_{k}$ |  $k^{\mathrm{th}}$ Lagrange polynomial
$p$ | solution polynomial order
$q_{h}$ | element shape factor
$\mathbf{Q}$ | spatial scheme matrix
$\mathbf{R}$ | update matrix
$u$ | primitive in real domain
Greek
$\gamma$ | grid geometric expansion factor
$\delta_{j}$ | mesh spacing
$\eta$ |  $2^{\mathrm{nd}}$ computational dimension variable
$\iota$ | VCJH scheme correction function variable
$\iota_{+}$ | variable $\iota$ for peak temporal stability
$\kappa(\mathbf{A})$ | condition number of matrix $\mathbf{A}$
$\xi$ |  $1^{\mathrm{st}}$ computational dimension variable
$\rho(\mathbf{A})$ | spectral radius of matrix $\mathbf{A}$
$\tau$ | time step
$\boldsymbol{\Omega}$ | solution domain
$\boldsymbol{\Omega}_{n}$ |  $n^{\mathrm{th}}$ solution sub-domain
$\hat{\boldsymbol{\Omega}}$ | standardised sub-domain
Subscript
$\mathrm{\bullet}_{B}$ | variable at bottom of cell
$\mathrm{\bullet}_{L}$ | variable at left of cell
$\mathrm{\bullet}_{R}$ | variable at right of cell
$\mathrm{\bullet}_{T}$ | variable at top of cell
Superscript
$\mathrm{\bullet}^{T}$ | transpose
$\mathrm{\bullet}^{\delta}$ | local polynomial fit of value
$\mathrm{\bullet}^{\delta C}$ | correction to value
$\mathrm{\bullet}^{\delta D}$ | localised discontinuous polynomial fit of value
$\mathrm{\bullet}^{\delta I}$ | common interface value based on local polynomial fit of value
$\hat{\mathrm{\bullet}}$ | transformed variable
$\overline{\mathrm{\bullet}}$ | locally fitted polynomial of variable
Other
$\hat{\nabla}$ | gradient operator in computational domain
$\Im{(z)}$ | imagine part of $z$ given $z\in\mathbb{C}$
$\Re{(z)}$ | real part of $z$ given $z\in\mathbb{C}$
$\mathbb{C}$ | set of complex number
$\mathbb{R}$ | set of real numbers