# FluTAS: A GPU-accelerated finite difference code for multiphase flows

Marco Crialesi-Esposito<EMAIL_ADDRESS>Nicolò Scapin<EMAIL_ADDRESS>Andreas D. Demou<EMAIL_ADDRESS>Marco Edoardo Rosti<EMAIL_ADDRESS>Pedro Costa<EMAIL_ADDRESS>Filippo Spiga<EMAIL_ADDRESS>Luca Brandt
<EMAIL_ADDRESS>Department of Engineering Mechanics, Royal Institute of
Technology (KTH), Stockholm, Sweden,  The Cyprus Institute, Nicosia, Cyprus,
Complex Fluids and Flows Unit, Okinawa Institute of Science and Technology
Graduate University (OIST), 1919-1 Tancha, Onna-son, Okinawa 904-0495, Japan,
Faculty of Industrial Engineering, Mechanical Engineering and Computer
Science, University of Iceland, Hjardarhagi 2-6, 107 Reykjavík, Iceland
NVIDIA Ltd, Cambridge (UK),  Department of Energy and Process Engineering,
Norwegian University of Science and Technology (NTNU), Trondheim, Norway.

###### Abstract

We present the Fluid Transport Accelerated Solver, FluTAS, a scalable GPU code
for multiphase flows with thermal effects. The code solves the incompressible
Navier-Stokes equation for two-fluid systems, with a direct FFT-based Poisson
solver for the pressure equation. The interface between the two fluids is
represented with the Volume of Fluid (VoF) method, which is mass conserving
and well suited for complex flows thanks to its capacity of handling
topological changes. The energy equation is explicitly solved and coupled with
the momentum equation through the Boussinesq approximation. The code is
conceived in a modular fashion so that different numerical methods can be used
independently, the existing routines can be modified, and new ones can be
included in a straightforward and sustainable manner. FluTAS is written in
modern Fortran and parallelized using hybrid MPI/OpenMP in the CPU-only
version and accelerated with OpenACC directives in the GPU implementation. We
present different benchmarks to validate the code, and two large-scale
simulations of fundamental interest in turbulent multiphase flows: isothermal
emulsions in HIT and two-layer Rayleigh-Bénard convection. FluTAS is
distributed through a MIT license and arises from a collaborative effort of
several scientists, aiming to become a flexible tool to study complex
multiphase flows.

###### keywords:

Multiphase flows, Volume-of-Fluid method, turbulence in multiphase flows,
High-performance computing, OpenACC directives.

††journal: Computer Physics Communications

Program summary

Program Title: Fluid Transport Accelerated Solver, FluTAS.
Developer’s repository link: https://github.com/Multiphysics-Flow-
Solvers/FluTAS.git.
Licensing provisions: MIT License.
Programming language: Fortran 90, parallelized using MPI and slab/pencil
decomposition, GPU accelerated using OpenACC directives.
External libraries/routines: FFTW, cuFFT.
Nature of problem: FluTAS is a GPU-accelerated numerical code tailored to
perform interface resolved simulations of incompressible multiphase flows,
optionally with heat transfer. The code combines a standard pressure
correction algorithm with an algebraic volume of fluid method, MTHINC [1].
Solution method: the code employs a second-order-finite difference
discretization and solves the two-fluid Navier-Stokes equation using a
projection method. It can be run both on CPU-architectures and GPU-
architectures.

##  1 Introduction

Multiphase flows are ubiquitous in many contexts, ranging from environmental
flows to industrial applications. The interaction between phases has a
prominent role in the formation and evolution of clouds [2], in sediment
transport [3, 4], oceanic sprays and bubbles generation [5] and more in
general it represents one of the grand challenges of environmental fluid
mechanics [6]. These flows are also crucial in several industrial
applications, such as pharmaceutical, transportation, food processing and
power generation [7]. From a theoretical point of view, the main difficulty
when analyzing multiphase flows relies on their multiscale nature, since the
length-scale of the interface is of the order of the mean-free path while in
most applications the typical length-scale is several orders of magnitude
larger ($\sim 10^{5}-10^{6}$). This huge separation of scales is magnified
when dealing with turbulent multiphase flows, which further broadens the
spectrum of length-scales, thus making unfeasible any attempt to bridge all of
them in a single and unique framework. For this reason, all the tools
developed so far, both of experimental and numerical nature, have focused on
only a portion of the scale spectrum while the remaining part is modelled or
neglected.

As regards multiphase turbulence, where most of our interests and applications
are, both experimental investigations and numerical simulations have been
extensively used in the last thirty years and have led to important
contributions in a variety of problems and configurations: to name a few,
particle laden flows and sediment transport discussed in [8, 4], bubbly and
droplet flows reviewed in [9, 10, 11], oceanic sprays and bubbles as
highlighted in [5]. Nevertheless, as already discussed in [10] and despite the
recent progress in the instantaneous measurement of bubble/droplet shape [12,
13, 14], there is still a lack of experimental data for the measurement of the
instantaneous velocity fields of both carried and dispersed phase as well as
for the turbulent kinetic energy and dissipation near the interface locations.
These limitations disappear when dealing with numerical simulations and,
therefore, in the last decades interface resolved simulations of multiphase
flows have become a central investigation tool. Despite the advantages,
numerical simulations are still limited to simple configurations and moderate
scale separation, and pose the challenge of the choice of the proper method to
fully resolve the two-phase interface. As discussed in [15], there is now
consensus that numerical methods suitable to perform interface-resolved
simulations of multiphase flow should have the following properties: i) be
able to enforce mass, momentum and kinetic energy conservation at discrete
level, ii) allow mismatches in the material properties, whose magnitude
depends on the application, and iii) handle complex and possibly arbitrary
topological changes. Among the four groups of numerical methods for multiphase
flows, Front-Tracking (FT) [16], Volume-of-Fluid (VoF) [17], Phase Field (PFM)
[18], Level-set (LS) [19], it exists at least a variant of each which
possesses the aforementioned numerical properties, giving some freedom to
researchers and scientists on the choice of their preferred numerical tool
(see [17, 20, 21] for a review).

Nevertheless, it is becoming ever more clear that another desirable property
of any numerical method is its straightforward adaptation to be able to run
massively parallel simulations, especially on accelerated architectures. With
the increase in the computing power driven by Graphics Processing Units (GPUs)
[22], several HPC centers are now shifting towards GPU-only and GPU-
accelerated architectures. This trend is making the GPU-parallelization of
numerical codes for fluid mechanics a mandatory requirement rather than a
simple advantage. This effort has been already taken for single-phase codes,
where at least three open-source codes for incompressible and fully
compressible simulations are able to run on accelerated architectures: AFiD
[23], STREAmS [24] and the accelerated version of CaNS [25]. Conversely, on
the multiphase counterpart, despite the large availability of CPU-based open
source codes, PARIS Simulator [26], TBFsolver [27], FS3D [28], NGA2 [29],
Basilisk [30] and MFC [31] to name few, limited effort has so far been devoted
to their adaptation to hybrid architectures.

In this work, we aim to fill this gap and present FluTAS (Fluid Transport
Accelerated Solver), a code for massive Direct Numerical Simulations on multi-
GPU and multi-CPU architectures targeting incompressible multiphase flows,
optionally with heat transfer. The numerical solution of these flows is
typically performed using finite-difference methods in a staggered variable
arrangement, and it involves the solution of a Poisson equation to enforce the
constraints on the velocity divergence. In this context, FluTAS uses as basis
the Navier-Stokes solver CaNS [32] and its GPU extension [25], whose key
feature is a general implementation incorporating all the possible homogeneous
pressure boundary conditions that can benefit of the FFT-based elliptic
solvers [33]. The single-phase Navier-Stokes solver is extended to a two-fluid
code using the one-domain formulation [17, 20] and coupled with the algebraic
VoF MTHINC [1] to capture the two-phase interface. This method combines the
exact mass conservation properties of certain geometric VoF methods with the
reduced number of local operations for the interface reconstruction of the
algebraic VoFs, making it a good candidate for properly exploiting hybrid and
accelerated architectures. The version available in our group has been
validated in [34] and extensively employed in a different variety of
multiphase configurations, both for laminar [35, 36, 37] and turbulent [38,
39, 40, 41] flows. Note that it has been extended to phase changing flows [42]
and also to handle weakly compressible multiphase flows (low-Mach
approximation) [43, 44].

This paper is organized as follows. In section 2, we introduce the governing
equations for the incompressible two-fluid system. The discretization details
of the VoF method, energy equation and Navier-Stokes solver are provided in
section 3, whereas the standard benchmarks for code validation are discussed
in section 4. Next, the parallelization for the GPU acceleration is presented
together with the scaling tests in section 5 and in section 6. The code
potentialities are shown in two demanding simulations of multiphase
turbulence: emulsions in homogeneous isotropic turbulence (HIT) and two-phase
thermal convection (see section 7). Finally, main conclusions and future
perspectives are summarized in section 8.

##  2 Governing equations

We consider a two-phase system of immiscible incompressible Newtonian fluids
(e.g., a gas-liquid system). The two phases are bounded by an infinitesimally
small interface, through which momentum and energy can be transferred. To
describe the system, we define a phase indicator function $H$ distinguishing
the two phases at position $\mathbf{x}$ and time $t$:

$H(\mathbf{x},t)=\begin{cases}1\hskip 14.22636pt\text{if
$\mathbf{x}\in\Omega_{1}$}\mathrm{,}\\\ 0\hskip 14.22636pt\text{if
$\mathbf{x}\in\Omega_{2}$}\mathrm{,}\end{cases}$ (1)

where $\Omega_{1}$ and $\Omega_{2}$ are the domains pertaining to phases $1$
and $2$. We can use $H$ to define the thermophysical properties in the whole
domain $\Omega=\Omega_{1}\cup\Omega_{2}$ as follows:

$\xi(\mathbf{x},t)=\xi_{1}H(\mathbf{x},t)+\xi_{2}(1-H(\mathbf{x},t))\mathrm{,}$
(2)

where $\xi_{i}$ ($i=1,2$) can be the mass density $\rho_{i}$, the dynamic
viscosity $\mu_{i}$, the thermal conductivity $k_{i}$ or the specific heat
capacity at constant pressure $c_{p,i}$. Hereafter, unless otherwise stated,
thermophysical quantities not specifically referring to one of the phases are
defined from eq. (2). The evolution of the indicator function is governed from
the following topological equation:

$\dfrac{\partial H}{\partial
t}+\nabla\cdot\left(\mathbf{u}_{\Gamma}H\right)=H\nabla\cdot\mathbf{u}_{\Gamma}\mathrm{,}$
(3)

where $\mathbf{u}_{\Gamma}$ is the interface velocity. In absence of phase
change, the one-fluid velocity $\mathbf{u}$ is continuous across the interface
and therefore, it can be employed as interface velocity in equation (3).

The equations governing the momentum and energy transport for the liquid and
gas phase are coupled through appropriate interfacial conditions [45],
reported below in the so-called one-fluid or whole-domain formulation, where
each transport equation is defined in $\Omega$ [20].

$\nabla\cdot{\mathbf{u}}=0\mathrm{,}$ (4)

$\rho\left[\dfrac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot\left(\mathbf{u}\mathbf{u}\right)\right]=-\nabla
p+\nabla\cdot\left[\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)\right]+\sigma\kappa\delta_{\Gamma}+\hat{\rho}\mathbf{g}\mathrm{,}$
(5)

$\rho c_{p}\left[\dfrac{\partial T}{\partial
t}+\nabla\cdot\left(\mathbf{u}T\right)\right]=\nabla\cdot\left(k\nabla
T\right)\mathrm{.}$ (6)

Here, $\mathbf{u}$ is the fluid velocity assumed to be continuous in $\Omega$,
$p$ is the hydrodynamic pressure, $T$ the temperature. In equation (5),
$\sigma$ is the surface tension, $\kappa$ the local interfacial curvature and
$\delta_{\Gamma}$ is a delta Dirac function, $\mathbf{g}$ is the gravity
acceleration and $\hat{\rho}$ is the volumetric density field modified to
account for the thermal effects in the gravity forces. Using the
Oberbeck–Boussinesq approximation, $\hat{\rho}$ reads as:

$\hat{\rho}=\rho_{1,r}\left[1-\beta_{l}\left(T-T_{r}\right)\right]H+\rho_{2,r}\left[1-\beta_{g}\left(T-T_{r}\right)\right](1-H)\mathrm{,}$
(7)

where $\rho_{i=1,2,r}$ are the reference phase densities and $\beta_{i=1,2}$
are the liquid and gas thermal expansion coefficients.

##  3 Numerical methodology

The numerical solution of the governing equations (3), (4), (5) and (6)
presented in section 2 is addressed on a fixed regular Cartesian grid with
uniform spacing $\Delta x$, $\Delta y$ and $\Delta z$ along each direction. A
marker-and-cell arrangement is employed for velocity and pressure points [46],
whereas all scalar fields are defined at the cell centers. Each time-step, the
governing equations are advanced in time by $\Delta t^{n+1}=t^{n+1}-t^{n}$,
with the previous time-step indicated with $\Delta t^{n}=t^{n}-t^{n-1}$.
Hereafter, we present the numerical discretization of the governing equations,
following the same order in which they are solved.

###  3.1 Volume of fluid: the MTHINC method

The first step of the time-marching algorithm consists in the interface
reconstruction and its subsequent advection. As previously mentioned, these
tasks are addressed within a fully Eulerian framework using a volume-of-fluid
(VoF) method. From a numerical point of view, this consists first in defining
the volume fraction $\phi$ in each cell of the computational domain as:

$\phi=\dfrac{1}{V_{c}}\int_{V_{c}}H(\mathbf{x},t)dV_{c}\mathrm{,}$ (8)

with $V_{c}=\Delta x\Delta y\Delta z$. Next, equation (3) is written in terms
of volume fraction as:

$\dfrac{\partial\phi}{\partial
t}+\nabla\cdot\left(H\mathbf{u}\right)=\phi\nabla\cdot\mathbf{u}\mathrm{.}$
(9)

The distinct feature of each class of VoF method lies in the way $H$ is
approximated. In this work we employ the algebraic volume-of-fluid method
based on the Multi-dimensional Tangent Hyperbola reconstruction, MTHINC [1],
whose central idea is to approximate $H$ with a hyperbolic tangent:

$H(\tilde{x},\tilde{y},\tilde{z})=\frac{1}{2}\left[1+\tanh(\beta_{th}\left(\mathcal{T}(\tilde{\mathbf{x}})+d_{th}\right)\right]\mathrm{,}$
(10)

where $\beta_{th}$, $d_{th}$ are the sharpness and the normalization
parameter, respectively, and $(\tilde{x},\tilde{y},\tilde{z})$ a local
coordinate system $\tilde{\mathbf{x}}=\left[\left(x-0.5\right)/\Delta
x,\left(y-0.5\right)/\Delta y,\left(z-0.5\right)/\Delta z\right]$. Employing
equation (10) has two distinct advantages with respect to a piecewise
approximation, commonly employed in the geometric VoF methods. First, the
phase indicator $H$ can be approximated with a reconstructing polynomial
$\mathcal{T}$ of arbitrary order in a straightforward manner. Next, once
$\mathcal{T}$ is known, the resulting interface at the two-phase boundary has
smooth but controlled thickness (with the parameter $\beta_{th}$), which also
allows computing accurately the normal vector $\mathbf{n}$ and curvature
tensor $\mathbf{K}$ directly from $\phi$. More details about the choice of
$\mathcal{T}$ and the calculations of $d_{th}$, $\mathbf{n}$ and $\mathbf{K}$
are found in the original paper by Ii et al. [1], but for completeness we
include them in the appendix A with the numerical implementation details.
After the reconstruction step, the interface is advected using a directional
splitting approach [47, 48], which consists in evaluating the numerical fluxes
sequentially in each direction using, for each split, the latest estimation of
VoF field. Accordingly, three provisional fields $\phi_{i,j,k}^{p}$ (with
$p=[x,y,z]$) are first computed:

$\phi_{i,j,k}^{p}=\dfrac{\phi_{i,j,k}^{s}-\dfrac{1}{\Delta
l^{p}}\left[f_{+}^{p}(\phi_{i,j,k}^{s})-f_{-}^{p}(\phi_{i,j,k}^{s})\right]}{1-\dfrac{\Delta
t^{n+1}}{\Delta l^{p}}\left(u_{+}^{p}-u_{-}^{p}\right)^{n}}\mathrm{,}$ (11)

where $s=[n,x,y]$, $[\Delta l^{x},\Delta l^{y},\Delta l^{z}]=[\Delta x,\Delta
y,\Delta z]$, $[u^{x},u^{y},u^{z}]=[u,v,w]$ with $u_{\pm}^{p}$ the $p$-th
velocity component. The calculation of the numerical fluxes $f_{\pm}$ in
equation (11) are evaluated using the hyperbolic tangent approximation of $H$
as detailed in appendix A. Next, the divergence correction step is applied in
order to impose the volume conservation of both phases at a discrete level:

$\phi_{i,j,k}^{n+1}=\phi_{i,j,k}^{z}-\sum_{p=x,y,z}\dfrac{\Delta
t^{n+1}}{\Delta l^{p}}\phi_{i,j,k}^{p}(u_{+}^{p}-u_{-}^{p})^{n}\mathrm{.}$
(12)

With the above approach, mass conservation is ensured up to the accuracy with
which the divergence free condition (4) is satisfied. Accordingly, if direct
methods are employed to solve the Poisson equation, the mass of each phase
results to be conserved up to machine precision. Another approach with a
similar property has been introduced in [49]. However, in that case the
dilatation terms at the denominator of equation (11) are treated in an
explicit manner, while here in an implicit strategy is employed. This comes at
a cost of the final correction step, given by equation (12), but with the
advantage of not introducing additional time-step restrictions (apart the
convective one) in the advection of the color function.

###  3.2 Thermal effects

The next step of the time-marching algorithm consists in advancing the
temperature field using an explicit second-order Adams-Bashforth method:

$T^{n+1}=T^{n}+\Delta
t^{n+1}\left(f_{t,1}\mathcal{M}_{T}^{n}-f_{t,2}\mathcal{M}_{T}^{n-1}\right)\mathrm{,}$
(13)

where $f_{t,1}=(1+0.5\Delta t^{n+1}/\Delta t^{n})$ and $f_{t,2}=0.5\Delta
t^{n+1}/\Delta t^{n}$ are the coefficients of the Adams-Bashforth scheme. In
equation (13), the operator $\mathcal{M}_{T}$ accounts for the advection and
diffusion contribution and it is provided below in a semi-discrete form:

$\mathcal{M}_{T}^{n}=-\nabla\cdot(\mathbf{u}^{n}T^{n})+\dfrac{1}{\rho^{n+1}c_{p}^{n+1}}\nabla\cdot(k^{n+1}\nabla
T^{n})\mathrm{.}$ (14)

All the spatial terms in equation (14) are discretized with second-order
central schemes, except for the temperature convection term. The
discretization of the latter is based on the 5th-order WENO5 scheme, as in
reference [50].

###  3.3 Pressure correction algorithm

Once the energy equation has been advanced, the momentum equation is solved
with a second-order pressure correction [51], reported below in a semi-
discrete form:

$\displaystyle\left(\dfrac{\mathbf{u}^{\star\star}-\mathbf{u}^{n}}{\Delta
t^{n+1}}\right)=f_{t,1}\mathcal{M}_{\mathbf{u}}^{n}-f_{t,2}\mathcal{M}_{\mathbf{u}}^{n-1}+\dfrac{\left(\sigma\kappa\delta_{\Gamma}+\hat{\rho}\mathbf{g}\right)^{n+1}}{\rho^{n+1}}\mathrm{,}$
(15) $\displaystyle\mathbf{u}^{\star}=\mathbf{u}^{\star\star}-\dfrac{\Delta
t^{n+1}}{\rho_{0}}\left[\left(1-\dfrac{\rho_{0}}{\rho^{n+1}}\nabla\hat{p}\right)+\nabla
p^{n}\right]\mathrm{,}$ (16)
$\displaystyle\nabla^{2}\psi^{n+1}=\dfrac{\rho_{0}}{\Delta
t^{n+1}}\nabla\cdot\mathbf{u}^{\star}\mathrm{,}$ (17)
$\displaystyle\mathbf{u}^{n+1}=\mathbf{u}^{\star}-\dfrac{\Delta
t^{n+1}}{\rho_{0}}\nabla\psi^{n+1}\mathrm{,}$ (18) $\displaystyle
p^{n+1}=p^{n}+\psi^{n+1}\mathrm{,}$ (19)

where the operator $\mathcal{M}_{\mathbf{u}}^{n}$ and
$\mathcal{M}_{\mathbf{u}}^{n-1}$ in equation (15) includes the convective and
diffusive terms computed at the current and previous time level, neglecting
the surface tension and gravity forces which are then included as source
terms. The spatial gradients in $\mathcal{M}_{\mathbf{u}}$ are discretized
with central schemes. The intermediate velocity $\mathbf{u}^{\star\star}$ is
then updated with the contribution from the terms due to the time-pressure
splitting, as in (16). Note that $\rho_{0}$ is the minimum value of the
density field in the computational domain and $\hat{p}$ represents the time-
extrapolated pressure between the current and the old time step, i.e.
$\hat{p}=(1+\Delta t^{n+1}/\Delta t^{n})p^{n}-(\Delta t^{n+1}/\Delta
t^{n})p^{n-1}$. Following [52] and contrary to [53, 54], the terms arising
from the pressure splittings are included in the prediction of the velocity
field (see eq. (16) before the imposition of the boundary conditions. This
approach has two distinct advantages. First, it represents an incremental
pressure projection which allows achieving an almost second-order accurate in
time pressure field [52]. Next, it ensures the consistency of the pressure
field near a solid boundary (i.e, $\mathbf{u}^{n+1}=\mathbf{u}^{\star}=0$),
where the pressure gradient component normal to the boundary (i.e.,
$\nabla_{\perp}\psi^{n+1}=0$) vanishes independently of the local density (see
eq. (18)).

Next, the constant coefficients Poisson equation (17) is solved with the
method of eigenexpansion technique that can be employed for different
combination of homogeneous pressure boundary conditions [33]. Finally, the
velocity field is corrected as in equation (18) in order to impose the
divergence constrain (i.e., solenoidal velocity field) and the pressure
updated as in equation (19).

####  3.3.1 Poisson solver

The code uses the FFT-based finite-difference direct solver developed and
implemented in the DNS code CaNS; see [32, 25]. The underlying numerical
approach dates back to the late 1970s [55, 33], and has regained popularity in
recent years, thanks to the improvements of hardware, and of and software
frameworks for collective data communications, provided by the MPI standard
and higher-level libraries like 2DECOMP&FFT. In a nutshell, the approach uses
Fourier-based expansions along two domain directions, which reduce the system
of equations resulting from the three-dimensional second-order finite-
difference Laplace operator (seven non-zero diagonals) to a simple,
tridiagonal system. These Fourier-based expansions depend on the boundary
conditions of the system, and can be computed using FFTs, some of them with
pre-/post-processing of the FFT input/output (see, e.g., [56]).
The FFT-based expansions are employed along directions $x$ and $y$, and the
resulting tridiagonal system along $z$ is then solved using Gauss elimination.
For calculations on CPUs, the method leverages the _guru_ interface of the
FFTW library [57], which allows for performing all possible combinations of
discrete transforms using the same syntax. On GPUs, the fast discrete cosine
and sine transforms have been implemented using real-to-complex/complex-to-
real FFTs from the cuFFT library, with pre- and post-processing of the input
and output signals to calculate the desired series expansion [56, 25]. We
refer to Refs. [32, 25] for more details on this method and its
implementation.
Concerning the parallelization of the method in a distributed-memory setting,
the FFT-based transforms and Gauss elimination steps require the data along
each direction to be local to each MPI task. The domain is decomposed using a
2D pencil decomposition, where collective _all-to-all_ communications are
required to transpose the orientation of the 2D data decomposition. These
transposes are performed using the 2DECOMP&FFT library [58], which was
modified to allow for GPU-GPU communication in [23, 25].
It is worth noting that, in line with the recent developments of CaNS, the
present method uses a default decomposition (i.e., “outside” the Poisson
solver) based on a partitioning along $y$ and $z$, resulting in $x$-aligned
pencils. This reduces the total number of data transposes to be performed
during the solution of the Poisson equation from $6$ to $4$. The approach has
been adopted for both CPUs and GPUs, and the required operations to solve the
Poisson equation are summarized as follows:

  1. 1.

perform forward FFT-based transforms along $x$;

  2. 2.

transpose $x$-to-$y$;

  3. 3.

perform forward FFT-based transforms along $y$;

  4. 4.

transpose $y$-to-$z$;

  5. 5.

solve tridiagonal system using Gauss elimination along $z$;

  6. 6.

transpose $z$-to-$y$;

  7. 7.

perform backward FFT-based transforms along $y$;

  8. 8.

transpose $y$-to-$x$;

  9. 9.

perform backward FFT-based transforms along $x$.

Moreover, for the GPU implementation, the solver explicitly reduces the number
of _all-to-all_ operations when the domain is not decomposed along $z$ (i.e.,
when a $x-y$ slab decomposition is prescribed). This effectively decreases the
number of collective operations from $4$ to $2$ (steps $2$ and $8$ above are
skipped). This is the approach adopted in the GPU runs presented here – due to
the higher memory bandwidth in GPUs, a slab decomposition suffices for
distributed-memory calculations with sufficiently small wall-clock time per
step. Explicitly skipping these two no-op resulted in a substantial reduction
in wall-clock time per step, and in an overall improvement in the parallel
scalability of the solver.

###  3.4 Complete solution algorithm

For clarity, a step by step description of the overall solution procedure is
presented in Algorithm 1.

Algorithm 1 Overall solution procedure

1:$\phi^{0},\ T^{0},\ \mathbf{u}^{0},\ p^{0}$ are initialized;

2:$\rho^{0}$, $\mu^{0}$, $k^{0}$ and $c_{p}^{0}$ are calculated using equation
(2) from $\phi^{0}$;

3:$n=0$ is set,

4:while ($t<t_{tot}$ $\parallel$ $n<N_{tot}$) do

5: Set $n=n+1$ and $\Delta t^{n+1}$;

6: $\phi^{n+1}$ is calculated from equation (11) and (12);

7: $\mathbf{n}^{n+1}$ and $\kappa^{n+1}$ are evaluated using the procedure
described in A;

8: $\rho^{n+1}$, $\mu^{n+1}$, $k^{n+1}$ and $c_{p}^{n+1}$ is calculated from
equation (2);

9: $T^{n+1}$ is calculated from Eq. (13);

10: $\mathbf{u}^{*}$ is calculated from Eq. (15) and Eq. (16);

11: $\psi^{n+1}$ is calculated from Eq. (17);

12: $\mathbf{u}^{n+1}$ is calculated from Eq. (18);

13: $p^{n+1}$ is computed from Eq. (19).

14:end while

15:End of simulation.

##  4 Validation

###  4.1 Two-dimensional Zalesak’s disk

The Zalesak problem represents a classical benchmark to assess the accuracy of
the interface capturing/tracking algorithm. It consists in the solid-body
rotation of a slotted disk immersed in an imposed two-dimensional velocity
field $\mathbf{u}=(0.5-y,x-0.5)$. The disk can be easily defined in a
Cartesian two-dimensional domain by setting the indicator function
$H_{i,j,k}^{0}$ equal to $1$ inside the following domain $\Omega_{H}$

$\Omega_{H}:\left[\left(x-0.5\right)^{2}+\left(y-0.75\right)^{2}\right]\leq
0.15^{2}\cap\left(|x-0.5|\geq 0.025\cup y\geq 0.85\right)\mathrm{.}$ (20)

The benchmark consists in comparing the deformation of the solid disk with
respect to the initial shape after one entire revolution. The VoF equation is
solved in a two-dimensional square domain $\Omega=[0,1]\times[0,1]$,
discretized with four different grid spacing $[\Delta x,\Delta
y]=[1/N_{x},1/N_{y}]$ with $N_{x}\times N_{y}=[32\times 32,64\times
64,128\times 128,256\times 256]$. Periodic boundary conditions are prescribed
in both directions. Simulations are conducted up to $t=2\pi$ (i.e., one
complete revolution of the slotted disk) using a constant time-step $\Delta
t=t/3200$. Note that this value has been chosen to ensure a stable time
integration for the highest grid resolutions cases (i.e., $256\times 256$) and
is employed for the coarser cases.

Figure 1: Deformation of the Zalesak’s slotted disk after $t=2\pi$ for
$\beta_{th}=2$ (left) and $\beta_{th}=3$ (right).

Figure 1 shows the final disk shape for different grid solutions and for two
sharpness parameters $\beta_{th}=2$ and $\beta_{th}=3$. Note that the highest
deviation from the initial shape are in the corner regions, where the high-
curvature regions are located. Moreover, the solution is weakly dependent on
the value of $\beta_{th}$ and deviations between the different employed
$\beta_{th}$ are visible only for the coarser simulations.
Finally, to assess the accuracy of the solution, we compute the $L_{1}$ norm
and the order of convergence as:

$L_{1}=\dfrac{1}{N_{x}N_{y}}\sum_{i=1}^{N_{x}}\sum_{y=1}^{N_{z}}|\phi(i,j)-\phi^{0}(i,j)|\mathrm{,}$
(21)
$n_{L1}=\dfrac{\log\left(\dfrac{L_{2,N}}{L_{1,N}}\right)}{\log(2)}\mathrm{,}$
(22)

where $L_{1,N}$ is the $L_{1}$-error using $N_{x}\times N_{y}$ grid points and
$L_{1,2N}$ is the $L_{1}$-error evaluated with $2N_{x}\times 2N_{y}$ grid
points.

Figure 2: Deformation of the Zalesak’s slotted disk after $t=2\pi$ for
$\beta_{th}=2$ (left) and $\beta_{th}=3$ (right).

Results are reported in figure 2, where an order of convergence between the
first and the second-order is achieved for $\phi$, almost independent of the
employed value of $\beta_{th}$.

###  4.2 Three-dimensional rising bubble

The rising bubble test case is a well-established numerical benchmark for
multiphase flows [59]. This test is presented here to showcase the ability of
the numerical tool to accurately capture the topological changes of a moving
interface. The flow is driven by the density difference between the two
phases, and is influenced by the viscosity difference and the surface tension.
The relevant dimensionless groups for this flow are the Reynolds number
Re$=\rho_{g}u_{r}l_{r}/\mu_{g}$, the Weber number
We$=\rho_{g}u_{r}^{2}l_{r}/\sigma$, the Froude number
Fr$=u_{r}/\sqrt{|\mathbf{g}|d_{0}}$, the density ratio
$\lambda_{\rho}=\rho_{l}/\rho_{g}$ and the viscosity ratio
$\lambda_{\mu}=\mu_{l}/\mu_{g}$. In these definitions, $l_{r}$ is a reference
length and $u_{r}$ the reference velocity. Moreover, $\sigma$ is the surface
tension coefficient, $\rho_{g}$ and $\rho_{l}$ the reference gas and liquid
densities, and $\mu_{g}$ and $\mu_{l}$ the reference gas and liquid dynamic
viscosity. Finally, $\mathbf{g}$ is the acceleration of gravity and $d_{0}$ is
the initial diameter of the spherical bubble.
Following the benchmark study [59], the values adopted for the dimensionless
groups are Re$=35$, We$=1$, Fr$=1$, $\lambda_{\rho}=10$ and
$\lambda_{\mu}=10$, setting $l_{r}=d_{0}$, $u_{r}=\sqrt{|\mathbf{g}|d_{0}}$
and the reference time $t_{r}=\sqrt{d_{0}/|\mathbf{g}|}$. The dimensions of
the computational domain are $l_{x}=l_{y}=2d_{0}$ and $l_{z}=4d_{0}$. The
acceleration of gravity acts along the $z$-direction. No-slip and no-
penetration boundary conditions are prescribed at the horizontal top and
bottom boundaries of the domain ($z$-normal) and periodic conditions are
prescribed at the vertical boundaries ($x$\- or $y$-normal). A uniform
Cartesian grid of 128$\times$128$\times$256 cells is used. Initially, stagnant
flow conditions are applied and the position of the center of mass of the
spherical bubble, denoted as $(x_{c}(t),y_{c}(t),z_{c}(t))$, is located at
$(d_{0},d_{0},d_{0})$. A constant time-step $\Delta t/t_{r}=2.8\times 10^{-4}$
is used to advance the solution in time.

Figure 3: Isosurfaces of $\phi=0.5$ at dimensionless times
$t\sqrt{|\mathbf{g}|/d_{0}}$ (a) 0, (b) 1.4, (c) 2.8 and (d) 4.2. Figure 4:
Bubble rise velocity as a function of time. Black solid line, present results;
red dashed line, reference results from [59].

Figure 3 shows the isosurfaces of $\phi=0.5$ at various time instances. It is
evident that as the initially spherical bubble rises, its surface topology
changes. To compare against the reference results from [59], figure 4 shows
the evolution of the bubble rising velocity $U_{c}$ in time. The bubble
velocity is defined as,

$U_{c}=\frac{\displaystyle{\int_{\Omega}\phi\ w\
d\Omega}}{\displaystyle{\int_{\Omega}\phi\ d\Omega}}\mathrm{,}$ (23)

| $t/t_{r}$ | $U_{c}/u_{r}$ | $A_{0}/A$
---|---|---|---
ref. [59] |  | 0.51013 | 0.97418
present | 1.4 | 0.51144 | 0.97892
dev. % |  | 0.26 | 0.49
ref. [59] |  | 0.49823 | 0.95925
present | 4.2 | 0.49512 | 0.96057
dev. % |  | 0.62 | 0.14
Table 1: Comparison of the dimensionless rise velocity $U_{c}/u_{r}$ and
bubble sphericity $A_{0}/A$ between reference and present results.

where $w$ is the vertical velocity component and $\Omega$ the volume of the
entire domain. After an initial period where the bubble accelerates, the rise
velocity reaches a maximum and then stabilizes. Figure 4 demonstrates an
excellent agreement between present and reference results. To further quantify
this agreement, table 1 presents benchmark quantities for comparison at
specific time instances. Besides the bubble velocity, the table shows the
bubble sphericity $A_{0}/A$, defined as the initial value of the bubble
surface area over the value at a later time. The deviation between reference
and present values is less than 1%.

###  4.3 Differentially heated cavity

To demonstrate the accuracy of the code in the presence of thermal effects,
this section considers the flow of air in a closed two-dimensional square
heated cavity. The cavity is heated and cooled by the vertical side walls
(y-normal), while the horizontal walls are adiabatic (z-normal). Within this
configuration, a circulation is formed and maintained by the ascending hot
fluid next to the heated wall and the descending cold fluid next to the cooled
wall.

(a)(b)

Figure 5: (a) Contour plot of the temperature field at $t/t_{r}=0.5$ (steady
state) for the differentially heated cavity test case, (b) Temporal evolution
of the wall-averaged Nusselt number on the heated wall. Black solid line,
present results; red dashed line, reference results from [60].

The flow is therefore purely thermally-driven and is characterized by the
Rayleigh number Ra$=|\mathbf{g}|\beta\Delta Tl_{r}^{3}/\left(\nu\alpha\right)$
and the Prandtl number Pr$=\nu/\alpha$. In these definitions, $\beta$ is the
fluid thermal expansion coefficient, $\nu$ is the fluid viscosity, $\alpha$ is
the fluid thermal diffusivity and $\Delta T=(T_{h}-T_{c})$ is the temperature
difference between the heated ($T_{h}$) and cooled ($T_{c}$) walls. Typically,
the height of the cavity is taken as the reference length ($l_{r}=L_{z}$),
while the reference velocity and time are defined as $u_{r}=\alpha/l_{r}$ and
$t_{r}=l_{r}^{2}/\alpha$. The case simulated here follows the setup presented
in several studies [61, 62, 60] with Ra$=10^{6}$ and Pr$=0.71$. The domain
boundaries are solid walls, and no-slip boundary conditions are applied. With
respect to the temperature field, constant temperature boundary conditions are
applied on the vertical walls and a zero temperature gradient along the normal
direction is applied on the horizontal walls. The domain is discretized in
space using a uniform Cartesian grid with $256\times 256$ cells. Initially,
the air in the cavity is stagnant and isothermal at a temperature
$T_{0}=T_{c}$. A constant time-step $\Delta t$ is used to advance the solution
in time, given by $\Delta t/t_{r}=5.0\times 10^{-7}$. Figure 5a) shows the
contour of the temperature field at $t/t_{r}=0.5$, at which point a steady
state has been reached. The temperature field is characterized by thin and
spatially developing thermal boundary layers next to the thermally active
vertical walls, and a stratified region at the central area of the cavity. The
heat transfer rate inside the cavity is expressed through the Nusselt number,
defined as,

$\mathrm{Nu}=\frac{hl_{r}}{k}=\frac{l_{r}}{\Delta
T}\mathbf{\nabla}T\big{|}_{w}\cdot\mathbf{n}_{w}\mathrm{,}$ (24)

where $h$ is the heat transfer coefficient, $k$ is the fluid thermal
conductivity, $\mathbf{\nabla}T\big{|}_{w}$ is the temperature gradient on any
of the thermally active vertical walls and $\mathbf{n}_{w}$ is the
corresponding unit normal vector on the wall. Figure 5b) shows the comparison
of the temporal evolution of the wall-averaged Nusselt number
$\left<Nu\right>_{z}$ on the heated wall between the present and reference
results from [60]. It is evident that the present results are in excellent
agreement with the reference solution for the entire duration of the
simulation. Furthermore, Table 2 presents the comparison of key benchmark
quantities at steady state, confirming the agreement between present and
reference results.

| $V_{max}/u_{r}$ | $W_{max}/u_{r}$ | $Nu_{max}$ | $Nu_{min}$ | $\left<Nu\right>_{z}$
---|---|---|---|---|---
Ref. [60] | 64.85 | 220.6 | 17.58 | 0.9794 | 8.830
Present | 64.86 | 220.3 | 17.67 | 0.9773 | 8.843
% dev. | 0.02 | 0.14 | 0.51 | 0.21 | 0.14
Table 2: Comparison of key benchmark quantities at steady state for the
differentially heated cavity test case. $V_{max}$ is the maximum horizontal
velocity along the vertical mid-plane ($y=0.5l_{r}$), $W_{max}$ is the maximum
vertical velocity along the horizontal mid-plane ($z=0.5l_{r}$), $Nu_{max}$
and $Nu_{min}$ are the maximum and minimum values of the Nusselt number on the
heated wall, and $\left<Nu\right>_{z}$ is the averaged Nusselt number value on
the heated wall.

##  5 Code parallelization and GPU acceleration

###  5.1 Domain decomposition

The code is designed to run both in multi-CPU and multi-GPU architectures. For
domain decomposition, both slab (1D) and pencil (2D) are allowed [32] through
the library 2DECOMP [58]. The type of decomposition can be implicitly set via
the 2-component array dims (e.g. $[1,n]$ for slabs and $[n,m]$ for pencils) in
one input file dns.in. The pencil/slab orientation can be arbitrarily chosen
as in CaNS, via the preprocessor flags -D_DECOMP_X, -D_DECOMP_Y and
-D_DECOMP_Z which set the direction over which the domain is not decomposed.
This flexibility allows improving the efficiency of both the CPU and the GPU
implementations. For CPU, using pencils allows increasing the number of
processes used per execution (i.e. up to $N^{2}$ for $n_{x}=n_{y}=n_{z}=N$),
hence reducing the time to solution. In the GPU implementation, only the
z-pencil and x-slabs decompositions are allowed. It is recommendable to use
x-slabs on GPU (i.e. compiling with -D_DECOMP_X and using dims$=[1,n]$) as
this implementation reduces the number of _all-to-all_ calls to the minimum,
hence reducing GPU-GPU communication and improving performances on multi-nodes
runs.

###  5.2 Code parallelization

The parallelization is performed using MPI. When GPU acceleration is enabled,
MPI allocates one rank for each GPU. The code assumes the chosen MPI library
is “CUDA-aware”, meaning GPU data is directly passed to MPI function call and
the MPI implementation takes care of moving the data in the most efficient
way. If available, GPU-to-GPU communication can leverage NVIDIA NVLink which
is a physical GPU-to-GPU interconnection known to have higher bandwidth (at
least one order of magnitude) than Infiniband. Throughout the code, all nested
for-loops, i.e. iterations over all the domain points, are accelerated on GPUs
using OpenACC [63], a portable standard directive-based programming model that
can execute code on multi-core CPUs as well as accelerators like NVIDIA GPU.
Such offload is not used for CPU-only compilation and execution. To execute
FluTAS, the platform needs to support NVIDIA Unified Memory, which has two
main advantages:

  * 1.

the ability of allocating and managing more GPU allocated memory than what is
physically present on the device;

  * 2.

the ability to avoid explicitly to handle data movements Host-to-Device and
Device-to-Host, leaving the runtime do the work for the developers.

Both features are used in the code and proved crucial for an efficient GPU
acceleration.

##  6 Code performance

We now present an analysis of the code performances on standard CPU-based and
accelerated GPU-based architectures. Tests on GPUs were performed on MeluXina
at LuxProvide (LXP, Luxembourg) [64] and Berzelius at National Supercomputer
Centre (NSC, Sweden) [65], while tests on CPUs were performed on Tetralith
also managed by NSC.

###  6.1 Weak and strong scaling

We first discuss the weak-scaling tests for a Rayleigh-Bénard problem with the
same set-up as it will be discussed in Section 7.2. For this test, we start
with a ”base” computational grid of $N_{x}\times N_{y}\times N_{z}=1024\times
512\times 256$ grid points on 2 GPUs. Then, while keeping fixed $N_{x}$ and
$N_{y}$, we increase $N_{z}$ proportionally to the number of GPUs, resembling
a procedure of spatial ensemble-average (i.e. more structures simulated to
improve the convergence of the large-scale statistics). As discussed in
Section 5, we adopt a slab parallelization along the $z$ direction using the
-D_DECOMP_X compilation option, which reduces to $2$ the number of _all-to-
all_ operations. It is worth noticing that, although both HPC machines are
equipped with A100-40GB NVIDIA cards, Berzelius has 8 GPUs/node while MeluXina
has 4 GPUs/node. Moreover, the interconnection between GPUs is handled through
NVLink, whereas node-to-node connection is performed through Infiniband (IB),
known for having lower bandwidth and for operating on a different protocol.
Hence, IB is less straightforward to handle as it requires a more careful
configuration from a hardware and software perspective. This implies choosing
the right MPI configurations and selecting compatible communication libraries,
resulting in an efficiency that may vary significantly among different HPC
centers. For these reasons, to prove the non IB-dependent scaling and maximize
the GPU-to-GPU communication throughput, Berzelius was used to perform weak-
scaling tests within a node. On the other hand, MeluXina was used for
multiple-node tests to assess the IB-dependent scaling.111Berzelius uses a HGX
8-way platform, 8 A100 GPU connected all together, designed primarily for
heavy AI workloads. MeluXina uses a HGX 4-way platform, 4 A100 GPU connected
all together, which is better suited for scale-out HPC workloads. The majority
of GPU-accelerated HPC clusters used primarily for simulation workloads in
various scientific fields adopt the HGX 4-way configuration (including many
European HPC systems funded by EuroHPC).

(a)(b)

Figure 6: For the two-layer Rayleigh-Bénard convection problem discussed in
Section 7.2: a) code performance on MeluXina and Berzelius, b) slowdown due to
transpose operation. For each set of data, we compute $t_{w,min}$ as the time
per-timestep at the minimum number of GPUs tested over $t_{w}$, i.e. the time-
per-timestep at the specific number of GPUs.

Figure 6a) shows that weak-scaling is linear when bounded by NVLink
communications (i.e. no IB communications), as clearly supported by tests on
Berzelius. When IB communications are required (i.e. node-to-node data
transfer) the code performances decrease. It is worth noticing that, while an
increasing communication overhead is provided by node-to-node communication on
IB network, additional slowdown is caused by the slab parallelization. By
increasing the number of elements along $z$, more data need to be transferred
during the x-to-z transposes, further increasing the communication load. This
is clearly shown in figure 6b), where the slowdown is found to increase
proportionally to the number of GPUs.
Results of the strong scaling tests are reported in figure 7. Here we use two
different grids, i.e. $1024\times 512\times 1024$ (grid-1) and $1024\times
1024\times 1024$ (grid-2) for the Rayleigh-Bénard problem discussed in Section
7.2. Tests are performed on Meluxina and Berzelius as for the weak scaling.
While keeping the problem sizes fixed, the number of GPUs is progressively
increased up to a maximum of $128$, starting from $N_{GPU}=16$ which
represents the minimum amount required to fit the two computational domains in
the available GPU memory. Despite a speed-up is always achieved, the code
shows a progressive loss in performance, i.e. a reduction of the benefits
derived by increasing the number of GPUs. Note, however, that a larger number
of grid point (e.g. grid-2) leads to lesser performance loss, as a higher GPU
occupancy can be obtained.
The decrease in performance observed in figure 7 is caused by two factors: the
increase in communication among GPUs, and the reduction in local problem size,
which does not leverage the full compute capacity of each GPU. While these
effects are present in a strong scaling test, weak scaling allows us to
isolate the effects of multi-GPU communication while keeping a higher GPU
saturation. Thus, we argue that weak scaling represents a better tool to
identify communication bottlenecks on multiple GPUs. Conversely, the strong
scaling is more useful to estimate how much a fixed domain can be partitioned
while keeping an efficient use of the computational resources.

Figure 7: Strong scaling test performed on Berzelius (black-dashed lines) and
MeluXina (red-dashed lines) clusters for two different grids: $1024\times
512\times 1024$ (grid-1) and $1024\times 1024\times 1024$ (grid-2). The black
continuous line indicates the ideal behavior desired for the strong-scaling
test.

Overall, the previous analysis suggests an important guideline for the user:
in presence of unbalanced compute-vs-network architectures (e.g., node-to-node
networking connection less efficient than the connection among GPUs within the
same node), the optimal number of GPUs to be employed should be chosen as
close as possible to the minimum amount required to fit the computational
domain in the available GPU memory. Indeed, this is not always the case with
older HPC architectures using previous generations of GPU hardware, where
NVLink connections across GPUs inside a node were typically missing. For a
fixed problem size, modern cards with high compute throughput will exhaust the
required computation faster, leaving the remaining part of the computation as
communication bound. In older GPU hardware, the acceleration is lower and
communication becomes the dominant component affecting scalability for a
larger number of GPUs. Hence, best practice dictates to use the least possible
number of GPUs; on modern units with 80 GByte of HBM memory, if possible, it
is therefore convenient to use one single 8-way GPU node (like a DGX A100)
where all GPU are also connected via NVLink reducing communication overhead
dramatically.

###  6.2 CPU-GPU comparison

In conclusion, we perform a comparison between the code performances on a CPU
and a GPU architecture. It is worth mentioning that such comparison is
notoriously not trivial. First, no exact and standard procedures to compare
the two systems are established. Next, code performances may exhibit large
variations among different architectures and using the same hybrid CPU-GPU
node to perform tests on both may be misleading. CPU-only nodes and CPU-GPU
nodes are intrinsically different in terms of network configuration and
GPU/CPU interconnection, hence an unbiased test may not be performed directly
on hybrid architectures (as CPU-GPU cluster would hardly be used to perform
CPU-only jobs). Therefore, the following analysis has to be taken as a first-
approximation estimate.

(a)(b)

Figure 8: Comparison of code-section load percentage on the total simulation
time for GPUs (panel a) and CPUs (panel b). The different ”slices” represent
different code sections: 1) VoF (i.e. interface reconstruction and advection,
update of the thermophysical properties), 2) RHS (i.e. discretization of the
governing equations), 3) Transposes (i.e. transpose operation in the solver),
4) Solver (i.e. only Gaussian elimination) and others (i.e. correction step,
divergence/time-step checks, output and post-processing routines).

Here, we repeat the weak-scaling simulation with $n_{GPU}=8$ GPUs on Berzelius
on $n_{CPU}=512$ CPUs on Tetralith, in both cases employing a slab
parallelization along $z$. The test shows that for GPUs the average wall-clock
time per-timestep is $t_{8,GPU}=0.191$ $s$, while for CPUs $t_{512,CPU}=1.075$
$s$. This results in an equivalent number of GPUs
$n_{eq}=(t_{512,CPU}n_{CPU})/(t_{8,GPU}n_{GPU})\approx$ $359$.
Finally, a comparison in terms of computing-load percentage for each code
section is displayed in figure 8. As previously anticipated, transposes during
GPU simulations (panel a) represents more than half of the computing load. The
remaining parts, mainly composed of stencil operations enclosed in for loops,
largely benefit from GPU-offload, while in CPUs (panel b) these account for
more than 70 % of the total wall-clock time per time-step.

##  7 Applications

###  7.1 Emulsions in HIT

In many multiphase flows, the dispersed phase interacts with the surrounding
turbulence induced through large-scale stirring mechanisms. While turbulence
introduces a vast range of scales, usually spanning over several order of
magnitudes, the presence of an interface introduces further complexity,
offering alternative paths for energy transmission through scales and
generating poly-dispersed droplet/bubble distributions. Resolving the
interplay between all these mechanisms leads to an extremely complex scenario
to simulate numerically. Even in simplified conditions, represented by
Homogeneous and Isotropic Turbulence (HIT), the number of grid points could
rapidly exceed $N\geq 1024^{3}$, quite challenging for multiphase flows.
Furthermore, variations of density and viscosity, and variations of the
surface tension coefficient may introduce smaller scales as lower viscosity in
the dispersed phase may accelerate vortices.

$Re_{\lambda}$ | $We_{\mathcal{L}}$ | $\mathcal{V}$ | $\mathcal{L}$ | $\rho_{1}/\rho_{2}$ | $\mu_{1}/\mu_{2}$
---|---|---|---|---|---
137 | 42.6 | 0.06 | $\pi$ | 1 | 1
Table 3: Physical dimensionless parameters of the present configuration:
$Re_{\lambda}=u_{rms}\lambda/\nu_{2}$ and
$We_{\mathcal{L}}=\rho_{1}u_{rms}^{2}\mathcal{L}/\sigma$ are the corresponding
Taylor microscale Reynolds number and the large scale Weber number,
$\mathcal{L}$ is the large-scale at which turbulence is forced, $\mathcal{V}$
is the volume fraction (ratio between the volume occupied by the disperse
phase and the total volume). The simulation is performed at matching density
and viscosity (i.e., $\rho_{1}/\rho_{2}=\mu_{1}/\mu_{2}=1$). Figure 9: Render
for the iso-contour of the volume fraction value $\phi=0.5$. Lateral planes
show the modulus of the vorticity field.

Fully developed turbulence is usually reached for Taylor-Reynolds number
$Re_{\lambda}=u_{rms}\lambda/\nu\gtrsim 200$, where $\lambda$ is the Taylor
scale and $u_{rms}$ is the root-mean-square fluctuation velocity. In
multiphase flow simulations, these intensities are rarely reached as typical
values are $Re_{\lambda}\lesssim 100$. In this example, we present a
simulation of a turbulent emulsion at $Re_{\lambda}\approx 137$ performed on a
grid of $512^{3}$ on a cube of side length $2\pi$ and with the main physical
parameters reported in table 3. Turbulence is sustained at large scale using
the Arnold-Beltrami-Childress forcing (see [66, 67]) throughout the whole
simulation.
Figure 9 shows a render at statistically-stationary state. Turbulence is first
simulated in single phase; when statistically-stationary conditions are
reached, the dispersed phase is introduced through a random distribution of
droplets and let to develop until convergence is reached. This is monitored in
terms of droplet-size-distribution and one-dimensional spectra, see figure 10.

(a)(b)

Figure 10: Results at statistical-stationary state for turbulent emulsions
simulation. Panel (a) shows the droplet size distribution, where the droplet
diameter $d$ is normalized by the single-phase Kolmogorov scale $\eta_{sp}$.
Also in black, the $-10/3$ and $-3/2$ power-laws are shown (see [68, 69]).
Panel (b) shows the one-dimensional turbulent kinetic energy spectra,
comparing the single-phase and multiphase (i.e. emulsion) cases. Here, the
$-5/3$ law for the inertial range is shown, which applies for almost a decade.

Despite the simplicity of this configuration, HIT offers a relevant framework
to study multiphase turbulence in complex configurations. To reach properly
convergent statistics, several large-eddies turnovers are required,
corresponding to $\sim 10^{7}$ time-steps. Hence, GPU acceleration will be
invaluable to reach fully-developed turbulent conditions in future studies.

###  7.2 Two-layer Rayleigh-Bénard convection

Rayleigh–Bénard convection is the flow developed inside a fluid layer that is
heated from below and cooled from above. It is driven by the density
differences that arise due to the temperature variation inside the fluid. Even
though it is a seemingly simple configuration, it encapsulates rich physics
that are encountered in a range of engineering applications and physical
phenomena. Beyond the classical setting, the study of the two-layer
Rayleigh–Bénard variant is crucial from both a fundamental and an applied
standpoint. First, regardless of the application, there is always some
dissolved gas in every liquid. Therefore, it is inevitable that a gaseous
phase will be formed in any realistic natural convection flow. Second,
physical phenomena such as the convection in the earth’s mantle [70] or
engineering applications such as the heat transfer inside magnetic confinement
systems in fusion reactors [71] are more accurately modelled as two-layer
convection, where the two fluid layers are dynamically coupled.

Figure 11 shows a schematic representation of the domain used for the
numerical simulations, as used previously in [72], in two dimensions. The
bottom and top walls are modeled as solid isothermal surfaces at constant
temperatures of 328 K and 318 K respectively. The x- and y-directions are
considered periodic, and the aspect ratio between the horizontal and vertical
dimensions of the cavity is $\Gamma=L_{x}/L_{z}=L_{y}/L_{z}=2$. The
dimensionless parameters adopted are shown in Table 4. The property ratios
between the two fluids are considered equal to 1 (kinematic viscosity $\nu$,
thermal diffusivity $\alpha$, specific heat $c_{p}$ and thermal expansion
coefficient $\beta$), except the density ratio
$\lambda_{\rho}=\rho_{2}/\rho_{1}=0.1$. This mismatch in densities is the
reason behind the arrangement of the fluids in a two-layer configuration.
Preliminary simulations revealed that a grid of $1024\times 1024\times 512$
cells and a CFL number of 0.50 were adequate for obtaining grid- and time
step-independent solutions.

Figure 11: Schematic representation of the geometry used for the two-layer Rayleigh–Bénard convection. The bottom heated surface is depicted in red and the top cooled surface in blue. $\lambda_{\rho}=\frac{\rho_{2}}{\rho_{1}}$ | Pr=$\frac{\nu_{1}}{\alpha_{1}}$ | Ra=$\frac{|\mathbf{g}|\beta_{1}\Delta TL^{3}}{\nu_{1}\alpha_{1}}$ | We= $\frac{\rho_{1}|\mathbf{g}|\beta_{1}\Delta TL^{2}}{\sigma}$ | Fr=$\sqrt{\beta_{1}\Delta T}$
---|---|---|---|---
0.1 | 1 | $10^{6},10^{7},10^{8}$ | $100$ | $1$
Table 4: Dimensionless parameters adopted for the study of two-layer
Rayleigh–Bénard convection. All other property ratios are equal to 1.

(a)(b)(c)

Figure 12: Instantaneous temperature fields in the x–y plane for the two-layer
Rayleigh–Bénard convection. (a) $Ra=10^{6}$, (b) $Ra=10^{7}$, (c) $Ra=10^{8}$.

Figure 12 shows instantaneous temperature fields in the x–y plane for the
three Rayleigh numbers considered. With increasing Rayleigh number, the
thermal structures in the cavity become finer, indicating increased turbulent
activity. This increased thermal agitation is not enough to induce any
significant interface movement, even at the highest Rayleigh number
considered. Furthermore, in all three cases, the temperature drop at the
bottom fluid layer is much smaller than the top layer. This is explained by
the fact that the two layers have the same thermal diffusivity and specific
heat but different densities, leading to a top layer with a smaller thermal
conductivity ($k=\alpha\rho c_{p}$) compared to the bottom layer. Therefore,
the top wall conducts heat less effectively than the bottom wall, which
explains the larger temperature gradients at the top layer. Focusing on the
highest Rayleigh number case, figure 13 shows instantaneous temperature
contours for $Ra=10^{8}$. As in classical Rayleigh–Bénard convection, hot and
cold plumes are ejected from the bottom and top walls. In the presence of a
single fluid layer, these plumes would typically get organized in large scale
circulation structures, extending from the bottom to the top wall. Here, the
existence of two fluid layers changes the classical picture; the interface
acts as a barrier, confining the thermal plumes in each fluid layer. The
interface also acts as a thermal conductor, promoting the exchange of heat
between the two layers. More specifically, the hot plumes ascending from the
bottom wall are cooled off when reaching the interface, forming colder plumes
that travel downwards. The opposite is true at the top layer, where the
descending cold plumes are heated by the interface and hotter plumes emerge
from the interface traveling upwards. Figure 13 clearly illustrates this
behavior, revealing the existence of regions dominated by ascending plumes and
regions dominated by descending plumes, hinting towards the organization of
the flow in three-dimensional large scale circulation structures in each
layer.

(a)(b)

Figure 13: Instantaneous temperature contours for the two-layer
Rayleigh–Bénard convection at $Ra=10^{8}$. (a) Bottom half of the cavity with
red plumes at 327 K and orange plumes at 326.5 K. (b) Top half of the cavity
with cyan plumes at 323 K and blue plumes at 321 K.

##  8 Conclusions and further developments

We present the code FluTAS, a numerical framework tailored for direct
numerical simulations of multiphase flows, with the option of heat transfer,
able to run efficiently on CPU-based standard architectures and on GPU-based
accelerated machines. The open source version, released under MIT license,
includes a pressure-correction algorithm for two-phase flows extended with an
algebraic Volume-of-Fluid method (MTHINC) for capturing the interface
dynamics.
We provide here a description of the employed numerical algorithm, with
details on the solution of the governing equations and of the advection of the
interface. After presenting different validation benchmarks both in single and
multiphase configurations, we discuss the code performance focusing on two
aspects: i) its current limitation when the communications among GPUs in
different nodes are considerable less efficient than the communication among
GPUs within the same node, ii) its advantages compared to CPUs in terms of
”time-to-solution”. Finally, we report results from two configurations of
fundamental interests in multiphase turbulence: emulsions in homogeneous
isotropic turbulence and the two-layer Rayleigh-Bénard convection.
In the future, we aim to improve the code maintainability and portability
(both on CPU and GPU) and to release additional modules under development,
e.g. weak compressibility and phase change [43, 44]. Further efforts will be
devoted to enhance the code performance on multiple GPUs nodes, reducing the
current communication bottlenecks. To this end, a promising strategy is the
one proposed in [73], i.e., implement the solution of the tridiagonal system
for the third direction on a distributed memory. The main advantage of this
approach is the elimination of _all-to-all_ operations in the Poisson solver.
This improvement combined with future enhancements in the software frameworks
for collective data communications will allow tackling several multiphase
problems while keeping an efficient use of the computational resources.

## Acknowledgements

We would also like to thank the staff at CINECA, notably Massimiliano Guarrasi
and Fabio Pitarri for their help in the development of the GPU functionalities
at an early stage of the work. We thank Francesco De Vita for the help with
the first implementation of the VoF MTHINC method.

## Funding

M.C-E., N.S., A.D. and L.B. acknowledge the support from the Swedish Research
Council via the multidisciplinary research environment INTERFACE, Hybrid
multiscale modelling of transport phenomena for energy efficient processes and
the Grant No. 2016-06119. M.E.R. was supported by the JSPS KAKENHI Grant
Number JP20K22402. P.C. was supported by the University of Iceland Recruitment
Fund grant No. 1515-151341, _TURBBLY_. The computer time was provided by SNIC
(Swedish National Infrastructure for Computing), by the National
Infrastructure for High Performance Computing and Data Storage in Norway
(project No. NN9561K) and by CINECA (Marconi100) in Italy. Scaling tests on
multiple GPUs has been performed on Berzelius (under the project No.
Berzelius-2022-29 and operated by SNIC) and MeluXina (under the project number
No. EHPC-BEN-2022B01-027 with the EuroHPC Benchmark Access and managed by
LXP). M.E.R. acknowledges the computer time provided by the Scientific
Computing section of Research Support Division at OIST.

##  Appendix A Implementation details of MTHINC method

###  A.1 Calculation of the reconstructing polynomial

In the current work, we consider a polynomial up to second order and,
therefore, $\mathcal{T}$ can be expressed using the following general
quadratic form:

$\mathcal{T}(\tilde{\mathbf{x}})=\mathbf{c}\cdot\left[\mathcal{Q}_{\mathcal{T},1}\cdot\left(\tilde{\mathbf{x}}\cdot\tilde{\mathbf{x}}^{T}\right)^{T}\right]+\mathbf{c}\cdot\mathbf{c}^{R}\cdot\mathcal{Q}_{\mathcal{T},2}\cdot\left(\tilde{\mathbf{x}}\cdot\tilde{\mathbf{x}}^{R}\right)^{T}+\mathcal{L}_{\mathcal{T}}\cdot\tilde{\mathbf{x}}^{T}\mathrm{.}$
(25)

where $\tilde{\mathbf{x}}^{R}=(\tilde{y},\tilde{z},\tilde{x})$ while the
components of the flag vectors (equal to $0$ or $1$) $\mathbf{c}$ and
$\mathbf{c}^{R}$ will be discussed later in A.1.2. Accordingly, to determine
$\mathcal{T}$, one needs just to compute the components of the vectors
$\mathcal{Q}_{\mathcal{T},i=1,2}$ and $\mathcal{L}_{\mathcal{T}}$. This can be
done by first imposing that the first-order and second-order gradient of
$\mathcal{T}$, evaluated for $\tilde{\mathbf{x}}=\tilde{\mathbf{x}}_{c}$, are
equal to the normal vector $\mathbf{n}$ and the curvature tensor $\mathbf{K}$:

$\begin{cases}\left.\nabla\mathcal{T}\right|_{\tilde{\mathbf{x}}=\tilde{\mathbf{x}}_{c}}&=\mathbf{n}\mathrm{,}\\\
\left.\nabla\left(\nabla\mathcal{T}\right)\right|_{\tilde{\mathbf{x}}=\tilde{\mathbf{x}}_{c}}&=\mathbf{K}\mathrm{.}\end{cases}$
(26)

Expanding equation 25 and the equations 26 lead after some manipulation to a
unique expression for the components of the vectors
$\mathcal{Q}_{\mathcal{T},i=1,2}$ and $\mathcal{L}_{\mathcal{T}}$:

$\mathcal{Q}_{\mathcal{T},1}=\left[a_{Q,1}^{x},a_{Q,1}^{y},a_{Q,1}^{z}\right]=\dfrac{1}{2}\left[K^{xx},K^{yy},K^{zz}\right]\mathrm{,}$
(27)
$\mathcal{Q}_{\mathcal{T},2}=\left[a_{Q,2}^{x},a_{Q,2}^{y},a_{Q,2}^{z}\right]=\left[K^{xy},K^{yz},K^{xz}\right]\mathrm{,}$
(28)
$\displaystyle\mathcal{L}_{\mathcal{T}}=\left[a_{L}^{x},a_{L}^{y},a_{L}^{z}\right]=\,\,$
$\displaystyle\left[n^{x}-\dfrac{c^{x}}{2}\left(K^{xx}+c^{y}K^{xy}+c^{z}K^{xz}\right)\mathrm{,}\right.$
$\displaystyle\hskip
5.5485pt\left.n^{y}-\dfrac{c^{y}}{2}\left(K^{yy}+c^{x}K^{xy}+c^{z}K^{xz}\right)\mathrm{,}\right.$
$\displaystyle\hskip
5.5485pt\left.n^{z}-\dfrac{c^{z}}{2}\left(K^{zz}+c^{x}K^{xz}+c^{y}K^{yz}\right)\right]\mathrm{.}$
(29)

Note that to recover a linear reconstruction, the curvature tensor is set
identically zero and thus also the components of the vector
$\mathcal{Q}_{\mathcal{T},i=1,2}$ in equations (27) and (28). Accordingly,
from equation (29), $\mathcal{L}_{\mathcal{T}}=[n^{x},n^{y},n^{z}]$ and
equation (25) reduces to
$\mathcal{T}(\tilde{\mathbf{x}})=\mathcal{L}_{\mathcal{T}}\cdot\tilde{\mathbf{x}}^{T}$
which represents the equation of a plane in the three-dimensional Cartesian
space.

####  A.1.1 Normal vector and curvature calculation

Given the smooth nature of the color function, the normal vector $\mathbf{n}$
and the curvature tensor $\mathbf{K}$ can be computed directly from the
corresponding geometrical definitions:

$\mathbf{n}=\dfrac{\nabla\phi}{|\nabla\phi|}\mathrm{,}$ (30)
$\mathbf{K}=-\nabla\mathbf{n}\mathrm{,}$ (31)

where $\mathbf{m}=\nabla\phi=(m^{x},m^{y},m^{z})$. Following the Youngs’
method [74, 75], the three components $m^{x}$, $m^{y}$ and $m^{z}$ (the
partial derivative of $\phi$ in each direction) are computed by first
evaluating the derivatives at the cell corners and then each corner-value is
averaged to find $\partial\phi/\partial x$, $\partial\phi/\partial y$ and
$\partial\phi/\partial z$. Once $\mathbf{n}$ is known, the curvature tensor
components are computed using (31) whereas the geometrical curvature is
derived from the sum of the three diagonal components of $\mathbf{K}$, i.e.,
$\kappa=-\left(K^{xx}+K^{yy}+K^{zz}\right)$.

####  A.1.2 Calculation of normalization parameter $d_{th}$

Once $\mathcal{T}$ is known, the last step to obtain $H(\tilde{\mathbf{x}})$
is to compute the normalization parameter $d_{th}$. To this purpose, we impose
the volume conservation at discrete level by setting that the volume integral
of $\hat{H}$ over local grid cell is equal to the VoF field in that cell.
Using the normalized Cartesian coordinate, this results in:

$\mathcal{I}:\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\hat{H}d\tilde{x}d\tilde{y}d\tilde{z}=\phi\mathrm{.}$
(32)

As remarked in [1], an exact integration of (32) is not possible while it
exists for the one dimensional integration. Taking for example the $\tilde{x}$
direction for an exact integration, we get:

$\int_{0}^{1}\hat{H}(\tilde{x})d\tilde{x}=\dfrac{1}{2}\log\left[\tilde{x}+\dfrac{\log\left(\cosh\left(\beta_{th}\left(\mathcal{T}\left(\tilde{x},\tilde{y},\tilde{z}\right)\right)\right)\right)}{\beta_{th}\left(\partial\mathcal{T}_{\tilde{x}}/\partial\tilde{x}\right)}\right]_{0}^{1}\mathrm{.}$
(33)

Note that in the integration of (33), we assume that the derivative of
$\mathcal{T}$ with respect to $\tilde{x}$ just depends on $\tilde{y}$ and
$\tilde{z}$. This is achieved by setting $\mathbf{c}=[0,1,1]$ and
$\mathbf{c}^{R}=[1,1,0]$. In the other two directions, numerical integration
should be performed and in this work we employ a two-point Gaussian quadrature
method. Therefore, the integral $\mathcal{I}$ in equation (32) results in:

$8\phi=\sum_{p=1}^{4}\left[\left(\tilde{x}+\dfrac{\log\left(\cosh\left(\beta_{th}\left(\mathcal{T}\left(\tilde{x},r_{p}(p),r_{m}(p)\right)+d_{th}\right)\right)\right)}{a_{L}^{x}}\right)\right]_{0}^{1}\mathbf{,}$
(34)

where $r_{p}=1+\sqrt{3}/2[-1,+1,-1,+1]$ and $r_{m}=1+\sqrt{3}/2[-1,-1,+1,+1]$.
In general, the direction along which the exact integration is performed
cannot be decided a-priori. On the other hand, a criterion for this choice is
based on the magnitude of the normal vector components and, therefore, three
cases are possible.

#### Case 1

If, $|n^{x}|\geq(|n^{y}|,|n^{z}|)$, the exact integration is performed only
along $x$ and we set $\mathbf{c}=[0,1,1]$ and $\mathbf{c}^{R}=[1,1,0]$ in
equation (25). Therefore $\mathcal{I}$ becomes:

$8\phi=\sum_{p=1}^{4}\left[\left(\tilde{x}+\dfrac{\log\left(\cosh\left(\beta_{th}\left(\mathcal{T}\left(\tilde{x},r_{p}(p),r_{m}(p)\right)+d_{th}\right)\right)\right)}{a_{L}^{x}}\right)\right]_{0}^{1}\mathbf{.}$
(35)

#### Case 2

On the other hand, if $|n^{y}|\geq(|n^{x}|,|n^{z}|)$, the exact integration is
performed only along $y$ and we set $\mathbf{c}=[1,0,1]$ and
$\mathbf{c}^{R}=[0,1,1]$ in equation (25). Therefore, $\mathcal{I}$ becomes:

$8\phi=\sum_{p=1}^{4}\left[\left(\tilde{y}+\dfrac{\log\left(\cosh\left(\beta_{th}\left(\mathcal{T}\left(r_{p}(p),\tilde{y},r_{m}(p)\right)+d_{th}\right)\right)\right)}{a_{L}^{y}}\right)\right]_{0}^{1}\mathbf{.}$
(36)

#### Case 3

Finally, if $|n^{z}|\geq(|n^{x}|,|n^{y}|)$, the exact integration is performed
only along $z$ and we set $\mathbf{c}=[0,0,1]$ and $\mathbf{c}^{R}=[1,0,1]$ in
equation (25). Therefore, $\mathcal{I}$ results:

$8\phi=\sum_{p=1}^{4}\left[\left(\tilde{z}+\dfrac{\log\left(\cosh\left(\beta_{th}\left(\mathcal{T}\left(r_{p}(p),r_{m}(p),\tilde{z}\right)+d_{th}\right)\right)\right)}{a_{L}^{z}}\right)\right]_{0}^{1}\mathbf{.}$
(37)

Depending on the three different cases, equations (35) or (36) and or (37) can
be solved for the only unknown $d_{th}$. To this purposed, they can be re-
written first as:

$4+\dfrac{1}{\beta_{th}a_{L}^{p}}\log\left(\dfrac{(AB^{-}D+1)(AB^{+}D+1)(AC^{-}D+1)(AC^{+}D+1)}{A^{2}(B^{-}D+1)(B^{+}D+1)(C^{-}D+1)(C^{+}D+1)}\right)=8\phi\mathrm{,}$
(38)

where $a_{L}^{p}$ is $a_{L}^{x}$, $a_{L}^{y}$ or $a_{L}^{z}$ according to the
three cases above. Finally, equation (38) can be further recast in a more
convenient quartic equation:

$\displaystyle 0$
$\displaystyle=\underbrace{A^{2}B^{-}B^{+}C^{-}C^{+}(A^{2}-Q)}_{\alpha_{4}}D^{4}\mathrm{,}$
$\displaystyle+\,\,\underbrace{A^{2}(B^{-}B^{+}(C^{-}+C^{+})+(B^{-}+B^{+})C^{-}C^{+})(A-Q)}_{\alpha_{3}}D^{3}\mathrm{,}$
$\displaystyle+\,\,\underbrace{A^{2}((B^{-}B^{+})(C^{-}+C^{+})+B^{-}B^{+}+C^{-}C^{+})(1-Q)}_{\alpha_{2}}D^{2}\mathrm{,}$
$\displaystyle+\,\,\underbrace{A(B^{-}+B^{+}+C^{-}+C^{+})(1-AQ)}_{\alpha_{1}}D\mathrm{,}$
$\displaystyle+\,\,\underbrace{(1-A^{2}Q)}_{\alpha_{0}}\mathrm{.}$ (39)

where $D=\exp(2\beta_{th}d_{th})$ while the constants $A$, $B^{\pm}$,
$C^{\pm}$ and $Q$ are given by:

$\begin{cases}A&=\exp\left(2\beta_{th}a_{L}^{p}\right)\mathrm{,}\\\
B^{\pm}&=\exp(2\beta_{th}\mathbf{a}_{v}\cdot\mathbf{r}_{p,B})\mathrm{,}\\\
C^{\pm}&=\exp(2\beta_{th}\mathbf{a}_{v}\cdot\mathbf{r}_{p,C})\mathrm{,}\\\
Q&=\exp(4\beta_{th}a_{L}^{p}(2\phi-1))\mathrm{.}\end{cases}$ (40)

The coefficients $a_{L}^{p}$ and the expressions for $\mathbf{a}_{v}$,
$\mathbf{r}_{p,B}$ and $\mathbf{r}_{p,C}$ are reported in tables 5 and 6,
differentiated among the three cases previously described. Once the
coefficients $\alpha_{i=1,4}$ of the quartic equation (39) are computed, a
solution for $D$ can be found. In the current work, we adopt the approach
proposed in [1]. Instead of computing all the four complex roots, we look for
the real and positive solution, fulfilling the constraint given by equation
32.

Case | $a_{L}^{p}$ | $\mathbf{a}_{v}$
---|---|---
1 | $a_{L}^{x}$ | $[a_{Q,2}^{z},a_{Q,2}^{z},a_{Q,1}^{yz},a_{L}^{y},a_{L}^{z}]^{T}$
2 | $a_{L}^{y}$ | $[a_{Q,2}^{x},a_{Q,2}^{z},a_{Q,1}^{xz},a_{L}^{x},a_{L}^{z}]^{T}$
3 | $a_{L}^{z}$ | $[a_{Q,2}^{x},a_{Q,2}^{y},a_{Q,1}^{xy},a_{L}^{x},a_{L}^{y}]^{T}$
Table 5: Coefficients $a_{L}^{p}$ and $\mathbf{a}_{v}$ (as given by equations (27), (28) and (29)) needed to compute $A$, $B^{\pm}$, $C^{\pm}$ and $D$ depending on the three cases. Case | $\mathbf{r}_{p,B}$ | $\mathbf{r}_{p,C}$
---|---|---
1 | $[(r_{p}^{\pm})^{2},(r_{p}^{-})^{2},(r_{p}^{\pm}r_{p}^{-}),r_{p}^{\pm},r_{p}^{-}]^{T}$ | $[(r_{p}^{-})^{2},(r_{p}^{\pm})^{2},(r_{p}^{-}r_{p}^{\pm}),r_{p}^{-},r_{p}^{\pm}]^{T}$
2 | $[(r_{p}^{\pm})^{2},(r_{p}^{-})^{2},(r_{p}^{\pm}r_{p}^{-}),r_{p}^{\pm},r_{p}^{-}]^{T}$ | $[(r_{p}^{-})^{2},(r_{p}^{\pm})^{2},(r_{p}^{-}r_{p}^{\pm}),r_{p}^{-},r_{p}^{\pm}]^{T}$
3 | $[(r_{p}^{\pm})^{2},(r_{p}^{-})^{2},(r_{p}^{\pm}r_{p}^{-}),r_{p}^{\pm},r_{p}^{-}]^{T}$ | $[(r_{p}^{-})^{2},(r_{p}^{\pm})^{2},(r_{p}^{-}r_{p}^{\pm}),r_{p}^{-},r_{p}^{\pm}]^{T}$
Table 6: Coefficients $\mathbf{r}_{p,B}$ and $\mathbf{r}_{p,C}$ needed to
compute $B^{\pm}$ and $C^{\pm}$ depending on the three cases. Note that
$r_{p}^{+}=1+\sqrt{3}/2$ and $r_{p}^{-}=1-\sqrt{3}/2$.

To this purpose, we first write equation (39) as:

$D^{4}+\gamma_{3}D^{3}+\gamma_{2}D+\gamma_{1}D+\gamma_{0}=0\mathrm{,}$ (41)

with $\gamma_{3}=a_{3}/a_{4}$, $\gamma_{2}=a_{2}/a_{4}$,
$\gamma_{1}=a_{1}/a_{4}$ and $\gamma_{0}=a_{0}/a_{4}$. Next, equation (41) is
recast in a quadratic form as:

$(D^{2}+\varepsilon_{1}x+\varepsilon_{2})^{2}-(\varepsilon_{3}D+\varepsilon_{4})^{2}=0\mathrm{,}$
(42)

where $2\gamma_{3}=2\varepsilon_{1}$,
$\gamma_{2}=\varepsilon_{1}^{2}+2\varepsilon_{2}-\varepsilon_{3}^{2}$,
$\gamma_{1}=2\varepsilon_{1}\varepsilon_{2}-2\varepsilon_{3}\varepsilon_{4}$
and $\gamma_{0}=\varepsilon^{2}-\varepsilon_{4}^{2}=\gamma_{0}$. Finally, by
introducing the variable $z=2\varepsilon_{2}$ and by comparing equation (41)
with equation (42), a cubic equation can be derived:

$z^{3}+\eta_{2}z^{2}+\eta_{1}z+\eta_{0}=0\mathrm{,}$ (43)

with $\eta_{2}=-\gamma_{2}$, $\eta_{1}=\gamma_{1}\gamma_{3}-4\gamma_{0}$ and
$\eta_{0}=\gamma_{0}(4\eta_{2}-\eta_{3}^{2})-\eta_{1}^{2}$. Equation (43) can
be easily solved with the Cardano’s formula. Excluding the complex solutions,
setting $\lambda=-\eta_{2}^{2}/9+\eta_{1}/3$,
$\mu=2\eta_{2}^{3}/27-\eta_{1}\eta_{2}/3+\eta_{0}$ and
$\Delta=\lambda^{2}+4\mu^{3}$, the real one $z_{r}$ is given by:

$z_{r}=\begin{cases}\left(\dfrac{-\lambda+\sqrt{\Delta}}{2}\right)^{1/3}-\left(\dfrac{+\lambda+\sqrt{\Delta}}{2}\right)^{1/3}-\dfrac{\eta_{2}}{3}\hskip
28.45274pt\text{if $\Delta\geq 0$}\mathrm{,}\\\
2\sqrt{-\mu}\cos\left[\dfrac{1}{3}\tan^{-1}\left(\dfrac{-\Delta}{-\lambda}\right)\right]-\dfrac{\eta_{2}}{3}\hskip
69.13998pt\text{if $\Delta<0$}\mathrm{.}\\\ \end{cases}$ (44)

Once $z_{r}$ is known, the coefficients $\varepsilon_{i=1,4}$ of equation (42)
can be found: $\varepsilon_{1}=\gamma_{1}/2$, $\varepsilon_{2}=z_{r}/2$,
$\varepsilon_{4}=\sqrt{\varepsilon_{2}^{2}-\gamma_{0}}$ and
$\varepsilon_{3}=(-\gamma_{3}/2+\varepsilon_{1}\varepsilon_{2})/\varepsilon_{4}$.
The last step requires the solution of equation (42). Once more, four
solutions are possible, but since $D>0$ by numerical constrains, only one
positive and real solution is acceptable. This can be easily computed as:

$D=\dfrac{-(\varepsilon_{1}-\varepsilon_{2})+\sqrt{(\varepsilon_{1}-\varepsilon_{2})^{2}-4(\varepsilon_{2}-\varepsilon_{4})}}{2}\mathrm{,}$
(45)

from which, the normalization parameter $d_{th}$ can be evaluated simply as:

$d_{th}=\dfrac{1}{2\beta_{th}}\log(D)\mathrm{.}$ (46)

###  A.2 Computation of the numerical flux

The approximate expression of $H$ is also used for the calculation of the
numerical fluxes in the interface advection step. These are evaluated as:

$f_{i\pm 1/2,j,k}^{x}=\dfrac{1}{\Delta y\Delta
z}\int_{t^{n}}^{t^{n+1}}\left[\int_{\Delta y}\int_{\Delta
z}\left.uH(\mathbf{x},t)\right|_{i\pm 1/2,j,k}^{n}\,\ dy\,\
dz\right]dt\mathrm{,}$ (47) $f_{i,j\pm 1/2,k}^{y}=\dfrac{1}{\Delta x\Delta
z}\int_{t^{n}}^{t^{n+1}}\left[\int_{\Delta x}\int_{\Delta
z}\left.vH(\mathbf{x},t)\right|_{i,j\pm 1/2,k}^{n}\,\ dx\,\
dz\right]dt\mathrm{,}$ (48) $f_{i,j,k\pm 1/2}^{z}=\dfrac{1}{\Delta x\Delta
y}\int_{t^{n}}^{t^{n+1}}\left[\int_{\Delta x}\int_{\Delta
y}\left.wH(\mathbf{x},t)\right|_{i,j,k\pm 1/2}^{n}\,\ dx\,\
dy\right]dt\mathrm{.}$ (49)

Note that the use of equations (47), (48) and (49) is impractical since they
contain both temporal and spatial integrations. For a more simple evaluation,
the time integral is replaced by a space integral performing a change of
variable. Taking as an example the integral along $x$ and using the cell-
centered coordinate system, we define

$\Delta\tilde{x}_{i\pm 1/2,j,k}=\begin{cases}\left[1-\dfrac{\Delta
t^{n+1}}{\Delta x}u_{i+1/2,j,k};1\right]&\text{$u_{i+1/2,j,k}\geq
0$}\mathrm{,}\\\ \left[0;-\dfrac{\Delta t^{n+1}}{\Delta
x}u_{i+1/2,j,k}\right]&\text{$u_{i+1/2,j,k}<0$}\mathrm{.}\\\ \end{cases}$ (50)

and we set $\Delta\tilde{y}=\Delta\tilde{z}=[0,1]$. Moreover, the indicator
function in equation (47) is approximated with the VOF function at the current
time step, $\phi^{n}$. Accordingly, $f_{i,j\pm 1/2,k}^{x}$ can be computed as:

$f_{i\pm
1/2,j,k}^{x}=\begin{cases}\displaystyle{+\Delta\tilde{x}\int_{\Delta\tilde{x}_{i+1/2}}\int_{\Delta\tilde{y}}\int_{\Delta\tilde{z}}\hat{H}_{i,j,k}^{x,n}(\phi^{n})d\tilde{V}\hskip
14.22636pt\text{$u_{i+1/2,j,k}\geq 0$}}\mathrm{,}\\\
\displaystyle{-\Delta\tilde{x}\int_{\Delta\tilde{x}_{i-1/2}}\int_{\Delta\tilde{y}}\int_{\Delta\tilde{z}}\hat{H}_{i,j,k}^{x,n}(\phi^{n})d\tilde{V}\hskip
14.22636pt\text{$u_{i-1/2,j,k}<0$}}\mathrm{.}\\\ \end{cases}$ (51)

where $d\tilde{V}=d\tilde{x}d\tilde{y}d\tilde{z}$. Likewise, $f_{i,j\pm
1/2,k}^{y,n}$ and $f_{i,j,k\pm 1/2}^{z,n}$ can be computed as:

$f_{i,j\pm
1/2,k}^{y}=\begin{cases}\displaystyle{+\Delta\tilde{y}\int_{\Delta\tilde{y}_{i+1/2}}\int_{\Delta\tilde{y}}\int_{\Delta\tilde{z}}\hat{H}_{i,j,k}^{y,n}(\phi^{x})d\tilde{V}\hskip
14.22636pt\text{$v_{i,j+1/2,k}\geq 0$}}\mathrm{,}\\\
\displaystyle{-\Delta\tilde{y}\int_{\Delta\tilde{y}_{i-1/2}}\int_{\Delta\tilde{y}}\int_{\Delta\tilde{z}}\hat{H}_{i,j,k}^{y,n}(\phi^{x})d\tilde{V}\hskip
14.22636pt\text{$v_{i,j-1/2,k}<0$}}\mathrm{.}\\\ \end{cases}$ (52)
$f_{i,j,k\pm
1/2}^{z}=\begin{cases}\displaystyle{+\Delta\tilde{z}\int_{\Delta\tilde{z}_{i+1/2}}\int_{\Delta\tilde{x}}\int_{\Delta\tilde{z}}\hat{H}_{i,j,k}^{z,n}(\phi^{y})d\tilde{V}\hskip
14.22636pt\text{$w_{i,j,k+1/2}\geq 0$}}\mathrm{,}\\\
\displaystyle{-\Delta\tilde{z}\int_{\Delta\tilde{z}_{i-1/2}}\int_{\Delta\tilde{y}}\int_{\Delta\tilde{z}}\hat{H}_{i,j,k}^{z,n}(\phi^{y})d\tilde{V}\hskip
14.22636pt\text{$w_{i,j,k-1/2}<0$}}\mathrm{.}\\\ \end{cases}$ (53)

with:

$\Delta\tilde{y}_{i,j\pm 1/2,k}=\begin{cases}\left[1-\dfrac{\Delta
t^{n+1}}{\Delta y}v_{i,j+1/2,k};1\right]&\text{$v_{i,j+1/2,k}\geq
0$}\mathrm{,}\\\ \left[0;-\dfrac{\Delta t^{n+1}}{\Delta
y}v_{i,j+1/2,k}\right]&\text{$v_{i,j+1/2,k}<0$}\mathrm{.}\\\ \end{cases}$ (54)
$\Delta\tilde{z}_{i,j,k\pm 1/2}=\begin{cases}\left[1-\dfrac{\Delta
t^{n+1}}{\Delta z}w_{i,j,k+1/2};1\right]&\text{$w_{i,j,k+1/2}\geq
0$}\mathrm{,}\\\ \left[0;-\dfrac{\Delta t^{n+1}}{\Delta
z}w_{i,j,k+1/2}\right]&\text{$w_{i,j,k+1/2}<0$}\mathrm{.}\\\ \end{cases}$ (55)

####  A.2.1 Overall VoF algorithm

Below, we report the overall VoF algorithm in the pseudocode 2. As a final
remark, note that the entire algorithm has been described assuming that the
first directional split is always oriented along $x$ (i.e., $x\rightarrow
y\rightarrow z$. Nevertheless, this solution proves to be only first-order
accurate in time. To improve the time accuracy of the solution, one
possibility is to alternate the splitting direction as suggested in [76].

Algorithm 2 Overall VoF algorithm

1:Set $(u^{x},u^{y},u^{z})=(u,v,w)$;

2:for $p=x,y,z$ do

3: Compute the numerical fluxes $f^{p}$;

4: Compute $\phi^{p}$;

5: Set the boundary conditions on $\phi^{p}$;

6: if $p!=z$ then

7: Using $\phi^{p}$, update $\mathbf{n}$ and $\kappa$ with the procedure
described in A.1.1 and $d_{th}$ with the procedure in A.1.2;

8: end if

9:end for

10:Compute $\phi^{n+1}$ using $\phi^{x}$, $\phi^{y}$ and $\phi^{z}$.

11:Using $\phi^{n+1}$, update $\mathbf{n}$ and $\kappa$ with the procedure
described in A.1.1 and $d_{th}$ with the procedure in A.1.2;

## References

  * [1] S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi, Y. Matsumoto, F. Xiao, An interface capturing method with a continuous function: the thinc method with multi-dimensional reconstruction, Journal of Computational Physics 231 (5) (2012) 2328–2358.
  * [2] W. W. Grabowski, L.-P. Wang, Growth of cloud droplets in a turbulent environment, Annual review of fluid mechanics 45 (2013) 293–324.
  * [3] G. Seminara, Fluvial sedimentary patterns, Annual Review of Fluid Mechanics 42 (2010) 43–66.
  * [4] L. Brandt, F. Coletti, Particle-laden turbulence: Progress and perspectives, Annual Review of Fluid Mechanics 54 (2021).
  * [5] F. Veron, Ocean spray, Annual Review of Fluid Mechanics 47 (2015) 507–538.
  * [6] T. Dauxois, T. Peacock, P. Bauer, C.-c. P. Caulfield, C. Cenedese, C. Gorlé, G. Haller, G. N. Ivey, P. F. Linden, E. Meiburg, et al., Confronting grand challenges in environmental fluid mechanics, Physical review fluids 6 (2) (2021) 020501.
  * [7] C. T. Crowe, Multiphase flow handbook, CRC press, 2005.
  * [8] G. A. Voth, A. Soldati, Anisotropic particles in turbulence, Annual Review of Fluid Mechanics 49 (2017) 249–276.
  * [9] F. Risso, Agitation, mixing, and transfers induced by bubbles, Annual Review of Fluid Mechanics 50 (2018) 25–48.
  * [10] S. Elghobashi, Direct numerical simulation of turbulent flows laden with droplets or bubbles, Annual Review of Fluid Mechanics 51 (2019) 217–244.
  * [11] V. Mathai, D. Lohse, C. Sun, Bubbly and buoyant particle–laden turbulent flows, Annual Review of Condensed Matter Physics 11 (2020) 529–559.
  * [12] A. U. M. Masuk, A. Salibindla, R. Ni, A robust virtual-camera 3d shape reconstruction of deforming bubbles/droplets with additional physical constraints, International Journal of Multiphase Flow 120 (2019) 103088.
  * [13] A. K. Salibindla, A. U. M. Masuk, S. Tan, R. Ni, Lift and drag coefficients of deformable bubbles in intense turbulence determined from bubble rise velocity, Journal of Fluid Mechanics 894 (2020).
  * [14] A. U. M. Masuk, A. K. Salibindla, R. Ni, Simultaneous measurements of deforming hinze-scale bubbles with surrounding turbulence, Journal of Fluid Mechanics 910 (2021).
  * [15] S. Mirjalili, S. S. Jain, M. Dodd, Interface-capturing methods for two-phase flows: An overview and recent developments, Center for Turbulence Research Annual Research Briefs 2017 (117-135) (2017) 13.
  * [16] S. O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, Journal of computational physics 100 (1) (1992) 25–37.
  * [17] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annual review of fluid mechanics 31 (1) (1999) 567–603.
  * [18] D. M. Anderson, G. B. McFadden, A. A. Wheeler, Diffuse-interface methods in fluid mechanics, Annual review of fluid mechanics 30 (1) (1998) 139–165.
  * [19] J. A. Sethian, P. Smereka, Level set methods for fluid interfaces, Annual review of fluid mechanics 35 (1) (2003) 341–372.
  * [20] A. Prosperetti, G. Tryggvason, Computational methods for multiphase flow, Cambridge university press, 2009.
  * [21] G. Soligo, A. Roccon, A. Soldati, Turbulent flows with drops and bubbles: What numerical simulations can tell us—freeman scholar lecture, Journal of Fluids Engineering 143 (8) (2021).
  * [22] A. Khan, H. Sim, S. S. Vazhkudai, A. R. Butt, Y. Kim, An analysis of system balance and architectural trends based on top500 supercomputers, in: The International Conference on High Performance Computing in Asia-Pacific Region, 2021, pp. 11–22.
  * [23] X. Zhu, E. Phillips, V. Spandan, J. Donners, G. Ruetsch, J. Romero, R. Ostilla-Mónico, Y. Yang, D. Lohse, R. Verzicco, et al., Afid-gpu: a versatile navier–stokes solver for wall-bounded turbulent flows on gpu clusters, Computer physics communications 229 (2018) 199–210.
  * [24] M. Bernardini, D. Modesti, F. Salvadore, S. Pirozzoli, Streams: A high-fidelity accelerated solver for direct numerical simulation of compressible turbulent flows, Computer Physics Communications 263 (2021) 107906.
  * [25] P. Costa, E. Phillips, L. Brandt, M. Fatica, Gpu acceleration of cans for massively-parallel direct numerical simulations of canonical fluid flows, Computers & Mathematics with Applications 81 (2021) 502–511.
  * [26] W. Aniszewski, T. Arrufat, M. Crialesi-Esposito, S. Dabiri, D. Fuster, Y. Ling, J. Lu, L. Malan, S. Pal, R. Scardovelli, et al., Parallel, robust, interface simulator (paris), Computer Physics Communications 263 (2021) 107849.
  * [27] P. Cifani, J. Kuerten, B. Geurts, Highly scalable dns solver for turbulent bubble-laden channel flow, Computers & Fluids 172 (2018) 67–83.
  * [28] K. Eisenschmidt, M. Ertl, H. Gomaa, C. Kieffer-Roth, C. Meister, P. Rauschenberger, M. Reitzle, K. Schlottke, B. Weigand, Direct numerical simulations for multiphase flows: An overview of the multiphase code fs3d, Applied Mathematics and Computation 272 (2016) 508–517.
  * [29] O. Desjardins, G. Blanquart, G. Balarac, H. Pitsch, High order conservative finite difference scheme for variable density low mach number turbulent flows, Journal of Computational Physics 227 (15) (2008) 7125–7159.
  * [30] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, Journal of Computational Physics 228 (16) (2009) 5838–5866.
  * [31] S. H. Bryngelson, K. Schmidmayer, V. Coralic, J. C. Meng, K. Maeda, T. Colonius, Mfc: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver, Computer Physics Communications (2020) 107396doi:10.1016/j.cpc.2020.107396.
  * [32] P. Costa, A fft-based finite-difference solver for massively-parallel direct numerical simulations of turbulent flows, Computers & Mathematics with Applications 76 (8) (2018) 1853–1862.
  * [33] U. Schumann, R. A. Sweet, Fast fourier transforms for direct solution of poisson’s equation with staggered boundary conditions, Journal of Computational Physics 75 (1) (1988) 123–137.
  * [34] M. E. Rosti, F. De Vita, L. Brandt, Numerical simulations of emulsions in shear flows, Acta Mechanica 230 (2) (2019) 667–682.
  * [35] F. De Vita, M. E. Rosti, S. Caserta, L. Brandt, On the effect of coalescence on the rheology of emulsions, Journal of Fluid Mechanics 880 (2019) 969–991.
  * [36] F. De Vita, M. E. Rosti, S. Caserta, L. Brandt, Numerical simulations of vorticity banding of emulsions in shear flows, Soft matter 16 (11) (2020) 2854–2863.
  * [37] M. E. Rosti, S. Takagi, Shear-thinning and shear-thickening emulsions in shear flows, Physics of Fluids 33 (8) (2021) 083319.
  * [38] M. E. Rosti, Z. Ge, S. S. Jain, M. S. Dodd, L. Brandt, Droplets in homogeneous shear turbulence, Journal of Fluid Mechanics 876 (2019) 962–984.
  * [39] M. Kozul, P. S. Costa, J. R. Dawson, L. Brandt, Aerodynamically driven rupture of a liquid film by turbulent shear flow, Physical Review Fluids 5 (12) (2020) 124302.
  * [40] M. Crialesi-Esposito, M. E. Rosti, S. Chibbaro, L. Brandt, Modulation of homogeneous and isotropic turbulence in emulsions, Journal of Fluid Mechanics 940 (2022).
  * [41] I. Cannon, D. Izbassarov, O. Tammisola, L. Brandt, M. E. Rosti, The effect of droplet coalescence on drag in turbulent channel flows, Physics of Fluids 33 (8) (2021) 085112.
  * [42] N. Scapin, P. Costa, L. Brandt, A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows, Journal of Computational Physics 407 (2020) 109251.
  * [43] F. Dalla Barba, N. Scapin, A. D. Demou, M. E. Rosti, F. Picano, L. Brandt, An interface capturing method for liquid-gas flows at low-mach number, Computers & Fluids 216 (2021) 104789.
  * [44] N. Scapin, F. Dalla Barba, G. Lupo, M. E. Rosti, C. Duwig, L. Brandt, Finite-size evaporating droplets in weakly compressible homogeneous shear turbulence, Journal of Fluid Mechanics 934 (2022).
  * [45] M. Ishii, T. Hibiki, Thermo-fluid dynamics of two-phase flow, Springer Science & Business Media, 2010.
  * [46] F. H. Harlow, J. E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, The physics of fluids 8 (12) (1965) 2182–2189.
  * [47] E. G. Puckett, A. S. Almgren, J. B. Bell, D. L. Marcus, W. J. Rider, A high-order projection method for tracking fluid interfaces in variable density incompressible flows, Journal of computational physics 130 (2) (1997) 269–282.
  * [48] E. Aulisa, S. Manservisi, R. Scardovelli, S. Zaleski, A geometrical area-preserving volume-of-fluid advection method, Journal of Computational Physics 192 (1) (2003) 355–364.
  * [49] G. D. Weymouth, D. K.-P. Yue, Conservative volume-of-fluid method for free-surface simulations on cartesian-grids, Journal of Computational Physics 229 (8) (2010) 2853–2865.
  * [50] M. Castro, B. Costa, W. S. Don, High order weighted essentially non-oscillatory weno-z schemes for hyperbolic conservation laws, Journal of Computational Physics 230 (5) (2011) 1766–1792.
  * [51] A. J. Chorin, Numerical solution of the navier-stokes equations, Mathematics of computation 22 (104) (1968) 745–762.
  * [52] C. Frantzis, D. G. Grigoriadis, An efficient method for two-fluid incompressible flows appropriate for the immersed boundary method, Journal of Computational Physics 376 (2019) 28–53.
  * [53] S. Dong, J. Shen, A time-stepping scheme involving constant coefficient matrices for phase-field simulations of two-phase incompressible flows with large density ratios, Journal of Computational Physics 231 (17) (2012) 5788–5804.
  * [54] M. S. Dodd, A. Ferrante, A fast pressure-correction method for incompressible two-fluid flows, Journal of Computational Physics 273 (2014) 416–434.
  * [55] P. N. Swarztrauber, The methods of cyclic reduction, fourier analysis and the facr algorithm for the discrete solution of poisson’s equation on a rectangle, Siam Review 19 (3) (1977) 490–501.
  * [56] J. Makhoul, A fast cosine transform in one and two dimensions, IEEE Transactions on Acoustics, Speech, and Signal Processing 28 (1) (1980) 27–34.
  * [57] M. Frigo, S. G. Johnson, Fftw: An adaptive software architecture for the fft, in: Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP’98 (Cat. No. 98CH36181), Vol. 3, IEEE, 1998, pp. 1381–1384.
  * [58] N. Li, S. Laizet, 2decomp & fft-a highly scalable 2d decomposition library and fft interface, in: Cray user group 2010 conference, 2010, pp. 1–13.
  * [59] S. Turek, O. Mierka, K. Bäumler, Numerical benchmarking for 3d multiphase flow: New results for a rising bubble, in: European Conference on Numerical Mathematics and Advanced Applications, Springer, 2017, pp. 593–601.
  * [60] J. Armengol, F. Bannwart, J. Xamán, R. Santos, Effects of variable air properties on transient natural convection for large temperature differences, International Journal of Thermal Sciences 120 (2017) 63–79.
  * [61] G. de Vahl Davis, I. Jones, Natural convection in a square cavity: a comparison exercise, International Journal for numerical methods in fluids 3 (3) (1983) 227–248.
  * [62] M. Leal, H. Machado, R. Cotta, Integral transform solutions of transient natural convection in enclosures with variable fluid properties, International Journal of Heat and Mass Transfer 43 (21) (2000) 3977–3990.
  * [63] Openacc.
URL https://www.openacc.org/sites/default/files/inline-files/OpenACC.2.7.pdf

  * [64] Meluxina.
URL https://luxprovide.lu/technical-structure/

  * [65] Berzelius.
URL https://www.nsc.liu.se/systems/berzelius/

  * [66] On the non-linear stability of the 1:1:1 ABC flow, Physica D: Nonlinear Phenomena 75 (4) (1994) 471–508.  doi:10.1016/0167-2789(94)00031-X.
  * [67] P. D. Mininni, A. Alexakis, A. Pouquet, Large-scale flow effects, energy transfer, and self-similarity on turbulence, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 74 (1) (2006) 1–13.  doi:10.1103/PhysRevE.74.016303.
  * [68] C. Garrett, M. Li, D. Farmer, The connection between bubble size spectra and energy dissipation rates in the upper ocean, Journal of Physical Oceanography 30 (9) (2000) 2163–2171.  doi:10.1175/1520-0485(2000)030<2163:TCBBSS>2.0.CO;2.
  * [69] G. B. Deane, M. D. Stokes, Scale dependence of bubble creation mechanisms in breaking waves, Nature 418 (6900) (2002) 839–844.  doi:10.1038/nature00967.
  * [70] F. H. Busse, On the aspect ratios of two-layer mantle convection, Physics of the Earth and Planetary Interiors 24 (4) (1981) 320–324.
  * [71] F. Wilczynski, D. W. Hughes, Stability of two-layer miscible convection, Physical Review Fluids 4 (10) (2019) 103502.
  * [72] H.-R. Liu, K. L. Chong, Q. Wang, C. S. Ng, R. Verzicco, D. Lohse, Two-layer thermally driven turbulence: mechanisms for interface breakup, Journal of Fluid Mechanics 913 (2021).
  * [73] S. Ha, J. Park, D. You, A multi-gpu method for adi-based fractional-step integration of incompressible navier-stokes equations, Computer Physics Communications 265 (2021) 107999.
  * [74] D. L. Youngs, Time-dependent multi-material flow with large fluid distortion, Numerical methods for fluid dynamics (1982).
  * [75] D. L. Youngs, An interface tracking method for a 3d eulerian hydrodynamics code, Atomic Weapons Research Establishment (AWRE) Technical Report 44 (92) (1984) 35.
  * [76] G. Strang, On the construction and comparison of difference schemes, SIAM journal on numerical analysis 5 (3) (1968) 506–517.