The growing availability of GPU-accelerated open-source solvers has boosted the capability of tackling complex single-phase and multiphase turbulent flows by means of direct and large-eddy simulations. GPU-accelerated solvers can leverage the heterogeneous computing architectures that are available in leading high-performance computing centers worldwide, taking advantage of the higher throughput and greater energy efficiency offered by GPUs as compared to CPUs. However, porting CPU-based numerical solvers to GPUs entails many outstanding challenges, such as parallelism exposure, inter-GPU communication, memory allocation constraints, and shared memory limitations. To overcome these challenges, GPU-friendly algorithms, performance portability strategies, and careful selection of computational paradigms and programming languages must be developed. Besides, adaptive mesh refinement and data compression may be integrated to mitigate I-O bottlenecks and enable simulations of more complex geometries on top of the existing requirements imposed by incompressible flows. When compressibility effects become significant, further considerations related to the adoption of high-performance preconditioners and multigrid solvers become crucial for tackling large, sparse linear systems and extending simulations to high-Mach flows. Finally, reduced-precision arithmetic can further enhance performance, energy efficiency, and scalability. In this work, we survey current applications of GPU-accelerated solvers in the broad area of fluid mechanics and turbulence simulations and discuss the main challenges and bottlenecks associated with code porting and optimization. We then conclude our analysis with an outlook on future perspectives for enabling efficient GPU-based exascale computing of turbulence.

In this work we propose a novel method to ensure important entropy inequalities are satisfied semi-discretely when constructing reduced order models (ROMs) on nonlinear reduced manifolds. We are in particular interested in ROMs of systems of nonlinear hyperbolic conservation laws. The so-called entropy stability property endows the semi-discrete ROMs with physically admissible behaviour. The method generalizes earlier results on entropy-stable ROMs constructed on linear spaces. The ROM works by evaluating the projected system on a well-chosen approximation of the state that ensures entropy stability. To ensure accuracy of the ROM after this approximation we locally enrich the tangent space of the reduced manifold with important quantities. Using numerical experiments on some well-known equations (the inviscid Burgers equation, shallow water equations and compressible Euler equations) we show the improved structure-preserving properties of our ROM compared to standard approaches and that our approximations have minimal impact on the accuracy of the ROM. We additionally generalize the recently proposed polynomial reduced manifolds to rational polynomial manifolds and show that this leads to an increase in accuracy for our experiments.

We introduce CaLES, a GPU-accelerated finite-difference solver designed for large-eddy simulations (LES) of incompressible wall-bounded flows in massively parallel environments. Built upon the existing direct numerical simulation (DNS) solver CaNS, CaLES relies on low-storage, third-order Runge-Kutta schemes for temporal discretization, with the option to treat viscous terms via an implicit Crank-Nicolson scheme in one or three directions. A fast direct solver, based on eigenfunction expansions, is used to solve the discretized Poisson/Helmholtz equations. For turbulence modeling, the classical Smagorinsky model with van Driest near-wall damping and the dynamic Smagorinsky model are implemented, along with a logarithmic law wall model. GPU acceleration is achieved through OpenACC directives, following CaNS-2.3.0. Performance assessments were conducted on the Leonardo cluster at CINECA, Italy. Each node is equipped with one Intel Xeon Platinum 8358 CPU (2.60 GHz, 32 cores) and four NVIDIA A100 GPUs (64 GB HBM2e), interconnected via NVLink 3.0 (200 GB/s). The inter-node communication bandwidth is 25 GB/s, supported by a DragonFly+ network architecture with NVIDIA Mellanox InfiniBand HDR. Results indicate that the computational speed on a single GPU is equivalent to approximately 15 CPU nodes, depending on the treatment of viscous terms and the subgrid-scale model, and that the solver efficiently scales across multiple GPUs. The predictive capability of CaLES has been tested using multiple flow cases, including decaying isotropic turbulence, turbulent channel flow, and turbulent duct flow. The high computational efficiency of the solver enables grid convergence studies on extremely fine grids, pinpointing non-monotonic grid convergence for wall-modeled LES.

In the last decades the urban population has increased a lot making cities objects of different studies. In this context, urban climate is an expanding research field and understanding the main features of the urban heat island (UHI) effect is one of the challenges. A discrete number of neighbourhoods has been object of study for this climate effect with growing interest in recent years, especially focusing on heat mitigation. Despite this, there is a lack of knowledge due to the complex nature of the problem given by the multi-physics involved, the multiple parameters that govern it and above all, the complexity of the city’s geometries that lack generality. In this respect, and to keep results applicable in a broader context, this work proposes an innovative approach to studying UHI effects, providing a unique framework for understanding the interaction between urban geometry and heat transport dynamics while addressing the complexities of urban configurations with a novel and methodological perspective.

We present a massively parallel GPU-accelerated solver for direct numerical simulations of transitional and turbulent flat-plate boundary layers and channel flows involving fluids in non-ideal thermodynamic states. While several high-fidelity solvers are currently available as open source, all of them are restricted to the ideal-gas region. In contrast, the CUBic Equation of state Navier-Stokes solver (CUBENS) can accurately model and simulate the non-ideal thermodynamics of single-phase compressible fluids in the vicinity of the vapor-liquid saturation line or the thermodynamic critical point. By employing high-order finite-difference schemes and convective terms in split, kinetic-energy-, and entropy-preserving form, the solver is numerically stable, and robust with minimal numerical dissipation, enabling it to capture the steep variations of non-ideal thermodynamic properties. For cost-effective high-fidelity simulations, in addition to MPI parallelization, CUBENS is GPU-accelerated using OpenACC directives for computation offloading, and asynchronous GPU-aware MPI for efficient GPU-GPU communication. Moreover, CUBENS is compatible with both NVIDIA and AMD GPU architectures, achieving significant performance results while ensuring energy-efficient simulations. For instance, using 64 NVIDIA A100 GPUs compared to 8192 CPUs at the same computational cost results in a speedup of approximately 130×. In multi-node and multi-GPU configurations ranging from 2 to 128 compute nodes (8 to 512 GPUs), a strong scaling efficiency of around 52% and a weak scaling efficiency of 0.88 with 1024³ points per GPU, corresponding to approximately 5 billion degrees of freedom, are achieved. The CUBENS solver is validated against selected cases from the literature, covering transitional to turbulent ideal and non-ideal flows up to the transonic regime. In particular, we demonstrate the solver's suitability and applicability for direct numerical simulations of transitional boundary layers with fluids at supercritical pressure and with buoyancy effects. The development of this high-fidelity solver offers the potential for future fundamental research in non-ideal compressible fluid dynamics. Program summary: Program Title: CUBic Equation of state Navier-Stokes (CUBENS) CPC Library link to program files: https://doi.org/10.17632/6jfy758gyv.1 Developer's repository link: https://github.com/pcboldini/CUBENS Licensing provisions: MIT Programming language: Fortran 90, OpenACC, MPI, Python, MATLAB Nature of problem: This code solves the three-dimensional Navier-Stokes equations for non-ideal gas flows in a Cartesian domain, applicable to boundary layers and channels. Solution method: This code uses high-order central finite-differences with split-convective form, preserving kinetic energy and entropy (KEEP) and pressure-equilibrium-preserving (PEP) property, for spatial discretization. The time advancement is performed with a third-order Total Variation Diminishing low-storage Runge-Kutta scheme. Flow non-ideality is accounted for by cubic equations of state and complex transport-properties models. Alongside MPI parallelization, the solver is GPU-accelerated using OpenACC for computation offloading and CPU-GPU data transfer, along with GPU-aware MPI for GPU-GPU communication.

We present a computational method for extreme-scale simulations of incompressible turbulent wall flows at high Reynolds numbers. The numerical algorithm extends a popular method for solving second-order finite differences Poisson/Helmholtz equations using a pencil-distributed parallel tridiagonal solver to improve computational performance at scale. The benefits of this approach were investigated for high-Reynolds-number turbulent channel flow simulations, with up to about 80 billion grid points and 1024 GPUs on the European flagship supercomputers Leonardo and LUMI. An additional GPU porting effort of the entire solver had to be undertaken for the latter. Our results confirm that, while 1D domain decompositions are favorable for smaller systems, they become inefficient or even impossible at large scales. This restriction is relaxed by adopting a pencil-distributed approach. The results show that, at scale, the revised Poisson solver is about twice as fast as the baseline approach with the full-transpose algorithm for 2D domain decompositions. Strong and weak scalability tests show that the performance gains are due to the lower communication footprint. Additionally, to secure high performance when solving for wall-normal implicit diffusion, we propose a reworked flavor of parallel cyclic reduction (PCR) that is split into pre-processing and runtime steps. During pre-processing, small sub-arrays with independent 1D coefficients are computed by parallel GPU threads, without any global GPU communication. Then, at runtime, the reworked PCR enables a fast solution of implicit 1D diffusion without computational overhead. Our results show that the entire numerical solver, coupled with the PCR algorithm, enables extreme-scale simulations with 2D pencil decompositions, which do not suffer performance losses even when compared to the best 1D slab configurations available for smaller systems.

The Monin–Obukhov similarity theory (MOST) is a cornerstone of atmospheric science for describing turbulence in stable boundary layers. Extending MOST to stably stratified turbulent channel flows, however, is non-trivial due to confinement by solid walls. In this study, we investigate the applicability of MOST in closed channels and identify where and to what extent the theory remains valid. A key finding is that the ratio of the half-channel height to the Obukhov length serves as a governing parameter for identifying distinct flow regions and determining their corresponding mean velocity scaling. Hence, we propose a relation to estimate this ratio directly from the governing input parameters: the friction Reynolds and friction Richardson numbers (Re_τ and Ri_τ). The framework is tested against a series of direct numerical simulations across a range of Re_τ and Ri_τ. The reconstructed velocity profiles enable accurate prediction of the skin-friction coefficient crucial for quantifying pressure losses in stratified flows in engineering applications.

In particle-laden turbulent wall flows, lift forces can influence the near-wall turbulence. This has been observed recently in particle-resolved simulations, which, however, are too expensive to be used in upscaled models. Instead, point-particle simulations have been the method of choice to simulate the dynamics of these flows during the last decades. While this approach is simpler, cheaper and physically sound for small inertial particles in turbulence, some issues remain. In the present work, we address challenges associated with lift force modelling in turbulent wall flows and the impact of lift forces in the near-wall flow. We performed direct numerical simulations of small inertial point particles in turbulent channel flow for fixed Stokes number and mass loading while varying the particle size. Our results show that the particle dynamics in the buffer region, causing the apparent particle-to-fluid slip velocity to vanish, raises major challenges for modelling lift forces accurately. While our results confirm that lift forces have little influence on particle dynamics for sufficiently small particle sizes, for inner-scaled diameters of order one and beyond, lift forces become quite important near the wall. The different particle dynamics under lift forces results in the modulation of streamwise momentum transport in the near-wall region. We analyse this lift-induced turbulence modulation for different lift force models, and the results indicate that realistic models are critical for particle-modelled simulations to correctly predict turbulence modulation by particles in the near-wall region.

Suspensions of finite-size solid particles in a turbulent pipe flow are found in many industrial and technical flows. Due to the ample parameter space consisting of particle size, concentration, density and Reynolds number, a complete picture of the particle–fluid interaction is still lacking. Pressure drop predictions are often made using viscosity models only considering the bulk solid volume fraction. For the case of turbulent pipe flow laden with neutrally buoyant spherical particles, we investigate the pressure drop and overall drag (friction factor), fluid velocity and particle distribution in the pipe. We use a combination of experimental (MRV) and numerical (DNS) techniques and a continuum flow model. We find that the particle size and the bulk flow rate influence the mean fluid velocity, velocity fluctuations and the particle distribution in the pipe for low flow rates. However, the effects of the added solid particles diminish as the flow rate increases. We created a master curve for drag change compared to single-phase flow for the particle-laden cases. This curve can be used to achieve more accurate friction factor predictions than the traditional modified viscosity approach that does not account for particle size.

We propose a numerical method tailored to perform interface-resolved simulations of evaporating multicomponent two-phase flows. The novelty of the method lies in the use of Robin boundary conditions to couple the transport equations for the vaporized species in the gas phase and the transport equations of the same species in the liquid phase. The Robin boundary condition is implemented with the cost-effective procedure proposed by Chai et al. [1] and consists of two steps: (1) calculating the normal derivative of the mass fraction fields in cells adjacent to the interface through the reconstruction of a linear polynomial system, and (2) extrapolating the normal derivative and the ghost value in the normal direction using a linear partial differential equation. This methodology yields a second-order accurate solution for the Poisson equation with a Robin boundary condition and a first-order accurate solution for the Stefan problem. The overall methodology is implemented in an efficient two-fluid solver, which includes a Volume-of-Fluid (VoF) approach for the interface representation, a divergence-free extension of the liquid velocity field onto the entire domain to transport the VoF, and the temperature equation to include thermal effects. We demonstrate the convergence of the numerical method to the analytical solution for multicomponent isothermal evaporation and observe good overall computational performance for simulating non-isothermal evaporating two-fluid flows in two and three dimensions.

The geometrically complex mechanisms of energy transfer in the compound space of scales and positions of wall turbulent flows are investigated in a temporally evolving boundary layer. The phenomena consist of spatially ascending reverse and forward cascades from the small production scales of the buffer layer to the small dissipative scales distributed among the entire boundary layer height. The observed qualitative behaviour conforms with previous results in turbulent channel flow, thus suggesting that the observed phenomenology is a robust statistical feature of wall turbulence in general. An interesting feature is the behaviour of energy transfer at the turbulent/non-turbulent interface, where forward energy cascade is found to be almost absent. In particular, the turbulent core is found to sustain a variety of large-scale wall-parallel motions at the turbulent interface through weak but persistent reverse energy cascades. This behaviour conforms with previous results in free shear flows, thus suggesting that the observed phenomenology is a robust statistical feature of turbulent shear flows featuring turbulent/non-turbulent interfaces in general.

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. Program summary: Program Title: : Fluid Transport Accelerated Solver, FluTAS. CPC Library link to program files: https://doi.org/10.17632/tp6k8wky8m.1 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.