# On the performance of a highly-scalable Computational Fluid Dynamics code on
AMD, ARM and Intel processors

Pablo Ouro<EMAIL_ADDRESS>School of Mechanical, Aerospace and
Civil Engineering, University of Manchester, Manchester, M13 9PL, UK  Hydro-
environmental Research Centre, School of Engineering, Cardiff University, CF24
3AA Cardiff, UK  Unai Lopez-Novoa<EMAIL_ADDRESS>Department of Computer
Languages and Systems, University of the Basque Country, Spain  Martyn F.
Guest<EMAIL_ADDRESS>Advanced Research Computing, Cardiff University,
Cardiff, UK

###### Abstract

No area of computing is hungrier for performance than High Performance
Computing (HPC), the demands of which continue to be a major driver for
processor performance and adoption of accelerators, and also advances in
memory, storage, and networking technologies. A key feature of the Intel
processor domination of the past decade has been the extensive adoption of
GPUs as coprocessors, whilst more recent developments have seen the increased
availability of a number of CPU processors, including the novel ARM-based
chips. This paper analyses the performance and scalability of a state-of-the-
art Computational Fluid Dynamics (CFD) code on three HPC cluster systems
equipped with AMD EPYC-Rome (EPYC, 4096 cores), ARM-based Marvell ThunderX2
(TX2, 8192 cores) and Intel Skylake (SKL, 8000 cores) processors. Three
benchmark cases are designed with increasing computation-to-communication
ratio and numerical complexity, namely lid-driven cavity flow, Taylor-Green
vortex and a travelling solitary wave using the level-set method, adopted with
$4^{th}$-order central-differences or a $5^{th}$-order WENO scheme. Our
results show that the EPYC cluster delivers the best code performance for all
the setups under consideration. In the first two benchmarks, the SKL cluster
demonstrates faster computing times than the TX2 system, whilst in the
solitary wave simulations, the TX2 cluster achieves good scalability and
similar performance to the EPYC system, both improving on that obtained with
the SKL cluster. These results suggest that while the Intel SKL cores deliver
the best strong scalability, the associated cluster performance is lower
compared to the EPYC system. The TX2 cluster performance is promising
considering its recent addition to the HPC portfolio.

###### keywords:

HPC , CFD , ARM ThunderX2 , AMD EPYC-Rome , Intel Skylake , Performance
Evaluation , Massively parallel implementation

††journal:

##  1 Introduction

The ability to perform large-scale simulations of turbulent flows has often
been restricted by the available capacity of High-Performance Computing (HPC)
facilities and the inherent limitations in the performance of the associated
processors. HPC-driven turbulence research has enabled new insights into
fundamental theory [1], and enabled engineers to build digital environments as
virtual representations of physical processes in fields such as hydraulics,
environmental turbulent flows, or offshore renewable energy [2, 3, 4].

For almost two decades, processor clock frequencies have been increasing at
very modest rates, favouring the adoption of multi-core CPUs with increasingly
complex memory hierarchies. Even if this approach has delivered double, four,
eight or more times the attainable Flop/s in a chip, rarely is this
accompanied by a similar improvement in memory bandwidth, in practice the main
bottleneck for Computational Fluid Dynamics (CFD) codes given their memory
bound characteristics. This was already foreseen in 2003, when the rate of
increase in CPU performance still followed Moore’s law, by Jiménez [5] who
pointed out that for CFD, memory bandwidth is as limiting a resource as
computing power. Following the widespread adoption of GPUs as compute
platforms over the past decade, an increasing variety of chip architectures
are now available for developers to improve the performance of massively
parallel codes [6]. However, improvements in code performance will of course
depend on the end-user application, e.g. a classical molecular dynamics code
may be compute-bound [7] whilst CFD codes are typically memory-bound [8]. This
calls for heterogeneous HPC systems that comprise a variety of chip
architectures, allowing users to maximise the performance of their specific
computing applications Gray et al. [9].

A new player has recently emerged in the HPC arena: ARM-based CPUs. The latter
have dominated the embedded and mobile market due to its low-power oriented
design, but manufacturers have now presented desktop and server level CPUs
based on this processor architecture [10]. Examples include Marvell with its
Thunder series of processors and Fujitsu with the A64FX.

The Mont-Blanc project, led by the Barcelona Supercomputing Centre (BSC),
pioneered the deployment of ARM chips in HPC environments with a focus on
studying energy efficiency. The first prototype in 2014 [11], Tibidabo,
comprised 128 nodes with dual-core ARM Cortex-A9 CPU, succeeded in 2018 by the
960 node system with ARM Cortex-A15 CPUs [12]. An increasing number of HPC
facilities are now building supercomputers with ARM-based nodes. The most
notable is Fugaku111Fugaku Supercomputer -
https://www.r-ccs.riken.jp/en/fugaku, a supercomputer hosted in the RIKEN
Center for Computational Science, Japan, which hosts 7.3 million ARM-based
Fujitsu A64FX cores. In the most recent Top500 list at the time of writing,
June 2020222Top500 list - https://www.top500.org, Fugaku entered at #1 by
running Linpack with 415.5 PFlop/s of performance and operating at 28.3 MW. In
the HPCG list of the same month333HPCG list -
https://www.top500.org/lists/hpcg, Fugaku also entered at #1 by running the
HPCG benchmark at 13.3 PFlop/s. This outstanding performance evidences the
potential of ARM chips for HPC.

The widespread increase in available computational resources has enabled
major, well documented advancements in the application of CFD to many
engineering fields, notably in hydraulics, atmospheric or environmental
turbulent flows. In 2015, Sotiropoulos [13] reviewed how supercomputers
enabled progress at the frontiers of hydraulics, including problems related to
sediment transport, free-surface flows, or renewable energy, all benefiting
from having high-resolution grids that appropriately resolve the turbulent
scales using Direct-Numerical Simulation (DNS) or Large-Eddy Simulation (LES)
[2, 4].

At present, DNS is restricted to relatively moderate Reynolds numbers (Re) for
classical homogeneous isotropic turbulence or wall-bounded channel flows [1],
as the computing demand depends on the required number of degrees of freedom,
in the order of Re${}^{11/4}_{\lambda}$ [14]. However, its use for engineering
applications remains impractical given the inhomogeneous boundary conditions
and high Reynolds numbers. For instance, Mazzuoli et al. [15] performed a DNS
study with over $2\times 10^{9}$ grid cells and 250,000 spherical particles at
a cost of $10^{7}$ CPU hours, running for approx. 480 days on 64 Intel Ivy
Bridge nodes. This makes it unfeasible to increase the flow regime or at least
it requires excessive computational power [15].

Alternatively, LES enables the resolution of complex multi-phase flows or
simulations in large numerical domains at a moderate-to-high computational
cost. Indeed, one of the most complex challenges in hydraulics and
environmental turbulent flows, the multi-phase flow simulation using
interface-capturing schemes, is amenable to study by LES [16, 17, 18]. In
free-surface open-channel flows, adopting the Level-Set Method (LSM) to
represent the air-water interface can provide more accurate results than using
the classic frictionless rigid-lid boundary condition, but the former comes at
five times higher computational expense (Krosronejad et al. [19]).

This paper studies the performance and scalability of an in-house CFD code on
three clusters featuring distinctive node architectures, namely dual-processor
AMD EPYC-Rome, ARMv8.1 ThunderX2, and Intel Skylake nodes, denoted hereafter
as EPYC, TX2 and SKL respectively. These have distinctive core-per-node
counts, underlying architecture and memory bandwidth, each impacting on the
performance of CFD codes. Initially, two well-known numerical benchmarks are
adopted for the computation of the incompressible Navier-Stokes equation: the
lid-driven cavity flow and Taylor-Green vortex. These tests are designed to
investigate code performance when using 4th-order central differences or 5th-
order WENO schemes, in application to problems of increasing sizes ranging
from four million to a billion grid cells, thereby enabling performance
analysis when increasing the computing workload and communication-to-
computation ratio. A third case involves the propagation of a solitary wave
using the level-set method as representative of a complex multi-phase flow
engineering application.

The structure of the paper is as follows: Section 2 introduces the in-house
code $\mathtt{Hydro3D}$, numerical schemes to compute the fluxes, and the
level-set method. The cluster systems are described in Section 3, while the
test cases used in the evaluation are outlined in Section 4. Performance and
scalability results are presented and discussed in Section 5, with the main
conclusions drawn in Section 6.

##  2 Hydro3D: an open-source CFD code

The software used in this research is $\mathtt{Hydro3D}$ [20, 21], an in-house
open-source code [22] for incompressible, viscous turbulent flows written in
$\mathtt{FORTRAN}$ and fully parallelised with MPI. $\mathtt{Hydro3D}$ allows
the simulation of moving bodies using the immersed boundary method [23, 24,
25, 26, 27], multi-phase flows with the Eulerian-Eulerian approach of the
level-set method [17, 28, 29] and Eulerian-Lagrangian framework of Lagrangian
particle tracking [30], and a local-mesh refinement method [20]. A hybrid
MPI/OpenMP version to compute the Immersed Boundary Method (IBM) module (Ouro
et al. [21]) showed that the MPI/OpenMP scheme can outperform pure MPI
computations when using IBM kernels with large stencils. An efficient multi-
grid solver is used to compute the Poisson pressure equation, whose details
are presented in Cevheri et al. [20].

$\mathtt{Hydro3D}$ is based on the Large-Eddy Simulation (LES) approach in
which the most energetic and largest flow structures are explicitly resolved,
whilst a sub-grid scale method is responsible for modelling the smallest flow
scales considered as those smaller than the filter size, in this case the grid
resolution $\Delta x_{i}$ [2]. The governing equations are the spatially
filtered Navier-Stokes equations for incompressible viscous flows, that read

$\displaystyle\frac{\partial\tilde{u}_{i}}{\partial x_{i}}=0$ (1)
$\displaystyle\frac{\partial\tilde{u}_{i}}{\partial
t}=-\frac{1}{\rho}\frac{\partial p}{\partial
x_{i}}-\tilde{u}_{j}\frac{\partial\tilde{u}_{i}}{\partial
x_{j}}+(\nu+\nu_{t})\frac{\partial^{2}\tilde{u}_{i}}{\partial
x_{j}^{2}}+S_{i}$ (2)

Here $\tilde{u}_{i}$ = $(u,v,w)^{T}$ is the spatially-filtered velocity vector
(for convenience the $(\tilde{\cdot})$ symbol is omitted hereafter), the
coordinates vector is $x_{i}$ = $(x,y,z)^{T}$, $\rho$ denotes the fluid
density, $p$ is the relative pressure, $\nu$ is the kinematic viscosity of the
fluid, and $S_{i}$ is a source term. The eddy-viscosity $\nu_{t}$ is
calculated using the Wall-Adapting Local Eddy-viscosity (WALE) sub-grid scale
model from Nicoud and Ducros [31] as,

$\displaystyle\nu_{t}=(C_{w}\Delta
x_{i})^{2}\frac{\left(S^{d}_{ij}S^{d}_{ij}\right)^{3/2}}{\left(\overline{S}_{ij}\overline{S}_{ij}\right)^{5/2}+\left(S^{d}_{ij}S^{d}_{ij}\right)^{5/4}}$
(3)

with $C_{w}$ being a constant equal to 0.46, $S^{d}_{ij}$ is the traceless
symmetric part of the square of the velocity gradient tensor, and
$\overline{S}_{ij}=0.5(\partial_{j}u_{i}+\partial_{i}u_{j})$ is the strain-
rate tensor of velocities.

The strategy for domain decomposition and MPI communication are described
below in Section 2.1, whilst the computation of fluxes is presented in Section
2.2. Section 2.3 presents the LSM, while the temporal advancement with the
fractional-step method is described in Section 2.4.

###  2.1 Domain partitioning and communication

$\mathtt{Hydro3D}$ adopts Cartesian meshes to discretise the computational
domain. Figure 1 presents the mesh distribution with uniform grid spacing in
every spatial direction with pressure computed on the cell centres and
velocities stored in a staggered fashion, i.e. on the cell faces. The grid
spacing is $\Delta x_{i}={x}_{i+1}-{x}_{i}$, and $\phi_{ij}$ represents the
LSM function that is equal to zero at the two-phase fluid interface, as
described later in Section 2.3.

Figure 1: Representation of the Cartesian grid and staggered velocity
arrangement.

Structured rectangular grids avoid the need for building and communicating
connectivity matrices required to identify the neighbours of every cell in
order to compute the fluxes. For instance, velocity derivative at a cell $i$
using 2nd\- or 4th-order centred finite differences use velocity values at
$\\{i-1,i+1\\}$ or $\\{i-2,i-1,i+1,i+2\\}$, respectively. These schemes allow
the adoption of a balanced domain decomposition of the numerical domain with
an even number of grid cells per sub-domain, in which $n_{i}\times n_{j}\times
n_{k}=n_{t}$ sub-domains are mapped onto $n_{p}$ processing units.

The communication protocol between sub-domains is performed using MPI in which
a number of halo or ghost cells, $n_{g}$ are used to exchange the required
information, such as velocities or pressure values, between adjacent,
overlapping sub-domains. The value of $n_{g}$ is set according to the stencil
of the discretisation scheme to calculate the fluxes, being equal to the
number of cells per direction plus one, i.e. $n_{g}$ = 2 and 3 for the 2nd\-
or 4th-order centred finite differences and 4 when using the 5th-order
Weighted Essentially Non-Oscillatory (WENO) scheme.

The communication overhead when exchanging information by overlapping sub-
domains is proportional to the number of ghost cells and sub-domain divisions
over the considered _ij_ direction, i.e. $n_{GC}$ = $n_{g}\times
n_{x_{i}}\times n_{x_{j}}$, as seen in Figure 2 which depicts the process of
data transfer between two adjacent sub-domains. Hence, in cases with a reduced
sub-domain division, there is a small number of large-size messages to be
communicated. In contrast, when using a large number of sub-domains,
communication is performed by a large number of small-sized messages. In this
version of $\mathtt{Hydro3D}$, point-to-point communications are achieved
using a sequence of non-blocking MPI_IRECV and blocking MPI_SEND and MPI_WAIT
primitives.

Figure 2: Communication between sub-domains using MPI.

###  2.2 Computation of fluxes

The computation of the convective fluxes is the most important component in
the discretisation of the numerical solution (Breuer [32]). The non-linear
convective flux, $\mathcal{C}=u_{j}\partial_{x_{j}}u_{i}$, in the r.h.s of Eq.
2 can be approximated in $\mathtt{Hydro3D}$ using 4th-order Central
Differences (CD) or 5th-order WENO schemes. Whilst CD schemes compute the
fluxes from a single, centred stencil, 5th-order WENO schemes perform a number
of operations to find the optimum stencils which, in turn, increase the
computational expense of the simulations. CD schemes are non-dissipative but
can lead to instabilities in the simulations, as in the case of shock-waves.
WENO schemes introduce numerical dissipation to provide a more stable solution
that needs to be controlled to avoid excessive damping. Due to the inherent
advantages of each scheme, hybrid CD-WENO schemes can be developed for
simulations with shock-turbulence interaction [33, 34], but these are not in
the scope of this paper.

In highly-accurate CD schemes a $p$-order in spatial accuracy is obtained by
approximating velocity derivatives using kernels of, at least, $p+1$ cells in
every spatial direction. Consequently, an increase in scheme accuracy leads to
larger computing times, which can be cumbersome in DNS or LES, and thus a
motivation to study the performance of $\mathtt{Hydro3D}$ when using each of
these.

In a 2nd-order CD scheme, the velocity gradient at a given cell results from
the linear interpolation of velocities from the immediately adjacent cells, as

$\displaystyle f_{i}=\frac{\partial\textit{$u_{i}$}}{\partial
x_{i}}=\frac{\textit{{u}}_{i+1}-\textit{{u}}_{i-1}}{2\Delta x_{i}}$ (4)

The velocity derivative computed with 4th-order CD fits at the $i$-cell a
third-order polynomial through the adjacent four cells, which reads:

$\displaystyle
f_{i}=\frac{-\textit{{u}}_{i+2}+9\textit{{u}}_{i+1}-9\textit{{u}}_{i-1}+\textit{{u}}_{i-2}}{16\Delta
x_{i}}$ (5)

WENO schemes are appropriate to deal with sharp gradients in the velocity
fields by introducing artificial dissipation as an upwind scheme [35]. This
results in a more suitable approach for the computation of free-surface flows
using LSM [28] as small discontinuities can appear at the interface. The basis
of the 5th-order WENO scheme is to build a high-order reconstruction of the
derivative of the flux at a grid cell $x_{i}$, which requires a flux splitting
to compute the positive and negative interface fluxes on the
$x_{i}$-direction, $f_{i\pm 1/2}$, at the cell interfaces, i.e. $\Delta
x_{i\pm 1/2}$, from three candidate stencils. In the calculation of the
velocity derivative $f_{i}$ with WENO, the optimal weights are: $c_{1}$ =
1/10, $c_{2}$ = 3/10 and $c_{3}$ = 6/10. The smoothness indicators,
$\beta_{k}$, serve to scale the optimal weights and are calculated as:

$\displaystyle\beta_{1}$
$\displaystyle=\frac{13}{12}(f_{i-2}-2f_{i-1}+f_{i})^{2}+\frac{1}{4}(f_{i-2}-4f_{i-1}+3f_{i})^{2}$
(6) $\displaystyle\beta_{2}$
$\displaystyle=\frac{13}{12}(f_{i-1}-2f_{i}+f_{i+1})^{2}+\frac{1}{4}(f_{i-1}+f_{i+1})^{2}$
(7) $\displaystyle\beta_{3}$
$\displaystyle=\frac{13}{12}(f_{i}-2f_{i+1}+f_{i+2})^{2}+\frac{1}{4}(3f_{i}-4f_{i+1}+f_{i+2})^{2}$
(8) $\displaystyle\alpha_{k}$
$\displaystyle=\frac{c_{k}}{(\beta_{k}+\varepsilon_{1})^{m}}\hskip
56.9055ptk\text{ = 1, 2, 3}$ (9)

Here $\varepsilon_{1}$ is added to avoid division by zero and set to a value
equal to 10-6, $m$ is set to 2, and weights $\alpha_{k}$ are normalised to
guarantee convexity. Thereafter, the non-linear weights $\omega_{k}$ are
obtained as,

$\displaystyle\omega_{k}=\frac{\alpha_{k}}{\sum_{k}\alpha_{k}}\hskip
56.9055ptk\text{ = 1, 2, 3}$ (10)

Finally, the velocity derivative at the cell $i$ is computed as:

$\displaystyle\begin{split}f_{i}&=\frac{1}{3}\omega_{1}f_{i-2}+\frac{1}{6}(7\omega_{1}+\omega_{2})f_{i-1}+\frac{1}{6}(11\omega_{1}+5\omega_{2}+2\omega_{3})f_{i}\\\
&+\frac{1}{6}(2\omega_{2}+5\omega_{3})f_{i+1}-\frac{1}{6}\omega_{3}f_{i+2}\end{split}$
(11)

Note that the latter flux reconstruction process is performed in the three
spatial directions, i.e. $f_{{i\pm 1/2},j,k}$, $f_{{i},{j\pm 1/2},k}$
$f_{{i},j,{k\pm 1/2}}$, for each of the quantities considered, e.g. $u_{i}$.

In $\mathtt{Hydro3D}$, the linear diffusive or viscous, fluxes,
$\mathcal{D}=\partial_{j}^{2}u_{i}$, are computed with a 4th-order CD scheme
independently to the scheme used for the convective terms as these are of
elliptic nature, which in its one-dimensional version is

$\displaystyle\mathcal{D}_{i}=$ $\displaystyle\nu\Delta
t\left(\frac{-\textit{{u}}_{i+2}+16\textit{{u}}_{i+1}-30\textit{{u}}_{i}+16\textit{{u}}_{i-1}-\textit{{u}}_{i-2}}{12\Delta
x_{i}^{2}}\right)$ (12)

###  2.3 Level-set method

$\mathtt{Hydro3D}$ adopts the Level-Set Method (LSM) to resolve Eulerian-
Eulerian multi-phase flows defining a continuous level-set function, $\phi$,
to determine the fluid density and viscosity fields across the computational
domain, which was validated in Kara et al. [17], Chua et al. [29] and McSherry
et al. [28]. The level-set function assigns positive or negative values
depending on which fluid occupies every grid cell. In the present study,
values of $\phi\geq 0$ denote water whilst $\phi\leq 0$ is air, with the air-
water interface being defined by $\phi=0$ [17] (see Figure 1), and a Heaviside
function $\mathcal{H}(\phi)$ is used to smooth the transition of the density
($\rho(\phi)$) and dynamic viscosity ($\mu(\phi)$) fields across a layer of
thickness $\varepsilon=1.5\Delta x_{i}$, as:

$\displaystyle\rho(\phi)=\rho_{a}+(\rho_{w}-\rho_{a})\mathcal{H}(\phi)$ (13)
$\displaystyle\mu(\phi)=\mu_{a}+(\mu_{w}-\mu_{a})\mathcal{H}(\phi)$ (14)

Here, sub-indices $a$ and $w$ denote values corresponding to the air and
water, respectively. The Heaviside function $\mathcal{H}(\phi)$ defined as:

$\mathcal{H}(\phi)=\begin{cases}0,&\text{\hskip 28.45274pt
$\phi<-\varepsilon$}\\\
\frac{1}{2}\left[1+\frac{\phi}{\varepsilon}+\frac{1}{\pi}\text{sin}\left(\frac{\pi\phi}{\varepsilon}\right)\right],&\text{\hskip
28.45274pt $\mathopen{|}\phi\mathclose{|}<\varepsilon$}\\\ 1,&\text{\hskip
28.45274pt $\phi>\varepsilon$}\end{cases}$

Free-surface flows resolved with LSM pose an additional computational overhead
in comparison to the standard rigid-lid approach in which the air-water
interface is fixed as the vertical upper boundary condition. LSM requires
resolving a non-linear hyperbolic advection equation in addition to the mass
and momentum conservation equations, to account for the transport of $\phi$:

$\displaystyle\frac{\partial\phi}{\partial
t}+\textit{u}_{i}\frac{\partial\phi}{\partial x_{i}}=0$ (15)

This equation is resolved computing the derivatives with the 5th-order WENO
scheme and advanced in time using a third-order Total-Variation Diminishing
Runge-Kutta (TVD-RK3) scheme, which is

$\displaystyle\phi^{k}=\phi(x_{i},t)$ (16)
$\displaystyle\phi^{k+1}=\phi^{k}-\Delta
t\left(u_{i}\frac{\partial\phi^{k}}{\partial x_{i}}\right)$ (17)
$\displaystyle\phi^{k+2}=\frac{3}{4}\phi^{k}+\frac{1}{2}\phi^{k+1}-\frac{1}{4}\Delta
t\left(u_{i}\frac{\partial\phi^{k+1}}{\partial x_{i}}\right)$ (18)
$\displaystyle\phi^{k+3}=\frac{1}{3}\phi^{k}+\frac{2}{3}\phi^{k+1}-\frac{2}{3}\Delta
t\left(u_{i}\frac{\partial\phi^{k+2}}{\partial x_{i}}\right)$ (19)

Note that flux splitting is applied to compute each gradient of $\phi^{k}$,
$\phi^{k+1}$ and $\phi^{k+2}$ at $\Delta x_{i\pm 1/2}$, which implies that for
Eq. 17 six WENO fluxes for $\phi^{k}$ are calculated, i.e. $f_{{i\pm
1/2},j,k}$, $f_{{i},{j\pm 1/2},k}$ $f_{{i},j,{k\pm 1/2}}$, which is then
repeated for $\phi^{k+1}$ and $\phi^{k+2}$.

Due to the intrinsic nature of the $\phi$-advection equation, mass
conservation is not guaranteed [36] as the required numerical stability
requirement of $|\nabla\phi|=1$ is not automatically fulfilled. This implies
that a re-initialisation technique is computed after Eq. 15 [36], which
requires resolving another advection equation, Eq. 20, for the signed-distance
function, $d_{0}$, solved with the TVD-RK3 scheme (Eq. 17 – 19) according to
an ”artificial time step” $\tau$ = CFLLSM$\cdot\texttt{MAX}(\Delta x_{i})$.
The smoothed signed function, $s(d_{0})$, in Eq. 20 is the calculated
according to Eq. 21 adopting the initial condition $d_{0}(x_{i},0)$ =
$\phi(x_{i},t)$.

$\displaystyle\frac{\partial d_{0}}{\partial\tau}=s(d_{0})(1-|\nabla d_{0}|)$
(20) $\displaystyle s(d_{0})=\frac{d_{0}}{\left(d_{0}^{2}+\left(|\nabla
d_{0}|\varepsilon_{r}\right)^{2}\right)^{0.5}}$ (21)

This re-initialisation process is iteratively solved until the condition of
$|\nabla\phi|=1$ is partially satisfied up to a residual $\varepsilon_{LS}$
set to 5 $\cdot$ 10-3. Additionally, a maximum of 15 iterations is set, since
the free-surface layer in turbulent flows will experience small high-frequency
oscillations due to the turbulent structures [28, 16]. These challenge fast
numerical convergence. Identically to the $\phi$-equation resolution, 18
stencil reconstructions are required at every re-initialisation step in
addition to the calculation of Eq. 20. This can notably increase the
computational cost if the number of iterations required during the re-
initialisation is large.

Alg. 1 presents the procedure to resolve the LSM, with TLS and
TLS_$\omega_{k}$ labelling the time spent on solving the whole LSM and that
only for the WENO weights in LSM, respectively.

Algorithm 1 Pseudocode for advection of $\phi$ function in LSM

1:$\phi^{t}=\phi(x_{i},t)$ $\triangleright$ TLS=MPI_WTIME

2:for k = 1, 3 do

3: WENO stencil for $\phi^{t}$ $\triangleright$ TLS $\omega_{k}$

4: Compute $\phi^{t+k}$ with Eqs. 17 \- 19

5:end for

6:Re-initialisation process

7:for m = 1, 15 do

8: $d_{0}=\phi^{t+1}$

9: while $|\nabla d_{0}|-1{>}\hskip 1.0pt\varepsilon_{LS}$ and m $\geq$ 2 do

10: for k = 1, 3 do

11: WENO stencil for $d_{0}^{k}$ $\triangleright$ TLS $\omega_{k}$

12: Compute $|\nabla d_{0}^{k}|$

13: Compute $s(d_{0}^{k})$

14: Compute $d^{k}$ with Eqs. 17 \- 18

15: end for

16: call MPI_ALLREDUCE($|\nabla d_{0}|$)

17: end while

18:end for

19:$\phi^{t+1}=d_{0}^{k}$

20:Heaviside function: $\mathcal{H}(\phi^{t+1})$

21:Update $\rho(\phi^{t+1})$ and $\mu(\phi^{t+1})$ fields $\triangleright$
TLS=TLS-MPI_WTIME

Note that the computational expense of using LSM in turbulent flows is
intrinsically high due to the fine grid resolution needed to capture the
instabilities at the air-water interface. Too coarse a grid resolution will
challenge numerical stability. For this reason, WENO is used to compute both
the LSM advection equations and convective fluxes in the momentum equation.
New insights into the computational cost of using LSM for free-surface flows
is of special relevance for the engineering research community, defining how
different HPC chip architectures impact on the performance of computing LSM.

###  2.4 Temporal advancement

The standard fractional-step method [37] is used to advance the equations in
time, as its step-wise nature allows the development of highly-scalable
parallelised codes. The adopted fractional-step method is fully explicit and
divides the computing sequence at each time step into three main stages as
depicted from Alg. 2. First, the predicted non-divergence-free velocity field,
$u_{i}^{*}$, is calculated from the pressure and velocities obtained at the
previous time step, $p^{t}$ and $u_{i}^{t}$, accomplished using an explicit
time advancement with a three-stage Runge-Kutta method. The second step
comprises the resolution of the Poisson pressure equation. This is iteratively
solved using a multi-grid method in order to correct the predicted velocity
field until the divergence-free condition, i.e. $\nabla\cdot u^{*}_{i}$, is
below a given residual $\varepsilon$ equal to $10^{-6}$. In the final step the
velocity and pressure are updated. If the LSM method is used, this is computed
before the velocity fluxes.

$\mathtt{Hydro3D}$ has an internal profiler to track the time spent in the
main subroutines at every time step, obtained with the master processor and
using the MPI_WTIME command. These are highlighted in Alg. 2 and account for:

  * 1.

TTT: Total time per time step or mean runtime.

  * 2.

TLS: Level-Set Method.

  * 3.

TCD: Convective and viscous fluxes computed over the three-step Runge-Kutta
(RK) to obtain the non-divergence free velocity.

  * 4.

TP: Resolution of the Poisson pressure equation.

  * 5.

Tup: Update of the velocity and pressure fields.

Note that subroutines consuming a reduced part of the overall computing time
per time step, e.g. less than 1.0%, are normally excluded from the analysis,
e.g. sub-grid scale model or time step calculation based on the CFL condition,
and thus not presented in the results.

In terms of collective communications, two main MPI_ALLREDUCE calls are
performed to determine velocity maxima that condition the time advancement of
the simulation. First, the maximum velocity from the previous time step,
$u_{i}^{t}$, is fetched to determine the time step $\Delta t$ according to the
adopted CFL value. Then the divergence free condition in the Poisson pressure
solver is checked with the non-solenoidal velocity $u_{i}^{*}$ after every
iteration of the solver (lines 17–19 in Alg. 2), with the code needing to
retrieve MAX$(\nabla\cdot u^{*}_{i})$ from every sub-domain, at every pressure
solver iteration.

Algorithm 2 Fractional-step method adopted in $\mathtt{Hydro3D}$

1:Variable allocation and initialisation

2:for t = $t_{0},t_{end}$ do $\triangleright$ TTT=MPI_WTIME

3: Store $u^{t}_{i}=u_{i}^{t-1}$

4: if Level-Set Method then $\triangleright$ TLS=MPI_WTIME

5: Alg. 1

6: end if$\triangleright$ TLS=TLS-MPI_WTIME

7: for $t_{RK}$ = 1, 3 do $\triangleright$ TCD=MPI_WTIME Convective fluxes,
$\mathcal{C}$ :

8: if $2^{nd}$-order CD then

9: Eq. 4

10: else if $4^{th}$-order CD then

11: Eq. 5

12: else WENO

13: Eqs. 6 \- 11

14: end if Viscous fluxes, $\mathcal{D}$: Eq. 12

15: end for$\triangleright$ TCD=TD-MPI_WTIME Non-divergence free velocity:

16: ${u_{i}}^{*}$ = ${u_{i}}^{t}$ \+ $\Delta t$ ($\mathcal{D}$ $({u_{i}}^{t}$)
+ $\mathcal{C}$ (${u_{i}}^{t}$) + $\nabla p^{t}$ \+ $S_{i}$) Solve Poisson
pressure equation: $\triangleright$ TP=MPI_WTIME

17: while $\nabla\cdot u_{i}^{*}{>}\hskip 1.0pt\varepsilon$ do

18: $\nabla^{2}\hat{p}=({\nabla\cdot u_{i}^{*}})/{\Delta t}$

19: end while$\triangleright$ TP=TP-MPI_WTIME Update velocity and pressure
fields:  $\triangleright$ Tup=MPI_WTIME

20: $u^{t+1}=u_{i}^{*}+\Delta t\nabla\hat{p}$

21: $p^{t+1}=p^{t}+\hat{p}-{\nu\Delta t\Delta\hat{p}}/{2}$ $\triangleright$
Tup=Tup-MPI_WTIME $\triangleright$ TTT=TTT-MPI_WTIME Time averaging and write
output files

22:end for

##  3 Cluster architectures

In this work, three HPC systems are used, each equipped with distinctive dual-
processor computing nodes. Again, for convenience, these systems are codenamed
as EPYC, SKL and TX2, which use AMD EPYC Rome, Intel Skylake and ARM Marvell
ThunderX2 processors respectively. Full system descriptions are provided in
Table 1.

Table 1: Description of the three HPC systems. System name | ThunderX2 | Skylake | EPYC-Rome
---|---|---|---
Computing Service | GW4 Isambard | Supercomputing | Supercomputing
|  | Wales (Hawk) | Wales (Hawk)
Number of nodes | 128 | 201 | 64
Core count/node |  2 CPUs$\times$32 = 64 |  2 CPUs$\times$20 = 40 |  2 CPUs$\times$32 = 64
Memory | 256 GB DDR4 | 192 GB DDR4 | 256 GB DDR4
per node | @2400MHz | @2666MHz | @2500MHz
Interconnect | Cray Aries | Infiniband EDR | Infiniband EDR
CPU spec: |  |  |
Architecture | ARM v8.1 | Intel Skylake | Rome (Zen 2)
Model | Marvell ThunderX2 | Intel Xeon 6148 Gold | EPYC-Rome 7502
Clock | 2.10 GHz | 2.40 GHz | 2.50 GHz
Memory channels | 8 | 6 | 8
TDP | 175 W | 150 W | 180 W

The TX2 system is the full ARM partition of Isambard, the first large-scale
ARM-based production level HPC system in the UK managed by the GW4 Alliance,
comprising the Universities of Bristol, Bath, Cardiff and Exeter, together
with Cray and hosted at the UK Met Office. Isambard is a Cray XC50 system
containing over 10,000 ARM cores and a full Cray software stack. More
information about the system, including the full specification, is available
at 444GW4 Isambard - https://gw4.ac.uk/isambard.

The SKL and EPYC systems are partitions of Hawk, an HPC system hosted by
Supercomputing Wales (SCW) at Cardiff University. SCW is a pan-Wales UK
project featuring the Universities of Cardiff, Swansea, Bangor and
Aberystwyth. At the time of writing, the SKL partition comprises 201 dual
processor nodes, totalling 8,040 cores and 46 TByte of memory, whilst the 64
node EPYC partition comprises 4096 EPYC-Rome 7502 cores with 50 TByte of
memory. More information about the system and the project is available at
555Supercomputing Wales - https://www.supercomputing.wales.

##  4 Benchmarks description

This section describes the three benchmarks selected to analyse the
performance of $\mathtt{Hydro3D}$, designed to deliver a range of
configurations that provide a varying balance of compute- and memory-bound
simulations.

###  4.1 3D lid-driven cavity flow

The first benchmark comprises a three-dimensional (3D) lid-driven cavity flow
simulated using 4th-order CD to compute the convective fluxes, which is the
fastest compute discretisation scheme using Cartesian grids. This minimises
the overhead in computing convective and diffusive fluxes. In order to analyse
the code performance, six problem sizes with increasing number of grid cells
are adopted.

The computational domain is a cube with dimensions
[0,1]$\times$[0,1]$\times$[0,1] and an inlet velocity of $U_{0}(x_{i})$ =
(0.0,0.0,1.0) set at the top boundary, i.e. at $z=1$, as a slip condition. No-
slip conditions are used at $x=0$ and 1 and at the cube bottom, $z=0$, whilst
periodic conditions are set at the transverse faces at $y=0$ and 1. The mesh
resolution is uniform in the three spatial directions and variable time step
is set with a CFL condition equal to 0.8. The Reynolds number is set to 400 as
adopted in previous scalability studies performed by Ouro et al. [21], who
previously validated the velocity field with $\mathtt{Hydro3D}$.

###  4.2 Taylor-Green vortex

The second benchmark comprises the simulation of the incompressible Taylor-
Green Vortex (TGV) adopting the 5th-order WENO scheme for the calculation of
the convective terms. This case is often selected to evaluate the numerical
properties of high-order discretisation schemes, e.g. scheme dissipation rate
[38, 39], such as the 5th-order WENO. The case simulated has a $Re$ =
$LU_{0}/\nu$ = 1,600 with $U_{0}$ equal to 1.0 and characteristic spatial
scale $L$ = 1 as the cube’s edge length. The spatial domain extends
$0<x,y,z<L$ with periodic boundary conditions on the three spatial directions.
A variable time step was set with a CFL condition of 0.3. The initial flow
field is defined according to the following velocity distribution:

$\displaystyle u(x,y,z)=U_{0}\text{ sin }\left(\frac{x}{L}\right)\text{ cos
}\left(\frac{y}{L}\right)\text{ cos }\left(\frac{z}{L}\right)$ (22)
$\displaystyle v(x,y,z)=-U_{0}\text{ cos }\left(\frac{x}{L}\right)\text{ sin
}\left(\frac{y}{L}\right)\text{ cos }\left(\frac{z}{L}\right)$ (23)
$\displaystyle w(x,y,z)=0$ (24)

Figure 3 presents iso-surfaces of $z$-vorticity
($\omega_{z}=\partial_{x}v-\partial_{y}u$) obtained at non-dimensional times
$t$ = 0 and 12 using a computational mesh with 32.8 million elements. These
depict how the initial large-scale vortices transition to smaller scale
structures as the simulation advances in time.

Figure 3: Iso-surfaces of $z$-vorticity at non-dimensional times $t$ = 0 and
12.

###  4.3 Solitary wave

In the third benchmark, a solitary wave propagating across a numerical wave
tank is simulated, which represents one of the current modelling challenges in
computational hydraulics [3]. The air-water interface is resolved using the
LSM which implies a higher computational load than the previous benchmarks, as
indicated in Section 2.3. The numerical wave tank measures 12.8 m $\times$
0.4m $\times$ 0.4m and is presented in Figure 4, with a water depth $d$ set to
0.2 m and solitary wave amplitude $H$ of 0.02 m. The density of the water is
1,000 kg m-3 and that of air 1.25 kg m-3, with viscosity values of
1$\cdot$10-6 m2 s-1 and 1.8$\cdot$10-5 m2 s-1, respectively. A constant time
step is set to $\Delta{t}$ = 0.001 s, whilst for the advection equation the
pseudo-time step ($\Delta\tau$) is variable with a CFL = 0.10 to guarantee
numerical stability [16, 40].

Figure 4: Schematic of the solitary wave benchmark case.

The solitary wave is generated using Boussinesq theory [41] with Eq. 25
defining the wave elevation, $\eta$, and Eqs. 26 and 27 denoting wavenumber,
$k_{w}$, and wave celerity, $c$, respectively.

$\displaystyle\eta(t)=H\text{cosh}^{-2}(k_{w}(x-ct))$ (25) $\displaystyle
k_{w}=\sqrt{{3h}/{4d^{3}}}$ (26) $\displaystyle c=\sqrt{g(H+d)}$ (27)

Here $g$ represents gravity acceleration equal to 9.81 m s-2 and $t$ is time.
At the start of the simulations, velocities are set to zero across the
computational domain. Thereafter, considering the wave elevation is normalised
as $\eta_{h}=\eta/H$ and $\varepsilon=H/d$, the velocities at the inflow
condition ($x_{in}=(0,y,z)$) are prescribed according to:

$\displaystyle\begin{split}&u(x_{in},t)=\varepsilon\sqrt{gd}\left[\eta_{h}-\frac{\varepsilon\eta_{h}^{2}}{4}+\frac{d^{2}}{3c^{2}}\left(1-\frac{3z^{2}}{2d^{2}}\right)\frac{\partial^{2}\eta_{h}}{\partial
t^{2}}\bigg{\rvert}_{t=0}\right]\end{split}$ (28) $\displaystyle
v(x_{in},t)=0$ (29)
$\displaystyle\begin{split}&w(x_{in},t)=z\frac{\varepsilon}{c}\sqrt{gd}\left[\left(1-\frac{\varepsilon\eta_{h}}{2}\right)\frac{\partial\eta_{h}}{\partial
t}+\frac{d^{2}}{3c^{2}}\left(1-\frac{z^{2}}{2d^{2}}\right)\frac{\partial^{3}\eta_{h}}{\partial
t^{3}}\bigg{\rvert}_{t=0}\right]\end{split}$ (30)

##  5 Results

This section presents the performance results of $\mathtt{Hydro3D}$ based on
an analysis of strong scaling behaviour with increasing problem size. The
comparison between the three HPC systems is reported on a node-to-node basis,
with each configuration designed to populate each node with a single MPI
process per core. Table 2 lists the configurations adopted in terms of node
count and sub-domain divisions which are equal to the number of cores used.
Note that the TX2 and EPYC nodes have 64 cores per node, whilst the SKL nodes
have 40. Given resource availability, the simulations with SKL and EPYC nodes
on the Hawk cluster are run on a maximum of 64 nodes, i.e. 2,560 and 4,096
cores respectively, whilst those with TX2 nodes on the Isambard cluster used
up to 125 nodes, i.e., 8,000 cores.

Table 2: Configurations of the simulations carried out on the ThunderX2 (TX2), Skylake (SKL) and EPYC-Rome (EPYC) systems for the three benchmark cases, including the number of nodes used (N) and sub-domains ($n_{t}=n_{x}\times n_{y}\times n_{z}$) that is equal to the number of cores. Benchmark 1: Cavity flow | Benchmark 2: TGV | Benchmark 3: solitary wave
---|---|---
TX2 | SKL | EPYC | TX2 | SKL | EPYC | TX2 | SKL | EPYC
N |  $n_{t}$ |  N |  $n_{t}$ |  N |  $n_{t}$ |  N |  $n_{t}$ |  N |  $n_{t}$ |  N |  $n_{t}$ |  N |  $n_{t}$ |  N |  $n_{t}$ |  N |  $n_{t}$
1 |  64 |  1 |  40 |  1 |  64 |  1 |  64 |  1 |  40 |  1 |  64 |  1 |  64 |  1 |  40 |  1 |  64
2 |  128 |  2 |  80 |  2 |  128 |  2 |  128 |  2 |  80 |  2 |  128 |  2 |  128 |  2 |  80 |  2 |  128
4 |  256 |  4 |  160 |  4 |  256 |  4 |  256 |  4 |  160 |  4 |  256 |  4 |  256 |  4 |  160 |  4 |  256
8 |  512 |  8 |  320 |  8 |  512 |  8 |  512 |  8 |  320 |  8 |  512 |  8 |  512 |  8 |  320 |  8 |  512
10 |  640 |  10 |  400 |  10 |  640 |  16 |  1000 |  16 |  640 |  16 |  1000 |  16 |  1024 |  16 |  640 |  16 |  1024
20 |  1280 |  16 |  640 |  16 |  1024 |  25 |  1600 |  25 |  1000 |  25 |  1600 |  32 |  2048 |  32 |  1280 |  32 |  2048
25 |  1600 |  25 |  1000 |  20 |  1280 |  32 |  2000 |  32 |  1280 |  32 |  2000 |  40 |  2560 |  64 |  2560 |  40 |  2560
32 |  2000 |  32 |  1280 |  25 |  1600 |  63 |  4000 |  63 |  2500 |  63 |  4000 |  64 |  4096 |  |  |  64 |  4096
63 |  4000 |  52 |  2048 |  32 |  2000 |  |  |  |  |  |  |  |  |  |  |  |
125 |  8000 |  64 |  2560 |  63 |  4000 |  |  |  |  |  |  |  |  |  |  |  |

The code was compiled with Intel Compiler v2017.4 using the -O2 flag for the
EPYC and SKL platforms, and executed using the Intel MPI Library v2017.4. On
the Isambard platform, the code was compiled using the Cray Compiling
Environment (CCE) v9.1.3 with the -O2 flag and executed with the Cray MPI
library using the ALPS system v6.6.

All benchmark simulations are run for 50 time steps, with mean computing times
calculated based on the MPI_WTIME command averaged over the last 40 time
steps, over which the standard deviation of the computing times is lower than
$0.1\%$. Note that during the domain partitioning, the sub-domains size is
kept uniform for all cases so as to avoid load imbalance between parallel
processes.

The runtime values presented from the internal code profiler includes both
computing and MPI operations. The breakdown of the main computing subroutines
are calculated using the profilers CrayPat on Isambard (TX2) and TAU v2.27
[42] on Supercomputing Wales (SKL and EPYC), of value in characterising how
systems distribute the computing workload for the same setup. The breakdown of
MPI directives is obtained on Hawk using the IPM profiler.

###  5.1 Node characterisation

In order to gain a sense of the raw computing power of the HPC clusters, a set
of stress benchmarks on individual nodes of each system were undertaken,
including Firestarter v1.7.4 to stress CPU cores and gather their peak
performance, and STREAM to measure memory bandwidth. These two benchmarks are
described in [43] and [44], respectively.

The Firestarter binary available in its official repository was used, while
STREAM was compiled with the Intel Compiler v2018.4 for the EPYC and SKL
systems. The STREAM array size was $2^{25}$, with the reported values
corresponding to the TRIAD test. The performance figures for TX2 are those
reported in Deakin et al. [45], as Firestarter did not run on the system due
to an incompatibility with ARMv8.

Table 3: Measured peak performance and memory bandwidth for individual nodes. Results from Firestarter and STREAM benchmarks. | ThunderX2 | Skylake | EPYC-Rome
---|---|---|---
TFLOP/s | 1.28 | 2.06 | 1.99
GB/s | 288.00 | 180.36 | 264.90

Results are shown in Table 3 for the three systems. SKL and EPYC clusters
deliver a similar level of TFLOP/s performance($\sim$3% difference between
them) whilst the EPYC system achieves 46.8% greater memory bandwidth. TX2
shows the lowest TFLOP/s performance figure but delivers the greatest memory
bandwidth amongst the three chips, i.e. 8.7% and 59.7% more than EPYC and SKL
processors, respectively.

###  5.2 Lid-driven cavity flow

Six grid resolutions are adopted for the first benchmark in order to perform a
set of tests that gradually evolve from a computing-bound scenario at low node
count to a communication-bound scenario at higher node count. Details of the
computational grids are provided in Table 4 with spatial resolution
($\Delta{x_{i}}$), number of elements per spatial direction ($n_{x_{i}}$) and
total number of elements comprising the entire problem ($N_{e}=n_{x_{i}}^{3}$)
in millions.

Table 4: Numerical details of the lid-driven cavity flow cases simulated. $\Delta{x_{i}}$ | $n_{x_{i}}$ | Ne ($\cdot 10^{6}$)
---|---|---
0.00625 | 160 | 4.1
0.00500 | 200 | 8.0
0.003125 | 320 | 32.8
0.00250 | 400 | 64.0
0.00125 | 800 | 512.0
0.00100 | 1,000 | 1,000.0

Figure 5 depicts the mean runtime, i.e. averaged total time per time step,
(TTT from Alg. 2) in log-log scale obtained with the three HPC systems for the
six problem sizes. The runtime on the SKL and EPYC systems follows a linear
decay for problem sizes Ne $\geq$ 64 million cells suggesting
$\mathtt{Hydro3D}$ scales well on these architectures. Conversely, when using
a large number of nodes for the three smallest problem sizes, i.e. Ne $\leq
32.8$, the average runtime plateaus as communication overheads start to
increase compared to a relatively lower computing demand given the small
number of grid cells per core. This is more noticeable in the TX2 system as,
for instance, with Ne = 8.0 the percentage of communications using two nodes
is approx. 11%, increasing to almost 90% for the highest node count.

Figure 5: Lid-driven Cavity flow: Results of averaged time per time step
values in seconds for the three HPC architectures and six grid sizes analysed.

For most problem sizes on the TX2 cluster, there is a steady decrease in
runtime with increasing number of processors although this follows a smaller
reduction rate than on the other systems, indicating that the TX2 processors
demonstrate less effective strong scaling in this case compared to the x86
processors. It is worth noting that for Ne $\leq 8$, TX2 fails to achieve any
reduction in the computation time when increasing the number of nodes beyond 8
due to communication overheads. For instance, for Ne = 4.1 using 2, 8 or 25
nodes, the proportion of runtime spent on MPI communications grows from 10%,
to 48%, and 83% respectively, whilst for Ne = 64 these figures are 16%, 16%
and 45%. Nevertheless, for relatively large problem sizes, i.e. Ne = 512 and
1,000, the performance and scalability of the TX2 cluster improves featuring a
linear decay in runtime similar to that found on the SKL and EPYC systems.
Overall, the results of Figure 5 show that the EPYC cluster provides the
fastest computations for all the lid-driven cavity flow simulations.

Figure 6: Lid-driven Cavity flow: Results of speed up for the three HPC
architectures and six grid sizes analysed. Ideal linear speedup is represented
with a straight black line.

Strong scaling results are presented in Figure 6 in terms of speed up
$S_{n}=T_{1}/T_{n}$, with $T_{n}$ denoting the averaged runtime per time step
when using $n$ processors, in log-log scale. It is observed that the speed up
on TX2 for small problem sizes does not increase linearly when the number of
nodes is over 8, again a consequence of the large communication overhead.
Conversely, increasing Ne improves the speed up on TX2 up to 125 nodes,
achieving comparable strong scalability to that found on the SKL and EPYC
systems. For these small size problems, scalability on EPYC is inferior to
that found on SKL as a result of the latter spending a lower percentage of
time in communications, as shown later in Figure 9.

Further insight into the code’s strong scalability is provided in terms of
parallel efficiency $E=\frac{PT_{p}}{QT_{q}}$, with $T_{p}$ and $T_{q}$ the
time per time step when using $P$ and $Q$ nodes, respectively. Ideal
scalability is achieved when the parallel efficiency is 1.0, i.e. the decrease
in runtime should be directly proportional to the increase in number of
processors. However, the time spent on communication and computational
processes, e.g. solving the Poisson equation with multi-grid [21], or even
caching effects, could lead to non-linear $E$ growth when increasing the
number of nodes.

Figure 7: Lid-driven Cavity flow: Results of parallel efficiency for the three
HPC architectures and six grid sizes analysed.

Figure 7 presents on a semi-log scale the values of $E$ for the lid-driven
cavity flow simulations. There is a good parallel efficiency on SKL with
increasing number of nodes for those problem sizes with Ne $\geq$ 32.8, whilst
for the two smallest problem sizes using a larger number of nodes impacts
negatively on the code’s scalability as $E\leq 1$. TX2 and EPYC results show
inferior parallel efficiency in problem sizes with Ne $\geq$ 64, but exhibit
improved performance when Ne $\leq$ 64 with values of $E$ near unity.

In order to better understand the performance results of Hydro3D, a breakdown
of the compute timestep is presented in Figure 8 when using two nodes and
problem sizes Ne = 8 and 64. In all cases, solving the iterative Poisson
pressure equation (TP) takes most of the computing time, even outweighing all
other subroutines combined for both SKL and EPYC systems. On TX2, the time
spent on convection (TC) is similar to TP and larger than that obtained on the
other two systems. Note that increasing the problem size has little effect on
the distribution of computing times.

Figure 8: Distribution of the main computing functions for the lid-driven
cavity flow benchmark when two nodes are used. Comparison between TX2 (left),
SKL (centre) and EPYC (right) for meshes Ne = 8 (a) and 64 (b). TP = Poisson
pressure solver, TD = diffusion, TC = convection, Tup = velocity update.

To complement the performance breakdown, Figure 9 (top) presents the
percentage of the total runtime spent in communication for the three cluster
systems, for all six problem sizes using 2, 8 and 25 nodes. For the smaller
node counts, communication accounts for less than 20% of the simulation time
in all cases, this being slightly higher on TX2. With 8 nodes, the Ne = 4.1
problem size shows an increased time in communication especially for TX2,
which decreases when increasing Ne as the computing workload increases.
However, with 25 nodes the relative time spent in communication with the x86
systems stays in the range of 20–25% while it reaches up to 70% for
$N_{e}\leq$ 32.8 on TX2.

Figure 9: Comparison for the lid-driven cavity flow benchmark of the
percentage of runtime spent on communication (top) when running on 2, 8 and 25
nodes with the three core architectures; and distribution of MPI directives
for the cases with SKL (red line) and EPYC (green line).

Figure 9 (bottom) provides a more detailed communication profile from
simulations on SKL and EPYC with the breakdown of the MPI time spent into
synchronisation with MPI_WAIT and MPI_BARRIER, collective communication
MPI_ALLREDUCE, and point-to-point communication MPI_SEND. These data were not
collected for TX2 as the IPM compiler used for this analysis was not
available. Synchronisation with MPI_WAIT arises as the major communication
overhead except for $N_{e}\geq$ 512 in which MPI_SEND increases considerably
as a result of communications involving many small messages distributed across
a large number of cores. The switch from MPI_WAIT to MPI_SEND consuming most
of the MPI communication overhead depends on the problem size and number of
nodes, and also varies slightly between EPYC and SKL systems. For this
benchmark case, collective communications have a reduced impact in the overall
communication time.

The results for this lid-driven cavity flow benchmark, in which 4th-order CD
is used to compute the convective fluxes, indicate that the EPYC system
delivers the best computational performance whilst the SKL cluster exhibits
the best strong scalability.

###  5.3 Taylor-Green vortex

For the Taylor-Green Vortex (TGV) benchmark, the computing overhead from the
convective fluxes increases from using a 5th-order WENO scheme compared to CD
schemes, as the former requires the computation of optimal weights [46].
Analogous to the previous case, four mesh sizes are adopted with details
provided in Table 5 regarding mesh size ($\Delta{x_{i}}$), number of elements
per spatial direction ($n_{x_{i}}$) and total number of grid cells in millions
(Ne). Note that this range of grid resolutions is similar to those available
in the literature [38].

Table 5: Numerical details of the meshes cases used for the TGV simulations. $\Delta{x_{i}}$ | $n_{x_{i}}$ | Ne ($\cdot 10^{6}$)
---|---|---
0.005 | 200 | 8.0
0.003125 | 320 | 32.8
0.0015625 | 640 | 262.1
0.00100 | 1,000 | 1,000

Results of the mean runtime per time step in seconds, speed up ($S_{n}$) and
parallel efficiency ($E$) are presented in Figure 10. For the smallest problem
size, EPYC and SKL systems achieve lower runtimes than TX2 for all core
counts, being up to 10$\times$ faster when 63 nodes are used. As with the lid-
driven cavity flow benchmark, the runtime achieved with TX2 for Ne = 8
plateaus when running on 16 or more nodes due to communication overheads. This
is also observed in the $S_{n}$ and $E$ plots, pointing to the need to improve
the strong scalability of $\mathtt{Hydro3D}$ on TX2 for relatively small
problem sizes, even though this is superior to that obtained for the lid-
driven cavity flow.

Figure 10: Results of runtime (top), speed up ($S_{n}$) (middle), and parallel
efficiency ($E$) (bottom) for the Taylor-Green vortex case.

With Ne = 32.8, $\mathtt{Hydro3D}$ achieves an almost linear $S_{n}$ increase
for each system. Taking a closer look at the relative performance between
processors for large problem sizes, for Ne = 262.1 SKL is on average approx.
1.4$\times$ faster than TX2 whilst EPYC is 2$\times$ and 2.6$\times$ faster
than SKL and TX2, respectively. Based on node-to-node comparison, with Ne =
1,000, SKL provides a 1.1$\times$ performance enhancement compared to TX2,
whose improved performance in large problem sizes is due to increased memory
bandwidth. On the other hand, EPYC is again the fastest system being
2.0$\times$ faster than SKL and 2.3$\times$ faster than TX2. Compared to the
lid-driven cavity flow case, the linearity observed in the SKL $S_{n}$ and $E$
curves is less apparent, resulting perhaps from the WENO weights becoming the
most expensive subroutine over the Poisson pressure solver.

Figure 11 shows the breakdown of the runtime spent on the main computing
subroutines, namely Poisson pressure solver (Tp), diffusion (TD), convection
(Tc), velocity update (Tup) and WENO weights (T$\omega_{k}$), for Ne = 8 and
262.1 and using two nodes. The computing workload is similarly distributed
across systems with the multi-grid pressure solver and WENO weights taking
approx. 18–20% and 41–47% of the computing time, except on TX2 with Ne = 8.
This is attributed to the lower computing capability of TX2, as seen in
Section 5.1, increasing the time required for iteratively solving the pressure
equation. Although not included here for sake of brevity, results for the
other two meshes and larger node count present a similar distribution,
pointing to the codes’ subroutines scaling well.

Figure 11: Breakdown of the main computing functions for the TGV benchmark
when two nodes are used. Comparison between TX2 (left), SKL (centre) and EPYC
(right) for meshes Ne = 8.0 (a) and 262.1 (b).

The communication profiles are presented in Figure 12, showing that for 2
nodes the simulations spend little time in communications. For 8 and 25 nodes,
the relative weight of data exchange increases, especially for small problem
sizes on TX2, with the percentage of communication notably increasing. In
relation to the MPI directives, MPI_WAIT and MPI_WAIT consume most of the
communication time with a balance similar to that observed in the lid-driven
cavity flow results when increasing the number of nodes.

Figure 12: Comparison of the percentage of runtime spent on communication
(top) when running on 2, 8 and 25 nodes for the TGV benchmark; and
distribution of MPI directives for the cases with Skylake (red line) and EPYC-
Rome (green line).

###  5.4 Solitary wave simulation with level-set method

The third benchmark presents the solitary wave simulation using LSM. The
previous TGV benchmark showcased the computing overhead from the 5th-order
WENO scheme, which is now expected to increase given the computation of both
convective fluxes and LSM advection and reinitialisation equations. In this
benchmark, the travelling wave is simulated using three grid resolutions with
the details of mesh resolution ($\Delta{x_{i}}$) that is uniform across the
domain, the number of divisions per spatial directions ($n_{x_{i}}$) and the
total number of grid cells (Ne, in millions) presented in Table 6.

Table 6: Numerical details of the meshes used for the solitary wave benchmark. $\Delta{x_{i}}$ | $n_{x}$ | $n_{y}$ | $n_{z}$ | Ne ($\cdot 10^{6}$)
---|---|---|---|---
0.0100 | 1,280 | 40 | 30 |  1.5
0.0050 | 2,560 | 80 | 60 | 12.3
0.0025 | 5,120 | 160 | 120 | 98.3

The mean time per time step obtained for the three problem sizes on each HPC
system is shown in Figure 13 (top) using a log-log scale. All systems show an
almost linear decrease in runtime values with increasing node count. Unlike
the former benchmarks, the percentage of communication in the solitary wave
simulation is less than 20% of the total runtime for configurations running on
the maximum number of nodes. $\mathtt{Hydro3D}$ achieves its best performance
on the EPYC cluster for all cases, although it is noticeable that simulation
times obtained with TX2 outperform those from SKL by approx. 1.4–1.6$\times$,
especially in the highest resolution case, in which TX2 and EPYC systems
achieve similar runtimes. In this benchmark, both systems spend a relatively
lower percentage of time in communications than does SKL, thereby enhancing
inter-node communications and proving advantageous when using LSM as
MPI_ALLREDUCE is called several times triggering the time in data exchange.

Figure 13: Results of averaged runtime per time step values in seconds (top),
speed up (mid) and parallel efficiency (bottom) for the three HPC
architectures and three grid sizes analysed in the solitary wave benchmark.

The resulting speed up ($S_{n}$) and parallel efficiency ($E$) for the
solitary wave case are shown in Figure 13 (mid and bottom). For Ne = 1.5,
$\mathtt{Hydro3D}$ features less than optimal strong scalability with $S_{n}$
being below the ideal linear performance increase and values of $E$ steadily
decreasing with increasing node count for all systems. The communication-to-
computation ratio grows almost linearly from 7% to 19% when increasing from 1
to 64 nodes irrespective of the problem size, indicating that communications
are not the major overhead. Instead, resolving the LSM with the present grid
resolutions requires 15 iterations in the $\phi$ re-initialisation which
triggers the overhead from computing WENO reconstructions [46], as observed in
Figure 14.

Increasing the resolution up to Ne = 12.3 provides an improvement in
scalability, as shown by the speed up curve, and results for Ne = 98.3 show
that SKL follows a linear increase in $S_{n}$, attaining values of $E$ above
one. TX2 and EPYC feature an almost linear $S_{n}$ increase below the ideal
threshold and steadily decrease in parallel efficiency with increasing node
count. Overall, these results suggest that Hydro3D demonstrates good
scalability depending on the overhead from the LSM re-initialisation that
appears as the most computationally demanding task in the solitary wave
benchmark.

The breakdown of the computing time (again excluding communications) is
presented in Figure 14 for the smaller and larger problem sizes using two
nodes, including time spent on the LSM advection equations (TLS) with the time
spent on computing the WENO weights for $\phi$ accounted separately as
TLS_$\omega_{k}$ (see Alg. 1). For the small problem size, Figure 14(a),
TLS_$\omega_{k}$ represents between 62–77% of the total computing time, with
the rest of the LSM workload, i.e. TLS, in the range of 11–16%. For Ne = 98.3,
the relative contribution from the pressure solver (TP) increases to 17–26% of
the total cost, in turn decreasing the percentage attributed to computing LSM.
These results outline the high additional computational cost of LSM for multi-
phase flows, which normally requires very fine grid resolutions and time steps
to guarantee numerical stability. Considering the workload distribution, the
code spends a similar proportion of time in the various subroutines when
running on the SKL and EPYC systems, whilst on TX2 the overall time spent in
convection, diffusion, velocity update and WENO weights for the convection
scheme are much reduced.

Figure 14: Breakdown of the main computing components for the solitary wave
benchmark case when running on two nodes. Comparison between TX2 (left), SKL
(centre) and EPYC (right) for meshes Ne = 1.5 (a) and 98.3 (b). TP = Poisson
pressure solver, T$\omega_{k}$ = WENO weights, TD = diffusion, TC =
convection, Tup = velocity update, TLS = LSM, and TLS_$\omega_{k}$ WENO
weights in LSM.

Profiles of communication overhead over the total runtime and MPI functions
for this benchmark are presented in Figure 15 for simulations with 2, 8 and 25
nodes and the three Ne adopted. The overall contribution of the communication
operations in the simulation is small for 2 and 8 nodes and slightly higher
for 25 nodes. Only for Ne = 1.5 does the code on TX2 experience an increase in
communication overheads analogous to the results obtained for the other
benchmark cases. The MPI directive consuming most of the communication time is
the collective MPI_ALLREDUCE required during the LSM, as indicated in Alg. 1.
The synchronisation primitive MPI_WAIT is the second MPI function consuming
most communication time, which represents approx 2–6% of the total runtime for
each case i.e., $\mathtt{Hydro3D}$ spends very little time constrained by
communication in this benchmark.

Figure 15: Comparison of the percentage of runtime spent on communication
(top) when running on 2, 8 and 25 nodes of the three systems; and distribution
of MPI functions for the cases with SKL (red line) and EPYC (green
line)systems. Solitary Wave benchmark case.

##  6 Conclusions

This paper compares the performance of three cluster systems characterised by
distinctive node architectures, namely Intel Skylake (SKL), ARMv8.1 ThunderX2
(TX2) and AMD EPYC-Rome (EPYC), in applications based on a state-of-the-art,
in-house Computational Fluid Dynamics (CFD) code, $\mathtt{Hydro3D}$. This
study investigated the code’s multi-node performance and scalability in three
benchmarks comprising several problem sizes and using up to 8,000 cores. The
benchmark cases considered were the lid-driven cavity flow, the Taylor-Green
vortex and the propagation of a solitary wave in a numerical tank, designed to
perform computations of increasing difficulty and expense. The computing
workload is varied changing the computation of the convective fluxes with 4th-
order central differences and 5th-order WENO schemes for the first two cases
respectively, whilst in the wave simulation the level-set method equations are
computed, also with the 5th-order WENO scheme.

The EPYC-based system provides the fastest computing times for all benchmark
cases irrespective of the problem size, attributed to the 64 cores per node in
comparison to the 40-core SKL nodes and larger clock frequency than the TX2
processors. TX2 performance for small problem sizes features sub-optimal
strong scalability as the communications overhead notably increases. However,
increasing the number of grid cells leads to reduced differences between TX2
and SKL with both demonstrating similar runtime values, although these were
greater than those obtained with EPYC processors. For the solitary wave case,
TX2 achieves better computational performance than SKL processors for all
problem sizes. Nonetheless, a notable feature in most cases is that
$\mathtt{Hydro3D}$ needs to improve its strong scalability on the TX2,
especially when the number of grid cells is relatively small, resulting from
an overhead in blocking MPI communications.

This study suggests that the major drawback when adopting high-order WENO
schemes is their larger computing overhead compared to central differences.
Results show that in simulations with the WENO scheme the code’s performance
on TX2 is similar to that on SKL processors. This suggests that CFD codes
using finite volumes or finite elements in unstructured grids relying on mesh
connectivity can demonstrate good performance when running on TX2 processors,
as such schemes are computationally more demanding than central differences.

Results shown here are of interest to the HPC and CFD communities, both for
those working on direct-numerical and large-eddy simulation applications whose
computational cost is extremely high, as well as those doing engineering-
oriented simulations. The present study is based on the performance of a
structured-grid CFD code, using central differences as widely used numerical
schemes to compute the viscous and convective fluxes in LES and DNS.

Considering the maturity of the ARM ecosystem for HPC, we note that
$\mathtt{Hydro3D}$ required no modification of the codebase to either compile
or run on the TX2 system. The Cray compilers and other performance tools used
in these tests e.g., profilers and debuggers, show a production-level
maturity. From the user perspective, they offer a similar interface to other
commercial HPC tools and produced no issue or error during their use.

Future work will develop in two directions: first, a study of the energy usage
of these systems. The present study has centred on raw computing times, but
energy efficiency is of increasing concern in the field of HPC and this will
be explored with the balance between FLOP/s and Watts. Second, an increased
focus on single-node performance. The current work has concentrated on the
scalability and bottlenecks at the cluster/multi-node level and the next step
will be to explore in greater detail the performance of the code at the node
level.

## Acknowledgements

The authors would like to acknowledge the support of the Supercomputing Wales
project, which is partially funded by the European Regional Development Fund
(ERDF) via the Welsh Government. The Isambard project is funded by UK’s EPSRC
(EP/P020224/1), the GW4 alliance, the Met Office, Cray and Arm. This work was
partially funded by the GW4 Data Science research project ”Evaluating the
performance of a large-eddy simulation code in 64-bit ARM processors”. The
authors acknowledge Prof. Simon McIntosh-Smith for his support and Drs. Ade
Fewings, Thomas Green and James Price for helping in the setup and debugging
of the code.

## References

  * Jiménez [2020] J. Jiménez,  Computers and turbulence,  European Journal of Mechanics, B/Fluids 79 (2020) 1–11.
  * Stoesser [2014] T. Stoesser,  Large-eddy simulation in hydraulics: Quo Vadis?,  Journal of Hydraulic Research 52 (2014) 441–452.
  * Rodi [2017] W. Rodi,  Turbulence Modeling and Simulation in Hydraulics : A Historical Review,  Journal of Hydraulic Engineering 143 (2017) 1–20.
  * Sotiropoulos [2019] F. Sotiropoulos,  Hydraulic Engineering in the Era of Big Data and Extreme Computing: Can Computers Simulate River Turbulence?,  Journal of Hydraulic Engineering 145 (2019) 1–13.
  * Jiménez [2003] J. Jiménez,  Computing high-Reynolds-number turbulence: Will simulations ever replace experiments?,  Journal of Turbulence 4 (2003) 37–41.
  * Lopez-Novoa et al. [2015] U. Lopez-Novoa, A. Mendiburu, J. Miguel-Alonso,  A survey of performance modeling and simulation techniques for accelerator-based computing,  IEEE Transactions on Parallel and Distributed Systems 26 (2015) 272–281.
  * Guest et al. [2019] M. Guest, A. Elena, A. Chalk,  DL_POLY - A performance overview analysing, understanding and exploiting available HPC technology,  Molecular Simulation (2019).
  * Celani [2007] A. Celani,  The frontiers of computing in turbulence: Challenges and perspectives,  Journal of Turbulence 8 (2007) 1–9.
  * Gray et al. [2012] A. Gray, I. Bethune, R. Kenway, L. Smith, M. Guest, C. Kitchen, P. Calleja, A. Korzynski, S. Rankin, M. Ashworth, A. Porter, I. Todorov, M. Plummer, E. Jones, L. Steenman-Clark, B. Ralston, C. Laughton,  Mapping application performance to HPC architecture,  Computer Physics Communications 183 (2012) 520–529.
  * McIntosh-Smith et al. [2019] S. McIntosh-Smith, J. Price, T. Deakin, A. Poenaru,  A performance analysis of the first generation of hpc-optimized arm processors,  Concurrency and Computation: Practice and Experience (2019) e5110.
  * Rajovic et al. [2014] N. Rajovic, A. Rico, N. Puzovic, C. Adeniyi-Jones, A. Ramirez,  Tibidabo: Making the case for an ARM-based HPC system,  Future Generation Computer Systems 36 (2014) 322–334.
  * Oyarzun et al. [2018] G. Oyarzun, R. Borrell, A. Gorobets, F. Mantovani, A. Oliva,  Efficient CFD code implementation for the ARM-based Mont-Blanc architecture,  Future Generation Computer Systems 79 (2018) 786–796.
  * Sotiropoulos [2015] F. Sotiropoulos,  Hydraulics in the era of exponentially growing computing power,  Journal of Hydraulic Research 53 (2015) 547–560.
  * Sagaut [2009] P. Sagaut, Large Eddy Simulation for Compressible Flows, Springer Verlag, Berlin, 2009\.
  * Mazzuoli et al. [2019] M. Mazzuoli, A. Kidanemariam, M. Uhlmann,  Direct numerical simulation of ripples in an oscillatory flow,  Journal of Fluid Mechanics 863 (2019) 572–600.
  * Kang and Sotiropoulos [2012] S. Kang, F. Sotiropoulos,  Numerical modeling of 3D turbulent free surface flow in natural waterways,  Advances in Water Resources 40 (2012) 23–36.
  * Kara et al. [2015] S. Kara, M. C. Kara, T. Stoesser, T. W. Sturm,  Free-Surface versus Rigid-Lid LES Computations for Bridge-Abutment Flow,  Journal of Hydraulic Engineering 141 (2015) 04015019.
  * Xie [????] Z. Xie,  A two-phase flow model for three-dimensional breaking waves over complex topography,  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences (????) 20150101.
  * Krosronejad et al. [2019] A. Krosronejad, M. Arabi, D. Angelidis, E. Bagherizadeh, K. Flora, A. Farhadzadeh,  Comparative Hydrodynamic Study of Rigid-Lid and Level-Set Methods for LES of Open-Channel Flow,  Journal of Hydraulic Engineering 145 (2019) 04018077.
  * Cevheri et al. [2016] M. Cevheri, R. McSherry, T. Stoesser,  A local mesh refinement approach for large-eddy simulations of turbulent flows,  International Journal for Numerical Methods in Fluids 82 (2016) 261–285.
  * Ouro et al. [2019] P. Ouro, B. Fraga, U. Lopez-Novoa, T. Stoesser,  Scalability of an Eulerian-Lagrangian large-eddy simulation solver with hybrid MPI/OpenMP parallelisation,  Computers & Fluids 179 (2019) 123–136.
  * Ouro et al. [2018] P. Ouro, T. Stoesser, B. Fraga, U. Lopez-Novoa, Hydro3D, 2018.
  * Kara et al. [2015] M. C. Kara, T. Stoesser, R. McSherry,  Calculation of fluid-structure interaction: methods, refinements, applications,  Proceedings of the ICE - Engineering and Computational Mechanics 168 (2015) 59–78.
  * Ouro and Stoesser [2017] P. Ouro, T. Stoesser,  An immersed boundary-based large-eddy simulation approach to predict the performance of vertical axis tidal turbines,  Computers & Fluids 152 (2017) 74–87.
  * Ouro et al. [2019a] P. Ouro, V. Muhawenimana, C. A. M. E. Wilson,  Asymmetric wake of a horizontal cylinder in close proximity to a solid boundary for Reynolds numbers in the subcritical turbulence regime,  Physical Review Fluids 4 (2019a) 104604.
  * Ouro et al. [2019b] P. Ouro, L. Ramírez, M. Harrold,  Analysis of array spacing on tidal stream turbine farm performance using Large-Eddy Simulation,  Journal of Fluids and Structures 91 (2019b) 102732\.
  * Liu et al. [2017] Y. Liu, T. Stoesser, H. Fang, A. Papanicolaou, A. G. Tsakiris,  Turbulent flow over an array of boulders placed on a rough, permeable bed,  Computers and Fluids 158 (2017) 120–132.
  * McSherry et al. [2018] R. McSherry, K. Chua, T. Stoesser, S. Mulahasan,  Free surface flow over square bars at intermediate relative submergence,  Journal of Hydraulic Research 56 (2018) 825–843.
  * Chua et al. [2019] K. V. Chua, B. Fraga, T. Stoesser, S. Ho Hong, T. Sturm,  Effect of Bridge Abutment Length on Turbulence Structure and Flow through the Opening,  Journal of Hydraulic Engineering 145 (2019) 04019024.
  * Fraga et al. [2016] B. Fraga, T. Stoesser, C. C. Lai, S. A. Socolofsky,  A LES-based Eulerianâ€“Lagrangian approach to predict the dynamics of bubble plumes,  Ocean Modelling 97 (2016) 27–36.
  * Nicoud and Ducros [1999] F. Nicoud, F. Ducros,  Subgrid-scale stress modelling based on the square of the velocity gradient tensor,  Flow, Turbulence and Combustion 62 (1999) 183–200.
  * Breuer [1998] M. Breuer,  Large eddy simulation of the subcritical flow past a circular cylinder: Numerical and modeling aspects,  International Journal for Numerical Methods in Fluids 28 (1998) 1281–1302.
  * Costa and Don [2007] B. Costa, W. Don,  High order Hybrid central—WENO finite difference scheme for conservation laws,  Journal of Computational and Applied Mathematics 204 (2007) 209–218.
  * Fernández-Fidalgo et al. [2018] J. Fernández-Fidalgo, X. Nogueira, L. Ramírez, I. Colominas,  Ana posteriori, efficient, high-spectral resolution hybridfinite-difference method for compressible flows,  Computer Methods in Applied Mechanics and Engineering 335 (2018) 91–127.
  * Shu and Osher [1989] C.-W. Shu, S. Osher,  Efficient implementation of essentially non-oscillatory shock-capturing schemes, II,  Journal of Computational Physics 83 (1989) 32–78.
  * Sussman [1994] M. Sussman,  A level set approach for computing solutions to incompressible two-phase flow,  Journal of Computational Physics (1994) 146–159.
  * Chorin [1968] A. Chorin,  Numerical solution of the Navier-Stokes equations,  Math Comput 22 (1968) 745–762.
  * Fehn et al. [2018] N. Fehn, W. Wall, M. Kronbichler,  Efficiency of high-performance discontinuous Galerkin spectral element methods for under-resolved turbulent incompressible flows,  Int J Numer Meth Fluids. 88 (2018) 32–54.
  * Fernández-Fidalgo et al. [2020] J. Fernández-Fidalgo, S. Clain, L. Ramírez, I. Colominas, X. Nogueira,  Very high-order method on immersed curved domains for finite difference schemes with regular Cartesian grids,  Computer Methods in Applied Mechanics and Engineering 360 (2020) 112782\.
  * Calderer et al. [2014] A. Calderer, S. Kang, F. Sotiropoulos,  Level set immersed boundary method for coupled simulation of air/water interaction with complex floating structures,  Journal of Computational Physics 277 (2014) 201–227.
  * Lee et al. [1982] J. Lee, J. Skjelbreia, F. Raichlen,  Measurements of velocities in solitary waves,  Journal of the Waterway, Port, Coastal and Ocean Division 108 (1982) 200–218.
  * Shende and Malony [2006] S. S. Shende, A. D. Malony,  The tau parallel performance system,  The International Journal of High Performance Computing Applications 20 (2006) 287–311.
  * Hackenberg et al. [2013] D. Hackenberg, R. Oldenburg, D. Molka, R. Schöne,  Introducing firestarter: A processor stress test utility,  in: 2013 International Green Computing Conference Proceedings, pp. 1–9.
  * McCalpin et al. [1995] J. D. McCalpin, et al.,  Memory bandwidth and machine balance in current high performance computers,  IEEE computer society technical committee on computer architecture (TCCA) newsletter 2 (1995).
  * Deakin et al. [2019] T. Deakin, S. McIntosh-Smith, J. Price, A. Poenaru, P. Atkinson, C. Popa, J. Salmon,  Performance portability across diverse computer architectures,  in: 2019 IEEE/ACM International Workshop on Performance, Portability and Productivity in HPC (P3HPC), pp. 1–13.
  * Tsoutsanis et al. [2018] P. Tsoutsanis, A. F. Antoniadis, K. W. Jenkins,  Improvement of the computational performance of a parallel unstructured WENO finite volume CFD code for Implicit Large Eddy Simulation,  Computers and Fluids 173 (2018) 157–170.