YEAR-NUMBER Conference Name, Date, and Location YEAR

# High Order Flux Reconstruction on Stretched and Warped Meshes

Will Trojak PhD Candidate, Department of Engineering, University of Cambridge,
AIAA Student Member. Rob Watson and Paul G. Tucker 0
Department of Engineering, University of Cambridge, Cambridge, UK, CB2 1PZ
Research Fellow, Department of Engineering, University of Cambridge, AIAA
MemberProfessor, Department of Engineering, University of Cambridge, AIAA
Assoc. Fellow

###### Abstract

High-order CFD is gathering a broadening interest as a future industrial tool,
with one such approach being Flux Reconstruction (FR). However, due to the
need to mesh complex geometries if FR is to displace current, lower order
methods, FR will likely have to be applied to stretched and warped meshes.
Therefore, it is proposed that the analytical and numerical behaviour of FR on
deformed meshes for both the 1D linear advection and the 2D Euler equations is
investigated. The analytical foundation of this work is based on a modified
von Neumann analysis for linearly deformed grids that is presented. The
temporal stability limits for linear advection on such grids are also explored
analytically and numerically, with CFL limits set out for several Runge-Kutta
schemes, with the primary trend being that contracting mesh regions give rise
to higher CFL limits whereas expansion leads to lower CFL limits. Lastly, the
benchmarks of FR are compared to finite difference and finite volumes schemes,
as are common in industry, with the comparison showing the increased wave
propagating ability on warped and stretched meshes, and hence, FR;s increased
resilience to mesh deformation.

## Nomenclature

Roman
---
$a$ | convective velocity
$c(k)$ | phase velocity at wavenumber $k$
$\mathbf{C}_{0}$ | centre cell FR matrix
$\mathbf{C}_{-1}$ | upwind cell FR matrix
$\mathbf{D}$ | first derivative matrix
$f$ | flux variable in physical domain
$\hat{G}(\hat{k})$ | computational filter kernel
$h_{l}\>\mathrm{\&}\>h_{r}$ | left and right correction functions
$J_{n}$ |  $n^{\mathrm{th}}$ cell Jacobian
$k$ | wavenumber
$k_{nq}$ | solution point Nyquist wavenumber, $(p+1)/\delta_{j}$
$\hat{k}$ |  $k_{nq}$ normalised wavenumber, $[0,\pi]$
$l_{n}$ |  $n^{\mathrm{th}}$ Lagrange basis function
$p$ | solution polynomial order
PPW | points per wavenumber at given error level
$\mathbf{Q}$ | spatial scheme matrix
$\mathbf{R}$ | update matrix
$u$ | primitive in physical domain
$\tilde{u}$ | Fourier reconstructed field
Greek
$\gamma$ | grid geometric expansion factor
$\delta_{j}$ | mesh spacing, $x_{j}-x_{j-1}$
$\epsilon$ | PPW error level
$\xi$ | transformed spatial variable
$\rho(\mathbf{A})$ | spectral radius of $\mathbf{A}$
$\tau$ | time step
$\boldsymbol{\Omega}$ | solution domain
$\boldsymbol{\Omega}_{n}$ |  $n^{\mathrm{th}}$ solution sub-domain
$\boldsymbol{\Omega}_{s}$ | standardised sub-domain
Subscript
$\mathrm{\bullet}_{l}$ | variable at left of cell
$\mathrm{\bullet}_{r}$ | variable at right of cell
Superscript
$\mathrm{\bullet}^{c}$ | common value at interface
$\mathrm{\bullet}^{T}$ | vector or matrix transpose
$\mathrm{\bullet}^{\delta}$ | discontinuous value
$\hat{\mathrm{\bullet}}$ | transformed variable
$\overline{\mathrm{\bullet}}$ | locally fitted polynomial of variable

##  1 Introduction

The potential of Large Eddy Simulation (LES) has been understood for some
time, however most current industrial LES and Reynolds-Averaged Navier Stokes
(RANS) implementations make use of second order spatial schemes. This kind of
lower spatial order for LES can be prohibitively expensive in some flow
regimes, with cost scaling, from Piomelli [1], with $\sim Re^{2.4}$ for the
innermost section of the boundary layer. Figure 1 aims to show that for a gas
turbine engine conventional LES has a much larger overhead than hybrid RANS-
LES, with full LES only really being currently feasible in the lower pressure
turbine stages. From metrics shown later, using FR could make conventional LES
far cheaper. The high cost of wall-resolved LES originates from all but 10-20%
of the vortical motions being directly simulated in both space and time,
Tucker [2]. The unresolved fraction of the vortical motion must, therefore, be
accounted for using a sub-grid scale model. As LES will make use of finite
discretisations on the underlying equation (be that through a Finite Volume
(FV), Discontinuous Galerkin (DG) or another method) the solution will incur
truncation and aliasing error. Therefore, for the sub-grid scale model to
correctly influence the flow, the sub-grid-scale error must be minimised. This
area was explored by Chow & Moin [3] and Ghosal [4], who found that for a
second order central FD LES the filter width was needed to be four times the
grid spacing. They also found that for an eighth order central FD LES a filter
width of twice the grid spacing was required. This effect occurs because for
second order schemes the numerical error and the sub-grid-scale forcing will
have the same order scaling with respect to the grid spacing, whereas a move
to higher order means the differencing error order increases and so the
filter-scale can decrease. Regardless, higher orders allow for coarser meshes
to capture the same flow regime, and can reduce the high computational barrier
to the use of LES. Hence high order methods have become of increasing interest
for industry. For example, studying the abstracts of AIAA SciTech Conferences
on Aerospace Sciences shows an increase from 62 to 99 papers that used or
studied high order methods between 2014 and 2015.

Figure 1: LES and hybrid RANS-LES mesh degrees of freedom requirement against
axial distance through a gas turbine engine. $Re$ is also plotted, showing
large variations in flow regimes occur. (Tucker et al. [5]) This shows that
with current methods LES is only routinely used for the low pressure turbines.
Data presented later can give an estimate of the 3D mesh requirement for
$4^{th}$ order Flux Reconstruction (FR).

Huynh [6] introduced a high order scheme which he termed Flux Reconstruction
(FR), and this has since evolved into a wider family of schemes explored by
Vincent et al. [7] using energy stability arguments to define a wide family of
correction functions. Extensions have been made to handle advection-diffusion
[8, 9], simplex and hypercube elements[10, 11, 12], and various optimisations
of the collocation points and correction functions [13, 14, 15]. The
quantitative evaluation performed analytically and numerically has largely
been confined to canonical meshes and the scenarios explored, although of
great importance, are often of less direct relevance to CFD practitioners.
This is due to their contrived nature leading to well presented problems, that
are useful as benchmarks, but do not always show the full picture. In essence,
it is common for the actual problems encountered by engineers to be highly
complex, making meshes with low levels of skew or low level inter-element
expansion very difficult to produce. Furthermore, FR is rooted in finite
element methods and consequently is an unstructured approach that
fundamentally enables far more irregularities within the mesh compared to
structured approaches. This highlights the importance of accurate
characterisation of the scheme on irregular meshes.

Previous investigations into the analytical behaviour of numerical methods on
irregular grids has been confined to schemes such as finite difference(FD),
where the mesh transformation can be more easily applied. For example Chung &
Tucker [16] in which the effect of hyperbolically transformed grids was
investigated for FD and compact FD schemes, clearly demonstrated the added
dissipation that transformation can cause. More recently You et al. [17]
demonstrated a more thorough approach, detailing the exact higher order terms
which lead to inaccuracy on hyperbolically skewed grids. Within the field of
FR, some numerical results have been presented pertaining to the accuracy of
FR applied to curvelinear grids, Mengaldo et al. [18], also showing some
interesting parallels between FR and Discontinuous Galerkin methods. This did,
however, limit the quantification to numerical experiments, which, although
very powerful, do not allow the bedrock of the scheme to be exposed. Further
investigation into the effect of curvilinear grids on the Jacobian was
presented by Kopriva [19]. Lastly, FR coupled to _r_ -type mesh adpatation was
performed by Sheshadri et al. [20] for the purposes of shock capture using a
’divide or conquer’ method that showed reasonable performance. These works
have laid a fine basis for the development of study into the effect of mesh
transformation on finite element schemes and in this case for FR.

The aims of this paper are broken down into several sections. An analytical
framework is constructed using the established von Neumann analysis, however
generality of the grid will be maintained such that non-uniform grids may be
investigated. The subsequent analytical investigations are concerned with the
dispersion and dissipation of various orders of FR for non-uniform 1D grids,
and with moving on to couple the spatial scheme to a temporal scheme to
establish theoretical CFL limits for various grids. A numerical methodology is
proposed such that the analytical results may validated, also enabling
evaluation of the fully discretised spatial-temporal scheme. Lastly the
numerical investigation is taken further, into 2D, using the Euler equations
to understand the behaviour of more relevant flows on meshes that have been
artificially degraded. This test case can also enable comparison to be made to
a more prevalent industrial FV scheme.

##  2 Flux Reconstruction

Flux Reconstruction [21, 6] (FR) applied to the linear advection equation will
form the basis of the initial investigation to be carried out, and, for the
readers’ convenience, an overview of the scheme is presented here. However,
for a more detailed understanding the reader should consult Castonguay [21] or
Huynh [6]. This 1D scheme can be readily converted to two dimensions (carried
out later) 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 [22] and
utilises the same subdivision of the domain into discontinuous sub-domains:

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

In the standardised sub-domain, $\boldsymbol{\Omega}_{s}\in\mathbb{R}^{d}$,
computational spatial variables are defined. When $d=1$,
$\boldsymbol{\Omega}_{s}=[-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 are placed at the edges of the sub-domain. (Figure
2(a) shows a 1D example). To transform from
$\boldsymbol{\Omega}_{n}\rightarrow\boldsymbol{\Omega}_{s}$ the Jacobian
$J_{n}$ is defined such that

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

and the solution point mapping in the physical domain used throughout this
investigation results in the Jacobian representing a linear mapping. With the
domain set up, we now proceed with defining the steps to construct a
continuous solution from discontinuous segments. The first stage is to define
a local solution polynomial in $\boldsymbol{\Omega}_{s}$ using Lagrange
interpolation.

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

Repeating the interpolation for the discontinuous flux in
$\boldsymbol{\Omega}_{s}$:

$\hat{f}^{\delta}=\hat{f}^{\delta}(\xi;t)=\sum_{i=1}^{p+1}{\hat{f}^{\delta}_{i}l_{i}(\xi)}$
(6)

Now using the Jacobian and the solution polynomials, the primitive and flux
values can be approximated in the physical domain $\boldsymbol{\Omega}_{n}$

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

This distinction is made due to the potentially approximate nature of $J_{n}$.
However for the case of linear transformations, as will be considered here,
$J_{n}$ fully captures the spatial transformation and thus
$\overline{u}^{\delta}$ will be a polynomial of order $p$. With a polynomial
formed in the reference domain the value of the primitive at the flux points,
see Fig. 2(a), can be defined as $\hat{u}^{\delta}_{l}=\hat{u}^{\delta}(-1)$
and $\hat{u}^{\delta}_{r}=\hat{u}^{\delta}(1)$, this can be repeated for the
flux values.

(a) Flux and solution point layout for $p=3$ in $\boldsymbol{\Omega}_{s}$,
with corresponding left and right Huynh, $g_{2}$, correction functions [6].

(b) Discontinuous primitive polynomials and the interpolated values at the
flux points for adjacent sub-domains.

Figure 2: Point layout in $\boldsymbol{\Omega}_{s}$ for $p=3$ and cell
interface topology.

Once the primitive values at the interface, $I$, have been interpolated, a
common interface flux can be calculated, $f^{c}_{I}$, in the physical domain.
For a general case this is done using a Riemann solver on the primitives at
the interface, such as: Roe [23]; flux-vector splitting [24]; or HLL [25]. In
order to get a spatially continuous solution over $\boldsymbol{\Omega}$ the
common interface flux must be incorporated into the solution. For FR this is
done by using a correction function to propagate the corrected flux gradient
into the $\boldsymbol{\Omega}_{n}$. Figure 2(a) shows the correction function
proposed by Huynh [6] and, in general, the left and right correction functions
are defined as $h_{l}(\xi)$ and $h_{r}(\xi)$. The procedure to apply the
correction to transformed discontinuous flux in $\boldsymbol{\Omega}_{s}$ is:

$\frac{\partial\hat{f}(\xi_{i})}{\partial\xi}=\sum_{j=1}^{p+1}{\hat{f}_{j}^{\delta}\frac{dl_{j}(\xi_{i})}{d\xi}}+(\hat{f}^{c}_{l}-\hat{f}^{\delta}_{l})\frac{dh_{l}(\xi_{i})}{d\xi}+(\hat{f}^{c}_{r}-\hat{f}^{\delta}_{r})\frac{dh_{r}(\xi_{i})}{d\xi}$
(8)

This then allows then transformed continuous equation can be written as:

$\frac{\partial\hat{u}(\xi_{i})}{\partial
t}=-\frac{\partial\hat{f}(\xi_{i})}{\partial\xi}$ (9)

and hence the solution can be advanced in time via some method of temporal
integration.

##  3 von Neumann Analysis

The main analysis to be carried out is a von Neumann analysis that follows the
work of Lele [26], Hesthaven et al. [27], Huynh [6], Vincent et al. [7], and
Asthana et al. [14]. The analysis shown here differs from that of previous
work in that it does not assume that the grid is uniform and instead retains
generality, the advantage being that the effect of grid stretching can be
investigated. For the analysis, consider the computational domain as before,
$\boldsymbol{\Omega}\in\mathbb{R}^{1}$ with $N-1$ sub-domains
$\boldsymbol{\Omega}_{n}$, with flux points of $\boldsymbol{\Omega}_{n}$
located at $\mathbf{x}_{j}\>\forall\>\\{j\in\mathbb{N},\>j\leqslant N\\}$. The
1D linear advection equation can be written as Eq. (10).

$\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0$ (10)

By projecting Eq. (10) onto the the space $\boldsymbol{\Omega}$, and combining
Eq. (8) & (9) result in:

$\frac{\partial\bar{\mathbf{u}}_{j}}{\partial
t}=-J_{j}^{-1}\Big{(}\mathbf{D}\mathbf{f}_{j}+\big{(}f^{c}_{l}-f_{j}(x_{j})\big{)}\mathbf{h_{l}}+\big{(}f^{c}_{r}-f_{j}(x_{j+1})\big{)}\mathbf{h_{r}}\Big{)}$
(11)

The notation used is compatible with that set out initially in [7], taking
$\mathbf{D}_{mn}$ as the contribution from the first derivative on the
$m^{th}$ Lagrange basis function to the $n^{th}$ solution point.
$\mathbf{h_{l}}$ and $\mathbf{h_{r}}$ are taken as
$dh_{l}(\boldsymbol{\xi})/d\xi$ and $dh_{r}(\boldsymbol{\xi})/d\xi$
respectively, where $\boldsymbol{\xi}$ are the reference coordinates of the
solution points. By setting $a=1$ and applying upwinding at the interfaces:

$\displaystyle J_{j}^{-1}\hat{f}^{c}_{l}$
$\displaystyle=J_{j-1}^{-1}\hat{u}_{j-1}(1)$ (12) $\displaystyle
J_{j}^{-1}\hat{f}^{c}_{r}$ $\displaystyle=J_{j}^{-1}\hat{u}_{j}(1)$ (13)

By substituting Eq. (12) & (13) into Eq. (11) and collecting the matrix
operators into $\mathbf{C}_{0}$ and $\mathbf{C}_{-1}$:

$\displaystyle\frac{\partial\bar{\mathbf{u}}_{j}}{\partial t}$
$\displaystyle=-J_{j}^{-1}\mathbf{C}_{0}\mathbf{u}_{j}-J_{j-1}^{-1}\mathbf{C}_{-1}\mathbf{u}_{j-1}$
(14) $\displaystyle\mathbf{C}_{0}$
$\displaystyle=\mathbf{D}-\mathbf{h_{l}}\mathbf{l_{l}}^{T}$ (15)
$\displaystyle\mathbf{C}_{-1}$
$\displaystyle=\mathbf{h_{l}}\mathbf{l_{r}}^{T}$ (16)

where $\mathbf{l_{l}}$ and $\mathbf{l_{r}}$ are again compatible with [7] and
defined such that $\mathbf{l_{l}}_{i}$ is the contribution of the
$i^{\mathrm{th}}$ Lagrange basis function evaluated at the left interface.
$\mathbf{l_{r}}$ is similarly defined.

Defining the continuous input as a Bloch wave, and projecting onto the
discrete solution domain:

$\displaystyle u(x,t)$ $\displaystyle=v\exp{\big{(}i(kx-\omega t)\big{)}}$
(17) $\displaystyle\mathbf{u}_{j}$
$\displaystyle=\mathbf{v}_{j}\exp{\big{(}ik(0.5(\xi+1)\delta_{j}+x_{j}-ct)\big{)}}$
(18)

Inputting this result into Eq. (14), and setting $\delta_{j}=x_{j}-x_{j-1}$
gives:

$c(k)\mathbf{v}=-\frac{i}{k}\Big{(}J_{j}^{-1}\mathbf{C}_{0}+J_{j-1}^{-1}\mathbf{C}_{-1}\exp{\big{(}-ik\delta_{j}\big{)}}\Big{)}\mathbf{v}$
(19)

Equation (19) shows that the modified phase velocity $c(k)$, is one of the
complex eigenvalues of a matrix describing the spatial transformation
performed by the scheme. For an FR scheme with order $p$ there will be $p$
eigenvalues to this problem, of which one is the physical result and the other
modes being phase shifted values to give an orthogonal set. The physical
interpretation of $c(k)$ is that a wave number’s dispersion factor is
$\Re(c(k))$ and its dissipation factor is $\Im(c(k))$.

The special case that $\delta_{j}=$  _const_ implies that $J^{-1}=$  _const_
which was the case investigated in [14, 7]. However, the more general form of
Eq. (19) allows von Neumann analysis to be performed on stretched meshes.
Importantly, Eq. (19) shows that the stencil of cells affecting the
dissipation and dispersion of an upwinded FR scheme is just the current cell
and its immediate neighbour and hence only the local expansion rate is
important for behaviour. This is clearly not the case for finite difference
schemes above second order. Repeating the analysis for an FD scheme will give
a basis of comparison and an example of the modified wavenumber for a
$4^{\mathrm{th}}$ order central difference scheme is given in Eq. (3):

$\displaystyle c(k)=\frac{i}{k}\Big{(}$ $\displaystyle
b_{-2,j}\exp{\big{(}-ik(\delta_{j}+\delta_{j-1})\big{)}}+b_{-1,j}\exp{\big{(}-ik(\delta_{j})\big{)}}+$
$\displaystyle
b_{2,j}\exp{\big{(}ik(\delta_{j+2}+\delta_{j+1})\big{)}}+b_{1,j}\exp{\big{(}ik(\delta_{j+1})\big{)}}\Big{)}$
(20)

where $b_{-2,j}$ is a weighting factor from the derivative of the Lagrange
polynomial basis function corresponding to the point $j-2$ evaluated at the
point $x_{j}$, and so on.

A further implication of Eq. (18) is that Eq. (11) can be rewritten in a form
called the update equation:

$\frac{\partial\bar{\mathbf{u}}_{j}}{\partial
t}=\mathbf{Q}\bar{\mathbf{u}}_{j}$ (21)

And for the case of pure upwinding at the interface it follows that
$\mathbf{Q}=-J_{j}^{-1}\mathbf{C}_{0}-J_{j-1}^{-1}\mathbf{C}_{-1}\exp{\big{(}-ik\delta_{j}\big{)}}$.
Putting the result into this form allows the analytical framework of Asthana
[14] to be used, hence:

$\displaystyle\bar{\mathbf{u}}_{j}(t+\tau)$
$\displaystyle=\mathbf{R}(\mathbf{Q})\bar{\mathbf{u}}_{j}(t)$ (22)
$\displaystyle\mathbf{R}_{33}$
$\displaystyle=\mathbf{I}+\frac{(\tau\mathbf{Q})^{1}}{1!}+\frac{(\tau\mathbf{Q})^{2}}{2!}+\frac{(\tau\mathbf{Q})^{3}}{3!}$
(23)

Equation (23) is an example definition of $\mathbf{R}$ for a $3^{rd}$-order
3-step low storage Runge-Kutta[28]. The form of Eq. (22) implies the von
Neumann condition of $\mathbf{R}$’s spectral radius,
$\rho(\mathbf{R})\leqslant 1\>\>\forall\>k\in\mathbb{R}$.

From the analytical form the main derived quantity to be considered is the
Points Per Wavelength (PPW) for dispersion $\mathrm{error}<1\%$ , defined as:

$\mathrm{PPW}=\frac{2\pi}{\big{\\{}\hat{k}|\inf{(|\Re{(c(\hat{k}))}-1|_{2},\epsilon)}\big{\\}}}$
(24)

where, $\epsilon$ is the error level and $c(k)$ is the convective velocity
from Eq.(19). This definition of PPW is based on the points being the solution
points in the case of FR. Consider the case of a awve whoser wavelength is the
length of an element. The normalised Nyquist wavenumber would then be
$2\pi/(p+1)$ and hence the use solution points is contained within the
definition. PPW is particularly interesting as it can be used to produce
minimum point requirements for a given region, if the scale of flow features
or the explicit filter width is known. A further derived quantity that will be
touched on briefly is the implicit filter kernel, Trefethen [29]:

$\hat{G}(\hat{k})\propto e^{t\Im(c(\hat{k}))}$ (25)

Although this is not the main subject of this paper it can be illuminating to
briefly look at the implicit filter.

To validate the analytical methods presented and gain insight into the fully
discretised scheme behaviour, a numerical approach is proposed similar to
Lele’s [26] analytical method. The methodology is to apply a scalar wave to a
one dimensional domain and allow the wave to be convected downstream. By
choosing a low CFL number of $0.01$, the spatial terms dominate the overall
error. (CFL numbers of 0.05 and 0.005 were also tested, and the results were
found to be largely unchanged, only affected by the discrete Fourier transform
as the sampling rate will vary with CFL. Hence a CFL number of 0.01 was used
to mitigate this error and give fast test turnover time). Fourier analysis is
then performed on the prescribed wave after convection through the grid. The
transform of the field $u$ is:

$\tilde{u}(x)=\sum_{k=-N/2}^{N/2}A_{k}\exp{\bigg{(}\frac{2\pi
ikx}{L}\bigg{)}}$ (26)

The method by Lele [26] would have, however, taken the spatial derivative of
the wave after convection to give:

$\tilde{u}^{\prime}(x)=\sum_{k=-N/2}^{N/2}A_{k}\bigg{(}\frac{2\pi
ik}{L}\bigg{)}\exp{\bigg{(}\frac{2\pi ikx}{L}\bigg{)}}$ (27)

By dividing Eq. (27) by Eq. (26) the modified wavenumber can be obtained as in
Eq. (28). In practice this approach has low throughput, due to need to
calculate the derivative and as the full wavenumber space has explored.
Furthermore, this method was found to be prone to the introduction of
additional sources of error caused by the method used to calculate the
derivative. Therefore, a direct comparison is made between the Fourier
Transform (FT) of the prescribed input wave and the FT of the prescribed wave
after convection through the grid. This is schematically represented in Fig.
3.

$k^{\prime}=\bigg{(}\frac{L}{2\pi
i}\bigg{)}\frac{\tilde{u}^{\prime}(x)}{\tilde{u}(x)}$ (28)

where $k^{\prime}$ is the modified wave number and equals $\omega/c(k)$. This
method is utilised for calculation of the PPW for numerical test cases in
accordance with Eq.(24).

Figure 3: Schematic of numerical scheme with incoming and outgoing wave
showing wavenumber transformation over mesh of length $L$.

##  4 Analytical Testing

###  4.1 Spatial Characteristics

The analytical method presented in Section 3 makes the implicit assumption
that there is a linear mapping between solution point placement in the real
and physical domain. The analysis did, however, carry through the ability for
the relative scaling of adjacent cells to be varied. Hence, this allows for
the characteristics of non-uniform grids to be investigated, with the
geometric expansion being one such linear transformation in common usage, for
example in the meshing of boundary layer. The geometric expansion is defined
as:

$x_{j+1}=x_{j}+\gamma(x_{j}-x_{j-1})$ (29)

where $\gamma$ is the grid expansion rate. The points that this transformation
defines are then used as the flux points for the element and linear
interpolation gives the solution points physical domain coordinates, using the
computational quadrature as weights to ensure a linear mapping. Proceeding to
analytically calculate the modified wave speed from Eq. (19), a preliminary
result that can be qualitatively informative is the filter kernel and is shown
in Fig. 4. In each case the convolution kernel is normalised independently by
the Nyquist wavenumber, $k_{nq}$. This is included as it clearly shows that
while there is merit in going to much higher orders ($p>4$), this does result
in a case of diminishing returns.

Figure 4: Analytical filter kernel, $\hat{G}(\hat{k})$, against Nyquist
normalised wavenumber, $\hat{k}$ for various order of upwinded FR with Huynh,
$g_{2}$, correction functions. The is for $t=100$ in Eq. (25).

Figure 5 shows normalised wavenumbers against normalised modified wavenumbers
for differing levels of mesh stretching. The results for the uniform grid case
can be compared to those found by Huynh [6] and found to be in agreement. What
the new results broadly exhibit is dispersion overshoot for expanding grids,
while contracting grids cause dispersion undershoot. For expanding meshes this
physically means that for a central band of wavenumbers, where the group
velocity ($\partial\omega/\partial k$) is not approximately constant, the
upwind group velocity is higher. This is caused by the change in the Nyquist
wavenumber between the smaller upwind and larger downwind elements in an
expanding mesh. Resultantly, as the solution advects downwind, this small
increase in the group velocity between elements means that a wave will be
advected into an element faster than it will exit. This gives rise to anti-
dissipation, seen in Fig 6. The opposite behaviour is exhibited in
contractions. At wavenumbers above this central band, Lagrange fitting becomes
ineffective at sufficiently projecting the prescribed wave into the functional
space, and so the dispersion relation goes zero regardless, and the
dissipation becomes high.

To highlight the practical impact of mesh deformation the PPW resulting in a
dispersion $\mathrm{error}<1\%$, Eq.(24) with $\epsilon=0.01$, is plotted
against expansion factor in Fig. 7. Over the range of expansion factors there
are some clear optimal PPW at varying polynomial orders. When Fig. 5 is
considered, it can be seen that the dispersion over- or under- shoot present
under some conditions can be counteracted by mesh warping. Therefore,
depending on mesh conditions, it may be beneficial to directionally vary the
spatial order, as this may reduce the point requirements locally. For example,
for an eddy passing through a complex mesh, fewer points would be needed while
passing through a contraction with $p=5$ compared to $p=4$ and vice versa. It
is proposed that this could be achieved by maintaining the order of the
polynomial interpolation, reflected in the number of solution points, but the
order of the correction function could be varied. Clearly this can only be
used as a means of dropping the order accuracy, and not as a means of
increasing order. From Vincent et al. [30] it can be seen that this method
results in a special case of the energy stability criterion and that the
correction functions proposed by Vincent et al. [30] with reduced order, will
still fulfil this criterion. A study was carried out to this effect, and it
was found that using a $p^{\mathrm{th}}$ order, as opposed to
$p+1^{\mathrm{th}}$ order, correction function on a $p^{\mathrm{th}}$ order
sub-domain results in a stable degradation of the spatial order, pointing to
the feasibility of this method.

(a) $p=2$

(b) $p=3$

(c) $p=4$

(d) $p=5$

Figure 5: Dispersion relations for FR, with Huynh, $g_{2}$, correction
functions, on various geometrically stretched meshes.
$--$ spectral performance, — $\gamma=0.4$,— $\gamma=0.6$, — $\gamma=0.8$, —
$\gamma=1.0$, — $\gamma=1.2$, — $\gamma=1.4$, — $\gamma=1.6$

(a) $p=2$

(b) $p=3$

(c) $p=4$

(d) $p=5$

Figure 6: Dissipation relations for FR, with Huynh, $g_{2}$, correction
functions, on various geometrically stretched meshes.
$--$ spectral performance, — $\gamma=0.4$,— $\gamma=0.6$, — $\gamma=0.8$, —
$\gamma=1.0$, — $\gamma=1.2$, — $\gamma=1.4$, — $\gamma=1.6$ Figure 7: Points
Per Wavelength (PPW) for $\mathrm{error}<1\%$ against expansion rate,
$\gamma$, for various spatial orders of upwinded FR with Huynh, $g_{2}$,
correction functions.

###  4.2 Coupled Spatial-Temporal Characteristics

The preceding analysis was predicated on using an analytical solution to the
linear advection equation that allowed the time derivative for the semi-
discrete linear system to be exactly calculated. This is of course a
simplification that cannot be used practically and Eq.(22) is far more
representative of a real implementation. The family of temporal integration
schemes to be coupled to FR here are low-storage Runge-Kutta [28], commonly
called RK33, RK44, etc.. The analysis performed is primarily in search of the
maximum stable CFL number, $\tau/\delta_{j}$ for $a=1$, obtained via varying
the time step, $\tau$, and calculation of the spectral radius of the update
matrix, $\mathbf{R}$. This is plotted in Fig. 9 for various geometric
expansion ratios. Initially focusing on contracting grids, the maximum stable
CFL number is shown to be higher than in the case of uniform grids. This could
have been expected from observation of Fig. 6, and furthermore it can be
reasoned that as a wave is swept from one cell to its smaller upwind
neighbour, the ability of the neighbouring element to resolve that wave
improves. This is due to the wave’s Nyquist normalised wavenumber, $\hat{k}$,
decreasing as it is advected through successively smaller elements.

Focusing on expanding grids, if third order is considered, $p=2$,
$\rho(\mathbf{R})\nless 1\>\>\forall\>\>\\{k\in\mathbb{R}:\gamma>1\\}$, this
means that while being strictly unstable, some wavenumbers are in practice
stable. This is displayed in Fig. 8(b) with $\rho(\mathbf{R}(k))$ being both
less than and greater than one. The practical implication is that a wave, $k$,
fed into the expanding grid can cause an instability if
$\rho(\mathbf{R}(k))>1$, however, as the wave advects the relative wavenumber
increases due to a decreasing $k_{nq}$. Hence a band of
$\rho(\mathbf{R}(k))<1$ will be encountered by the wave and the instability
will be attenuated. This procedure would be expected to repeat until the wave
is beyond the grid resolution. 111Also of note from Fig. 8 is that the
spectral radius is a periodic function which depends on the element Nyquist
wavenumber, rather than the solution point Nyquist wavenumber.

For $p>2$, different expanding mesh characteristics are seen with
$\rho(\mathbf{R})\geqslant 1\>\>\forall\>\>\\{k\in\mathbb{R}:\gamma>1\\}$,
Fig. 8(c). Later numerical tests will show that for $p>2$ instability is also
encountered, but this is likewise attenuated as the wave advects. This result
can be seen analytically be observing the dissipation relation in Fig. 6 and
using the same logic as before. The full impact on stability that this implies
will be discussed later alongside numerical findings. However, for simplicity,
the scheme stability limit (as shown in Table. 4.2) will be taken as the
higher point between either the point of the sharp increase in the spectral
radius or the point at which the spectral radius increases above 1, (see Fig.
9).

(a) $\gamma=0.9$, $p=2$

(b) $\gamma=1.1$, $p=2$

(c) $\gamma=1.1$, $p=3$

Figure 8: Spectral radius of $\mathbf{R}$ for RK44, upwinded FR with Huynh,
$g_{2}$, correction functions, against Nyquist normalised wavenumber for
$\gamma=0.9$ and $\gamma=1.1$. with various time steps $\tau$. $\tau$
increasing is shown as a decrease in colour intensity.

(a) RK33

(b) RK44

Figure 9: Spectral radius of time scheme specific update matrix, $\mathbf{R}$,
for $p=3$ upwinded FR with Huynh, $g_{2}$, corrections, against CFL number on
various grid.

Table 1: Analytical CFL limit of FR for various grid expansion factors and
temporal integration schemes. Using the Huynh correction function.

|  | CFL
---|---|---
Time Scheme | Spatial Order | $\gamma=0.7$ | $0.8$ | $0.9$ | $1.0$ | $1.1$ | $1.2$ | $1.3$
RK33 | 3 | 0.519 | 0.482 | 0.463 | 0.448 | 0.442 | 0.436 | 0.424
4 | 0.284 | 0.269 | 0.261 | 0.254 | 0.250 | 0.245 | 0.239
5 | 0.183 | 0.177 | 0.172 | 0.167 | 0.164 | 0.161 | 0.159
RK44 | 3 | 0.592 | 0.547 | 0.531 | 0.513 | 0.505 | 0.495 | 0.507
4 | 0.318 | 0.307 | 0.297 | 0.288 | 0.282 | 0.278 | 0.270
5 | 0.218 | 0.199 | 0.194 | 0.189 | 0.186 | 0.182 | 0.179
RK55 | 3 | 0.702 | 0.634 | 0.611 | 0.590 | 0.579 | 0.567 | 0.558
4 | 0.353 | 0.352 | 0.342 | 0.332 | 0.326 | 0.320 | 0.311
5 | 0.246 | 0.230 | 0.224 | 0.217 | 0.214 | 0.210 | 0.204

##  5 Numerical Results

###  5.1 Grid Stretching for Linear Advection

The analytical procedures set out up to this point have been semi-discrete and
idealised. For CFD practitioners the comparative performance and implemented
performance of FR is highly important. From this the process of mesh
generation can be informed as well as greater understanding of the expected
results gained. To this end, numerical tests are performed for wavenumbers
$0\leqslant k\leqslant k_{nq}$, where $k_{nq}$ is the Nyquist wavenumber for a
uniform mesh of unit length. For the purposes of comparison, $k$ and
$k^{\prime}$ (the modified wavenumber) are normalized by the mesh averaged
Nyquist wavenumber. This gives $0\leqslant\hat{k}=\pi k/k_{nq}\leqslant\pi$
and $0\leqslant\Re(\hat{k}^{\prime}=k^{\prime}/k_{nq})\leqslant\pi$.

Finite Difference (FD) schemes are used to provide a comparison akin to high
quality industrial codes. At higher orders, central difference schemes begin
to become unstable as, for unstretched grids, second order and greater central
difference schemes offer no dissipation. Therefore, the only sources of
numerical error are from dispersion and temporal integration. Because of this,
and as increasing order leads to better dispersion performance, the cell
Reynolds number increases. The result is that FD schemes at intermediate wave
numbers become unable to damp out disturbances. (It must be noted that due to
the very low CFL number, the temporal scheme introduces negligible numerical
error). This is combated by adding smoothing to the solution, a common
practice in industry, and is here accomplished by adding a small amount
($0.5-2\%$) of Lax-Friedrichs differencing.

Two key points are highlighted by Fig. 10(a). The first of these is that FR
requires fewer PPW than FD schemes at equivalent order. Importantly, this
means that coarser meshes can be utilised by FR for similar wave resolving
chFaracteristics to FD schemes at the same order. It is believed that the
increased accuracy of FR originates from the polynomial reconstruction in a
reference sub-domain, hence the propagation of information in FR is largely
controlled by the correction function, which can lead to superior performance.
Whereas, FD methods use a stencil, for which information can freely propagate
through, hence a less coherent solution is produced due to each point
effectively producing its own polynomial fit of the solution.

(a) Points per wavelength (PPW) for $1\%$ error against spatial order of FD
and FR schemes. FD schemes are central differencing except order 3, which is
upwinding. $\mathrm{DoF}=\mathrm{const.}=180$ and $\tau=1\times 10^{-4}\>s$.
The analytical results shown are found from Eq. (19) and similar analysis.

(b) PPW against mesh stretching factor for $4^{\mathrm{th}}$ order FD and FR
schemes. $\mathrm{DoF}=180$ and $\tau=1\times 10^{-5}\>s$.

Figure 10:

Secondly, Fig. 10(a) shows a discrepancy between the theoretical and
analytical results for both FR and FD. The origin of the error in both of
these schemes is the numerical diffusion. For FD, this is due to the scheme’s
subtle instability, meaning that for useful implementation, some diffusion
must be added to ensure the survival of the solution. FR is also affected by
numerical diffusion, but this is caused by its own intrinsic dissipation,
apparent in Fig. 6. When numerical tests are run, the dispersion and
dissipation are inseparable and although a dispersion relationship of sorts
can be found, it is impacted by the dissipation of the scheme damping out
higher wavenumbers. This is the reality of any application of a scheme and so
it can be informative to run both analytical and numerical tests as the
numerically derived PPW shown in Fig. 10(a) are those that an end user will
experience. The same effect can be seen in Fig. 10(b): as the mesh becomes
stretched the PPW rises more quickly than the analytical results would
predict. This is caused by the onset of dissipation at lower wavenumbers for
deformed meshes, but, importantly, a wave passed through a multi-element mesh
will have the transfer function applied multiple times, thereby causing
greater attenuation.

Results of great significance displayed both numerically and analytically show
that FR has the ability to resolve waves better than FD schemes and that FR is
more numerically robust when applied to geometrically stretched meshes, with
FR requiring $33\%$ of the mesh points compared to FD in 1D for severely
stretched meshes ($\gamma=1.2$ at fourth order). Moving to two or three
dimensions, this result, in the most extreme examples, can be $11\%$ or $4\%$
respectively. The increased ability of FR to handle stretched meshes is again
because of the localised fitting within sub-domains, and here the linear
transformation caused by the stretching of the elements is exactly captured in
the Jacobian. For this case the impact of adjacent cell stretching is felt
only through convection of the solution through one interface. However, for a
fully compressible Euler or Navier-Stokes implementation the effect of
adjacent cells could be increased as Riemann solvers at all interfaces would
be necessary and will give rise to more inter-cell communication.

To further understand the stability of the full numerical scheme it is
necessary to consider the spatio-temporal coupling. In the previous section
the effect of this coupling was considered and it was said that for
$p\geqslant 3$ on expanding grids, the behaviour is slightly different, _i.e._
$\rho(\mathbf{R})\geqslant 1\>\forall\>k$. The implication this has for the
stability of higher order expanding grids is not clear from the spectral
radius, however, because as the wave moves through the expanding grid
$\hat{k}$ will increase, so the scheme dissipation will add a stabilising
effect. To show this, a similar numerical method is used, however now taking a
spatial slice for various orders, wavenumbers, and grids. Figure 11 shows two
such slices. Initially, the fed wave shows some instability but advection
through the grid means dissipation from the spatial scheme will begin to
cancel some of the negative dissipation of the time scheme. The result is that
after an initial band of instability both orders show recovery of the solution
before beginning to decay on the sparser cells.

(a) $p=2$

(b) $p=4$

Figure 11: For 3rd order and 5th order FR on an expanding grid with
$\gamma=1.1$ using RK44 time integration a spatial slice is shown. With the
convected parameter $u$, the CFL number and the Nyquist normalised
($\hat{k}=\pi k/k_{nq}$).

This also illustrated a limitation of the analytical approach adopted here -
that taking a solution of the form of Eq. (17) means that the solution is
static, _i.e._ evolutions of the solution from far upstream are not permitted.
Again, this emphasises the importance of running numerical tests alongside
analytical ones.

###  5.2 Grid Warping for the Euler Equations

Few problems confronted in engineering are ever sufficiently simple that they
can solved with sufficient accuracy by 1D methods, making extension to higher
dimensionality crucial. Two dimensions also allows for a greater range of
geometrical deformations to occur, even while maintaining a linear
transformation of elements. Included within this, each element has a higher
number of degrees of freedom, revealing a potential mechanism for inarrucacies
to enter the solutions. To evaluate the effect of higher dimensionality the
isentropic convecting vortex (ICV) test case is used, as it has a known
analytical solution, so numerical error can be straightforwardly calculated
for the Euler equations. A mixture of mesh qualities are to be tested, so mesh
quality was artificially reduced by stochastically jittering corner nodes of a
uniform grid via time seeded random numbers. The degree of jitter is
controlled by a multiplying factor and the mesh quality and warp is then
characterised by a skew angle. This is defined as the mesh average absolute
angle by which the element cross diagonals deviate from square (Fig. 12) and
encompasses both the skewness and aspect ratio of a mesh. Some sample meshes
are shown in Fig. 13(a) \- 13(c).

Figure 12: Cross diagonal angle definition. $\alpha=\beta-90^{\circ}$.

(a) Average skew angle: $\alpha=\>1.0^{\circ}$

(b) Average skew angle: $\alpha=\>6.1^{\circ}$

(c) Average skew angle: $\alpha=\>15.0^{\circ}$

Figure 13: $19\times 19$ quadrilateral meshes showing differing degrees of
node jittered mesh warp.

A finite volume (FV) scheme (a simplified version of T-block), with the same
explicit time integration as used in the FR calculation was used for
performance comparison, as it is representative of a family of schemes widely
used in industry. To evaluate the spatial error, the temporal error has to be
minimised by use of an appropriately small time step, which in this case
corresponds to $\mathrm{CFL}=0.01$. ($\mathrm{CFL}=0.05\>\&\>0.005$ were also
tested and the error was found to be independent of the temporal scheme at
this level). By comparing the exact solution, $u$, and computed solution,
$u^{\prime}$, the error, $\theta$, can be calculated and the spatial order of
accuracy (OOA) can be obtained:

$\theta=u-u^{\prime}=\mathcal{O}(\delta^{p+1})$ (30)

Recovery of the spatial OOA is shown in Fig. 14, _i.e_ , $\mathrm{OOA}=p+1$,
and a comparison can be made with the $5^{th}$ order test gradient to see
this. The plotted results for the FV scheme display an $\mathrm{OOA}\approx
2$, as well as a large increase in the cell averaged $l_{2}$ error compared to
FR for the same number of Degrees of Freedom (DoF). Also shown in Fig. 14 are
the results of moderate mesh warping. In FR’s case, a move towards
$\mathrm{OOA}\approx 4$ occurs, and the FV scheme result becomes aphysical,
_i.e_ $\mathrm{OOA}=0$.

Figure 14: Point averaged $l_{2}$ norm of
error,$\frac{1}{N}\sum_{n=1}^{N}{|\theta_{n}|_{2}}$, against degrees of
freedom (DoF) for $p=4$ FR and nominally second order FV scheme. $CFL=0.01$
for $3000$ time steps.

A more detailed investigation into the deterioration of the spatial order is
performed via variation of the degree of node jittering. Tests were run with
$\mathrm{CFL}=0.01$ for 500 time steps (again, results were found to be
independent of CFL number at this level). The results of numerical tests and a
predictive procedure are shown in Fig. 15. The predictive procedure uses the
error data from a mesh of fewer degrees of freedom together with the desired
OOA to make a prediction of the error at a higher quantity of degrees of
freedom.

(a) Magnitude of mesh warp against error for various mesh densities using an
FR scheme $p=4$ ($\bullet$). Also plotted are the predicted error values for
$4^{th}$ order recovery ($-$).

(b) Magnitude of warp against error for various mesh densities using an FV
scheme ($\bullet$). Also plotted is the predicted error based on the order
recovered from Fig. 14.(-). The DoF are chosen to be comparable to the DoF of
the FR tests.

Figure 15:

Figure 15 shows that, consistent with Fig. 14, at low skew angles, FR’s
$\mathrm{OOA}$ is unaffected by mesh quality: $\mathrm{OOA}\rightarrow p+1$ as
$\alpha\rightarrow 0^{\circ}$. Figure 15(a) then confirms that on poorer
quality meshes the $\mathrm{OOA}$ drops by approximately one, here shown as a
drop from $\mathrm{OOA}=5$ to $\mathrm{OOA}\approx 4$, this is shown by
comparison to the $4^{\mathrm{th}}$ order predicted error line. This warp
induced change occurs at $\alpha\approx 1.5^{\circ}$, and as mesh skewness
increases to its maximum, FR is still able to give an accurate solution. By
comparison the FV scheme (Fig. 15(b)) undergoes warp induced error change at
$\alpha\approx 0.55^{\circ}$, transitioning from $\mathrm{OOA}\approx 2$ to
$\mathrm{OOA}=0$. Although the FV scheme is much simpler, this shows that the
performance of such schemes widely used in industry can be rapidly eroded by a
lack of mesh quality.

The root cause of the loss of order accuracy of FR on warped grids is not
fully revealed by these tests. The Jacobian that maps between the physical and
computational domain in this case provides an exact mapping due to the node
jittering providing a linear transformation. Hence, error is not introduced to
the convective velocity through the Jacobian. The additional error has two
potential sources. Firstly the linear component of FR could introduce error
via ill-conditioning of the projection to the functional space when waves are
advected at an angle. Secondly, Jameson et al. [31] derived the aliasing error
for 1D non-linear problems in the FR framework, yet in higher dimensions, on
arbitrary grids, cross multiplication of projection terms will be present.
Therefore, as the mesh become skewed aliasing has the scope to introduce
larger quantities of error. Both of these topics are quite expansive and are
left for further investigation.

As a final aside to illustrate the increase in performance that FR offers
compared to a typical FV method, the ICV test case on a uniform mesh was used.
For this test the number of degrees of freedom was varied such that the grid
averaged $l_{2}$ error was comparable - this was done with a
$\mathrm{CFL}=0.01$ for $100$ time steps. (A small variation of the CFL number
was made and the results were found to be invariant with CFL number). The test
was carried out on a single core of an Intel® Xeon® L5630 which was otherwise
idle. The results of testing are shown in Table 2, and show that FR requires
$\sim 2.6$ orders of magnitude less wall time for the same error in 2D as the
FV method used here.

Table 2: Comparable errors in a 2D ICV test for FR and FV schemes.

| Flux Reconstruction | Finite Volume
---|---|---
Point Averaged Error | $4.7\times 10^{-5}$ | $4.01\times 10^{-5}$
Wall Time | $1.6\>s$ | $639\>s$
Cells | $64,\>p=4$ | $4,000,000$

##  6 Conclusions

The use of FR on warped meshes is important for the likely future applications
of the scheme. It has been shown that FR is more resilient to distorted meshes
than some FD and FV families of schemes which are currently widely used in
industry. A detailed look at the PPW of FR on stretched grids with varying
order, as well as study of the FR stability criterion, shows that, depending
on geometry, the order of the scheme can be varied to increased performance.
In particular, within a given cell the correction function order can be
directionally varied to increase wave resolving ability. It is also shown from
the linear advection equation that the CFL limit with non-regular grids is
dependent on the dominant wave direction, with contracting grids providing a
stabilising effect. This is a feature of FR that will impact boundary layer
meshes. It was proposed that the ill-conditioning of the functional projection
of FR causes the degradation in accuracy for warped meshes. A more complete
study of this is left as future work. Lastly, in some runtime comparisons, FR
was found to require $\sim 2.6$ orders of magnitude less wall time for the
same error in two dimensional test, compared to a widely used FV method.

## Acknowledgements

The support of the Engineering and Physical Sciences Research Council of the
United Kingdom is gratefully acknowledged.

## References

  * [1] Piomelli, U., “Wall-layer models for large-eddy simulations,” Progress in Aerospace Sciences, Vol. 44, No. 6, 2008, pp. 437–446.
  * [2] Tucker, P., Unsteady Computational Fluid Dynamics in Aeronautics, Vol. 104 of Fluid Mechanics and Its Applications, Springer Netherlands, Dordrecht, 1st ed., 2014.
  * [3] Chow, F. K. and Moin, P., “A Further Study of Numerical Errors in Large-Eddy Simulations,” Journal of Computational Physics, Vol. 184, No. 2, 2003, pp. 366–380.
  * [4] Ghosal, S., “An Analysis of Numerical Errors in Large-Eddy Simulations of Turbulence,” Journal of Computational Physics, Vol. 125, No. 1, 1996, pp. 187–206.
  * [5] Tucker, P., Eastwood, S., Klostermeier, C., Jefferson-Loveday, R., Tyacke, J., and Liu, Y., “Hybrid LES Approach for Practical Turbomachinery Flows: Part 1 - Hierarchy and Example Simulations,” ASME Turbo Expo 2010: Power for Land, Sea and Air, Vol. 134, Glasgow, 2010, p. 13.
  * [6] Huynh, H. T., “A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods,” 18th AIAA Computational Fluid Dynamics Conference, Vol. 2007-4079, 2007, pp. 1–42.
  * [7] Vincent, P. E., Castonguay, P., and Jameson, A., “Insights From von Neumann Analysis of High-Order Flux Reconstruction Schemes,” Journal of Computational Physics, Vol. 230, No. 22, 2011, pp. 8134–8154.
  * [8] Huynh, H. T., “A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin for Diffusion,” 47th AIAA Aerospace Science Meeting, No. January in Fluid Dynamics and Co-located Conferences, American Institute of Aeronautics and Astronautics, jun 2009, pp. 1–34.
  * [9] Castonguay, P., Williams, D. M., Vincent, P. E., and Jameson, A., “Energy Stable Flux Reconstruction Schemes for Advection-Diffusion Problems,” Computer Methods in Applied Mechanics and Engineering, Vol. 267, No. 1, 2013, pp. 400–417.
  * [10] Huynh, H. T., “High-Order Methods Including Discontinuous Galerkin by Reconstructions on Triangular Meshes,” 49th AIAA Aerospace Sciences Meeting, No. January, 2011, pp. 1–28.
  * [11] Williams, D. M. and Jameson, A., “Energy stable flux reconstruction schemes for advection-diffusion problems on tetrahedra,” Journal of Scientific Computing, Vol. 59, No. 3, 2014, pp. 721–759.
  * [12] Sheshadri, A., An Analysis of Stability of the Flux Reconstruction Formulation With Applications to Shock Capturing, Phd thesis, Stanford University, 2016.
  * [13] Vincent, P. E., Farrington, A. M., Witherden, F. D., and Jameson, A., “An Extended Range of Stable-Symmetric-Conservative Flux Reconstruction Correction Functions,” Computer Methods in Applied Mechanics and Engineering, Vol. 296, 2015, pp. 248–272.
  * [14] Asthana, K. and Jameson, A., “High-Order Flux Reconstruction Schemes with Minimal Dispersion and Dissipation,” Journal of Scientific Computing, Vol. 62, No. 3, 2014, pp. 913–944.
  * [15] Asthana, K., López-Morales, M. R., and Jameson, A., “Non-Linear Stabilization of High-Order Flux Reconstruction Schemes via Fourier-Spectral Filtering,” Journal of Computational Physics, Vol. 303, 2015, pp. 269–294.
  * [16] Chung, Y. M. and Tucker, P. G., “Accuracy of higher-order finite difference schemes on nonuniform grids,” AIAA journal, Vol. 41, No. 8, 2003, pp. 1–3.
  * [17] You, D., Mittal, R., Wang, M., and Moin, P., “Analysis of Stability and Accuracy of Finite-Difference Schemes on a Skewed Mesh,” Journal of Computational Physics, Vol. 213, No. 1, 2006, pp. 184–204.
  * [18] Mengaldo, G., De Grazia, D., Vincent, P. E., and Sherwin, S. J., “On the Connections Between Discontinuous Galerkin and Flux Reconstruction Schemes: Extension to Curvilinear Meshes,” Journal of Scientific Computing, Vol. 67, No. 3, 2016, pp. 1272–1292.
  * [19] Kopriva, D. A., “Metric identities and the discontinuous spectral element method on curvilinear meshes,” Journal of Scientific Computing, Vol. 26, No. 3, 2006, pp. 301–327.
  * [20] Sheshadri, A., Crabill, J., and Jameson, A., “Mesh Deformation and Shock Capturing Techniques for High-Order Simulation of Unsteady Compressible Flows on Dynamic Meshes,” SciTech 2015, , No. January, 2015, pp. 1–15.
  * [21] Castonguay, P., High-order Energy Stable Flux Reconstruction Schemes For Fluid Flow Simulations on Unstructured Grids, Ph.D. thesis, Stanford University, 2012.
  * [22] Reed, W. H. and Hill, T. R., “Triangular Mesh Methods for the Neutron Transport Equation,” Los Alamos Report LA-UR-73-479, , No. 836, 1973, pp. 10.
  * [23] Roe, P. L., “Approximate Riemann Solvers , Parameter Vectors, and Difference Schemes,” Journal of Computational Physics, Vol. 43, 1981, pp. 357–372.
  * [24] van Leer, B., “Fourteenth International Conference on Numerical Methods in Fluid Dynamics,” Lecture Notes In Physics, , No. February, 1982, pp. 11–15.
  * [25] Harten, A., Lax, P. D., van Leer, B., and Vanleer, B., “On upstream differencing and Godunov-type schemes for hyperbolic conservation laws,” SIAM Review, Vol. 25, No. 1, 1983, pp. 35–61.
  * [26] Lele, S. K., “Compact Finite Difference Schemes With Spectral-Like Resolution,” Journal of Computational Physics, Vol. 103, No. 1, 1992, pp. 16–42.
  * [27] Hesthaven, J. S. and Warburton, T., Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Vol. 54 of Texts in Applied Mathematics, Springer New York, New York, NY, 1st ed., 2008.
  * [28] Kennedy, C. A., Carpenter, M. H., and Lewis, R. M., “Low-Storage, Explicit Runge-Kutta Schemes for the Compressible Navier-Stokes Equations,” Applied Numerical Mathematics, Vol. 35, No. 1, 2000, pp. 177–219.
  * [29] Trefethen, L. N., Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations, Ithaca, New York, 1st ed., 1994.
  * [30] Vincent, P. E., Castonguay, P., and Jameson, A., “A New Class of High-Order Energy Stable Flux Reconstruction Schemes,” Journal of Scientific Computing, Vol. 47, No. 1, 2010, pp. 50–72.
  * [31] Jameson, A., Vincent, P. E., and Castonguay, P., “On The Non-Linear Stability of Flux Reconstruction Schemes,” Journal of Scientific Computing, Vol. 50, No. 2, 2012, pp. 434–445.


* [171] L. Spangher, A. Tawade, A. Devonport, C. Spanos, Engineering vs. ambient type visualizations: Quantifying effects of different data visualizations on energy consumption, in: Proceedings of the 1st ACM International Workshop on Urban Building Energy Sensing, Controls, Big Data Analysis, and Visualization, UrbSys’19, Association for Computing Machinery, New York, NY, USA, 2019, p. 14–22.  doi:10.1145/3363459.3363527.
URL https://doi.org/10.1145/3363459.3363527

  * [172] H. P. Das, I. C. Konstantakopoulos, A. B. Manasawala, T. Veeravalli, H. Liu, C. J. Spanos, A novel graphical lasso based approach towards segmentation analysis in energy game-theoretic frameworks, in: 2019 18th IEEE International Conference On Machine Learning And Applications (ICMLA), IEEE, 2019, pp. 1702–1709.
  * [173] H. P. Das, I. Konstantakopoulos, A. B. Manasawala, T. Veeravalli, H. Liu, C. J. Spanos, Do occupants in a building exhibit patterns in energy consumption? analyzing clusters in energy social games, in: Workshop on Tackling Climate Change with Machine Learning, Conference on Neural Information Processing Systems (NeurIPS), 2020.
  * [174] L. Spangher, A. Gokul, M. Khattar, J. Palakapilly, A. Tawade, A. Bouyamourn, A. Devonport, C. Spanos, Prospective experiment for reinforcement learning on demand response in a social game framework, in: Proceedings of the Eleventh ACM International Conference on Future Energy Systems, 2020, pp. 438–444.
  * [175] B. Orland, N. Ram, D. Lang, K. Houser, N. Kling, M. Coccia, Saving energy in an office environment: A serious game intervention, Energy and Buildings 74 (2014) 43–52.
  * [176] D. Geelen, D. Keyson, S. Boess, H. Brezet, Exploring the use of a game to stimulate energy saving in households, Journal of Design Research 14 10 (1-2) (2012) 102–120.
  * [177] M. Casals, M. Gangolells, M. Macarulla, N. Forcada, A. Fuertes, R. V. Jones, Assessing the effectiveness of gamification in reducing domestic energy consumption: Lessons learned from the energaware project, Energy and Buildings 210 (2020) 109753.
  * [178] D. S. Watson, S. Kiliccote, N. Motegi, M. A. Piette, Strategies for demand response in commercial buildings (2006).
  * [179] P. Alstone, J. Potter, M. A. Piette, P. Schwartz, M. A. Berger, L. N. Dunn, S. J. Smith, M. D. Sohn, A. Aghajanzadeh, S. Stensson, et al., 2025 california demand response potential study-charting california’s demand response future. final report on phase 2 results, Tech. rep., Lawrence Berkeley National Lab.(LBNL), Berkeley, CA (United States) (2017).
  * [180] M. H. Albadi, E. F. El-Saadany, Demand response in electricity markets: An overview, in: 2007 IEEE power engineering society general meeting, IEEE, 2007, pp. 1–5.
  * [181] Federal Energy Regulatory Commission, Energy primer: a handbook for energy market basics.
URL https://www.ferc.gov/sites/default/files/2020-06/energy-primer-2020_0.pdf

  * [182] C. Deb, F. Zhang, J. Yang, S. E. Lee, K. W. Shah, A review on time series forecasting techniques for building energy consumption, Renewable and Sustainable Energy Reviews 74 (2017) 902–924.
  * [183] J. Jazaeri, T. Alpcan, R. Gordon, M. Brandao, T. Hoban, C. Seeling, Baseline methodologies for small scale residential demand response, in: 2016 IEEE Innovative Smart Grid Technologies-Asia (ISGT-Asia), IEEE, 2016, pp. 747–752.
  * [184] L. Xu, S. Wang, R. Tang, Probabilistic load forecasting for buildings considering weather forecasting uncertainty and uncertain peak load, Applied energy 237 (2019) 180–195.
  * [185] Y. Kim, H.-g. Son, S. Kim, Short term electricity load forecasting for institutional buildings, Energy Reports 5 (2019) 1270–1280.
  * [186] T. X. Nghiem, C. N. Jones, Data-driven demand response modeling and control of buildings with gaussian processes, in: 2017 American Control Conference (ACC), IEEE, 2017, pp. 2919–2924.
  * [187] J. L. Mathieu, P. N. Price, S. Kiliccote, M. A. Piette, Quantifying changes in building electricity use, with application to demand response, IEEE Transactions on Smart Grid 2 (3) (2011) 507–518.
  * [188] D. Jung, V. B. Krishna, W. G. Temple, D. K. Yau, Data-driven evaluation of building demand response capacity, in: 2014 IEEE International Conference on Smart Grid Communications (SmartGridComm), IEEE, 2014, pp. 541–547.
  * [189] S. Burger, J. P. Chaves-Avila, C. Batlle, I. Perez-Arriaga, The value of aggregators in electricity systems (2016).
  * [190] U. Agwan, C. J. Spanos, K. Poolla, Asset participation and aggregation in incentive-based demand response programs, in: 2021 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm), IEEE, 2021, pp. 89–94.
  * [191] P. Shinde, M. Amelin, A literature review of intraday electricity markets and prices, 2019 IEEE Milan PowerTech (2019) 1–6.
  * [192] A. Brusaferri, M. Matteucci, P. Portolani, A. Vitali, Bayesian deep learning based method for probabilistic forecast of day-ahead electricity prices, Applied Energy 250 (2019) 1158–1175.
  * [193] U. Agwan, Optimal prosumer aggregations: Design and modeling (2020).
  * [194] W. Saad, Z. Han, H. V. Poor, Coalitional game theory for cooperative micro-grid distribution networks, in: 2011 IEEE international conference on communications workshops (ICC), IEEE, 2011, pp. 1–5.
  * [195] D. Wu, X. Ma, S. Huang, T. Fu, P. Balducci, Stochastic optimal sizing of distributed energy resources for a cost-effective and resilient microgrid, Energy 198 (2020) 117284.
  * [196] M. Quashie, C. Marnay, F. Bouffard, G. Joós, Optimal planning of microgrid power and operating reserve capacity, Applied Energy 210 (2018) 1229–1236.
  * [197] U. Agwan, K. Poolla, C. J. Spanos, Optimal composition of prosumer aggregations, in: 2021 IEEE PES Innovative Smart Grid Technologies Europe (ISGT Europe), IEEE, 2021, pp. 1–5.
  * [198] S.-J. Kim, G. B. Giannakis, An online convex optimization approach to real-time energy pricing for demand response, IEEE Transactions on Smart Grid 8 (6) (2016) 2784–2793.
  * [199] N. Liu, X. Yu, C. Wang, C. Li, L. Ma, J. Lei, Energy-sharing model with price-based demand response for microgrids of peer-to-peer prosumers, IEEE Transactions on Power Systems 32 (5) (2017) 3569–3583.
  * [200] L. Jia, Q. Zhao, L. Tong, Retail pricing for stochastic demand with unknown parameters: An online machine learning approach, in: 2013 51st Annual Allerton Conference on Communication, Control, and Computing (Allerton), IEEE, 2013, pp. 1353–1358.
  * [201] U. Agwan, L. Spangher, W. Arnold, T. Srivastava, K. Poolla, C. J. Spanos, Pricing in prosumer aggregations using reinforcement learning, in: Proceedings of the Twelfth ACM International Conference on Future Energy Systems, 2021, pp. 220–224.
  * [202] B. Dong, Y. Liu, W. Mu, Z. Jiang, P. Pandey, T. Hong, B. Olesen, T. Lawrence, Z. O’Neil, C. Andrews, et al., A global building occupant behavior database, Scientific data 9 (1) (2022) 1–15.
  * [203] G. Fierro, M. Pritoni, M. AbdelBaky, D. Lengyel, J. Leyden, A. Prakash, P. Gupta, P. Raftery, T. Peffer, G. Thomson, et al., Mortar: an open testbed for portable building analytics, ACM Transactions on Sensor Networks (TOSN) 16 (1) (2019) 1–31.
  * [204] C. Miller, A. Kathirgamanathan, B. Picchetti, P. Arjunan, J. Y. Park, Z. Nagy, P. Raftery, B. W. Hobson, Z. Shi, F. Meggers, The building data genome project 2, energy meter data from the ASHRAE great energy predictor iii competition, Scientific data 7 (1) (2020) 1–13.
  * [205] B. Balaji, A. Bhattacharya, G. Fierro, J. Gao, J. Gluck, D. Hong, A. Johansen, J. Koh, J. Ploennigs, Y. Agarwal, et al., Brick: Metadata schema for portable smart building applications, Applied energy 226 (2018) 1273–1292.