# Spectral/hp element methods: recent developments, applications, and
perspectives

Hui Xu  Chris D. Cantwell  Carlos Monteserin  Claes Eskilsson  Allan P.
Engsig-Karup  Spencer J. Sherwin  Department of Aeronautics, Imperial College
London, London SW7 2AZ, United Kingdom  Department of Applied Mathematics and
Computer Science, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark
Center for Energy Resources Engineering (CERE), Technical University of
Denmark, 2800 Kgs. Lyngby, Denmark  Department of Civil Engineering, Aalborg
University, DK-9220 Aalborg Ø, Denmark  Division Safety and Transport,
Research Institutes of Sweden (RISE), SE-50115 Borås, Sweden

###### Abstract

The spectral/hp element method combines the geometric flexibility of the
classical h-type finite element technique with the desirable numerical
properties of spectral methods, employing high-degree piecewise polynomial
basis functions on coarse finite element-type meshes. The spatial
approximation is based upon orthogonal polynomials, such as Legendre or
Chebychev polynomials, modified to accommodate $\mathcal{C}^{0}$-continuous
expansions. Computationally and theoretically, by increasing the polynomial
order p, high-precision solutions and fast convergence can be obtained and, in
particular, under certain regularity assumptions an exponential reduction in
approximation error between numerical and exact solutions can be achieved.
This method has now been applied in many simulation studies of both
fundamental and practical engineering flows. This paper briefly describes the
formulation of the spectral/hp element method and provides an overview of its
application to computational fluid dynamics. In particular, it focuses on the
use the spectral/hp element method in transitional flows and ocean
engineering. Finally, some of the major challenges to be overcome in order to
use the spectral/hp element method in more complex science and engineering
applications are discussed.

###### keywords:

High-precision spectral/hp elements, Continuous Galerkin method, Discontinuous
Galerkin method, Implicit large eddy simulation

††journal: Journal of Hydrodynamics

## Project Leader

Spencer Sherwin is Professor of Computational Fluid Mechanics and the Head of
Aerodynamics Section in the Department of Aeronautics at Imperial College
London. He received his MSE and PhD from the Department of Mechanical and
Aerospace Engineering Department at Princeton University in 1995. Prior to
this he received his BEng from the Department of Aeronautics at Imperial
College London in 1990. In 1995 he joined the Department of Aeronautics at
Imperial College as a lecturer and subsequently became a full professor in
2005. Over the past 27 years he has specialised in the development and
application of advanced parallel spectral/hp element methods for flow around
complex geometries with a particular emphasis on vortical and bluff body
flows, biomedical modelling of the cardiovascular system and more recently in
industrial practice through partnerships with McLaren Racing and Rolls Royce.
Professor Sherwin’s research group also develops and distributes the open-
source spectral/hp element package Nektar++ (www.nektar.info) which has been
applied to direct numerical simulation, large eddy simulation and stability
analysis to a range of applications including vortex flows of relevance to
offshore engineering and vehicle aerodynamics as well as biomedical flows
associated with arterial atherosclerosis. He has published numerous peer-
reviewed papers in international journals covering topics from numerical
analysis to applied and fundamental fluid mechanics and co-authored a highly
cited book on the spectral/hp element method. Since 2014 Professor Sherwin has
served as an associate editor of the Journal of Fluid Mechanics. He is a
Fellow of the Royal Aeronautical Society, a Fellow of the American Physical
Society and in 2017 he was elected a Fellow of the Royal Academy of
Engineering.

##  1 Introduction

Over the past few decades, computational fluid dynamics (CFD) has become
increasingly powerful and has therefore been seen as the natural starting
point to investigate a variety of mathematical and physical problems in
science and engineering. Traditionally, “low-order methods” with up to second-
order spatial accuracy have been widely adopted as the default implementation
for the simulation of fluid flows, often based around the Reynolds Averaged
Navier-Stokes (RANS) equations. This approach has achieved a great deal of
success in many applications due to its well-established robustness and
efficiency. However, as the demands on the accuracy of CFD outputs have
increased, such as requiring the solution of the unsteady flow equations in
complex geometries, low-order methods are less able to provide the necessary
level of precision in capturing transient dynamics, as compared to higher
accuracy schemes for a given computational cost. Therefore there is currently
a great interest in the development and application of high-order methods such
as the spectral/hp element discretisation.

Finite element methods are widely used across a broad range of engineering and
scientific disciplines. They can be categorised into three classes [1, 2]: (1)
the classical h-version finite element method; (2) the spectral element
method, or p-version finite element method; (3) the hp-version finite element
method, called the spectral/hp element method. Once the computational domain
is partitioned into a non-overlapping element set, the spectral/hp element
method employs a “spectral-like” approach in each element, representing the
solution through a basis of polynomials. Therefore, the spectral/hp element
method combines the advantages of the spectral element method, in terms of the
properties of accuracy and rapid convergence, with those of the classical
h-version finite element method, that allows complex geometries to be
effectively captured. Compared with traditional low-order finite element
schemes, the spectral/hp element method can provide an arbitrary order of
spatial-approximation accuracy under the assumption of sufficient smoothness
of the solution. It therefore combines the advantages of the low-order finite
element method family, whilst also providing an additional attractive higher-
precision, approximation for solving partial differential equations [3, 4].
Recent studies also indicate that the compact nature of the spectral/hp
element method is well-placed to take advantage of modern multi-core and many-
core computing hardware [5, 6]. All of the above properties have positioned
the spectral/hp element method as an attractive numerical strategy for
obtaining high-precision solutions with a relatively low computational cost.

The spectral/hp element method is gaining increasing traction in the field of
CFD [3], and it has achieved great success in both the direct and large eddy
simulation (LES) of complex flows [7, 8, 9, 10, 11]. It has also been
successfully applied to a broad range of other applications in various
research fields, for instance, cardiac electrophysiology [12], solid mechanics
[13, 14], porous media [15] and oceanographic modelling [16].

One of the challenges of these methods is that they are more challenging to
implement than low-order methods. There are now, however, a number of open-
source packages which encapsulate the complexities of the method. One such
package that the authors have been developing is Nektar++, a cross-platform
spectral/hp element framework (http://www.nektar.info) which has made high-
order finite element methods more accessible to the broader community and can
be used to solve a range of the emerging challenges in high-fidelity
scientific computations [4]. It enables the construction of high-order
discretisations for solving a wide range of partial differential equations,
supporting hybrid meshes, for example, triangles and quadrilaterals or
tetrahedra, prisms, pyramids and hexa-hedra. The current software package
includes a number of pre-written scalable solvers, including for
incompressible flow, compressible flow, shallow water equations, acoustic
perturbation equations and others.

In this paper, we present a review of the state-of-the-art of the spectral/hp
element method, and its applications. In §2, we briefly describe the
spectral/hp element method and discuss the numerical advantages and robustness
of the method in solving partial differential equations. In §3, we highlight
some of the implementational challenges of the method in order to achieve
computational efficiency and discuss potential solutions through our efforts
with Nektar++. In §4, we then provide a survey of the applications in
transitional flows and waves (wave propagation and wave-body interaction) in
ocean engineering to which the spectral/hp element method has been applied.
Finally, in §5, we conclude with further discussions and some perspectives on
future directions.

##  2 Spectral/hp element method

In this section we provide an overview of the mathematical foundations of the
spectral/hp element method. A more detailed derivation of the mathematical
theory can be found in [3] but is beyond the scope of this paper.

Figure 1: Illustration of the construction of a two-dimensional fourth-order
$\mathcal{C}^{0}$-continuous modal triangular expansion basis using a
generalised tensor-product procedure.

###  2.1 Overview of method

Spectral/hp element discretisations are the underpinning approximation for
both continuous and discontinuous Galerkin formulations and can be constructed
in 1-D, 2-D and 3-D. To obtain a general view of the method we can consider
the example in Fig. 1. At the bottom of this figure we observe a triangulation
of the British Isles. In general, the computational mesh can comprise of a mix
of different shaped elements which could be triangles and quadrilaterals in
two dimensions, and tetrahedra, pyramids, prisms and hexahedra in three
dimensions. Within each element we develop a polynomial expansion, as
represented in the lower-right of Fig. 1, where we observe all of the
expansion modes to represent a fourth-order polynomial using a modal or
hierarchical expansion basis. For Continuous Galerkin expansions, the design
of these elements typically has the property of being decomposed into boundary
and interior modes so that $\mathcal{C}^{0}$-continuity can be achieved by
matching boundary expansion modes that have support on the edges of the
element. The interior modes are necessarily defined to be zero on the
boundaries.

There are a number of ways of developing the polynomial expansion within an
element $K_{i}$. In the example shown in Figure 1 we illustrate the popular
construction using a tensor product of one-dimensional expansions which are
defined in the regular regions $K$ and $K^{\prime}$ which can then be mapped
to the physical elemental regions using the mapping $\vec{\chi}(\vec{\xi})$.
Tensor-based expansions are very popular in quadrilateral and hexahedral
regions but require some modification for triangular or tetrahedral regions
using a Duffy [17] or collapsed-coordinate system (see [3] for more details).
It is possible to define polynomial expansions that are hierarchical/modal in
construction, much like a small-scale spectral expansion, or alternatively one
can use a nodal basis where the expansion is defined at a specific set of
points (often related to Gaussian quadrature) and only one expansion mode has
a unit value, whilst all other expansions basis are zero, at this point. This
provides a Kronecker delta property to the expansion which can be useful when
evaluating non-linear products. Again, for triangular elements the
construction of a nodal basis is more involved and typically requires the use
of a generalised Vandemonde matrix to relate them to a known expansion which
is often of a modal construction (see [18]).

Once a polynomial expansion has been defined in each element, the approach
used to “bolt together” these individual expansions defines the approximation
method. For example, if one enforces a $\mathcal{C}^{0}$-continuous expansion
by ensuring the polynomial expansion is continuous across elemental boundaries
we obtain the classical continuous Galerkin method. Alternatively, if one does
not directly enforce the expansions to be $\mathcal{C}^{0}$-continuous but
ensures appropriate flux quantities are continuous between elements, one can
construct a so-called Discontinuous Galerkin scheme [3, 18, 19].

###  2.2 Outline of mathematical formulation

In general we consider the numerical solution of partial differential
equations (PDEs) of the form $\mathcal{L}u=0$ on a domain $\Omega$, which may
be geometrically complex, for some solution $u$. Practically, $\Omega$ takes
the form of a $d$-dimensional finite element mesh consisting of elements
$K_{i}$, embedded in a space of dimension $\hat{d}$, such that
$d\leq\hat{d}\leq 3$, with $\Omega=\cup_{i}K_{i}$ and $K_{i}\cap K_{j}$ is an
empty set or an interface between elements of dimension $\bar{d}<d$. We solve
the PDE problem in the weak sense and, in general, $u|_{K_{i}}$ must be smooth
with at least a first-order derivative; we therefore require that $u|_{K_{i}}$
is in the Sobolev space $W^{1,2}(K_{i})\equiv H^{1}(K_{i})$ [20]. For a
continuous discretisation, we additionally impose continuity along element
interfaces. For illustrative purposes, we assume that it can be expressed as
follows:

$\mbox{find $u\in H^{1}(\Omega)$ such that }a(u,v)=l(v)\ \forall v\in
H^{1}(\Omega).$ (1)

where $a(\cdot,\cdot)$ is a bilinear form, $l(v)$ is a linear form and
$H^{1}(\Omega)$ is formally defined as

$H^{1}(\Omega):=\\{v\in L^{2}(\Omega)|D^{\alpha}u\in L^{2}(\Omega)\ \forall\
|\alpha|\leq 1\\}.$ (2)

To solve this problem numerically, we consider solutions in a finite-
dimensional subspace $V^{\delta}\in H^{1}(\Omega)$ and cast our problem as:

$\mbox{find $u^{\delta}\in H^{1}(\Omega)$ such that
}a(u^{\delta},v^{\delta})=l(v^{\delta})\ \forall\ v^{\delta}\in V^{\delta},$
(3)

augmented with appropriate boundary conditions. For a projection which
enforces continuity across elements, we impose the additional constraint that
$V^{\delta}\in\mathcal{C}^{0}$. We assume the solution can be represented as
$u^{\delta}({\bf x})=\sum_{n}\hat{u}_{n}\Phi_{n}({\bf x})$, a weighted sum of
$N$ trial functions $\Phi_{n}({\bf x})$ defined on $\Omega$ and our problem
becomes that of finding the coefficients $\hat{u}_{n}$. The approximation
$u^{\delta}$ does not directly give rise to unique choices for the
coefficients $\hat{u}_{n}$. To achieve this we place a restriction on the
residual $R={\mathcal{L}}u^{\delta}$ that its $L^{2}$ inner product, with
respect to the test functions $\Psi_{n}({\bf x})$, is zero. For a Galerkin
projection we choose the test functions to be the same as the trial functions,
that is $\Psi_{n}=\Phi_{n}$. As outlined previously, to construct the global
basis $\Phi_{n}$ we first consider the contributions from each element in the
domain. Each $K_{i}$ is mapped from a standard reference space
$K\in[−1,1]^{d}$ by a parametric mapping ${\chi}_{e}:K\rightarrow K_{i}$ given
by ${\rm x}={\chi}_{e}({\xi})$, where $K$ is one of the supported region
shapes, and ${\xi}$ are $d$-dimensional coordinates representing positions in
a reference element, distinguishing them from ${\rm x}$ which are
$\hat{d}$-dimensional coordinates in the Cartesian coordinate space. The
mapping need not necessarily exhibit a constant Jacobian, supporting deformed
and curved elements through an isoparametric mapping. For triangular,
tetrahedral, prismatic and pyramid elements, the Duffy transform [17] is used
and operations, such as calculating derivatives, map the coordinate system to
the non-collapsed coordinate system.

A local polynomial basis is constructed on each reference element with which
to represent solutions. A one-dimensional order-$P$ basis is a set of
polynomials $\phi_{p}(\xi),0\leq p\leq P$, defined on the reference segment.
The particular choice of basis is usually made based on its mathematical or
numerical properties and may be modal or nodal in nature. For two- and three-
dimensional regions, a tensorial basis may be used, where the polynomial space
is constructed as the tensor-product of one-dimensional bases on segments,
quadrilaterals or hexahedral regions. In particular, a common choice is to use
a modified hierarchical Legendre basis, given as a function of one variable by

$\phi_{p}(\xi)\left\\{\begin{array}[]{ll}\frac{1-\xi}{2},&p=0,\\\
\left(\frac{1-\xi}{2}\right)\left(\frac{1+\xi}{2}\right)P^{1,1}_{p-1}(\xi),&0<p<P,\\\
\frac{1+\xi}{2},&p=P,\end{array}\right.$ (4)

which supports the aforementioned boundary-interior decomposition and
therefore improves numerical efficiency when solving the globally assembled
system. Equivalently, $\phi_{p}$ could be defined by the Lagrange polynomials
through the Gauss-Lobatto-Legendre quadrature points which would lead to a
traditional spectral element method. On a physical element $K_{i}$ the
discrete approximation $u^{\delta}$ to the solution $u$ may be expressed as

$u^{\delta}({\rm x})=\sum_{n}\hat{u}_{n}\phi_{n}([{\chi}_{e}]^{-1}({\rm x}))$
(5)

where $\hat{u}_{n}$ are the coefficients from Eq. (1), obtained through
projecting $u$ onto the discrete space. Therefore, we restrict our solution
space to

$V\coloneqq\\{u\in H^{1}(\Omega)\ |\ u|_{K_{i}}\in\mathcal{P}_{p}(K_{i})\\}$
(6)

where $\mathcal{P}_{p}(K_{i})$ is the space of order-$p$ polynomials on
$K_{i}$. Elemental contributions to the solution may be assembled to form a
global solution through an assembly operator. In a continuous Galerkin
setting, the assembly operator sums contributions from neighbouring elements
to enforce the $\mathcal{C}^{0}$-continuity requirement. In a discontinuous
Galerkin formulation, such mappings transfer flux values from the element
interfaces into global solution vector.

Figure 2: Propogation of a cone initial condition using a DG formulation of
the advection equation. After one rotation we see the solution using a $P$=2,
$P$=6 and $P$=8 order approximations.

###  2.3 What are the beneficial properties of using a spectral/hp element
method?

For smooth problems one of the benefits of a spectral/hp formulation is often
quoted as the exponential rate of convergence one observes when increasing the
polynomial order of the approximation, providing the solution has sufficient
regularity. In practice, the good approximation properties are also realised
in terms of the better diffusion and dispersion properties one observes as you
increase the polynomial order. This is exemplified in Fig. 2 where we have
taken an initial condition of a Gaussian and advected it under a rotational
velocity field for one revolution. We then observe the solution of a DG
approximation of the linear advection problem using a $P$=1, $P$=3 and $P$=8
order approximation with $128,32$ and $8$ elements, respectively. We see from
this simple test that the shape is far less diffused in the $P$=8 problem when
compared to the $P$=2 discretisation, even though the overall number of
degrees of freedom is comparable.

Recently, a comparison between high-order flux reconstruction and the
industrial-standard solver STAR-CCM+ was made, in which the accuracy and
computational cost of flux reconstruction (FR) methods in PyFR and STAR-CCM+
for a range of test cases including scale-resolving simulations of turbulent
flow were investigated [21]. The results from both PyFR and STAR-CCM+ show
that third-order schemes can be more accurate than second-order schemes for a
given cost. Moreover, advancing to higher-order schemes on GPUs with PyFR was
found to offer even further accuracy vs cost benefits relative to industry-
standard tools. These demonstrate the potential utility of high-order methods
on GPUs for scale-resolving simulations of unsteady turbulent flows.

###  2.4 Using the spectral/hp element method in marginally resolved problems

In the previous sections we outlined the underlying formulation of the
spectral/hp element methodology and we now describe some of the state-of-the-
art developments necessary for simulating problems which are only marginally
resolved. These tend to arise at higher Reynolds numbers, which generally
promote greater levels of turbulence that cannot realistically be captured,
even on the latest parallel computers. This leads to the build-up of energy in
the larger resolved scales, due to the inability of the discretisation to
fully capture the smaller scales of the flow. This is addressed with two
techniques. The first involves performing the numerical integration of non-
linear terms consistently. Even with this, non-linear interactions may still
lead to an artificial build-up of energy in the smaller resolved scales. A
second technique is then to apply a spectral filtering, known as spectral
vanishing viscosity, to dampen such energies.

####  2.4.1 Dealiasing

Dealiasing strategies have been observed to be effective in enhancing the
numerical stability when solving problems using the spectral/hp element method
[22, 23, 24, 25, 26]. The errors are typically caused by insufficient
quadrature employed in the Galerkin discretisation of nonlinear terms [3].
When the simulation is not sufficiently well-resolved, this leads to energy in
the shorter length scales to be transferred back into the longer length
scales. When simulations are well-resolved, the numerical errors created by
insufficient quadrature are negligible [23, 24]. However, in marginally
resolved or under-resolved simulations, aliasing errors may significantly
pollute the solution [23]. As discussed in [26], there are three kinds of
aliasing sources:

  * •

PDE aliasing, which relates to quasi-linear and non-linear terms (see [23] for
a simple example of the Burgers equation);

  * •

geometrical aliasing, which arises due to deformed/curved elements.

  * •

interface-flux aliasing, in the case of discontinuous methods (see [26] for
more details).

We start by giving an example from incompressible viscous flow past a NACA
0012 wing tip at a Reynolds number of $Re=1.2\times 10^{6}$, originally
discussed in [26], to demonstrate the aliasing error in the near surface
region. Figure 3 illustrates the dynamics of the flow by showing the iso-
contours of helicity with a $P$=3 approximation. In the boundary region,
prismatic elements are employed while tetrahedral elements are used in the
rest of the domain. Figure 4, shows the first two kinds of aliasing errors. In
Figure 4(a), 30% aliasing error with respect to the magnitude of the nonlinear
terms is observed and in Figure 4(b) we show a close up of regions with high
geometrical aliasing error near the wing surface where 600% error is observed.
The aliasing errors in this numerical calculation significantly pollute the
solution and cause it to very rapidly become numerically unstable.

Figure 3: Helicity for an incompressible viscous flow over a NACA 0012 wing
(from the reference [26])

Figure 4: PDE aliasing errors at 30% of the magnitude of the nonlinear terms
(a) and close up regions near the wing surface showing geometric aliasing
errors (b) (from the reference [26])

We can categorise dealiasing techniques into two strategies.

  * •

Local dealiasing: PDE-dealiasing through consistent integration of the
nonlinear terms only [23, 26];

  * •

Global dealiasing: PDE- and geometrical-dealiasing through consistent
integration of all the terms of the discretisation [26].

The phrase consistent integration refers to polynomial nonlinearities, where
it is effectively possible to know a priori the number of quadrature points
necessary for the integration to be exact [23, 27]. Therefore, the non-
linearities can be consistently integrated. For non-polynomial functions, as
typically arise in the Euler equations, the concept of consistent integration
is not well defined since it is impossible to fully control the quadrature
error. We refer to the reader to [26] for a more complete discussion of these
strategies as well as a discussion of interface aliasing.

####  2.4.2 Spectral vanishing viscosity

As previously discussed, it is well known that spectral/hp element
discretisations generally lead to approximations that have low dissipation and
low dispersion [3] per degree of freedom when compared to lower-order methods.
In solving advection-diffusion equations and nonlinear partial differential
equations such as advection-dominated flows, at marginal resolutions,
oscillations appear that may render the computation unstable [3]. Artificial
viscosity has been used in may discretisation methods to suppress wiggles
associated with high wavenumbers, for example hyperviscosity has been broadly
and effectively used in simulations using the Fourier method. A related
concept is the so-called spectral vanishing viscosity (SVV), which was
originally proposed based on a second-order diffusion (convolution) operator
for spectral Fourier methods [28]. The SVV concept was originally motivated
through the inviscid Burgers equation where an additional diffusion term was
added, i.e.

$\frac{\partial}{\partial t}u(x,t)+\frac{\partial}{\partial
x}\left(\frac{u^{2}(x,t)}{2}\right)=\epsilon\frac{\partial}{\partial
x}\left(\mathscr{D}_{\epsilon}\frac{\partial u(x,t)}{\partial x}\right),$ (7)

where $\epsilon(\rightarrow 0)$ is a viscosity amplitude and
$\mathscr{D}_{\epsilon}$ is a viscosity kernel, which may be nonlinear, can be
a function of $x$ and is only activated for high wave numbers. With a small
amount of controlled dissipation, spectral accuracy can be retained. This
approach was extended to the Navier-Stokes equations by Karamanos &
Karniadakis [29] and Kirby & Sherwin [30] and has also been investigated and
applied to Large Eddy Simulation by Pasquetti [31, 32].

To provide some intuition about the influence of SVV, we show in Figure 5 an
example from [30] of applying SVV to the parabolic equation

$\frac{\partial u}{\partial t}=\nu\nabla^{2}u+S_{VV}(u)$ (8)

on $\Omega=[0,2]\times[0,2]$ with $\nu=10^{-5}$ and periodic boundary
conditions. The initial condition is given by $u(x,y,t=0)=\sin(\pi x)\sin(\pi
y)$, in which a perturbation is introduced. A total of sixteen quadrilateral
elements (four per direction) with $15th$ order polynomials per element
direction were used. The following “classical” one-dimensional SVV kernel is
employed,

$\mathscr{D}=\left\\{\begin{array}[]{ll}0&p\leq P_{\rm cut},\\\
\exp\left(-\frac{(p-P)^{2}}{(p-P_{\rm cut})^{2}}\right)&p>P_{\rm
cut},\end{array}\right.$ (9)

where $P$ is the total number of modes employed and $P_{\rm cut}$ is the
cutoff polynomial order. SVV with the kernel function $\mathscr{D}$ can be
regarded as a low-pass filter. We see that the SVV dissipation added to the
high mode numbers with respect to the spectral element discretisation does
indeed yield dissipation at the global high wavenumber scales of the solution
in Fig. 5.

Figure 5: A demonstration of stabilisation by SVV: (left) standard diffusion
at $t$=0.1; (centre) standard diffusion with SVV $P_{\rm cut}$=7,
$\epsilon$=0.1; (right) standard diffusion with SVV $P_{\rm cut}$=3,
$\epsilon$=0.1 ( from the reference [30]).

In applications, Karamanos and Karniadakis applied SVV to high Reynolds number
turbulent flows and the method was viewed as an alternative LES approach [29].
Kirby and Karniadakis computed the artificial viscosity by incorporating the
local strain and the Panton function and the methodology was called SVV-LES
[33]. The SVV approach has been reasonably widely used to simulate turbulent
flow and vortical flows [7, 30, 31, 32, 34, 35, 36]. SVV has been explicitly
regarded as a turbulent model of implementing implicit large eddy simulation
(iLES) under the assumption “The action of subgrid scales on the resolved
scales is equivalent to strictly dissipative action” [37], even though SVV is
not explicitly designed as a subgrid-scale model.

Recently, in order to achieve a better performance of SVV, the linear
advection equation augmented with spectral vanishing viscosity (SVV) has been
analysed [38]. Dispersion and diffusion characteristics of the spectral/hp
continuous Galerkin (CG) formulation with SVV-based stabilization are verified
to display similar non-smooth features in flow regions where convection is
much stronger than dissipation or vice-versa, owing to a dependency of the
standard SVV operator on a local Péclet number $Pe=ch/\epsilon$ where $a$ is
the advection speed, $h$ is the mesh spacing and $\epsilon$ is the base SVV
magnitude [38]. A modification was proposed [38] to the traditional SVV
scaling, which enforces a globally constant Péclet number. In this approach
the artificial dissipation strength parameter $\epsilon$ is defined by
$\epsilon=\epsilon_{0}ch/P$, where $\epsilon_{0}$ is a globally fixed
parameter used to adjust the dissipation strength, $P$ would be the polynomial
order used in each element and $c$ is a measure of the advection velocity
magnitude on each local element. Furthermore, the following “power” kernel
function is proposed [38]

$\mathscr{D}=(p/P)^{P_{SVV}},$ (10)

where $p$ denotes mode index and $P_{SVV}$ is now not an activation threshold,
but a similar effect in terms of computing the viscous effects on the highest
modes. The mesh Péclet number is now held constant globally. In addition, the
“power kernel” function has been devised for the advocated SVV operator to
provide a consistent increase in resolution power (per degree of freedom) when
the polynomial order is increased – a feature not naturally achieved through
the widely used “exponential kernel” introduced in [39].

In turbulence simulations using Discontinuous Galerkin methods, there were
early discussions [40, 41], highlighting that the scheme successfully
predicted low-order statistics with fewer degrees of freedom (DoFs) than
traditional numerical methods. More extensive assessments [42, 43] indicated
that DG can predict high-order statistics with accuracy comparable to that of
spectral methods for an equivalent number of DoFs. However, there is still
little understanding of why and how one can use the spectral/hp element method
for under-resolved turbulence simulations when either spectral vanishing
viscosity (as sometime employed in Continuous Galerin methods) or upwind
fluxes (as naturally arise in DG methods) are used to provide dissipation at
under-resolved scales. Recently, exploratory studies, focusing on the DG
formulation, were undertaken to address how numerical dissipation affects the
under-resolved scales [44]. By exploring DG’s propagation characteristics
directly from the dispersion-diffusion curves of a linear advection problem, a
simple criterion (named ‘the 1% rule’) was proposed to estimate the effective
resolution of the DG scheme for a given hp approximation space.

Numerical experiments on Burger turbulence helped to clarify why the DG
formulation can be suitable for under-resolved simulations of turbulence: the
numerical discretisation is capable of resolving scales up to $k_{1\%}$ with
good accuracy while dissipation is provided at the end of the energy spectrum
in the form of numerical diffusion. A further study for improving SVV by the
“1%” rule was conducted with a new SVV kernel function[38]. The “1%” rule
enables the SVV to mimic the property of the DG method with upwind numerical
fluxes. A CG formulation based on the suggested SVV operator and kernel
function was investigated. Being able to emulate the upwind properties of a DG
scheme within a CG approximation using SVV appears to have relatively robust
properties which are attractive for the marginally resolved simulations and
this approach is referred to as a ‘DGKernel’ SVV within the Nektar++ framework
discussed in the next section.

##  3 Implementation of spectral/hp element method

One of the challenges of enabling the spectral/hp element method to become
more pervasive for industrial and environmental problems is the
implementational complexity of the methods. One of the overarching aims of the
Nektar++ project is to provide an efficient framework upon which a broad range
of physical processes can be modelled using these approaches. With this in
mind, we outline some of efficiency challenges associated with implementing
these methods for practical use and detail some suitable solutions in the
context of Nektar++.

###  3.1 Nektar++ overview

Figure 6: Diagrammatic representation of the libraries in Nektar++ (from the
reference [4]).

The design of Nektar++ is intended to mirror the mathematical formulation of
the spectral/hp element method, from one-dimensional basis functions up to
multi-dimensional, multi-element discretisations of complex geometries. The
package consists of a tiered collection of libraries which implement different
aspects of the formulation, as illustrated in Figure 6. On top of this,
domain-specific application solvers can be readily developed which specify the
physical processes to be modelled, while leveraging the implementation of the
discretisation provided by the libraries in a relatively transparent manner.

At the lowest level the LibUtilities implements one-dimensional bases of
orthogonal polynomials which are used to construct the expansion bases on
individual elements. Their analytical derivatives and the distributions of
points used for performing Gaussian quadrature are also captured. On top of
this, the reference elements are defined for tessellating one-, two- and
three-dimensional domains in the StdRegions library. Nektar++ supports a range
of elemental regions to enable maximum flexibility in defining meshes which
can effectively capture complex geometries, including segments in 1D,
quadrilaterals and triangles in 2D, and hexahedra, prisms, pyramids and
tetrahedra in 3D. For multi-dimensional regions expansion bases may be defined
using natively two-dimensional bases (e.g. Lagrange polynomials at Fekete
points), or though a tensor-product construction of 1D polynomials. For each
of these regions, core finite element operators, including backward transform,
inner product and derivative operators are defined.

Constituent elements of the computational mesh are represented in the
LocalRegions library through a mapping from the reference region to the
physical coordinates of the element. These mappings are captured in the
SpatialDomains library and can capture high-order geometric curvature in the
element. These geometric factors or used to extend the finite element
operators defined for the standard reference regions to the physical elements.
The MultiRegions encapsulates the construction of multi-element global
expansions through the construction of assembly operators which identify
corresponding degrees of freedom in the local and global representations.
Through this mapping, global finite element operators are constructed. Global
linear systems can then be solved using a number of direct and iterative
linear algebra techniques. Iterative approaches, such as the preconditioned
conjugate gradient method, solve a system through the repeated action of the
operator, providing scope for a range of performance optimisations, discussed
in Section 3.3.

While the central concept captured in Nektar++ are high-order spectral/hp
element spatial discretisations, high-order time integration algorithms are
also implemented to allow for highly accurate transient simulations. These
include fully explicit, implicit and mixed implicit-explicit (IMEX) schemes.

###  3.2 Combined spectral/hp-Fourier discretisation

While the spectral/hp element method exhibits some of the numerical
characteristics of a pure spectral method, the use of the latter is still
preferred in the case of regular geometries such as rectangles and cuboidal
domains. This situation occurs widely in many, particularly fundamental,
hydrodynamics problems. In Nektar++, we leverage this by allowing for a mixed
spectral/hp-Fourier discretisation. In this configuration, geometrically
homogeneous coordinate directions are represented with a truncated Fourier
expansion, while coordinate directions with geometric complexity use the
spectral/hp element discretisation. This has the effect of decoupling the
individual Fourier modes, allowing the spectral/hp operator to be applied
independently within each component of each Fourier mode.

The Fourier decomposition also enables parallel decomposition of the modes to
different processes [45]. Combined with the parallelisation of the spectral/hp
element discretisation this leads to a hybrid parallelism. An open challenge
here is identifying the optimal parallelisation of such mixed-discretisation
domains to achieve optimal performance.

The GlobalMappings library further extends this discretisation strategy by
supporting geometries in which a coordinate direction can be transformed to a
homogeneous representation by an analytical function. This enables more
complex geometries to be discretised using this efficient formulation.

###  3.3 From h to p efficiently

One of the challenges in implementing the spectral/hp method is achieving
computational performance across a broad range of polynomial orders. The
classic approach to implementing linear finite element operators is to
assemble the global form of the operator and act directly on the global
degrees of freedom. This is performant since mesh vertex has a high valency
and the assembly operation dramatically reduces the size of the global linear
system used within the iterative solver. With increasing polynomial order,
this reduction is less and the memory indirection associated with the global
operator leads to reduced through-put and degraded performance.

A number of studies [46, 47, 48, 49] have concluded that for higher polynomial
orders, one should instead scatter the global degrees of freedom back onto
their local elemental representation and perform the action of the operator in
an element-wise fashion using a matrix-vector operation. This enables a more
compact representation of the operator in memory, greater cache-locality in
applying the operator and a reduction in overall memory usage. At very high
polynomial orders, the size of the local operator grows as
$\mathcal{O}(P^{4})$. A more efficient approach in these cases is to take
advantage of the tensor-product construction of the operators to decompose the
elemental operator into a sequence of one-dimensional operators, known as sum-
factorisation. Algebraically, this corresponds to a sequence of matrix-matrix
operations which, besides the reduction in floating operation count, can be
generally more efficiently implemented.

Recalling the construction of the local elemental operators as the reference
region operator, transformed under a geometric mapping, one can exploit the
decomposition of the elemental operators into these constituent parts. The
Collections library in Nektar++ exploits this decomposition on groups of
similar elements (same shape and expansion basis) by applying first the
geometric factors associated with the operator, before acting with the
reference region operator on all of the elements in the group simultaneously
[50]. This effectively replaces a sequence of matrix-vector operations with a
single matrix-matrix operation, reducing data movement and improving floating-
point performance. This decomposition can be applied both to the local
elemental matrix representation as well as the sum-factorisation approach.

###  3.4 High-order mesh generation

Figure 7: Splitting a reference prismatic element and applying the mapping to
obtain a high-order layer of prisms from the macro-element (from the reference
[51])

Figure 8: Surface mesh of a NACA 0012 wing geometry at polynomial order $P$=
5. The coarse prismatic boundary layer is highlighted in blue (left) and the
refined boundary layer which uses three layers and a progression ratio $r$= 3
(right) (from the reference [51])

For complex three-dimensional geometries, the generation of high-order curved
meshes along solid surfaces is a challenging topic [51, 52, 53, 54, 55].
Approaches for generating unstructured high-order meshes, extending coarse
linear grids to incorporate high-order curvature, are a relatively new
development[52, 56, 57]. The main challenge is robustness, since near-wall
curvature must be introduced in such a way as to prevent the generation of
self-intersecting elements [51]. Recently, several methods have been proposed
to cope with high-order curvilinear boundary layer meshing, namely (1) the
isoparametric approach [51]; (2) the variational approach [53], and; (3) the
thermo-elastic analogy [54]. In the isoparametric approach, a reference
element is mapped to a physical element with the supplied curvature
information. The generation of a boundary layer mesh is achieved by splitting
the reference element along the wall normal coordinate and then applying the
mapping from a reference element to a physical element yielding the desired
boundary layer mesh as illustrated in Fig. 7. In Fig. 8, a mesh for the
NACA0012 aerofoil profile is shown which is generated by the isoparametric
approach (see [51] for details). In the variational approach, a functional
which defines a measure of energy over mesh and takes the mesh displacement
and its derivatives as its arguments, is minimised using a nonlinear
optimisation strategy [53]. It is demonstrated that the variational framework
is efficient for both mesh quality optimisation and untangling of invalid
meshes. In order to circumvent the drawback of elastic-model based methods
which are used to tackle element self-intersection, a thermo-elastic analogy,
as an extension of the elastic formulation, is proposed to ‘heat’ or ‘cool’
elements. The thermo-elastic formulation leads to an additional degree of
robustness, which shows the potential to significantly improve meshing [54].
NekMesh, which is a mesh generation and manipulation utility bundled with
Nektar++, supports these strategies for high-order mesh generation.

##  4 Applications

In this section, we illustrate some typical applications of spectral/hp
element methods in fluid mechanics, typically in hydrodynamics. Some examples
will be given through use of the pre-written solvers in Nektar++.

###  4.1 The incompressible Navier-Stokes equations

We first introduce the incompressible Navier-Stokes equations on a bounded
domain $\Omega$ which allows one to solve the governing equations for viscous
Newtonian fluids governed by

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}$ $\displaystyle=-\nabla
p+\nu\nabla^{2}\mathbf{u},$ (11a) $\displaystyle\nabla\cdot\mathbf{u}$
$\displaystyle=0,$ (11b)

with appropriate boundary conditions In Nektar++, a velocity correction scheme
is employed, which uses a splitting/projection method where the velocity and
the pressure are decoupled [58]. In the original approach, a stiffly stable
time integration scheme was proposed. Briefly, high-order splitting schemes
comprise of three steps involving explicit eval-uation of the non-linear
terms, followed by the implicit solution of the pressure Poisson system and
finally solving a series of Helmholtz problems to enforce the viscous terms
and velocity boundary conditions. Time integration is handled using a generic
time-stepping framework [59], utilising one of a number of implicit-explicit
(IMEX) time-integration schemes.

###  4.2 Transitional flows

In the past few decades, the spectral/hp element method has been applied to
laminar and transitional flows (e.g. [58, 60, 61, 62, 63, 64]). As discussed
earlier, the advantages of the spectral/hp element method are its low
numerical dissipation and dispersion errors, which allow individual flow
structures to be accurately captured over long time periods. Combined with the
development of effective stabilisation techniques (e.g., SVV or filter-based
stabilisation [65]), the spectral/hp element method has in most recent years
proved effective in modelling turbulent flow [33, 66, 67, 68, 69, 70].

Transitional problems are those in which turbulence dominates the flow domain,
or in which the transition to turbulence is a fundamental aspect of the
simulation. Here we provide a survey of the range of studies undertaken using
the spectral/hp element method in relation to transitional flows, focusing on
(1) turbulence; (2) separated flows; (3) hydrodynamic stability, and; (4)
vortical flows and high Reynolds number wingtip vortices.

Figure 9: Turbulent channel flow (left) and turbulent pipe flow (right)

.

####  4.2.1 Turbulent flows

As noted in [4], highly resolved turbulent simulations were traditionally
undertaken using spectral methods, imposing strong geometrical restrictions.
Spectral/hp element methods alleviate this although, as discussed in §3, when
the domain of interest has a geometrically homogeneous direction, a
combination of the spectral/hp element method and the traditional spectral
method is still particularly advantageous [45, 71].

A benchmark simulation of turbulent flow over a periodic hill has been
performed using Nektar++. It is challenging to resolve accurately due to the
detachment of the fluid from the smooth surface and the generation of a
recirculation region downstream of the hill [4]. A two-dimensional mesh of
3626 quadrilateral elements with $P=6$ in the streamwise direction was
constructed, and a Fourier pseudo-spectral method with 160 collocation points
in the spanwise direction was employed to perform the simulation at a Reynolds
number 2800. Excellent agreement with the benchmark statistics was obtained.
With the quasi-3D solver, turbulent channel and pipe flows have also been
accurately simulated and validated [45]; these are illustrated in Fig. 9 and
are used as benchmark cases in Nektar++.

Figure 10: An illustration of the computational domain for the turbulent
simulation in a boundary layer (left) and laminar-turbulent transition
triggered by TS instability in a boundary layer. The iso-surfaces indicate the
different pressure levels (from the reference [11])

.

Recently, a global mapping technique has been developed which allows the
simulation of fluid flows with geometrically periodic variation along
homogeneous directions [8, 11]. A typical example is to simulate turbulent
flows induced by a localised surface deformation in the boundary layer [11].
The computational domain is illustrated in Fig. 10(left) and the turbulent
flow which forms dowstream of the surface deformation is shown in Fig.
10(right). The freestream unit Reynolds number is $1.2\times 10^{6}$ and the
reference free-stream velocity is 18 m/s. The laminar-turbulent transition is
initialised by the nonlinear process of upstream Tollmien-Schlichting (TS)
disturbance breakdown to turbulence downstream of the surface indentation. The
disturbance is generated by a disturbance strip as illustrated in Fig.
10(left), which has a vibration frequency of 172Hz and is used to mimic the TS
disturbance by the receptivity mechanism. The numerical calculations of this
case are validated with an experimental study for which excellent agreement is
obtained.

####  4.2.2 Separated flows

Flow separation can be one of the most important topics in fluid mechanics,
due to its relevance to aerodynamic performance in many engineering
applications. There are typically two kinds of separations – geometrical,
where the flow separates from a sharp obstacle in the flow; and pressure-
induced, depending on the pressure gradient over a smooth surface. For both
kinds of separation, the spectral/hp element method has been used to
accurately capture detachment and reattachment points, e.g. laminar flow in a
channel expansion [72], turbulent separation in a three-dimensional diffuser
[25], the flow around a wall-mounted square cylinder [73], flows over periodic
hills [74], small separation bubbles induced by localised imperfections [11,
75]. Practically, the studies of separated flows are always accompanied by
turbulent flows due to the occurrence of sensitive instabilities induced by
the geometrical discontinuity.

Recently, with the global mapping technique, direct numerical simulations of
the flow around wings with spanwise waviness were investigated to explore its
effect on the wing performance [9, 76]. The geometries are based on a NACA0012
profile with a small modification to obtain zero-thickness trailing edge to
overcome meshing challenges. The wavy geometries were obtained by applying the
following coordinate transformation to the straight infinite wing

$\bar{x}=x-\frac{h}{2}\cos\left(\frac{2\pi}{\lambda}z\right),$ (12)

(a) Perspective view (b) Planform .

Figure 11: The geometry of a wavy wing with $h/c$=0.1 and $\lambda/d$=0.5
(from the reference [9])

(a) Recirculation zones (b) Contours of spanwise velocity .

Figure 12: Flow properties around the wing with an attack angle $12^{\circ}$,
$h/c$=0.1 and $\lambda/d$=0.5 (from the reference [9])

where $h$ is the waviness peak-to-peak amplitude, $\lambda$ is its wavelength,
$\bar{x}$ is the physical coordinate in the cord direction, and $x$ and $z$
are the chordwise and spanwise coordinates in the computational domain in Fig.
11. At a very low Reynolds number 1000, for moderate angles of attack, the
waviness leads to decreases in both drag and lift forces, which lead to a
decrease in the lift- to-drag ratio and a suppression of the fluctuating lift
coefficient [9]. The physical mechanism behind these were explained. In Fig.
12(a), the recirculation regions are shown and the visualisation was made
considering the regions where streamwise velocity is negative. In Fig. 12(b),
the coloured contours represent spanwise velocity at the plane $x$=-0.03,
which is close to the location of the leading edge in the peak of the waviness
at $x$=-0.05, depicting how the flow moves away from the waviness peak in the
lower portion of the wing and it moves towards it in the upper part. A further
study was undertaken at Reynolds numbers 10,000 and 50,000 for different
attack angles through highly resolved direct numerical simulations [76], which
provides a better understanding of wing performance with the use of spanwise
waviness.

####  4.2.3 Hydrodynamic stability

In addition to solving the fully non-linear incompressible Naiver-Stokes
equations in Nektar++, the solution of the linearised incompressible Navier-
Stokes equations is also supported to enable global flow stability analyses to
be performed with respect to a steady or a time-periodic base flow. This
process identifies whether such flows are susceptible to a fundamental change
of state when an infinitesimal disturbance is introduced. The linearised
incompressible Navier-Stokes equations take the form

$\displaystyle\frac{\partial\mathbf{u^{\prime}}}{\partial
t}+\mathbf{U}\cdot\nabla\mathbf{u^{\prime}}+\mathbf{u^{\prime}}\cdot\nabla\mathbf{U}$
$\displaystyle=-\nabla p+\nu\nabla^{2}\mathbf{u^{\prime}}+\mathbf{f},$ (13a)
$\displaystyle\nabla\cdot\mathbf{u^{\prime}}$ $\displaystyle=0,$ (13b)

where $\mathbf{u}^{\prime}$ is the perturbation and $\mathbf{U}$ denotes the
base flow or the time-dependent periodic flow sampled at regular intervals and
interpolated. For stability analysis, suitable boundary conditions should be
imposed. The linear evolution of a perturbation under Eq. (13) can be
expressed as

$\frac{\partial\mathbf{u}^{\prime}(t)}{\partial
t}=\mathscr{A}(t)\mathbf{u}^{\prime}(t),$ (14)

with an initial condition $\mathbf{u}^{\prime}(0)$. If the base flow
$\mathbf{U}$ is steady, the perturbation $\mathbf{u}^{\prime}(t)$ be expressed
by eigenmode solutions of $\mathscr{A}(t)$ as follows

$\mathbf{u}^{\prime}(\vec{x},t)=\exp(\lambda_{j}t)\tilde{\mathbf{u}}_{j}+{\rm
c}.{\rm c}.$ (15)

The dominant eigenvalues and eigenmodes of the operator $\mathscr{A}(t)$ are
defined by solutions to the following equation

$\mathscr{A}(t)\tilde{\mathbf{u}}_{j}=\lambda_{j}\tilde{\mathbf{u}}_{j},$ (16)

The leading eigenvalue $\lambda_{j}$ is used to detect the global stability of
the flow. Using a similar approach, the operator $\mathscr{A}^{*}(t)$ for the
adjoint form of the linearised Navier-Stokes evolution operator can be used to
examine the receptivity of the flow and, in combination with the direct mode,
identify the sensitivity to base flow modification and local feedback. The
direct and adjoint methods can be combined to explore convective instabilities
over different time horizons $\tau$ in a flow by computing the leading
eigenmodes of $(\mathscr{A}^{*}\mathscr{A})(\tau)$ known as transient growth
analysis. The eigenvalues and eigenmodes of $\mathscr{A}(\tau)$ and
$(\mathscr{A}^{*}\mathscr{A})(\tau)$ can be solved directly if the dimensions
of the matrices are small. Practically, time-stepper-based methodologies based
on Arnoldi iteration methods (which are now supported in Nektar++) [77], can
be more efficient for complex problems [78, 79, 80, 81] than direct
calculation of the eigenvalue spectrum from the linearised Navier-Stokes
operator.

(a) $u$ (b) $v$

(c) $\tilde{u}$ (d) $\tilde{v}$ .

Figure 13: Direct stability analysis of two-dimensional flow past a circular
cylinder at $Re$=50: (a) the x-component of the base flow $\mathbf{u}$; (b)
the y-component of the base flow $\mathbf{u}$; (c) the x-component of the
dominant direct mode;(d) the y-component of the dominant direct mode.

To illustrate the direct stability analysis with a classic example, we
calculated the dominant mode of the two-dimensional flow past a circular
cylinder at $Re=50$. We show the two components of the base flow in Figs. 13
(a) and (b) and the two components of the dominant direct mode in Figs. 13 (c)
and (d). The leading eigenmode is characterised by the asymmetry in the
streamwise component and symmetry in the cross-stream component. We also
observe the spatial distribution of the modes, with the leading direct modes
extending far downstream of the cylinder. Recently, the instability methods
using the spectral/hp element method have been employed for investigating the
stability analysis of vortical flows and controlling wakes of flows past bluff
bodies [82, 83, 84]. With the linearised Navier-Stokes equations, the
interaction between instability waves and surface distortion in a 2-D and 3-D
boundary layer have been precisely investigated [75, 10, 11].

####  4.2.4 Vortical flows and wingtip vortex

Simulating and understanding vortical flows are vital in hydrodynamics and
aerodynamics. In many applications, manipulation of vortices in the vicinity
of flow boundaries is crucial for improving performance in engineering
practice [7, 85, 86, 87, 88, 89]. Developing a better understanding of the
near wake of the vortex, lying within one chord length of the trailing edge,
is essential in understanding the complex flow structure interactions [7],
which is also crucial in understanding cavitation in vortical structures
generated by hydro propellers [86]. Computationally, in high Reynolds number
flows, accurate numerical simulations of these kinds of vortical structures
are challenging for the traditional numerical methods due to numerical
dissipation. A recent study in which wingtip vortices were simulated
demonstrated that the adjustable and controllable low-dissipation properties
of the spectral/hp element method were beneficial for modelling and simulating
vortical flows [7, 44, 90].

To illustrate the ability of the spectral/hp element method to accurately
capture vortical flows at a Reynolds number $1.2\times 10^{6}$, a rectangular
wing with a NACA0012 profile is investigated with a rounded wing cap
(consequently, a longer semispan where the wing is thickest) and a blunt
trailing edge [7]. The results showed better correlation with experimental
results than previous numerical results, both in terms of the static pressure
distribution, prediction of the jetting velocity, vortex spanwise location,
and the ability to resolve the secondary vortex interaction with the main
wingtip vortex. The iLES method based on SVV has been shown to be a compelling
alternative for computing complex vortex-dominated flows, such as the wingtip
vortex, motivating its use for complex industrially relevant cases where high-
fidelity computational fluid dynamics can become an enabling technology.

Figure 14: Tip vortex of an elliptic hydrofoil

.

Figure 15: Vorticity components $\omega_{x}$, $\omega_{y}$ and $\omega_{z}$
along axial (left), tangential (mid) and spanwise (right) directions at the
section $x/c=1.06$

.

Built on the confidence of simulating the NACA0012 wingtip vortex, recently,
an initial simulation with a Reynolds number $1.217\times 10^{6}$ is currently
being conducted to simulate the tip vortex of an elliptic hydrofoil where the
incident flow attack angle is $7^{\circ}$. The tip vortex structure is
illustrated by the iso-surfaces of the total pressure in Fig. 14. The contour
figures of the vorticity components are shown at the section $x/c=1.06$ in
Fig. 15 and the origin is located at the tip of the hydrofoil. In Fig. 15
(left), the vortex core is shown in axial vorticity $\omega_{x}$ and
tangential and spanwise vorticity components ($\omega_{y}$ and $\omega_{z}$)
are shown in dipoles in Fig. 15 (mid) and (right). This calculation will be
used to provide vortex core pressure for studying cavitation in a vortex
structure. The numerical results will be going to be compared with
experimental data.

###  4.3 Waves in ocean engineering

This section is devoted to outlining the recent progress in spectral/hp
element modelling of wave propagation and wave-body interaction. For water
waves propagation the fully nonlinear potential flow (FNPF) equations is the
fundamental governing equation and can be derived from the Navier-Stokes
equations by assuming inviscid and irrotational flow. One of the main
challenges in the last decade have focused on use of spectral/hp element
methods for coastal engineering [91], and the development of robust and
efficient solvers with support for unstructured meshes to capture realistic
shore lines, geometric features, and adapt meshes to relevant features of the
solution [92].

####  4.3.1 Fully nonlinear potential flow

Finite elements are widely used for solving the FNPF equations [93, 94, 95,
96], but the use of spectral/hp elements remains scarse. The first attempt to
solve the FNPF equations using spectral/hp elements is due to Robertson and
Sherwin [97] using an Arbitrary Lagrangian Eulerian (ALE) approach for the
free water surface in a two-dimensional setting. Solving the FNPF equations is
non-trivial, due to the need to evolve a set of highly nonlinear free surface
boundary conditions in a robust way, together with the efficiency solution of
a Laplace problems – possibly – in large ocean areas. For spectral/hp element
models, Robertson and Sherwin identified a mesh asymmetry problem associated
with triangulation of the fluid domain, which gives rise to real-valued
eigenvalues and thus represented a severe instability problem. This was
handled by adding a diffusive term in the kinematic free surface boundary
condition proportional to the mesh asymmetry, but at the cost of reduced
convergence rates. The problem of mesh asymmetry which destroys the dispersion
relation was circumvented by Engsig-Karup et al. [98, 99] by solving the FNPF
equations in a $\sigma$-transformed domain [94] using a single layer of
quadrilaterals in 2D and prisms in 3D. The FNPF model was shown to exhibit
exponential convergence even for steep nonlinear stream function waves which
most simpler wave models cannot handle. To stabilise the nonlinear
simulations, it was found that the quartic nonlinear terms in the free surface
equation need to be integrated exactly and for the steep waves a mild (1%)
modal filter could be used to remove high frequency noise that can stem from
aliasing in nonlinear terms for marginal resolution, and the gradient recovery
steps that reconstruct the gradients of the solution (e.g. velocities) from
the $\mathcal{C}^{0}$ approximations.

Letting $\bar{\nabla}$ denote the horizontal gradient operator, the velocity
potential $\phi$ satisfy the Laplace equation

$\bar{\nabla}\phi+\frac{\partial^{2}\phi}{\partial z^{2}}=0\,,$ (17)

with time-dependent free surface boundary conditions on $\Gamma_{z=\eta}$
given in the Zakahrov form [100]

$\displaystyle\frac{\partial\eta}{\partial
t}+\bar{\nabla}\cdot\bar{\nabla}\tilde{w}-\tilde{w}\left(1+\bar{\nabla}\eta\cdot\bar{\nabla}\eta\right)=0\,,$
(18) $\displaystyle\frac{\partial\tilde{\phi}}{\partial
t}+g\eta+\frac{1}{2}\left(\bar{\nabla}\tilde{\phi}\cdot\bar{\nabla}\tilde{\phi}+\tilde{w}^{2}\left(1+\bar{\nabla}\eta\cdot\bar{\nabla}\eta\right)\right)=0\,,$
(19)

and zero normal flow on the bottom $\Gamma_{z=-h}$ and on structures
$\Gamma_{b}$

$\frac{\partial\phi}{\partial n}=0\,.$ (20)

Here $\eta$ is the free surface elevation, $w$ the vertical velocity, $h$ the
still water depth and $g$ the acceleration of gravity. The
$\sigma$-transformed FNPF solver as presented in [99] has been implemented in
the Nektar++ framework in order to allow for larger-scale modelling [101]. The
solution is advanced explicitly in time by solving the free-surface equations.
The computed $\tilde{\phi}^{n+1}$ are used as Dirichlet condition for the
Laplace equation at $\Gamma_{z=\eta}$. After solving for the velocity
potential, $\tilde{w}^{n+1}$ is obtained by a $C^{0}$ gradient recovery as
$\tilde{w}^{n+1}$ is needed to advance the free surface conditions another
time step. Explicit time-stepping schemes are effective for FNPF as the time-
step restriction is not dependent on the horizontal mesh size, but only on the
vertical resolution and water depth, see [98]. A $\sigma$-transformed FNPF
solver is well-suited for large-scale wave propagation simulations, since re-
meshing is avoided. Also, an intrinsic property of FNPF solvers is the ability
to get accurate kinematics in all of the fluid domain as a part of the
solutions, which is relevant for wave-induced load predictions on marine
structures. The $\sigma$-transformed FNPF solver can be used for computing
e.g. short wave disturbances in ports, see Fig. 16, and wave scattering from
bottom mounted vertical cylinders, see Fig. 17.

Figure 16: Snap shot of free surface elevation in the Port of Visby using
fifth order method.

(a)

(b)

(c)

Figure 17: Wave scattering of regular waves by a cylinder. Wave propagating
from left to right. (a) Maximum wave run-up on front-side of cylinder, (b)
wave passing cylinder and (c) wave run-up on lee-side

The use of $\sigma$-transformed domains excludes any truncated bodies in the
domain. In order to handle arbitrarily shaped bodies in the domain the mixed-
Eulerian-Lagrangian (MEL) approach should be used and there are ongoing
efforts to implement a SEM based on the MEL approach. In [102] a work-around
on the mesh asymmetry problem associated with the MEL was presented. By using
hybrid meshes consisting of a single layer of vertically aligned quads (in 2D)
at the free surface, the eigenvalues were shown to be purely imaginary for any
triangulation of the inner field. When using higher-order local polynomial
expansions as done in the spectral/hp element method, a local re-meshing
(movement of quadrature points inside the local elements) as well as global
re-meshing (movement of vertices) are needed in the MEL approach. An example
of a spectral/$hp$ element simulation using MEL for wave-body interaction are
illustrated in Fig. 18. Here a solitary wave is passing over a horizontal
submerged cylinder.

(a)

(b)

(c)

Figure 18: Solitary wave passing a submerged horizontal cylinder. (a) zoom in
on the hybrid mesh, (b) horizontal velocity and (c) vertical velocoity (from
the reference [102])

####  4.3.2 Shallow water wave equations

While the application of spectral/hp element method for FNPF equations is
emerging, in the field of depth-integrated shallow water long wave equations
there is a wider use of spectral/hp elements methods. These equations can be
derived from the FNPF equations by expanding the velocity potential in terms
of the vertical coordinate and integrating the Laplace problem over the fluid
depth. This removes the vertical dependence of the problem. Non-hydrostatic
wave equations (such as Boussinesq-, Serre- or Green-Naghdi-type equations)
are used for wave propagation and transformation in nearshore regions. The
dispersive effects are included in the equations through higher-order mixed
and spatial derivatives. Eskilsson and Sherwin presented the first spectral/hp
model for 1D Boussinesq-type equations [103], and later in 2D [104]. In the
latter study was illustrated that the use of high-order schemes is
computationally efficient compared to low-order methods for long-time
integration. This is general knowledge, mathematically proven in [105], but
the importance is accenturated in wave equation where the numerical dispersion
error must be kept small compared to the physical dispersion terms. Additional
spectral/hp studies of weakly dispersive Boussinesq equations are [106, 107].
Engsig-Karup et al. [108, 109, 110] solved a set highly dispersive and
nonlinear Boussinesq-type equations. More recently the interest has turned
toward Green-Naghdi type equations due to their description of nonlinearity
and a number of spectral/hp element methods have been put forward [111, 112,
113]. We note that all the above models are based on high-order discontinuous
Galerkin (DG) methods.

The classical Boussinesq equations due to Peregrine [114] read

$\displaystyle\frac{\partial\eta}{\partial
t}+\bar{\nabla}\cdot\left(\left(h+\eta\right)\bar{\mathbf{u}}\right)=0\,,$
(21) $\displaystyle\frac{\partial\bar{\mathbf{u}}}{\partial
t}+g\bar{\nabla}\eta+\bar{\mathbf{u}}\cdot\bar{\nabla}\bar{\mathbf{u}}+\frac{h^{2}}{6}\bar{\nabla}\left(\bar{\nabla}\cdot\frac{\partial\mathbf{u}}{\partial
t}\right)-\frac{h^{2}}{2}\bar{\nabla}\left(\bar{\nabla}\cdot\frac{\partial(h\mathbf{u})}{\partial
t}\right)=0\,,$ (22)

where $\bar{\mathbf{u}}$ is the horizontal velocities. The exist a pre-written
solver for the Peregrine equations in Nektar++. The equations are solved in
terms of conservative variables and the solution scheme is using the wave-
continuity step as introduced in [104]. First the advection part is
approximated and used in an intermediate step. In the intermediate step the
dispersive terms are reformulated into a Helmholtz equation in terms of the
auxiliary variable $q=\bar{\nabla}\cdot\partial_{t}(d\mathbf{u})$, in which
$d=h+\eta$. Subsequently, a projection back to conservative variables is
performed. The Boussinesq equations can be used for computing the generation
of higher harmonics. A classical example is the semi-circular shoal case, see
Fig. 19.

(a)

(b)

Figure 19: Higher harmonics generation on a semi-circular shoal. (a) 3D
surface view and (b) computed and measured harmonics through the centerline
(from the reference [115])

With regard to non-dispersive, hydrostatic wave equations the shallow water
equations (SWE) are typically used in marine applications to predict tidal
forcing but can also be used for other long wave phenomena such as tsunami and
storm-surge predictions. For tsunami modelling hp-adaptive DG models have been
presented in order to resolve the shock wave propagation in detail [116, 117,
118]. However, the main bulk of work on spectral/hp element methods for the
SWE are linked to numerical weather prediction or ocean circulation where the
SWE is used as a stepping stone towards solving the 3D primitive equations,
see e.g. [119, 120, 121, 122, 123, 124, 125].

The SWE is solved in Nektar++ as a sub-model for the dispersive Boussinesq
solver simply by ignoring the dispersive step. Recently a version of the SWE
model capable of solving the flow on curved surfaces using Nektar++ was
presented in [126]. The model uses the so-called method of moving frames [127]
to express the SWE on curved and rotating geometries. The model was shown to
exhibit exponential convergence and handled the standard test cases for
spherical SWE, see Fig. 20 for the flow around an isolated mountain.

(a)

(b)

(c)

Figure 20: Spherical shallow water flow over an isolated mountain. (a) mesh,
(b) initial surface elevation and (c) surface elevation after 15 days. (from
the reference [126])

##  5 Future directions and perspectives

Significant progress has been made in theoretical and applied aspects of the
spectral/hp element method over the past two decades which have enabled its
application to a wide range of challenging hydrodynamics and broader
industrial problems. However, there still remain a number of challenges [55].
In order to enable the spectral/hp element method to be used to address
difficult, large-scale, high Reynolds number flow problems in complex
geometries, we see a number of areas which still require improvement:

  1. 1.

High-order mesh generation: Robustly generating high-order meshes with CAD
conforming representation of the underlying geometry remains one of the key
challenges to industrial application.

  2. 2.

One of the major benefits of the spectral/hp element method is the
opportunities presented for optimising the computational discretisation in
accordance with the physics being modelled through hp adaptivity. p-adaptivity
and r-adaptivity have recently been explored [128, 129] but further effort is
needed in identifying the optimal refinement algorithms to choose to optimise
the computational cost while also improving the numerical accuracy of
solutions.

  3. 3.

The stabilisation techniques discussed in §2, such as SVV and projection-based
filtering techniques, require further improvement through further calibration
for marginal/under-resolved computations. As mentioned by Sagaut [37], “some
specific numerical stabilisation procedures can be defined, which tune the
numerical dissipation in such a way that the results remain sensitive to
subgrid modelling”. In terms of the similarity between artificial dissipation
and the direct energy cascade model, some numerical stabilisation techniques
can be used to perform “no-model” large eddy simulations (iLES). However,
Sagaut stated that “certain studies have shown that, for coarse grids, i.e.
high values of the numerical cutoff length, increasing the order of accuracy
of the upwind scheme can lead to a degradation of the results [130]”. The
proposed ‘DGKernel’, which emulates the upwind properties of a DG scheme,
still requires more detailed calibration for a range of high- Reynolds fluid
flows.

  4. 4.

Verification and Validation of spectral/hp element simulations for under-
resolved computations: The use of stabilisation techniques, such as SVV,
allows us to interpret the under-resolved numerical calculations as implicit
large eddy simulation. However, we need to continue to investigate how the
stabilisation methods and numerical properties interact to ensure the benefits
of high-precision are not destroyed through artificial/numerical pollution.

In the context of hydrodynamics applications, there is substantial opportunity
for gaining an improved understanding of a wide range of physical mechanisms
through the use of the spectral/hp element method. In particular:

  1. 1.

Tip vortex modelling: As indicated in [7], without the accurate modelling of
the three dimensional boundary layer, the developing vortex remained
challenging to compute accurately, even for the advanced RANS models
correcting for the high degree of curvature in the flow. In addition,
accurately capturing the low-pressure region within the vortex core and
sustaining this low pressure even just one chord length downstream of the
trailing edge is particularly challenging. As demonstrated [7], the SVV-based
iLES method has been shown to be potential alternative for computing complex
unsteady vortex dominated flows. Simulating vortical flows using the
spectral/hp element method with the concept of the SVV-based iLES can be
regarded as an interesting alternative method in hydrodynamics.

  2. 2.

Implicit large eddy simulation: A notable effort has been made to refine the
SVV-based strategy for under-resolved simulations of high Reynolds fluid flows
[38, 44]. This strategy has great potential to provide high fidelity iLES
computations for a broad range of applications in hydrodynamics.

  3. 3.

Fluid-structure interaction: As indicated in [131, 132, 133, 134, 135],
spectral elements offer superior wave propagation capabilities, which has been
used for simulations of fluid-structure interaction (e.g.[133, 134, 136]).
Interaction subject to surface waves can be attractive

  4. 4.

Cavitating flow: Recently, discontinuous spectral/hp element method has been
implemented for the efficient and accurate simulation of multiphase flows
[137]. The compressible Navier-Stokes equations, coupled with a highly
realistic equation of state, are solved by a discontinuous Galerkin spectral
element method for flows with cavitation. In Nektar++, the discontinuous
Galerkin method and the flux reconstruction approach have been used to solve
the compressible Navier-Stokes equations. This therefore provides a basis to
apply the spectral/hp element method in applications of cavitating flows.

  5. 5.

Free surface flows: As is well known free surface flows can be modelled as a
depth average “shallow-water” approach or using ALE moving mesh techniques.
The spectral/hp element method has been used to treat these two kinds of free
surface flows [97, 138, 139, 140, 141]). For the second free surface problem,
discontinuous spectral/hp element methods have a great potential in
applications (e.g. [139, 140, 142, 143, 144]), especially when h/p-adaptivity
is employed.

##  6 Acknowledgments

The authors kindly thank the Executive Editor-in-Chief Prof. Liandi Zhou for
the invitation to contribute this review article and Dr. Wei Zhang for his
contribution to conducting the tip-vortex simulation. H.X. and S.J.S would
like to acknowledge support under EPSRC grant EP/L000407/1.

## References

## References

  * [1] I. Babuška, M. Suri, The p- and h-p versions of the finite element method, an overview, Computer Methods in Applied Mechanics and Engineering 80 (1) (1990) 5–26.
  * [2] I. Babuška, M. Suri, The p and h-p versions of the finite element method, basic principles and properties, SIAM Rev. 36 (4) (1994) 578–632.
  * [3] G. Karniadakis, S. Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics, second ed., Oxford University Press, 2005.
  * [4] C. D. Cantwell, D. Moxey, A. Comerford, A. Bolis, G. Rocco, G. Mengaldo, D. D. Grazia, S. Yakovlev, J.-E. Lombard, D. Ekelschot, B. Jordi, H. Xu, Y. Mohamied, C. Eskilsson, B. Nelson, P. Vos, C. Biotto, R. M. Kirby, S. J. Sherwin, Nektar++: An open-source spectral/hp element framework, Computer Physics Communications 192 (2015) 205–219.
  * [5] S. Yakovlev, D. Moxey, R. M. Kirby, S. J. Sherwin, To CG or to HDG: A comparative study in 3d, Journal of Scientific Computing 67 (1) (2016) 192–220.
  * [6] J. King, S. Yakovlev, Z. S. Fu, R. M. Kirby, S. J. Sherwin, Exploiting batch processing on streaming architectures to solve 2d elliptic finite element problems: A hybridized discontinuous galerkin (hdg) case study, Journal of Scientific Computing 60 (2) (2014) 457–482.
  * [7] J. E. W. Lombard, D. Moxey, S. J. Sherwin, J. F. A. Hoessler, S. Dhandapani, M. J. Taylor, Implicit large-eddy simulation of a wingtip vortex, AIAA Journal 54 (2) (2016) 506–518.
  * [8] D. Serson, J. R. Meneghini, S. J. Sherwin, Velocity-correction schemes for the incompressible navier-stokes equations in general coordinate systems, Journal of Computational Physics 316 (2016) 243–254.
  * [9] D. Serson, J. R. Meneghini, S. J. Sherwin, Direct numerical simulations of the flow around wings with spanwise waviness at a very low reynolds number, Computers & Fluids 146 (2017) 117–124.
  * [10] H. Xu, J.-E. W. Lombard, S. J. Sherwin, Influence of localised smooth steps on the instability of a boundary layer, Journal of Fluid Mechanics 817 (2017) 138–170.
  * [11] H. Xu, S. M. Mughal, E. R. Gowree, C. J. Atkin, S. J. Sherwin, Destabilisation and modification of tollmien–schlichting disturbances by a three-dimensional surface indentation, Journal of Fluid Mechanics 819 (2017) 592–620.
  * [12] C. D. Cantwell, S. Yakovlev, R. M. Kirby, N. S. Peters, S. J. Sherwin, High-order spectral/hp element discretisation for reaction-diffusion problems on surfaces: Application to cardiac electrophysiology, Journal of Computational Physics 257 (Part A) (2014) 813–829.
  * [13] B. Szabó, I. Babuka, Finite element analysis, Wiley, 1991.
  * [14] D. Moxey, D. Ekelschot, Ü. Keskin, S. J. Sherwin, J. Peiró, A thermo-elastic analogy for high-order curvilinear meshing with control of mesh validity and quality, Procedia Engineering 82 (2014) 127–135, 23rd International Meshing Roundtable (IMR23).
  * [15] A. Comerford, K. Y. Chooi, M. Nowak, P. D. Weinberg, S. J. Sherwin, A combined numerical and experimental framework for determining permeability properties of the arterial media, Biomechanics and Modeling in Mechanobiology 14 (2) (2015) 297–313.
  * [16] C. Eskilsson, S. J. Sherwin, Spectral/hp discontinuous galerkin methods for modelling 2d boussinesq equations, Journal of Computational Physics 212 (2) (2006) 566 – 589.
  * [17] M. G. Duffy, Quadrature over a pyramid or cube of integrands with a singularity at a vertex, SIAM Journal on Numerical Analysis 19 (6) (1982) 1260 – 1262.
  * [18] J. S. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods, Springer-Verlag New York, 2008.
  * [19] W. H. Reed, T. R. Hill, Triangular mesh methods for the neutron transport equation, Tech. rep., Los Alamos National Laboratory, lA-UR-73-479 (1973).
  * [20] R. A. Adams, Sobolev spaces, Academic Press, New York, 1975.
  * [21] B. Vermeire, F. Witherden, P. Vincent, On the utility of gpu accelerated high-order methods for unsteady flow simulations: A comparison with industry-standard tools, Journal of Computational Physics 334 (2017) 497 – 521\.
  * [22] D. Gottlieb, S. A. Orszag, Numerical analysis of spectral methods: theory and applications, Society for Industrial and Applied Mathematics, 1977.
  * [23] R. M. Kirby, G. E. Karniadakis, De-aliasing on non-uniform grids: algorithms and applications, Journal of Computational Physics 191 (1) (2003) 249 – 264.
  * [24] R. M. Kirby, S. J. Sherwin, Aliasing errors due to quadratic nonlinearities on triangular spectral /hp element discretisations, Journal of Engineering Mathematics 56 (3) (2006) 273–288.
  * [25] J. Malm, P. Schlatter, D. S. Henningson, Coherent structures and dominant frequencies in a turbulent three-dimensional diffuser, Journal of Fluid Mechanics 699 (2012) 320–351.
  * [26] G. Mengaldo, D. D. Grazia, D. Moxey, P. E. Vincent, S. J. Sherwin, Dealiasing techniques for high-order spectral element methods on regular and irregular grids, Journal of Computational Physics 299 (2015) 56 – 81.
  * [27] B. Cockburn, S. C. Hou, C. W. Shu, The runge-kutta local projection discontinuous galerkin finite element method for conservation laws. iv. the multidimensional case, Mathematics of Computation 54 (1990) 545 – 581.
  * [28] E. Tadmor, Convergence of spectral methods for nonlinear conservation laws, SIAM Journal on Numerical Analysis 26 (1) (1989) 30 – 44.
  * [29] G. S. Karamanos, G. E. Karniadakis, A spectral vanishing viscosity method for large-eddy simulations, Journal of Computational Physics 163 (1) (2000) 22 – 50\.
  * [30] R. M. Kirby, S. J. Sherwin, Stabilisation of spectral/hp element methods through spectral vanishing viscosity: Application to fluid mechanics modelling, Computer Methods in Applied Mechanics and Engineering 195 (23) (2006) 3128 – 3144.
  * [31] R. Pasquetti, Spectral vanishing viscosity method for les: sensitivity to the svv control parameters, Journal of Turbulence 6 (12) (2005) 1 – 14.
  * [32] R. Pasquetti, Spectral vanishing viscosity method for large-eddy simulation of turbulent flows, Journal of Scientific Computing 27 (1) (2006) 365–375.
  * [33] R. M. Kirby, G. E. Karniadakis, Coarse resolution turbulence simulations with spectral vanishing viscosity—large-eddy simulations (SVV-LES), Journal of Fluids Engineering 124 (4) (2002) 886 – 891.
  * [34] C. Xu, Stabilization methods for spectral element computations of incompressible flows, Journal of Scientific Computing 27 (1) (2006) 495 – 505\.
  * [35] E. Severac, E. Serre, A spectral vanishing viscosity for the les of turbulent flows within rotating cavities, Journal of Computational Physics 226 (2) (2007) 1234 – 1255.
  * [36] K. Koal, J. Stiller, H. Blackburn, Adapting the spectral vanishing viscosity method for large-eddy simulations in cylindrical configurations, Journal of Computational Physics 231 (8) (2012) 3389 – 3405.
  * [37] P. Sagaut, Large eddy simulation for incompressible flows, Springer-Verlag Berlin Heidelberg, 2001.
  * [38] R. C. Moura, S. J. Sherwin, J. Peiró, Eigensolution analysis of spectral/hp continuous galerkin approximations to advection-diffusion problems: Insights into spectral vanishing viscosity, Journal of Computational Physics 307 (2016) 401 – 422.
  * [39] Y. Maday, S. M. O. Kaber, E. Tadmor, Legendre pseudospectral viscosity method for nonlinear conservation laws, SIAM Journal on Numerical Analysis 30 (2) (1993) 321 – 342.
  * [40] S. Collis, Discontinuous galerkin methods for turbulence simulation, in: Proceedings of the 2002 Center for Turbulence Research Summer Program, 2002, pp. 155 – 167.
  * [41] S. Collis, K. Ghayour, Discontinuous galerkin methods for compressible dns, in: Proceedings of the 4th ASME/JSME Joint Fluids Summer Engineering Conference, 2003, pp. 1777 – 1786.
  * [42] L. Wei, A. Pollard, Direct numerical simulation of compressible turbulent channel flows using the discontinuous galerkin method, Computers & Fluids 47 (1) (2011) 85 – 100.
  * [43] J.-B. Chapelier, M. de la Llave Plata, F. Renac, E. Lamballais, Evaluation of a high-order discontinuous galerkin method for the dns of turbulent flows, Computers & Fluids 95 (2014) 210 – 226.
  * [44] R. C. Moura, S. J. Sherwin, J. Peiró, Linear dispersion-diffusion analysis and its application to under-resolved turbulence simulations using discontinuous galerkin spectral/hp methods, Journal of Computational Physics 298 (2015) 695 – 710.
  * [45] A. Bolis, C. D. Cantwell, D. Moxey, D. Serson, S. J. Sherwin, An adaptable parallel algorithm for the direct numerical simulation of incompressible turbulent flows using a fourier spectral/hp element method and mpi virtual topologies, Computer Physics Communications 206 (2016) 17–25.
  * [46] P. E. Vos, S. J. Sherwin, R. M. Kirby, From h to p efficiently: Implementing finite and spectral/hp element methods to achieve optimal performance for low- and high-order discretisations, Journal of Computational Physics 229 (13) (2010) 5161 – 5181.
  * [47] C. D. Cantwell, S. J. Sherwin, R. M. Kirby, P. H. J. Kelly, From h to p efficiently: Strategy selection for operator evaluation on hexahedral and tetrahedral elements, Computers & Fluids 43 (1) (2011) 23–28.
  * [48] C. D. Cantwell, S. J. Sherwin, R. M. Kirby, P. H. J. Kelly, From h to p efficiently: selecting the optimal spectral/hp discretisation in three dimensions, Mathematical Modelling of Natural Phenomena 6 (3) (2011) 84–96.
  * [49] A. Bolis, C. D. Cantwell, R. M. Kirby, S. J. Sherwin, From h to p efficiently: optimal implementation strategies for explicit time-dependent problems using the spectral/hp element method, International journal for numerical methods in fluids 75 (8) (2014) 591–607.
  * [50] D. Moxey, C. D. Cantwell, R. M. Kirby, S. J. Sherwin, Optimising the performance of the spectral/hp element method with collective linear algebra operations, Computer Methods in Applied Mechanics and Engineering 310 (Supplement C) (2016) 628–645.
  * [51] D. Moxey, M. D. Green, S. J. Sherwin, J. Peiró, An isoparametric approach to high-order curvilinear boundary-layer meshing, Computer Methods in Applied Mechanics and Engineering 283 (2015) 636 – 650.
  * [52] S. J. Sherwin, J. Peiró, Mesh generation in curvilinear domains using high-order elements, International Journal for Numerical Methods in Engineering 53 (1) (2002) 207–223.
  * [53] M. Turner, J. Peiró, D. Moxey, A variational framework for high-order mesh generation, Procedia Engineering 163 (2016) 340 – 352, 25th International Meshing Roundtable.
  * [54] D. Moxey, D. Ekelschot, Ü. Keskin, S. J. Sherwin, J. Peiró, High-order curvilinear meshing using a thermo-elastic analogy, Computer-Aided Design 72 (2016) 130 – 139, 23rd International Meshing Roundtable Special Issue: Advances in Mesh Generation.
  * [55] Z. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert, H. Huynh, N. Kroll, G. May, P.-O. Persson, B. van Leer, M. Visbal, High-order cfd methods: current status and perspective, International Journal for Numerical Methods in Fluids 72 (8).
  * [56] S. Dey, R. M. O’Bara, M. S. Shephard, Curvilinear mesh generation in 3d, in: Proceedings of the 8th International Meshing Roundtable, Lake Tahoe, California, 1999.
  * [57] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331.
  * [58] G. E. Karniadakis, M. Israeli, S. A. Orszag, High-order splitting methods for the incompressible navier-stokes equations, Journal of Computational Physics 97 (2) (1991) 414 – 443.
  * [59] P. E. Vos, C. Eskilsson, A. Bolis, S. Chun, R. M. Kirby, S. J. Sherwin, A generic framework for time-stepping partial differential equations (pdes): general linear methods, object-oriented implementation and application to fluid problems, International Journal of Computational Fluid Dynamics 25 (3) (2011) 107–125.
  * [60] P. F. Fischer, E. M. Rønquist, Spectral element methods for large scale parallel navier-stokes calculations, Computer Methods in Applied Mechanics and Engineering 116 (1) (1994) 69 – 76.
  * [61] S. J. Sherwin, G. E. Karniadakis, A triangular spectral element method; applications to the incompressible navier-stokes equations, Computer Methods in Applied Mechanics and Engineering 123 (1) (1995) 189 – 229.
  * [62] S. J. Sherwin, G. E. Karniadakis, Tetrahedralhpfinite elements: Algorithms and flow simulations, Journal of Computational Physics 124 (1) (1996) 14 – 45.
  * [63] P. F. Fischer, An overlapping schwarz method for spectral element solution of the incompressible navier-stokes equations, Journal of Computational Physics 133 (1) (1997) 84 – 101.
  * [64] A. G. Tomboulides, S. A. Orszag, Numerical investigation of transitional and weak turbulent flow past a sphere, Journal of Fluid Mechanics 416 (2000) 45 – 73.
  * [65] P. Fischer, J. Mullen, Filter-based stabilization of spectral element methods, Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 332 (3) (2001) 265 – 270.
  * [66] C. E. Wasberg, T. Gjesdal, B. A. P. Reif, Ø. Andreassen, Variational multiscale turbulence modelling in a high order spectral element method, Journal of Computational Physics 228 (19) (2009) 7333 – 7356.
  * [67] H. Blackburn, S. Schmidt, Spectral element filtering techniques for large eddy simulation with dynamic estimation, Journal of Computational Physics 186 (2) (2003) 610 – 629.
  * [68] T. Iliescu, P. F. Fischer, Large eddy simulation of turbulent channel flows by the rational large eddy simulation model, Physics of Fluids 15 (10) (2003) 3036–3047.
  * [69] J. Ohlsson, P. Schlatter, P. F. Fischer, D. S. Henningson, Stabilization of the Spectral-Element Method in Turbulent Flow Simulations, Springer Berlin Heidelberg, Berlin, Heidelberg, 2011, pp. 449–458.
  * [70] A. D. Beck, T. Bolemann, D. Flad, H. Frank, G. J. Gassner, F. Hindenlang, C.-D. Munz, High-order discontinuous galerkin spectral element methods for transitional and turbulent flow simulations, International Journal for Numerical Methods in Fluids 76 (8) (2014) 522–548.
  * [71] H. M. Blackburn, S. J. Sherwin, Formulation of a galerkin spectral element-fourier method for three-dimensional incompressible flows in cylindrical geometries, Journal of Computational Physics 197 (2) (2004) 759 – 778.
  * [72] A. T. Patera, A spectral element method for fluid dynamics: Laminar flow in a channel expansion, Journal of Computational Physics 54 (3) (1984) 468 – 488.
  * [73] R. Vinuesa, P. Schlatter, J. Malm, C. Mavriplis, D. S. Henningson, Direct numerical simulation of the flow around a wall-mounted square cylinder under various inflow conditions, Journal of Turbulence 16 (6) (2015) 555–587.
  * [74] L. T. Diosady, S. M. Murman, DNS of Flows over Periodic Hills using a Discontinuous Galerkin Spectral-Element Method, American Institute of Aeronautics and Astronautics, 2014.
  * [75] H. Xu, S. J. Sherwin, P. Hall, X. S. Wu, The behaviour of tollmien-schlichting waves undergoing small-scale localised distortions, J. Fluid Mech. 792 (2016) 499–525.
  * [76] D. Serson, J. R. Meneghini, S. J. Sherwin, Direct numerical simulations of the flow around wings with spanwise waviness, Journal of Fluid Mechanics 826 (2017) 714 – 731.
  * [77] D. Barkley, H. M. Blackburn, S. J. Sherwin, Direct optimal growth analysis for timesteppers, International Journal for Numerical Methods in Fluids 57 (9) (2008) 1435 – 1458.
  * [78] H. Blackburn, D. Barkley, S. J. Sherwin, Convective instability and transient growth in flow over a backward-facing step, Journal of Fluid Mechanics 603 (2008) 271–304.
  * [79] H. Blackburn, S. J. Sherwin, D. Barkley, Convective instability and transient growth in steady and pulsatile stenotic flows, Journal of Fluid Mechanics 607 (2008) 267–277.
  * [80] C. D. Cantwell, D. Barkley, H. M. Blackburn, Transient growth analysis of flow through a sudden expansion in a circular pipe, Physics of Fluids 22 (3) (2010) 034101.
  * [81] C. D. Cantwell, D. Barkley, Computational study of subcritical response in flow past a circular cylinder, Physical review E 82 (2) (2010) 026315.
  * [82] M. S. Broadhurst, V. Theofilis, S. J. Sherwin, Spectral element stability analysis of vortical flows, Springer Netherlands, Dordrecht, 2006, pp. 153 – 158\.
  * [83] G. Rocco, Advanced instability methods using spectral/hp discretisations and their applications to complex geometries, Ph.D. thesis, Department of Aeronautics, Imperial College London (2014).
  * [84] B. Jordi, Steady-state solvers for stability analysis of vortex dominated flows, Ph.D. thesis, Department of Aeronautics, Imperial College London (2015).
  * [85] W. J. Devenport, M. C. Rife, S. I. Liapis, G. J. Follin, The structure and development of a wing-tip vortex, Journal of Fluid Mechanics 312 (1996) 67 – 106\.
  * [86] R. E. Arndt, Cavitation in vortical flows, Annual Review of Fluid Mechanics 34 (1) (2002) 143 – 175.
  * [87] C. W. Pitt Ford, H. Babinsky, Lift and the leading-edge vortex, Journal of Fluid Mechanics (2013) 280 – 313.
  * [88] T. Leweke, S. L. Dizès, C. H. Williamson, Dynamics and instabilities of vortex pairs, Annual Review of Fluid Mechanics 48 (1) (2016) 507 – 541.
  * [89] J. Feys, S. A. Maslowe, Elliptical instability of the moore–saffman model for a trailing wingtip vortex, Journal of Fluid Mechanics 803 (2016) 556–590.
  * [90] R. C. Moura, G. Mengaldo, J. Peiró, S. J. Sherwin, On the eddy-resolving capability of high-order discontinuous galerkin approaches to implicit les / under-resolved dns of euler turbulence, Journal of Computational Physics 330 (2017) 615 – 623.
  * [91] C. Eskilsson, A. P. Engsig-Karup, S. J. Sherwin, J. S. Hesthaven, L. Bergdahl, The next step in coastal numerical models: spectral/hp element methods?, in: Proceedings of the WAVES2005 Conference, Madrid, 2005.
  * [92] M. Brocchini, A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics, Proc. Roy. Soc. London Ser. A 469 (2160) (2013) 1–27.
  * [93] G. X. Wu, R. E. Taylor, Finite element analysis of two-dimensional non-linear transient water waves, Applied Ocean Research.
  * [94] X. Cai, H. P. Langtangen, B. F. Nielsen, A. Tveito, A finite element method for fully nonlinear water waves, J. Comput. Phys. 143 (1998) 544–568.
  * [95] Q. W. Ma, G. X.Wu, R. E. Taylor, Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves. part 1: methodology and numerical procedure., International Journal for Numerical Methods in Fluids 36 (3) (2001) 265–285.
  * [96] Q. W. Ma, S. Yan., Quasi ALE finite element method for nonlinear water waves, Journal of Computational Physics 212 (1) (2006) 52 –72.
  * [97] I. Robertson, S. J. Sherwin, Free-surface flow simulation using $hp$/spectral elements, Journal of Computational Physics 155 (1999) 26–53.
  * [98] A. P. Engsig-Karup, C. Eskilsson, D. Bigoni, A stabilised nodal spectral element method for fully nonlinear water waves, Journal of Computational Physics 318 (2016) 1–21.
  * [99] A. P. Engsig-Karup, C. Eskilsson, D. Bigoni, Unstructured spectral element model for dispersive and nonlinear wave propagation, in: Proceedings of The Twenty-sixth (2016) International Ocean and Polar Engineering Conference, Vol. 3, 2016, pp. 661–669.
  * [100] V. E. Zakharov, Stability of periodic waves of finite amplitude on the surface of a deep fluid, Journal of Applied Mechanics and Technical Physics.
  * [101] A. Mieritz, A. Engsig-Karup, C. Eskilsson, D. Moxely, S. J. Sherwin, Large-scale water wave simulations using a spectral/$hp$ element fully nonlinear potential flow model, manuscript in preparation.
  * [102] A. P. Engsig-Karup, C. Monteserin, C. Eskilsson, A stabilised nodal spectral element method for fully nonlinear water waves, Part 2: Wave-body interaction, arXiv:1703.09697.
  * [103] C. Eskilsson, S. J. Sherwin, A discontinuous spectral element model for Boussinesq-type equations, Journal of Scientific Computing 17 (1-4) (2002) 143–152.
  * [104] C. Eskilsson, S. J. Sherwin, Discontinuous Galerkin spectral/hp element modelling of dispersive shallow water systems, Journal of Scientific Computing 22 (2005) 269–288.
  * [105] H.-O. Kreiss, J. Oliger, Comparison of accurate methods for the integration of hyperbolic equations, Tellus 24 (1972) 199–215.
  * [106] C. Eskilsson, S. J. Sherwin, L. Bergdahl, An unstructured spectral/hp element model for enhanced Boussinesq-type equations, Coastal Engineering.
  * [107] M. Dumbser, M. Facchini, A space-time discontinuous Galerkin method for Boussinesq-type equations, Applied Mathematics and Computation 272 (2016) 336–346.
  * [108] A. P. Engsig-Karup, Unstructured nodal DG-FEM solution of high-order boussinesq-type equations, Ph.D. thesis, PhD. Thesis. Department of Mechanical Engineering, Technical University of Denmark (2006).
  * [109] A. P. Engsig-Karup, J. S. Hesthaven, H. B. Bingham, P. Madsen, Nodal DG-FEM solutions of high-order Boussinesq-type equations, Journal of Engineering Mathematics 56 (2006) 351–370.
  * [110] A. P. Engsig-Karup, J. S. Hesthaven, H. B. Bingham, T. Warburton, DG-FEM solution for nonlinear wave-structure interaction using boussinesq-type equations, Coastal Engineering 55 (2008) 197–208.
  * [111] N. Panda, C. Dawson, Y. Zhang, A. Kennedy, J. Westerink, A. Donahue, Discontinuous Galerkin methods for solving Boussinesq–Green–Naghdi equations in resolving non-linear and dispersive surface water waves, Journal of Computational Physics 273 (2014) 572–588.
  * [112] A. Duran, F. Marche, Discontinuous-Galerkin discretization of a new class of Green-Naghdi equations, Communications in Computational Physics 17 (3) (2015) 721–760.
  * [113] A. Samii, C. Dawson, An explicit hybridized discontinuous Galerkin method for Serre-Green-Naghdi wave model, Computer Methods in Applied Mechanics and Engineering 330 (2018) 447–470.
  * [114] D. Peregrine, Long waves on a beach, Journal of Fluid Mechanics 27 (1967) 815–827.
  * [115] C. Eskilsson, S. J. Sherwin, Spectral/hp discontinuous Galerkin methods for modelling 2D Boussinesq equations, Journal of Computational Physics 212 (2006) 566–589.
  * [116] S. Blaise, A. St-Cyr, A dynamic hp-adaptive discontinuous Galerkin method for shallow-water flows on the sphere with application to a global tsunami simulation, Monthly Weather Review 140 (2012) 978–996.
  * [117] S. Blaise, A. St-Cyr, D. Mavriplis, B. Lockwood, Discontinuous Galerkin unsteady discrete adjoint method for real-time efficient tsunami simulations, Journal of Computational Physics 232 (2013) 416–430.
  * [118] B. Bonev, J. Hesthaven, F. Giraldo, M. Kopera, Discontinuous Galerkin scheme for the spherical shallow water equations with applications to tsunami modeling and prediction, EPFL-ARTICLE-232449.
  * [119] H. Ma, A spectral element basin model for the shallow water equations, Journal of Computational Physics 109 (1993) 133–149.
  * [120] M. Iskandarani, D. Haidvogel, J. Boyd, A staggered spectral element model with application to the oceanic shallow water equations, International Journal for Numerical Methods in Fluids 20 (1995) 393–414.
  * [121] M. Taylor, J. Tribbia, M. Iskandarani, The spectral element method for the shallow water equations on the sphere, Journal of Computational Physics 130 (1997) 92–108.
  * [122] F. X. Giraldo, J. S. Hesthaven, T. Warburton, Nodal high-order Discontinuous galerkin methods for the spherical shallow water equations, Journal of Computational Physics 181 (2002) 499–525.
  * [123] R. Nair, S. Thomas, R. Loft, A discontinuous Galerkin global shallow water model, Montly Weather Review 133 (2005) 876–888.
  * [124] M. Lauter, F. Giraldo, D. Handorf, K. Dethloff, A discontinuous galerkin method for the shallow water equations in spherical triangular coordinates, Journal of Computational Physics 227 (2008) 10226–10242.
  * [125] L. Bao, R. Nair, H. Tufo, A mass and momentum flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere, Journal of Computational Physics 271 (2014) 224–243.
  * [126] S. Chun, C. Eskilsson, Method of moving frames to solve the shallow water equations on arbitrary rotating curved surfaces, Journal of Computational Physics 333 (2017) 1–23.
  * [127] S. Chun, Method of moving frames to solve conservation laws on curved surfaces, Journal of Scientific Computing 53 (2) (2012) 268–294.
  * [128] D. Moxey, C. D. Cantwell, G. Mengaldo, D. Serson, D. Ekelschot, J. Peiró, S. J. Sherwin, R. M. Kirby, Towards p-adaptive spectral/hp element methods for modelling industrial flows, in: Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2016, Springer, 2017, pp. 63–79.
  * [129] D. Ekelschot, D. Moxey, S. J. Sherwin, J. Peiró, A p-adaptation method for compressible flow problems using a goal-based error indicator, Computers & Structures 181 (2017) 55–69.
  * [130] D. Tafti, Comparison of some upwind-biased high-order formulations with a second-order central-difference scheme for time integration of the incompressible navier-stokes equations, Computers & Fluids 25 (7) (1996) 647 – 665.
  * [131] E. Priolo, J. M. Carcione, G. Seriani, Numerical simulation of interface waves by high-order spectral modeling techniques, The Journal of the Acoustical Society of America 95 (2) (1994) 681–693.
  * [132] M. A. Sprague, T. L. Geers, A spectral-element method for modelling cavitation in transient fluid–structure interaction, International Journal for Numerical Methods in Engineering 60 (15) (2004) 2467–2499.
  * [133] N. Bodard, M. O. Deville, Fluid-structure interaction by the spectral element method, Journal of Scientific Computing 27 (1) (2006) 123–136.
  * [134] M. A. Sprague, T. L. Geers, A spectral-element/finite-element analysis of a ship-like structure subjected to an underwater explosion, Computer Methods in Applied Mechanics and Engineering 195 (17) (2006) 2149 – 2167.
  * [135] M. A. Sprague, T. L. Geers, Spectral elements and field separation for an acoustic fluid subject to cavitation, Journal of Computational Physics 184 (1) (2003) 149 – 162.
  * [136] R. van Loon, P. D. Anderson, F. N. van de Vosse, S. J. Sherwin, Comparison of various fluid-structure interaction methods for deformable bodies, Computers & Structures 85 (11) (2007) 833 – 843, fourth MIT Conference on Computational Fluid and Solid Mechanics.
  * [137] M. Hoffmann, S. Boblest, P. Offenhäuser, C.-D. Munz, A robust high-order discontinuous galerkin solver for fluid flow with cavitation, in: Ninth International Conference on ICCFD9-196 Computational Fluid Dynamics (ICCFD9), Istanbul, Turkey, 2016.
  * [138] H. Lee-Wing, A. T. Patera, A legendre spectral element method for simulation of unsteady incompressible viscous free-surface flows, Computer Methods in Applied Mechanics and Engineering 80 (1) (1990) 355 – 366.
  * [139] M. Sussman, M. Y. Hussaini, A discontinuous spectral element method for the level set equation, Journal of Scientific Computing 19 (1) (2003) 479–500.
  * [140] C.-L. Lin, H. Lee, T. Lee, L. J. Weber, A level set characteristic galerkin finite element method for free surface flows, International Journal for Numerical Methods in Fluids 49 (5) (2005) 521–547.
  * [141] R. Bouffanais, M. O. Deville, Mesh update techniques for free-surface flow solvers using spectral element method, Journal of Scientific Computing 27 (1) (2006) 137–149.
  * [142] J. Grooss, J. Hesthaven, A level set discontinuous galerkin method for free surface flows, Computer Methods in Applied Mechanics and Engineering 195 (25) (2006) 3406 – 3429, discontinuous Galerkin Methods.
  * [143] E. Marchandise, J.-F. Remacle, A stabilized finite element method using a discontinuous level set approach for solving two phase incompressible flows, Journal of Computational Physics 219 (2) (2006) 780 – 800.
  * [144] V.-T. Nguyen, J. Peraire, B. C. Khoo, P.-O. Persson, A discontinuous galerkin front tracking method for two-phase flows with surface tension, Computers & Fluids 39 (1) (2010) 1 – 14.