Source-linked AI summary
AFiD-GPU: a versatile Navier-Stokes Solver for Wall-Bounded Turbulent Flows on GPU Clusters
Xiaojue Zhu, Everett Phillips, Vamsi Spandan, John Donners, Gregory Ruetsch, Josh Romero, Rodolfo Ostilla-Mónico, Yantao Yang, Detlef Lohse, Roberto Verzicco, Massimiliano Fatica, Richard J. A. M. Stevens
TL;DR
Turbulence demands large meshes and immense computing power to resolve its many scales, while extreme Rayleigh-number convection remains computationally inaccessible. This paper ports AFiD to GPUs and improves its efficiency, reporting an order-of-magnitude wall-clock improvement over CPU and simulations of two wall-bounded flows.
Problem
Direct numerical simulations must resolve turbulent scales with large meshes and immense computing power, while the ultimate regime of Rayleigh–Bénard flow has not been reached computationally because of its extreme Rayleigh and Reynolds numbers.
Source: §1 · pp. 2–3Method
The paper ports AFiD to GPU clusters and introduces efficiency improvements for large fluid-dynamics simulations.
Source: §1 · p. 4Results
The GPU code improves wall-clock time by an order of magnitude over the CPU code and is shown working with Rayleigh–Bénard convection and plane Couette flow.
Source: §6 · p. 25Takeaways & Limitations
The authors conclude that the GPU code can work with both tested flows, while higher-Rayleigh and higher-Reynolds-number simulations remain more demanding than those tested.
Source: §6 · p. 25For Rayleigh–Bénard convection, the Boussinesq approximation makes buoyancy depend linearly on temperature and ignores temperature dependence in density, viscosity, and thermal diffusivity.
Source: §2.1 · p. 5
Abstract
from arXiv · showhide
The AFiD code, an open source solver for the incompressible Navier-Stokes equations ({\color{blue}\burl{http://www.afid.eu}}), has been ported to GPU clusters to tackle large-scale wall-bounded turbulent flow simulations. The GPU porting has been carried out in CUDA Fortran with the extensive use of kernel loop directives (CUF kernels) in order to have a source code as close as possible to the original CPU version; just a few routines have been manually rewritten. A new transpose scheme, which is not limited to the GPU version only and can be generally applied to any CFD code that uses pencil distributed parallelization, has been devised to improve the scaling of the Poisson solver, the main bottleneck of incompressible solvers. The GPU version can reduce the wall clock time by an order of magnitude compared to the CPU version for large meshes. Due to the increased performance and efficient use of memory, the GPU version of AFiD can perform simulations in parameter ranges that are unprecedented in thermally-driven wall-bounded turbulence. To verify the accuracy of the code, turbulent Rayleigh-B\'enard convection and plane Couette flow are simulated and the results are in good agreement with the experimental and computational data that published in previous literatures.
PROGRAM SUMMARY
AFiD-GPU is an open-source, GPU-accelerated solver for three-dimensional Navier–Stokes equations coupled to a scalar field. It uses finite differences, Runge–Kutta and Crank–Nicolson time advancement, and pencil-distributed MPI parallelization.
Program overview
- AFiD-GPU solves three-dimensional Navier–Stokes equations coupled with a scalar field in a cubic box bounded by two walls and four periodic boundaries.The program is written in Fortran 90, CUDA Fortran, and MPI, and is licensed under GPLv3.Source: PROGRAM SUMMARY · p. 2
Solution method
- The solver uses second-order finite differences for spatial discretization and third-order Runge–Kutta with Crank–Nicolson time advancement.Its parallelization uses two-dimensional pencil-distributed MPI, with GPU-accelerated routines.Source: PROGRAM SUMMARY · p. 2
Software
- The open-source code is supported and updated through the AFiD website.External routines include PGI, CUDA Toolkit, FFTW3, and HDF5.Source: PROGRAM SUMMARY · p. 2
1. Introduction
The paper targets the computational demands of resolving turbulent flows across their many scales, including high-Rayleigh-number convection and high-Reynolds-number wall-bounded flows. It ports AFiD to GPU clusters and develops efficiency improvements for large simulations.
Computational motivation
- Direct numerical simulations must resolve the scales of turbulent motion, requiring large meshes and immense computational power.As fluid velocity rises, energy transfers to smaller scales and the range of motion scales increases.Source: §1 · p. 2
Example flows
- Rayleigh–Bénard convection and plane Couette flow serve as paradigmatic examples of convection and wall-bounded shear flow.The paper describes these classical systems as settings for studying concepts including nonlinear dynamics, pattern formation, and turbulence.Source: §1 · pp. 2–3
Rayleigh–Bénard convection
- High-Rayleigh-number convection remains computationally challenging because the ultimate regime—with turbulence in both bulk and boundary layers—has not been reached computationally.The paper identifies the extreme Rayleigh and Reynolds numbers of this regime as the computational challenge.Source: §1 · p. 3
GPU motivation and contribution
- The AFiD GPU port and efficiency improvements are intended to support unprecedentedly large fluid-dynamics simulations.The paper notes that GPUs suit data-parallel work and offer high memory bandwidth, useful for low-order finite-difference CFD with little data reuse.Source: §1 · pp. 3–4
2. AFiD code
AFiD solves three-dimensional flow equations on Cartesian meshes and advances them with a staggered-grid finite-difference scheme. Its pencil decomposition supports local wall-normal solves, while pressure correction requires Fourier transforms, data transposes, and a tridiagonal solve.
2. AFiD code
- AFiD solves three-dimensional Navier–Stokes equations with a temperature equation on a Cartesian mesh, using periodic wall-parallel directions and a bounded wall-normal direction.Temperature is active for Rayleigh–Bénard flow and disabled for plane Couette flow; other features are shared except boundary conditions.Source: §2.1 · p. 4
2.1. Numerical scheme
- The solver uses conservative, central, second-order finite differences on a staggered grid, with pressure at cell centers and temperature collocated with vertical velocity for thermal convection.The arrangement avoids interpolation error in the buoyancy term and is energy-conserving as the time step tends to zero.Source: §2.1 · p. 5, Figure 1 · p. 5
- Time advancement combines fractional steps, low-storage third-order Runge–Kutta, and Crank–Nicolson for implicit terms; pressure correction enforces continuity and yields divergence-free velocity.The pressure correction is computed from a Poisson equation; simulations use a maximum CFL number of approximately 1.3.Source: §2.1 · p. 6
2.2. Parallelization strategy
- The two-dimensional pencil decomposition aligns data in the wall-normal direction, enabling local implicit solves there while avoiding all-to-all communication for those calculations.Horizontal viscous terms are integrated explicitly; halo updates are still needed for intermediate velocity calculations.Source: §2.2 · p. 6
- The Poisson solver applies Fourier transforms in both periodic directions, transposes data among decompositions, and solves the resulting wall-normal systems with a tridiagonal Thomas solver.The described workflow computes divergence, performs forward transforms and transposes, solves the tridiagonal systems, then reverses the transforms and transposes to return the correction in real space.Source: §2.2 · p. 7
- The original CPU Poisson algorithm transposes one three-dimensional array and uses all possible transpose combinations, including non-neighbor x-to-z and z-to-x transfers.Those two transpose types use MPI_Neighbor_alltoallw in the CPU implementation; the GPU version uses a different scheme described later.Source: §2.2 · p. 8, Figure 2 · p. 8
3. GPU implementation
AFiD’s GPU implementation uses CUDA Fortran and CUF kernels to keep the code close to its CPU version, while adding memory reuse and a revised transpose scheme. The implementation also instruments profiler traces and manages multi-GPU data transfers.
3. GPU implementation
- AFiD’s GPU port uses CUDA Fortran and CUF kernels extensively, with explicit kernels retained for selected routines.The implementation aims to keep the source close to the CPU version and allows explicit control of data movement.Source: §3 · pp. 8–9
3.1. CUDA & CUDA Fortran
- CUDA Fortran extends regular Fortran, and CUF directives let the compiler generate GPU kernels from nested loops while leaving loop contents usually unchanged.The implementation relies on the PGI compiler because the other compiler described does not support CUF kernels.Source: §3.1 · pp. 9–10
3.2. CUF kernels
- CUF kernels handle nested-loop parallelization and scalar reductions, while manually written routines cover vector statistics, batched tridiagonal solves, and local transposes.The described CUF reductions work only on scalars, so the statistics routine is manually implemented.Source: §3.2 · p. 12
3.3. Reducing the memory footprint
- Memory reuse reduces the GPU count needed for a 1024^3 grid from 48 K20x to 25 K20x.The implementation shares storage for temporary real and complex work arrays and reuses FFT workspaces across plans.Source: §3.3 · pp. 12–13
3.4. Multi GPU implementation
- Each MPI rank maps to one GPU, and asynchronous copies combined with IRECV/ISEND overlap host-device transfers with communication.For this initial version, standard MPI calls require data to be in CPU memory.Source: §3.4 · p. 13
3.5. Efficient data transposes
- The revised transpose adds a z-to-x rotation so processors exchange data within row sub-communicators, supporting the IRECV/ISEND communication strategy.The tridiagonal solver can use a decomposition different from the one used to compute the Poisson equation’s right-hand side.Source: §3.5 · pp. 13–14
3.6. Profiling using NVTX
- NVTX annotations support profiler traces that show CPU and GPU activity; an example trace shows host-device copies overlapping MPI calls during TranYZ.The example is from one rank of a 1 × 16 processor-grid run on a 1024^3 mesh.Source: §3.6 · pp. 14–15, Figure 3 · p. 14
4. Code performance
AFiD shows strong scaling on both CPU and GPU, while the GPU version continues scaling on the 2048^3 grid after CPU efficiency plateaus. The optimal process-grid shape also shifts with processor count as transpose and halo-exchange communication change in importance.
4.1. Optimal configuration
- The code tests process-grid configurations using transpose and halo-exchange timings; the reported best shape shifts from probably 1×18 at low count to probably 32×32 at higher count.The reported timings are 0.5901 s for 18×1 in the low-count sweep and 0.03014 s for 32×32 in the high-count sweep.Source: §4.1 · pp. 15–17
- Four of six Poisson-solver transposes involve the first processor-grid dimension, favoring 1×N layouts until halo-exchange time becomes more dominant.At higher processor counts, the preferred two-dimensional decomposition minimizes halo-exchange communications; decomposition shape affects strong scaling.Source: §4.1 · p. 17
4.2. Performance comparison
- On a 2048^3 grid, 128 P100 GPUs reach performance close to about 6K CPU cores, and the GPU code continues scaling after CPU efficiency plateaus.Figure 4 presents the scaling comparison; the authors also report that upgraded P100 GPUs nearly double performance over K20X GPUs, with slightly worse strong scaling.Source: Figure 4 · p. 16, §4.2 · p. 18
- On the 2048×3072×3072 mesh, P100 GPUs’ 16 GB memory enables the simulation to run on 180 nodes and permits finer meshes.Source: §4.2 · p. 18
5. Validation
The GPU and CPU versions produce identical results up to machine precision, and AFiD’s Rayleigh–Bénard convection and plane Couette results agree with prior experimental or computational data. The section also reports flow-structure and profile features that frame these validations.
5. Validation
- GPU results are identical to CPU results up to machine precision across the validation cases.This indicates that the two code versions are consistent.Source: §5 · p. 20
5.1. Rayleigh-Bénard convection
- The Rayleigh-Bénard cases use an aspect-ratio-one cell, and Nusselt number is evaluated as the non-dimensional heat transport between the plates.The reported flow structures become smaller as Rayleigh number increases, motivating powerful simulations at high Rayleigh numbers.Source: §5.1 · pp. 20–21
- Experiments, simulations, and Grossmann-Lohse theory agree very well up to Ra = 1011, while experiments differ in the reported transition ranges.The Grenoble facilities show an increased Nusselt number around Ra = 5 × 1011; Göttingen experiments report transitions around Ra∗ 1 ≈2×1013 and Ra∗ 2 ≈7×1013, with no settled explanation for the discrepancy.Source: Figure 6 · p. 21, §5.1 · p. 21
- The Rayleigh-Bénard simulations use precisely assigned temperature boundaries, smooth surfaces, and the Boussinesq approximation, allowing control parameters to be adjusted and flow quantities to be accessed.The authors present these simulation features as useful for investigating the ultimate Rayleigh-Bénard transition.Source: §5.1 · pp. 21–22
5.2. Plane Couette flow
- At bulk Reynolds number Rec = 3000, the plane Couette simulation gives Reτ = 171, in excellent agreement with the prior study; its mean-velocity profiles also agree excellently with Ref..The profile has a viscous layer for y+ < 5 and a logarithmic layer for 50 < y+ < 171.Source: §5.2 · p. 22, Figure 7 · p. 23
- Plane Couette flow contains coherent high- and low-speed streaks spanning the channel, with a reported large-scale spanwise width of (4-5)h.The authors note that capturing the largest structures requires an extremely large box; the largest cited plane Couette DNS reaches Reτ ≈1000, versus Reτ ≈5200 for channel flow.Source: §5.2 · p. 23, Figure 8 · p. 24
6. Conclusions and future plans
The paper concludes that AFiD-GPU offers an order-of-magnitude wall-clock improvement over the CPU code and supports wall-bounded-flow studies. Future work includes hybrid CPU-GPU execution, asynchronous output, and multiple-resolution simulations.
6. Conclusions and future plans
- The GPU code achieves an order-of-magnitude wall-clock improvement over the CPU code, while the tested Rayleigh-Bénard and plane Couette cases validate its operation.The authors describe higher-Rayleigh- and higher-Reynolds-number simulations as more demanding than the tested cases.Source: §6 · p. 25
- The GPU solver improves AFiD’s capabilities for wall-bounded flows, including Rayleigh-Bénard convection and plane Couette flow, and the authors report initial GPU simulations of Rayleigh-Bénard flow with external shearing.They also note possible use for other high-Reynolds-number wall-turbulent flows, such as channel flow.Source: §6 · pp. 24–25
- A planned hybrid version will use otherwise idle CPU cores alongside GPUs in the implicit solver, splitting the vertical domain between the processors.A preliminary version is being tested; the split requires additional local CPU-GPU halo transfers.Source: §6 · pp. 25–26
- The authors are considering asynchronous HDF5 output so the GPU can advance the solution while CPUs finish writing copied data to disk.Source: §6 · p. 26
- A multiple-resolution approach uses a finer temperature grid than momentum grid to reduce computational time and memory occupancy.The authors report that about 80%–90% of resources are spent solving the momentum equations and use a separate refined time-step procedure to stabilize temperature integration.Source: §6 · p. 26
- The GPU code is planned for open-source release at afid.eu.Source: §6 · p. 26