# A Python/Fortran Implementation of the Lattice-Boltzmann Kernel on Multiple GPU Using the OpenACC Framework

Carlos Junqueira-junior, Erwan Zamora Medina, Noureddine Taibi, Simon Marié

► To cite this version:

Carlos Junqueira-junior, Erwan Zamora Medina, Noureddine Taibi, Simon Marié. A Python/Fortran Implementation of the Lattice-Boltzmann Kernel on Multiple GPU Using the OpenACC Framework. Concurrency and Computation: Practice and Experience, 2026, 38 (1), 10.1002/cpe.70518. hal-05448760

HAL Id: hal-05448760

https://hal.science/hal-05448760v1

Submitted on 8 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

OPEN ACCESS

# A Python/Fortran Implementation of the Lattice-Boltzmann Kernel on Multiple GPU Using the OpenACC Framework

Carlos Junqueira-Junior<fnref n="1" /> ID | Erwan Zamora Medina² ID | Noureddine Taibi² | Simon Marié²

**Correspondence:** Carlos Junqueira-Junior (junior.junqueira@ensam.eu)

**Received:** 1 July 2025 | **Revised:** 2 December 2025 | **Accepted:** 10 December 2025

**Keywords:** computational fluid dynamics | GPU | high-performance computing | lattice Boltzmann | OpenACC

## ABSTRACT

The increasing availability of GPU accelerated architectures for high-performance computing presents new opportunities for scientific software but also challenges due to the complexity of porting legacy codes to accelerator platforms. Directive-based programming models such as OpenACC offer a minimally intrusive pathway to exploit GPU acceleration without requiring extensive rewriting of existing codes. The current work presents a comprehensive performance and portability study of a LatticeBoltzmann Method solver (PyLB) originally written in Python, Mpi4Py, and Fortran for CPU architectures, which is ported to GPUs using OpenACC directives applied to the Fortran routines. The performance of the solver is evaluated on NVIDIA V100, A100, and H100 GPUs available on the Jean Zay supercomputer from Institute for Development and Resources in Intensive Scientific Computing (IDRIS) in France. Roofline analysis and extensive strong and weak scalability tests are conducted, showing that the GPU-enabled version of PyLB scales efficiently across multiple GPUs. The solver achieves performance on the H100 GPU equivalent to thousands of CPU cores and shows strong energy and carbon efficiency advantages over traditional CPU-based simulations. The implementation is validated using classical benchmarks, including the decaying Taylor-Green vortex and the flow over a 3-D sphere. The results confirm the physical accuracy of the GPU port while highlighting its computational and environmental advantages.

## 1 | Introduction

Graphical processing units (GPUs) are becoming increasingly prevalent in the high-performance computing (HPC) landscape, particularly in current exascale computing initiatives. French computing centers are actively embracing this trend by significantly integrating GPU architectures into their resources. Since the installation of the Joliot-Curie computer [1] from the TGCC in 2017, which predominantly comprised scalar nodes, subsequent supercomputers like Jean Zay [2], from IDRIS, and Adastra [3], from CINES, have been added to

the prestigious Top500 list [4] of supercomputers in 2018 and 2022, respectively. Notably, these systems incorporate approximately 30% and 40% GPU nodes, respectively. Recognizing the necessity of harnessing the computational capabilities offered by these advanced architectures, as detailed in references [5, 6], the research group has undertaken efforts to port a solver based on the Lattice-Boltzmann Method (LBM) to GPU platforms using OpenACC directives, a minimally intrusive approach that enables GPU support without requiring extensive code rewriting or significant changes to existing data structures [7–9].

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2026 The Author(s). Concurrency and Computation: Practice and Experience published by John Wiley & Sons Ltd.

<footnote n="1">¹Laboratoire DynFluid, Ecole Nationale Supérieur d'Arts et Métiers, Paris, France | ²Laboratoire DynFluid, Conservatoire National des Arts et Métiers, Paris, France</footnote>
The Lattice Boltzmann Method [10, 11] is nowadays recognized as a fast and reliable algorithm to numerically solve the Boltzmann equation. Moreover, it has already been shown to be efficient running on accelerator systems [12]. The physics of this kind of model is governed by a mesoscopic description of the collision between particles. Hence, if a given collision operator is chosen with a reliable equilibrium state, a wide variety of physical modeling could be obtained, from turbulent to relativistic flows [13–15]. In order to describe fluid dynamics, governed by the Navier-Stokes equation, the Bhatnagar-Gross-Krook (BGK) collision operator, based on a relaxation towards the equilibrium, has been shown to be an efficient mesoscopic description [16]. The LBM then performs a discretization of the velocity space in which the fluid particles are allowed to displace. This discretization has to be highly connected to the mesh and induces strong constraints in the choice of the lattice velocities. These constraints are often coupled to the algorithmic advection which basically relies on a collision and a propagation step. The propagation step is led by the mesh and the collision step depends on the description of the relaxation process. The traditional way to describe this step is to assign a relaxation parameter to the main statistical moments when they relax to their equilibrium state. This model, also referred to as the Multiple Relaxation Times (MRT) model [17] has been shown to recover the behavior of the weakly compressible Navier-Stokes equation. It has been shown [18, 19] that this kind of method has a lower dissipation error compared to traditional finite-difference schemes. As a counterpart, LBM suffers from numerical instabilities when the Reynolds number becomes high. The origins of these instabilities have been actively studied and remain an open subject [20–22]. Some modern collision models have now emerged and could improve those aspects by changing the moments definition [23], by re-normalizing the post-collision step [24], or by enforcing energy conservation [25–27]. A theoretical comparison of a wide variety of collision models is proposed in References [22, 28]. In the present study, the LBM method will be restricted to the standard BGK and MRT models with optimized relaxation times as defined in Reference [17].

The LBM solver used in the present work is the PyLB code [18, 29–31], an in-house code developed at the DynFluid laboratory in France. The code is designed to offer a fast and user-friendly hands-on approach, allowing undergraduate students and postdoctoral researchers who do not have extensive high-performance computing expertise to begin using PyLB for LBM calculations early in their research. The PyLB source code can be made available to the scientific community upon request to the authors. The PyLB framework solves the LBM equations and can be classified in the same category as other scientific solvers available in the literature, such as Palabos [32], waLBerla [33], and Sailfish [34].

PyLB is structured using a two-layer programming approach in which the high-level layer is implemented in Python, while the low-level layer is written in Fortran to execute the computational routines for solving the LBM equations. In the literature, programming languages dedicated to GPU computing, such as CUDA [35–37] and OpenCL [38], have been shown to achieve high performance by enabling optimization of low-level GPU instructions. However, this approach can be highly intrusive and often requires extensive code rewriting, which may alter

the solver's structure and hinder the hands-on philosophy when adapting PyLB to GPU architectures. Instead, GPU directives like OpenACC [7, 39] and OpenMP [40, 41] provide attractive alternatives to CUDA and OpenCL, as they allow code loops to be targeted for GPU execution without significant code modifications, while still delivering high computational performance. In this work, the OpenACC approach is adopted because its user community is more mature compared to that of OpenMP for GPUs.

The LBM algorithm has been extensively studied and optimized [42–44] to be run on a wide range of computing architectures from the early GPU implementations [45–47] to the AMD-GPUs [48] and FPGAs [49, 50]. The rise of powerful GPU devices for computational physics often comes with long software rewriting and optimization to fit the hardware specific language constraints (CUDA or Open-CL) needed to reach the newest hardware performances. When the development time is not a problem, an emerging strategy is to produce automatic code generation to fit the hardware performance [51]. However, these software optimizations could be demanding in the context of massive and historical codes with specific physical features. This study then proposes to demonstrate the performances obtained from a low-development strategy with GPU directives introduced in low-level existing codes. All the LBM simulations exposed in this work are based on this approach, which is implemented in a proper in-house parallel scientific Python/Fortran code called PyLB.

This article is structured as follows: the first section presents the theoretical and numerical formulation of the Lattice Boltzmann Method (LBM), followed by an overview of the solver algorithm and its structure. The third section presents the validation of the code using canonical flow configurations. The next section describes the computing platform used to evaluate the performance of the LBM framework on GPUs. The subsequent section details the implementation steps taken to port the code to OpenACC and analyzes its Roofline performance on a single GPU. This is followed by a discussion of the strong and weak scalability of the solver when using multiple graphics processing units. The last results section compares the performance, energy consumption, and carbon footprint of the solver on CPU and GPU architectures. Finally, the last section summarizes the conclusions drawn from this work.

## 2 | Numerical Formulation

### 2.1 | Lattice Boltzmann Method

The Lattice Boltzmann method [10], used to perform fluid flow simulations, is not directly based on the resolution of the Navier-Stokes equations but is a particular discretization of the BGK Boltzmann equation [16], which describes the dynamics of a gas with the BGK collision operator. The LBM equation is then obtained by restricting the velocities to a particular set of discrete velocities, which gives the discrete Boltzmann equation:

$$ \frac{\partial f_{\alpha}(\mathbf{x}, t)}{\partial t} + c_{\alpha,i} \frac{\partial f_{\alpha}(\mathbf{x}, t)}{\partial x_i} = -\frac{1}{\tau} (f_{\alpha} - f_{\alpha}^{eq}), \quad (1) $$

where $f_{\alpha}(\mathbf{x}, t)$ is the distribution of particles with a given velocity $c_{\alpha}$ at a given position $\mathbf{x}$ at time $t$. The collision between particles
is given by the BGK collision operator, which describes an average collision effect through the relaxation to a local equilibrium, $f^{eq}$, with a relaxation parameter, $\tau$. The usual 19-velocity lattice (D3Q19) [13] allows the recovery of dynamics described by the 3-D isothermal Navier-Stokes equations for small Mach numbers. From this lattice, the equilibrium function is written as:

$$f_{\alpha}^{eq}(\mathbf{x}, t) = \rho\omega_{\alpha}\left(1 + \frac{\mathbf{u}\cdot\mathbf{c}_{\alpha}}{\tilde{c}_{0}^{2}} + \frac{(\mathbf{u}\cdot\mathbf{c}_{\alpha})^{2}}{2\tilde{c}_{0}^{4}} - \frac{|\mathbf{u}|^{2}}{2\tilde{c}_{0}^{2}}\right), \quad (2)$$

in which the coefficients $\omega_{\alpha}$ and $\tilde{c}_{0}$ are defined by:

$$\left\{ \begin{array}{l} \omega_{\alpha} = \frac{1}{3}, \quad \frac{1}{18}, \quad \frac{1}{36} \quad \text{for } \alpha = 0, 1 \le \alpha \le 6, 7 \le \alpha \le 18, \\ \tilde{c}_{0}^{2} = \frac{1}{3}. \end{array} \right. \quad (3)$$

The conserved moments are then defined as:

$$\rho = \sum_{\alpha} f_{\alpha}, \quad \text{and} \quad \rho \mathbf{u} = \sum_{\alpha} \mathbf{c}_{\alpha} f_{\alpha}, \quad (4)$$

in which $\rho$ stands for the fluid density and $\mathbf{u} = (u_x, u_y, u_z)$ is the Cartesian flow velocity vector. The final step to obtain the LBM algorithm is to perform a space and time discretization. This is achieved by using the advective properties of the left-hand side of equation (1), which can be integrated along the characteristic $\mathbf{c}_{\alpha}$ to obtain the following LBM algorithm:

$$\left\{ \begin{array}{l} g_{\alpha}^{\text{coll}}(\mathbf{x}, t) = g_{\alpha}(\mathbf{x}, t) - \frac{dt}{\tau_g}(g_{\alpha}(\mathbf{x}, t) - g_{\alpha}^{\text{eq}}(\mathbf{x}, t)), \\ g_{\alpha}(\mathbf{x}, t) = g_{\alpha}^{\text{coll}}(\mathbf{x} - \mathbf{c}_{\alpha} dt, t - dt), \end{array} \right. \quad (5)$$

where the $g_{\alpha}$ distribution function comes from the integration step to obtain an explicit formulation and is related to the distribution $f_{\alpha}$ by the equation $g_{\alpha} = f_{\alpha} + \frac{dt}{2\tau}(f_{\alpha} - f_{\alpha}^{eq})$, which implies $g_{\alpha}^{eq} = f_{\alpha}^{eq}$ and $\tau_g = \tau + \frac{dt}{2}$. From this, the algorithm imposes $dt = dx = 1$ in order to be consistent with a uniform grid size imposed by the lattice. The physical grid step is obtained by discretizing a reference length scale, $L$, with a given number of points, $N$. Based on these parameters, the LBM algorithm can recover the Navier-Stokes dynamics with a second-order accuracy in space and time.

In the BGK collision operator, the distribution functions relax toward the equilibrium according to a single relaxation time. A more sophisticated idea is to relax each moment according to a proper relaxation time. This method is called multiple relaxation time (MRT) [17] and is known to alleviate some stability issues encountered with the BGK operator. The implementation of the MRT model is based on the modification of the collision step of equation (5), which is done in the momentum space:

$$\mathbf{m}^{\text{coll}}(\mathbf{x}, t) = \mathbf{m}(\mathbf{x}, t) - S(\mathbf{m}(\mathbf{x}, t) - \mathbf{m}^{\text{eq}}(\mathbf{x}, t)), \quad (6)$$

$$\mathbf{g}(\mathbf{x}, t) = \mathcal{M}^{-1} \mathbf{m}^{\text{coll}}(\mathbf{x} - \mathbf{c}_{\alpha} dt, t - dt), \quad (7)$$

where the matrix $\mathcal{M}$, transforms the $\mathbf{g}$ distribution functions into $\mathbf{m}$ moments:

$$\mathbf{m} = \mathcal{M} \mathbf{g}, \quad (8)$$

in which $\mathcal{M}$ is a square transformation matrix. If the first row of this matrix is filled only by 1, then the first moment is the density. The inverse transformation from the moments to distribution functions is simply $\mathbf{g} = \mathcal{M}^{-1}\mathbf{m}$. Further details about the construction of the matrix $\mathcal{M}$ can be found in [17].

The equilibrium moments are obtained from $\mathbf{m}^{eq} = \mathcal{M}\mathbf{g}^{eq}$. The diagonal of the $\mathcal{S}$ matrix corresponds to the inverse of the relaxation time, also called the relaxation rate, associated with each moment:

$$S = \operatorname{diag}(0, s_1, s_2, 0, s_3, 0, s_3, 0, s_3, s_v, s_2, s_v, s_2, s_v, s_v, s_v, s_4, s_4, s_4), \quad (9)$$

where $s_v$ is related to the fluid viscosity, $\nu$:

$$\frac{1}{s_v} = 3\nu - \frac{1}{2}. \quad (10)$$

The other relaxation rates, $s_1, s_2, s_3$, and $s_4$ do not appear in the macroscopic equations and are chosen according to stability optimization [17, 52], leading to $s_1 = 1.19, s_2 = 1.4, s_3 = 1.2$, and $s_4 = 1.98$. The BGK operator is recovered if all the relaxation rates are the same.

Variations of the collision operator can be found in the literature [23, 24, 53–57], as well as high-order discretizations of Equation (2), which can enhance numerical stability [58] and enable simulations at higher Mach numbers [59]. However, a detailed comparison of these alternative formulations lies beyond the scope of the present work and will not be addressed in this article. For a comprehensive theoretical comparison of various collision models and their impact on physical accuracy, the reader is referred to [28] and [22]. In this study, the widely used D3Q19 lattice with the BGK and MRT collision operators serves as the reference LBM model for all subsequent comparisons.

## 3 | Algorithm and Code Structure

The LBM algorithm used in this study relies on Equations (2) and (4) for the BGK model and on Equations (6) and (7) for the MRT model, with the physical parameters defined in Equation (10). The PyLB solver was originally developed for CPU architectures and is structured in two layers: a Python interface that manages reading the simulation setup, initializing conditions, memory allocation, and communication via the Mpi4Py library [60, 61]; and a Fortran layer for the core computational routines and boundary conditions. Figure 1 presents the flowchart of the algorithm, with the Python and Fortran layers highlighted in gray and blue, respectively. The use of the Mpi4Py library is represented in light green.

To enable GPU acceleration, the research group opted to port the LBM solver using OpenACC directives applied to the Fortran routines. This programming model is minimally intrusive, allowing GPU support without extensive code refactoring or modification of existing data structures [7]. The introduction of this paradigm has been described in the literature as relatively straightforward [7, 37, 39]. However, due to the two-layer design of the PyLB solver, integrating the required compiler tool-chain and associated compilation flags to support Python and Fortran code, OpenACC directives, and MPI communication posed a challenge.
<figure><img src="image_3.png" /><figcaption>FIGURE 1 | Structure of the CPU and OpenACC implementations of the PyLB solver.</figcaption></figure>

Resolving these issues involved support from the technical staff of the Jean Zay supercomputer, as this process required specific knowledge of the system architecture, compiler configuration, operating system, and batch scheduler setup.

Additionally, the original CPU version used the Mpi4Py library, which is not GPU-aware and does not allow explicit access to intra-node MPI task indices [62], which imposes implementation challenges on multi-GPU intra-node communications. As a result, the research group wrote all MPI exchanges into the Fortran layer to ensure compatibility with the OpenACC implementation. The Python Mpi4Py library remains in the code solely to perform domain decomposition during the initialization stage. This avoids extensive code rewriting and has no impact on performance, since this operation is only executed once and is not part of the compute-intensive routines. Overall, setting up the compilation environment and rewriting MPI-related code proved to be more complex than inserting the OpenACC directives into the Fortran routines themselves.

The flowchart of the OpenACC-enabled version of PyLB is also shown in Figure 1, where the Python and Fortran layers are again indicated in gray and blue, respectively. The MPI routines implemented in Fortran are marked in dark green, and the OpenACC directives are shown in beige.

## 4 | Numerical Validation

The present section validates the OpenACC implementation of the PyLB solver through the simulation of two well-documented physical configurations: the 3-D decaying Taylor-Green vortex and the flow over a 3-D sphere. The former is used for the Roofline, strong scalability, and energy & carbon footprint analyses, while the latter serves in the weak scalability study. Both test cases are run on the GPU partitions of the Jean Zay supercomputer, using the BGK and MRT collision models, and the results are compared against theoretical, numerical, and experimental data available in the literature. This validation ensures that the modifications do not compromise the ability of the solver to accurately handle realistic physical flow configurations.

## 4.1 | Taylor-Green Vortex Decay

The decaying Taylor-Green Vortex (TGV) test case is a well-established benchmark for evaluating vortex stretching and the generation of small-scale eddies. This test case provides insights into the transition to turbulence and has been extensively applied to analyze dissipation errors in numerical schemes, as well as the influence of collision operators in the lattice Boltzmann method (LBM) [63]. The CPU version of the PyLB solver has already been validated using the TGV test case as a comparison [29, 64, 65].

In the current study, simulations are conducted within a 3-D periodic cubic domain $\Omega$ defined as $0 \le x, y, z \le 2\pi$, with a Reynolds number set to $Re = 1600$. The initial conditions for the Taylor-Green vortex are prescribed by defining the velocity and pressure fields, $p$, as follows:

$$ 
\left\{ 
\begin{aligned}
p &= p_{\infty} + \frac{\rho_{\infty} U_{\infty}^2}{16} [\cos(2z) + 2][\cos(2x) + \cos(2y)] \\
u_x &= U_{\infty} \sin(x) \cos(y) \cos(z) \\
u_y &= -U_{\infty} \cos(x) \sin(y) \cos(z) \\
u_z &= 0
\end{aligned} 
\right. 
\quad (11) 
$$

in which the $\infty$ subscripts represent reference values. All quantities are normalized using a characteristic length scale $L$ and velocity scale $U_\infty$, ensuring a dimensionless formulation of the governing equations. In this context, the dimensionless time is defined as $t^* = tU_\infty/L$, while the Reynolds number is expressed as $Re = U_\infty L/v$. The initialization of the distribution functions follows the methodology proposed by Skordos [66] and Mei et al. [67], which accounts for the presence of significant initial velocity gradients. To validate the GPU code, the entrophy, $\mathcal{E}$, in the domain is computed as follows:

$$ \omega_{ij} = \left( \frac{\partial u_i}{\partial x_j} - \frac{\partial u_j}{\partial x_i} \right) \quad (12a) $$

$$ \mathcal{E} = \frac{1}{\rho_{\infty} \Omega} \int_{\Omega} \rho \frac{\omega \cdot \omega}{2} d\Omega \quad (12b) $$
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Comparison of enstrophy between PyLB for BGK and MRT models at two different resolutions and the spectral reference [68].</figcaption></figure>

The discrete computation of the vorticity field $\omega$ uses a second-order estimation of velocity gradients, which can then be integrated into the enstrophy $\mathcal{E}$.

The validation study is conducted using two domain sizes, 256³ and 512³ grid points, uniformly distributed in all three spatial directions, and both collision models available in the OpenACC version of the PyLB code. Simulations are executed using different numbers of GPUs and compute nodes to assess the impact of OpenACC acceleration and MPI communication on the physical results. All calculations produced identical enstrophy values up to the 14th significant digit when using different numbers of resources for the same test case. Figure 2 compares the evolution of $\mathcal{E}$ enstrophy computed by PyLB against spectral results from Wang et al. [68]. It can be observed that the numerical calculations for $N = 256^3$ slightly underestimate the enstrophy peak at $t^* \approx 9$, which is primarily attributed to the numerical dissipation inherent to the second-order LBM scheme. However, results indicate good agreement for the enstrophy calculation, when spectral and LBM resolutions match, as in the case of $N = 512^3$.

## 4.2 | Flow Over a 3-D Sphere Test Case

The LBM solver is also validated using the 3-D flow over a sphere to evaluate boundary conditions other than the periodic ones imposed in the TGV configuration. This second test case involves a sphere placed in a uniform flow at $Re = 300$, where $Re = u_\infty D/v$, with $u_\infty = 0.1c_s$, $D = 50$ lattice units, and kinematic viscosity $\nu$ tuned to match the Reynolds number. The numerical setup includes a simulation box spanning $12D \times 6D \times 6D$, ensuring minimal blockage effects while resolving wake dynamics. The sphere is located inside the simulation box at $(x_c, y_c, z_c) = (2D, 3D, 3D)$. Boundary conditions are implemented as follows: no-slip walls based on the bounce-back (BB) scheme [69] at the sphere, an equilibrium-based Dirichlet condition at the inlet enforcing $(\rho_\infty, u_\infty, v_\infty, w_\infty) = (1.0, 0.1c_s, 0.0, 0.0)$, and periodic boundary conditions on the other four domain surfaces. The sphere surface is resolved with approximately 50 lattice units across its diameter.

Figure 3 presents a snapshot of the iso-contour of the $\lambda_2$ criterion in the wake behind the sphere, colored by the stream-wise velocity $u_x$. This visualization serves as a qualitative representation of the vortex structures in the flow over the sphere [70], as calculated using the PyLB solver, and illustrates the characteristic hairpin vortex structures observed in the literature [71, 72] for this flow configuration. Table 1 presents the values of the drag coefficient ($C_d$), vertical lift coefficient ($C_l$), and Strouhal number ($St$) obtained using the present method and compares them with numerical and experimental reference data. The definitions of these coefficients and the Strouhal number can be found in Reference [73]. The results produced by the PyLB solver with the OpenACC implementation for the flow over a sphere at $Re = 300$, show good agreement with values reported in the literature, supporting the validity of the code enhancements introduced in this work.

# 5 | Computing Architecture Description

For the present work, the computing performance of the solver is evaluated on the Jean Zay system [2], a state-of-the-art supercomputer that features both an HPE SGI 8600 section and an Eviden BullSequana XH3000 section, comprising a total of five partitions. The latter was added during the upgrade of Jean Zay in July 2024. These partitions include one with CPU-only nodes, three with GPU-accelerated nodes, and one dedicated to pre- and post-processing as well as visualization, each designed to meet a range of computational needs. The HPE SGI 8600 compute nodes are interconnected via a 100 Gb/s Intel Omni-Path (OPA) network [80], while the newer Eviden BullSequana XH3000 nodes use a 400 Gb/s InfiniBand network [81, 82]. All nodes have access to a high-bandwidth parallel file system. As of the mentioned upgrade, the Jean Zay system's cumulative peak performance has reached 125.9 petaflops.

The scalar partition includes 720 compute nodes, each equipped with two Intel Cascade Lake (CSL) 6248 processors, providing a total of 40 cores per node, with 192 GB of memory per node. The accelerated partitions of Jean Zay provide multiple configurations
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | Iso-contour of $\lambda_2$ criterion on the wake past the sphere, colored by stream-wise velocity $u_x$.</figcaption></figure>

<table><caption>**TABLE 1** | Comparison of drag and lift coefficients, and Strouhal number for flow past a sphere at $Re = 300$. The superscript (*) denotes experimental results.</caption><thead><tr><th rowspan="2">Authors</th><th colspan="3">Re = 300</th></tr><tr><th>C<sub>d</sub></th><th>C<sub>l</sub></th><th>St</th></tr></thead><tbody><tr><td>Roos and Willmarth* [74]</td><td>0.629</td><td>—</td><td>—</td></tr><tr><td>Johnson and Patel [72]</td><td>0.656</td><td>-0.069</td><td>0.137</td></tr><tr><td>Tomboulides and Orszag [75]</td><td>0.671</td><td>—</td><td>0.136</td></tr><tr><td>Constantinescu and Squires [76]</td><td>0.655</td><td>-0.065</td><td>0.136</td></tr><tr><td>Kim and Choi [77]</td><td>0.657</td><td>-0.067</td><td>0.134</td></tr><tr><td>Ploumhans et al. [78]</td><td>0.683</td><td>-0.061</td><td>0.135</td></tr><tr><td>Mimeau et al. [79]</td><td>0.673</td><td>-0.066</td><td>0.133</td></tr><tr><td>PyLB-OpenACC</td><td>0.695</td><td>-0.065</td><td>0.139</td></tr></tbody></table>

for varied computational demands. The primary GPU partition, comprising 396 four-GPU nodes, features dual CSL 6248 processors, with 40 cores and 192 GB of memory per node. Among these nodes, 126 are equipped with four Nvidia Tesla Volta (V100) SXM2 16GB GPUs while the remaining 270 nodes include four Nvidia V100 SXM2 32GB GPUs. Additionally, the system includes 31 compute nodes, each outfitted with eight Nvidia V100 SXM2 32GB GPUs. These nodes feature dual Intel Cascade Lake 6226 processors, with 24 cores per node. Of these nodes, 20 provide 384 GB of memory, while the remaining 11 offer 768 GB of memory.

Another GPU partition includes 52 nodes, each equipped with eight Nvidia Ampere (A100) SXM4 80GB GPUs. These nodes are powered by dual AMD Milan EPYC 7543 processors, with a total of 64 cores and 512 GB of memory per node. The latest expansion, the Eviden Bull Sequana XH300 section, includes a new GPU partition with 364 nodes. Each node is equipped with four Nvidia Hopper (H100) SXM5 80GB GPUs and dual Intel Xeon Platinum 8468 processors, providing a total of 96 cores per node, with 512 GB of memory. Additionally, Jean Zay offers a high-performance parallel file system, including 2.5 PB of SSD storage and over 30 PB of disk capacity, to support data-intensive computations. The system operates under Red Hat version 8.6 with Slurm version 23.02.6 as its workload manager.

The NVIDIA HPC Software Development Kit (SDK) [83] is used in this work to compile the Fortran layer of the PyLB solver, specifically, version 23.1 for the V100 and A100 architectures, and version 24.3 for the H100 GPU. An open-source implementation of the Message Passing Interface (MPI), OpenMPI [84], is employed to handle data exchange among MPI tasks. For the V100 and A100 architectures, OpenMPI version 4.1.1 is used in conjunction with CUDA toolkit version 12.0.0. For the H100 GPU, OpenMPI version 4.1.5 is used with CUDA toolkit version 12.4.1. The Python layer of the LBM solver is executed using Python version 3.10.4 across all three accelerators considered in this study. The choice of compiler and software versions is dictated by the software environment provided on the Jean Zay supercomputer.

## 6 | Acceleration of PyLB on a Single GPU Using OpenACC

The LBM solver was originally developed with a Python interface that handles memory allocation and communication via the Mpi4Py library, while the core computational routines are implemented in a lower-level Fortran layer, as illustrated in Figure 1. As an initial step toward GPU acceleration, the research group began porting the solver using OpenACC directives applied to the Fortran routines responsible for computing the BGK collision model. The OpenACC paradigm is known to be minimally intrusive, avoiding the need for extensive code rewriting or modifications to existing data structures when porting frameworks to GPU architectures [7, 37, 39]. In this study, the GPU implementation is based on the use of *acc kernels* directives, which define regions of code to be offloaded to the accelerator, while delegating architecture-specific optimizations to the compiler.

Algorithm 1 illustrates the first version of PyLB implementing the BGK collision model, ported to GPUs using the OpenACC standard. This version is referred to as BGK4 in the present work. GPU directives are applied to the four main computational blocks that together perform a complete iteration of the solver. The first block, spanning lines 1 to 7, computes four 3-D arrays: the density $\rho$ and the three velocity components $u_i$. The second block, from lines 8 to 12, calculates the 4-D array representing the equilibrium distribution function $g^{\text{eq}}$, as defined in Equation (2). The third block, from lines 13 to 16, computes the 4-D array $g^{\text{coll}}$. The
1: $acc kernels

2: **for** k, j, i **do**

3: **for** p **do**

4: $\rho_{i,j,k}^+ = g_{i,j,k,p}$

5: **end for**

6: **end for**

7: **for** k, j, i **do**

8: $u_{i,j,k}^{l=1,2,3} = c_p^l * (g_{i,j,k,p}) / \rho_{i,j,k}$

9: **end for**

10: $acc end kernels

11: $acc kernels

12: **for** k, j, i **do**

13: **for** p **do**

14: $g_{i,j,k,p}^{eq} = \text{Eq. (2)}$

15: **end for**

16: **end for**

17: $acc end kernels

18: $acc kernels

19: **for** p, k, j, i **do**

20: $g_{i,j,k,p}^{coll} = g_{i,j,k,p} - \tau_g * (g_{i,j,k,p} - g_{i,j,k,p}^{eq})$

21: **end for**

22: $acc end kernels

23: $acc kernels

24: **for** p, k, j, i **do**

25: $g_{i,j,k,p} = g_{i,j,k,p}^{coll} * (g_{i,j,k,p}^{eq}) * (g_{i,j,k,p}^{coll})$

26: **end for**

27: $acc end kernels

final block, covering lines 17 to 20, performs the propagation step, updating the 4-D distribution function g. The loop indices i, j, and k denote the three spatial dimensions, while p indexes the 19 discrete lattice velocities.

The NVIDIA Nsight Systems performance analysis tool [85] is employed to evaluate the performance of the LBM solver when simulating the 3-D decaying TGV test case. Three domain sizes are considered for this evaluation: 128, 256, and 512 points in each of the three directions. The calculations are run for three iterations and the Nsight analysis toolkit repeats each one several times to collect detailed statistical data that is post-processed within the context of the Roofline model. This model is a performance analysis framework that plots achieved floating-point performance (in FLOP/s) against operational intensity (FLOPs per byte of memory traffic), illustrating whether an application is limited by memory bandwidth or compute capacity [86]. To determine the hardware-specific ceilings required for this analysis, the Empirical Roofline Tool (ERT) [87–89] is used to benchmark the characteristics of the accelerated nodes of the Jean Zay supercomputer equipped with V100 (32GB), A100, and H100 GPUs. The benchmark collects the attainable memory bandwidth (in TB/s) and computational throughput (in TFLOP/s) of these nodes.

The post-processed Roofline performance data corresponding to Alg. 1 are shown in Figure 4, where blue, red, and green colors denote the V100, A100, and H100 GPU architectures, respectively. Square, diamond, and circle markers indicate the Taylor–Green vortex test cases using domain sizes of 128³, 256³,

1: $acc kernels

2: **for** k, j, i **do**

3: $\rho = g_{i,j,k,1} + g_{i,j,k,2} + \cdots + g_{i,j,k,19}$

4: $u_{i,j,k}^{l=1,2,3} = c_p^l * (g_{i,j,k,p}) / \rho$

5: $g_1^{eq} = \text{Eq. (2)}$ **for** p = 1

6: $g_2^{eq} = \text{Eq. (2)}$ **for** p = 2

7: $\vdots$

8: $g_{19}^{eq} = \text{Eq. (2)}$ **for** p = 19

9: $g_{i,j,k,1}^{coll} = g_{i,j,k,1} - \tau_g * (g_{i,j,k,1} - g_1^{eq})$

10: $g_{i,j,k,2}^{coll} = g_{i,j,k,2} - \tau_g * (g_{i,j,k,2} - g_2^{eq})$

11: $\vdots$

12: $g_{i,j,k,19}^{coll} = g_{i,j,k,19} - \tau_g * (g_{i,j,k,19} - g_{19}^{eq})$

13: **end for**

14: $acc end kernels

15: $acc kernels

16: **for** k, j, i **do**

17: $g_{i,j,k,1} = g_{i,j,k,1}^{coll} * (g_{i-j-c_1^1}, (j-c_1^2), (k-c_1^3), 1)$

18: $\vdots$

19: $g_{i,j,k,19} = g_{i-j-c_{19}^1} * (g_{i-j-c_{19}^2}, (j-c_{19}^3), (k-c_{19}^4), 19)$

20: **end for**

21: $acc end kernels

and 512³ points, respectively. Filled markers represent the original BGK4 implementation in the PyLB code. Among the three GPUs, only the H100 was capable of executing the 512³ case within the GPU memory and run-time constraints imposed by the Jean Zay batch scheduler. The V100 was unable to run this case due to its limited 32 GB of memory, compared to the 80 GB available on both the A100 and H100. The A100, while capable of handling the memory load, exceeded the maximum job time allowed by the scheduler when executing three TGV iterations with the Nsight profiling tool.

The BGK implementation from Alg. 1 exhibits low arithmetic intensity relative to the peak TFLOP/s ceilings of the tested GPUs, indicating that the solver is primarily memory-bound. Interestingly, the measured arithmetic intensity increases with domain size, an unexpected trend given that the algorithm itself remains unchanged. This behavior stems from the original code structure, which was designed for readability rather than performance. As a result, Alg. 1 includes several performance-limiting factors by construction, such as multiple acc kernels spread across Fortran routines and non-contiguous or data-dependent loop operations. The former directives introduce overhead from kernel launches, while the latter inhibit effective vectorization. Notably, the inner loops on lines 3 and 10 of Alg. 1 are compiled as sequential due to non-contiguous memory access and data dependencies. These issues collectively contribute to the observed growth in arithmetic intensity, even for a fixed algorithm (BGK4), when increasing problem size.

To better exploit the characteristics of accelerator architectures, the PyLB-BGK solver is improved as shown in Alg. 2, hereafter referred to as BGK2. In this version, the four 3-D arrays ($\rho$ and $u_l$), as well as the 4-D equilibrium distribution function, $g^{eq}$, are replaced with scalar variables since they are temporary variables, and all computations involving the lattice velocity index p
<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | Roofline performance model of PyLB using the BGK2, BGK4, and MRT collision algorithms for the decaying Taylor-Green vortex test case on different GPU architectures available on the Jean Zay supercomputer. Blue (—), red (—), and green (—) colors represent the V100, A100, and H100 GPUs, respectively. The □, ◦, and ○ markers correspond to domain sizes of 128³, 256³, and 512³ points, respectively. Filled markers indicate the BGK4 algorithm, empty markers represent the BGK2 implementation, and lightly filled markers denote the MRT algorithm.</figcaption></figure>

are explicitly unrolled. These modifications eliminate loops over $p$, thereby avoiding loop-carried data dependencies and improving memory coalescence. Such changes enhance vectorization, which is a key performance factor on accelerator architectures. Additionally, there are no longer any sequential loops reported during the code compilation. Another notable improvement is that the collision routine is now implemented as a single major *acc kernel*, eliminating the overheads observed in the BGK4 version.

The optimizations introduced in the BGK2 algorithm lead to a significant performance improvement: approximately a ten-fold increase on the V100 and A100 GPUs, and a 25% gain on the H100, in both arithmetic intensity and computational throughput, compared to the legacy BGK4 implementation. The relatively smaller improvement observed on the H100 can be attributed to its advanced architecture, which already delivers high baseline performance even for less optimized code. With significantly higher memory bandwidth and more sophisticated latency-hiding mechanisms, the H100 is more resilient to the inefficiencies of the BGK4 implementation. Consequently, the impact of the BGK2 optimizations is less pronounced compared to older GPUs. Furthermore, the compiler toolchain for the H100 may apply more effective low-level optimizations, and the optimized BGK2 version may approach architectural resource limits, such as register usage or shared memory, which can constrain additional performance gains. Algorithm 2 achieves approximately 0.8, 20, and 40 TFLOP/s on the V100, A100, and H100 GPUs, respectively. Moreover, the arithmetic intensity is no longer dependent on the problem size.

The MRT collision model is also ported to GPUs using the OpenACC standard, and its Roofline performance is reported in

Figure 4 using lightly filled markers. Its pseudo-code is presented in Alg. 3. The code modifications follow a similar approach to those in the BGK2 algorithm, where $\rho$ and $u_l$ are stored as scalar variables, along with all temporary variables listed between lines 7 and 10. The computation and update of the 4-D distribution functions $g$ and $g^{\text{coll}}$ are unrolled, eliminating the need for a fourth loop over the lattice velocity index $p$. When comparing the MRT and BGK2 algorithms, the former exhibits slightly higher arithmetic intensity, which consequently results in slightly better arithmetic performance.

# 7 | Multiple GPU Development: Scalability Study

Computational fluid dynamics can require a significant amount of computational resources to produce high-quality results. Therefore, it is important to ensure that a computational solver is capable of using parallel resources efficiently. Strong and weak scalability have been widely used in the CFD community for decades as performance metric for parallel algorithms [90–95]. The former is based on measuring the computing time when increasing the number of computational resources for fixed problem size while the latter measures the cost of simulations when increasing both the problem domain size and the computing resources proportionally. The same metrics are applied in the current study to evaluate the performance of the enhanced version of the Lattice Boltzmann solver using OpenACC directives.

## 7.1 | Strong Scaling

The strong scalability study [90] is performed to evaluate the performance of the solver as computational resources increase
1: $acc \text{ kernels}
2: \textbf{for } k, j, i \textbf{ do}
3:     $\rho = g_{i,j,k,1} + g_{i,j,k,2} + \cdots + g_{i,j,k,19}$
4:     $u^{1,2,3} = c_p^l * (g_{i,j,k,p}) / \rho$
5:     // Update distribution functions $g^{\text{coll}}$
6:     // based on Eq. 6 using scalar temp.
7:     vars.
8:     $g_{i,j,k,1}^{\text{coll}} = \text{Eq. (6) for } p = 1$
9:     $g_{i,j,k,2}^{\text{coll}} = \text{Eq. (6) for } p = 2$
10:    $\vdots$
11:    $g_{i,j,k,19}^{\text{coll}} = \text{Eq. (6) for } p = 19$
12: \textbf{end for}
13: \textbf{end for}
14: \textbf{for } k, j, i \textbf{ do}
15:     $g_{i,j,k,1} = g_{i-c_1^1}^{\text{coll}} \cdot (j-c_1^2) \cdot (k-c_1^3), 1$
16:     $g_{i,j,k,2} = g_{i-c_2^1}^{\text{coll}} \cdot (j-c_2^2) \cdot (k-c_2^3), 2$
17:     $\vdots$
18:     $g_{i,j,k,19} = g_{i-c_1^1}^{\text{coll}} \cdot (j-c_1^2) \cdot (k-c_1^3), 19$
19: \textbf{end for}
20: \textbf{end for}

for a fixed computational domain. The study is conducted on NVIDIA V100, A100, and H100 GPU architectures available on the Jean Zay supercomputer. The MRT and BGK models are also evaluated in this study. Moreover, various problem sizes are considered, defined by the number of grid points per spatial direction, denoted as $N_p$, which remains constant in all three Cartesian directions, such that $N_p = N_x = N_y = N_z$. Consequently, the total number of grid points in the domain is given by $N_p^3$. In the present work, $N_p$ starts at 128 and is progressively doubled up to 1024, that is, $N_p = 128, 256, 512$, and $1024$, resulting in an eight-fold increase in the total domain size. For each fixed domain size, the number of allocated GPUs, denoted as $N_g$, ranges from 1 and is progressively doubled up to 1024 GPUs, depending on the memory capacity of the architecture in use.

Each calculation runs 1000 iterations, $N_{it}$, to measure the total computational time of the solver, denoted as $T$. Once the data are collected, the speedup, $S$, the performance in Lattice Updates Per Second (LUPS), and the strong scaling efficiency, $\eta^{\text{stg}}$, are computed. The speedup for a given problem size $N_p$ running on $N_g$ is defined as:

$$S_{N_g}^{N_p} = \frac{T_0^{N_p}}{T_{N_g}^{N_p}}, \qquad (13)$$

where $T_0^{N_p}$ and $T_{N_g}^{N_p}$ represent the computational time required to complete 1000 iterations using $N_g$ graphics cards and the baseline configuration, defined as the starting point for scalability analysis, respectively. The LUPS metric, widely used in the Lattice Boltzmann community [65, 96, 97], quantifies the number of lattice updates per second. In the present work, it is reported in billions of lattice updates per second (GLUPS) and defined as:

$$GLUPS_{N_g}^{N_p} = \frac{N_p^3 \times N_{it}}{10^9 \times T_{N_g}^{N_p}}. \qquad (14)$$

<table><caption>TABLE 2 | Scaling starting point for different workloads and GPU architecture.</caption><thead><tr><th rowspan="2">N<sub>p</sub><sup>3</sup></th><th colspan="3">N<sub>g0</sub><sup>N<sub>p</sub></sup></th></tr><tr><th>V100</th><th>A100</th><th>H100</th></tr></thead><tbody><tr><td>128<sup>3</sup></td><td>1</td><td>1</td><td>1</td></tr><tr><td>256<sup>3</sup></td><td>1</td><td>1</td><td>1</td></tr><tr><td>512<sup>3</sup></td><td>2</td><td>1</td><td>1</td></tr><tr><td>1024<sup>3</sup></td><td>16</td><td>8</td><td>8 (BGK) &amp; 4 (MRT)</td></tr></tbody></table>

The strong scaling efficiency $\eta^{\text{stg}}$ for a problem size of $N_p$ grid points per spatial direction and using $N_g$ is given by:

$$\eta_{N_g}^{N_p, \text{stg}} = \frac{S_{N_g}^{N_p} \times N_{g0}^{N_p}}{N_g}, \qquad (15)$$

in which $N_{g0}^{N_p}$ is the number of GPUs used as the scaling starting point for the given problem size $N_p$. Ideally, $N_{g0}^{N_p}$ corresponds to a single GPU, and $T_0^{N_p}$ represents the execution time measured with one GPU for the $N_p$ problem size. However, depending on the problem size and the memory constraints of the GPU architecture, this baseline configuration may not always be feasible. Table 2 indicates the starting point for each workload when using different GPUs. The V100 architecture can accommodate up to $N_p^3 = 256^3$ grid points per GPU for the TGV test case in the PyLB solver, as it has the smallest total memory among the three architectures compared in this study, with 32 GB of HBM2 memory. In contrast, the A100 and H100 architectures, each equipped with 80 GB of HBM2e and HBM3 memory, respectively, can handle up to $512^3$ grid points for the 3-D TGV configuration in PyLB.

Figure 5 presents the GLUPS performance of both MRT and BGK models implemented in the PyLB solver across the three GPU architectures available on the Jean Zay supercomputer. The results indicate that the solver scales efficiently on all three NVIDIA GPU nodes tested in this study. Additionally, large workloads contribute to high GLUPS performance as the number of GPUs operating in parallel increases. Moreover, the ratio of $N_p^3$ to memory capacity can be a key factor when selecting an optimal GPU configuration since GPU architectures tend to achieve higher computational efficiency when fully utilized, as they primarily rely on SIMT parallelism. This effect is evident in Figure 5, where the V100 GPUs reach the same GLUPS performance as the A100 GPUs under specific conditions: at $N_p^3 = 128^3$ for $16 \le N_g \le 32$; at $N_p^3 = 256^3$ for $64 \le N_g \le 128$; and at $N_p^3 = 512^3$ and $N_p^3 = 1024^3$ for $N_g = 128$.

The LBM solver scales efficiently on both the V100 and A100 nodes of Jean Zay when employing CUDA-aware MPI exchanges, which leverage GPU-Direct technology to enable direct communication between GPU memory buffers. For all evaluated workloads, $128^3 \le N_p^3 \le 1024^3$, once the number of GPUs reaches the per-node limit ($N_g = 4$ for the V100 architecture and $N_g = 8$ for the A100 GPUs) MPI communications must be performed over the InfiniBand network. This introduces higher communication overhead and results in a decline in performance. However, this limitation is less pronounced for PyLB on the H100 computing
<figure><img src="image_7.png" /><figcaption>**FIGURE 5** | Performance in GLUPS from the strong scalability study of the TGV test case using different problem sizes, $N_p^3 = 128^3, 256^3, 512^3$, and $1024^3$, on the three GPU architectures available on the Jean Zay supercomputer. The violet (—), green (—), and yellow (—) colors represent the V100, A100, and H100 GPU architectures, respectively. The square markers (■—) represent the MRT model, while the dashed line with circular markers (○—) corresponds to the BGK model. The ○ marker highlights the V100's peak GLUPS result and the minimum number of A100 or H100 accelerators required to outperform it. (a) $N_p^3 = 128^3$, (b) $N_p^3 = 256^3$, (c) $N_p^3 = 512^3$, (d) $N_p^3 = 1024^3$.</figcaption></figure>

<table><caption>**TABLE 3** | Maximum V100 GLUPS performance equivalence on A100 and H100 architectures for both BGK and MRT models.</caption><thead><tr><th rowspan="2">BGK</th><th colspan="2">V100</th><th colspan="2">A100</th><th colspan="2">H100</th></tr><tr><th>N<sub>p</sub></th><th>GLUPS</th><th>N<sub>g</sub></th><th>GLUPS</th><th>N<sub>g</sub></th><th>GLUPS</th><th>N<sub>g</sub></th></tr></thead><tbody><tr><td>128</td><td>1.31</td><td>64</td><td>1.56</td><td>4</td><td>1.33</td><td>1</td></tr><tr><td>256</td><td>5.29</td><td>128</td><td>6.85</td><td>8</td><td>5.82</td><td>4</td></tr><tr><td>512</td><td>11.38</td><td>128</td><td>11.67</td><td>8</td><td>15.88</td><td>8</td></tr><tr><td>1024</td><td>22.12</td><td>128</td><td>22.64</td><td>32</td><td>36.18</td><td>16</td></tr></tbody></table><table><thead><tr><th rowspan="2">MRT</th><th colspan="2">V100</th><th colspan="2">A100</th><th colspan="2">H100</th></tr><tr><th>N<sub>p</sub></th><th>GLUPS</th><th>N<sub>g</sub></th><th>GLUPS</th><th>N<sub>g</sub></th><th>GLUPS</th><th>N<sub>g</sub></th></tr></thead><tbody><tr><td>128</td><td>1.31</td><td>64</td><td>1.56</td><td>4</td><td>1.60</td><td>2</td></tr><tr><td>256</td><td>5.22</td><td>128</td><td>6.84</td><td>8</td><td>8.20</td><td>8</td></tr><tr><td>512</td><td>10.42</td><td>128</td><td>11.73</td><td>32</td><td>13.20</td><td>8</td></tr><tr><td>1024</td><td>17.48</td><td>128</td><td>20.98</td><td>32</td><td>21.55</td><td>16</td></tr></tbody></table>

nodes analyzed in this study, as the solver demonstrated scalability up to 128 GPUs for all workloads. This improvement can be attributed to the fact that the Jean Zay nodes equipped with H100 GPUs are connected via a 400 Gb/s InfiniBand network, as opposed to the 100 Gb/s Intel Omni-Path interconnect used by the V100 and A100 nodes of the IDRIS supercomputer. Moreover, the H100 GPUs use third-generation high-bandwidth memory (HBM3), which offers a substantial increase in bandwidth, approximately 3 TB/s, compared to the A100's HBM2e (around 2 TB/s) and the V100's HBM2 (around 900 GB/s). This combination of enhanced memory bandwidth and improved network connectivity significantly accelerates data transfers, benefiting large-scale HPC applications.

As expected, newer GPU generations provide increased computing power, as reflected in the performance progression from the V100 to the A100 and finally to the H100 architecture. Table 3 presents, for each workload evaluated in this study and for both LBM models, the best performance obtained by the PyLB solver on $N_g$ V100 GPUs from Jean Zay, as well as the number of A100
<figure><img src="image_8.png" /><figcaption>**FIGURE 6** | Speedup and efficiency from the strong scalability study of the TGV test case with the LBM-BGK model using different problem sizes, $N_p = 128, 25, 512$, and $1024$, on the three GPU architectures available on the Jean Zay supercomputer. The violet (—), green (—), and yellow (—) colors represent the V100, A100, and H100 GPUs architectures, respectively. The solid line with the square marker (■—) represents the speedup, while the dashed line with the circle marker (○—) stands for the efficiency. The ○ marker highlights the V100's peak GLUPS result and the minimum number of A100 or H100 accelerators required to outperform it, while the △ marker indicates the largest $N_p$ that can be allocated while maintaining a strong scalability efficiency above 50%.
(a) $N_p = 128$, (b) $N_p = 256$, (c) $N_p = 512$, (d) $N_p = 1024$.</figcaption></figure>

and H100 GPUs required to surpass this reference performance. These values are also highlighted in Figures 5–7, where circle markers indicate the maximum GLUPS achieved on the V100 and the corresponding number of A100 and H100 GPUs needed to exceed it in the strong-scaling evaluation. For the largest problem size, 1024³ grid points, only 32 A100 GPUs and 16 H100 GPUs are required to outperform the scalability results obtained with 128 V100 GPUs, corresponding to speedup factors of approximately 4 and 8, respectively.

Both BGK and MRT Lattice Boltzmann collision models exhibit similar GLUPS performance for workloads of $N_p^3 = 128^3$ and $256^3$. However, for larger workloads of $N_p^3 = 512^3$ and $1024^3$, the BGK model demonstrates slightly better performance than the MRT model while maintaining similar scaling trends. Although both the A100 and H100 architectures feature 80 GB of memory, the H100 can handle slightly larger LBM workloads due to its significantly higher memory bandwidth (HBM3 vs. HBM2e), improved memory compression, larger and more optimized L2 cache, and more efficient data movement strategies. These advancements reduce memory overhead and enhance computational efficiency, enabling the H100 to support larger grid sizes

compared to the A100. This advantage is reflected in the scaling starting point for 1024 grid points per direction, $N_{g0}^{1024}$, as shown in Table 2. The A100 architecture can handle up to 512 points per direction for both the MRT and BGK LBM models. While the H100 maintains the same capacity for the BGK model, it is more tolerant of memory-efficient models such as MRT, resulting in a lower starting point of $N_{g0}^{1024} = 4$ compared to the A100 and H100 with the BGK model starting point.

Figures 6 and 7 illustrate the speedup and strong scalability efficiency of the PyLB solver on the V100, A100, and H100 nodes of the Jean Zay supercomputer when using the BGK and MRT models, respectively. Table 4 summarizes the best strong-scaling performance in GLUPS and the corresponding $\eta^{\text{stg}}$ for each problem size, collision model, and GPU architecture. These performance indicators help guide users in selecting the appropriate GPU architecture and number of GPUs for a given problem size. Consistent with the previously presented GLUPS results, the PyLB solver demonstrates increasing speedup as more GPUs are allocated in parallel. However, the strong-scalability efficiency can drop significantly when inter-node communication becomes dominant and/or when GPU resources are not fully exploited.
<figure><img src="image_9.png" /><figcaption>**FIGURE 7** | Speedup and efficiency from the strong scalability study of the TGV test case with the LBM-MRT model using different problem sizes, $N_p = 128, 256, 512$, and $1024$, on the three GPU architectures available on the Jean Zay supercomputer. The violet (—), green (—), and yellow (—) colors represent the V100, A100, and H100 GPUs architectures, respectively. The solid line with the square marker (—) represents the speedup, while the dashed line with the circle marker (○) stands for the efficiency. The ○ marker highlights the V100's peak GLUPS result and the minimum number of A100 or H100 accelerators required to outperform it, while the △ marker indicates the largest $N_p$ that can be allocated while maintaining a strong scalability efficiency above 50%. (a) $N_p = 128$, (b) $N_p = 256$, (c) $N_p = 512$, (d) $N_p = 1024$.</figcaption></figure>

<table><caption>**TABLE 4** | Maximum GLUPS performance and corresponding strong-scaling efficiency for each GPU architecture and collision model.</caption><thead><tr><th rowspan="2">Arch.</th><th rowspan="2">Model</th><th colspan="3">Np = 128</th><th colspan="3">Np = 256</th><th colspan="3">Np = 512</th><th colspan="3">Np = 1024</th></tr><tr><th>GLUPS</th><th>N<sub>g</sub></th><th>η<sup>stg</sup> (%)</th><th>GLUPS</th><th>N<sub>g</sub></th><th>η<sup>stg</sup> (%)</th><th>GLUPS</th><th>N<sub>g</sub></th><th>η<sup>stg</sup> (%)</th><th>GLUPS</th><th>N<sub>g</sub></th><th>η<sup>stg</sup> (%)</th></tr></thead><tbody><tr><td rowspan="2">V100</td><td>BGK</td><td>1.31</td><td>64</td><td>3.98</td><td>5.29</td><td>128</td><td>6.75</td><td>11.38</td><td>128</td><td>15.86</td><td>22.12</td><td>128</td><td>49.27</td></tr><tr><td>MRT</td><td>1.31</td><td>64</td><td>3.93</td><td>5.22</td><td>128</td><td>6.91</td><td>10.42</td><td>128</td><td>16.28</td><td>17.49</td><td>128</td><td>56.73</td></tr><tr><td rowspan="2">A100</td><td>BGK</td><td>1.72</td><td>8</td><td>23.10</td><td>6.85</td><td>8</td><td>62.75</td><td>14.56</td><td>32</td><td>29.43</td><td>28.27</td><td>128</td><td>14.41</td></tr><tr><td>MRT</td><td>1.79</td><td>8</td><td>24.23</td><td>6.35</td><td>32</td><td>67.01</td><td>15.86</td><td>64</td><td>18.81</td><td>27.05</td><td>128</td><td>19.51</td></tr><tr><td rowspan="2">H100</td><td>BGK</td><td>2.17</td><td>32</td><td>5.09</td><td>15.25</td><td>128</td><td>6.21</td><td>80.57</td><td>128</td><td>27.71</td><td>184.63</td><td>128</td><td>61.55</td></tr><tr><td>MRT</td><td>2.22</td><td>32</td><td>6.00</td><td>15.49</td><td>128</td><td>7.61</td><td>67.00</td><td>128</td><td>28.58</td><td>138.57</td><td>128</td><td>55.77</td></tr></tbody></table>

The authors defined a minimum strong scalability efficiency threshold of $\eta^{stg} \ge 50\%$ for all evaluated $N_p$ values across the three GPU generations to ensure an effective operational point. Table 5 presents, for each GPU architecture and $N_p$, the maximum number of GPUs that can be allocated while maintaining a strong scalability efficiency above the 50% threshold in the TGV test case. These values are also highlighted in Figures 6 and 7

by the triangle markers. Data indicate that the LBM-BGK collision model running on the V100 architecture can only exceed the minimum threshold at $N_p^3 \ge 512^3$ when using multiple nodes of Jean Zay. The same communication latency is not observed for the LBM-MRT model, which can start using multiple V100 computing nodes for $N_p^3 \ge 256^3$ while maintaining $\eta^{stg} \ge 50\%$. When the code runs on the H100 nodes of Jean Zay, which
<table><caption>**TABLE 5** | Maximum $N_g$ that respects a minimum threshold of $\eta^{stg} \ge 50\%$ for both BGK and MRT models and three generations of NVIDIA GPUs.</caption><thead><tr><th>BGK</th><th colspan="2">V100</th><th colspan="2">A100</th><th colspan="2">H100</th></tr><tr><th>N<sub>p</sub></th><th>N<sub>g</sub></th><th>η<sup>stg</sup></th><th>N<sub>g</sub></th><th>η<sup>stg</sup></th><th>N<sub>g</sub></th><th>η<sup>stg</sup></th></tr></thead><tbody><tr><td>128</td><td>4</td><td>50%</td><td>2</td><td>68%</td><td>2</td><td>65%</td></tr><tr><td>256</td><td>4</td><td>94%</td><td>8</td><td>64%</td><td>8</td><td>53%</td></tr><tr><td>512</td><td>8</td><td>80%</td><td>8</td><td>94%</td><td>32</td><td>50%</td></tr><tr><td>1024</td><td>64</td><td>59%</td><td>16</td><td>65%</td><td>128</td><td>61%</td></tr></tbody><thead><tr><th>MRT</th><th colspan="2">V100</th><th colspan="2">A100</th><th colspan="2">H100</th></tr><tr><th>N<sub>p</sub></th><th>N<sub>g</sub></th><th>η<sup>stg</sup></th><th>N<sub>g</sub></th><th>η<sup>stg</sup></th><th>N<sub>g</sub></th><th>η<sup>stg</sup></th></tr></thead><tbody><tr><td>128</td><td>4</td><td>53%</td><td>2</td><td>68%</td><td>2</td><td>69%</td></tr><tr><td>256</td><td>8</td><td>65%</td><td>8</td><td>67%</td><td>8</td><td>64%</td></tr><tr><td>512</td><td>8</td><td>63%</td><td>8</td><td>95%</td><td>32</td><td>57%</td></tr><tr><td>1024</td><td>128</td><td>57%</td><td>32</td><td>60%</td><td>128</td><td>55%</td></tr></tbody></table>

consist of 4 GPUs per node as well, both collision models implemented in PyLB exceed the threshold when using multiple nodes for problem sizes larger than 256³ lattice points. For the A100 nodes of the IDRIS supercomputer, the PyLB solver only compensates for inter-node communication latency at the largest workload evaluated in this study, $N_p^3 = 1024^3$. This performance is attributed to the fact that the A100 nodes of Jean Zay are composed of 8 GPUs, allowing them to efficiently allocate significantly larger problems compared to the V100 and H100 nodes of Jean Zay. It is worth mentioning that the threshold defined here by the authors is merely a recommendation for efficient computing. Nevertheless, the PyLB solver can still achieve performance improvements even with lower strong scaling efficiency, as previously indicated in Figure 5.

## 7.2 | Weak Scaling

The weak scaling of the PyLB solver is evaluated for a streamwise periodic flow over a 3-D sphere, using the three generations of NVIDIA GPUs available on the Jean Zay supercomputer: Volta, Ampere, and Hopper. This evaluation is performed across different computational loads, denoted as $C_g$, to assess the ability of the solver to maintain steady performance as both the problem size in the streamwise direction, $N_x$, and the number of GPUs, $N_g$, are increased proportionally.

Unlike the strong scaling analysis in the previous section, weak scalability involves scaling only the number of points in the streamwise direction, $N_x$, proportionally to $N_g$ while keeping the computational load per GPU, $C_g$, constant. The computational load is defined as:

$$C_g = \frac{N_x \times N_y \times N_z}{N_g} \qquad (16)$$

The initial domain in this study consists of a single 3-D wall-bounded sphere with an equal number of points in each direction, $N_x = N_y = N_z$. A periodic boundary condition is

imposed in the streamwise direction. As the number of GPUs increases proportionally, the domain is extended in the x Cartesian direction by concatenating identical subdomains, each containing a wall-bounded sphere, to maintain a constant computational load per GPU.

The evaluated computational loads are 64³, 128³, 256³, and 512³ points per GPU. The number of GPUs used ranges from 1 to 128. Similar to the strong scalability tests, calculations are executed for 1000 iterations to measure the total computing time, T. The weak scaling efficiency, for a given computational load $C_g$ and number of GPUs $N_g$, is calculated as:

$$\eta_{N_g}^{C_g, \text{weak}} = \frac{T_0^{C_g}}{T_{N_g}^{C_g}}, \qquad (17)$$

where $T_0^{C_g}$ and $T_{N_g}^{C_g}$ represent, for a given $C_g$ computing charge, the computing time required to perform 1000 iterations using one and $N_g$ GPUs, respectively. Unlike from the strong scaling study, the starting point for the weak scalability is defined as $N_g = 1$.

Figures 8 and 9 present the weak scaling efficiency, $\eta^{\text{weak}}$, and the computational performance measured in GLUPS for the PyLB solver using the BGK and MRT collision models, respectively. Both collision models exhibit similar trends when increasing the problem size proportionally to the number of GPUs. The code demonstrates scalability across the three computational loads and GPU architectures evaluated, while the weak scalability efficiency gradually decreases, eventually converging to a minimum value. This minimal point can reach high efficiency values for large $C_g$. This stabilization suggests that computational costs associated with increasing the problem size, $C_g$, such as network communication overhead, synchronization bottlenecks, cache inefficiencies, and memory limitations, become significant beyond a certain scale. Consequently, a saturation point is reached, beyond which additional GPUs no longer yield linear performance improvements.

Table 6 presents the saturation point, $N_g$, at which $\eta^{\text{weak}}$ converges to its minimum value during the weak scaling test on the V100, A100, and H100 GPUs for both collision models. The LBM solver exhibits a high efficiency for $C_g \ge 256^3$, where $70 \le \eta^{\text{weak}} \le 95\%$. For a computing charge of $256^3$, both BGK and MRT collision models on the V100 and A100 nodes show a weak scaling efficiency stagnation around 80%, occurring at $N_g = 16$ and $N_g = 4$, respectively. Under the same conditions, the Jean Zay nodes equipped with H100 GPUs reach a saturation efficiency of 70% at $N_g = 4$.

The case of $C_g = 512^3$ is evaluated only for the A100 and H100 architectures, as the V100 nodes cannot allocate such a problem size within the PyLB solver. For this computing charge, the A100 architecture shoes shows weak scaling efficiency stagnating at 85% for the BGK model, for $N_g = 2$, and at 75% for the MRT model, for $N_g = 16$. Meanwhile, the H100 architecture exhibits weak scaling efficiency stagnation of 95% at $N_g = 4$ for the BGK model and at $N_g = 16$ for the MRT model.

The weak scaling efficiency of PyLB for $C_g = 128^3$ on the V100 and A100 GPUs, for both collision models, is lower, $\approx 50\%$,
<figure><img src="image_10.png" /><figcaption>**FIGURE 8** | Weak scalability efficiency, $\eta^{\text{weak}}$, and GLUPS for the flow over a sphere test case with the LBM-BGK model using different computing charges, $C_g = 128^3$, $256^3$, and $512^3$ on the three GPU architectures available on the Jean Zay supercomputer. The violet (—), green (—), and yellow (—) colors represent the V100, A100, and H100 GPUs architectures, respectively. The solid line with the square marker (■) represents the performance in GLUPS, while the dashed line with the circle marker (●) stands for the weak scaling efficiency. (a) $C_g = 128^3$, (b) $C_g = 256^3$, (c) $C_g = 512^3$.</figcaption></figure>

compared to the efficiency observed for higher computing charges. The $\eta^{\text{weak}}$ values for the same conditions are even lower when running on the H100 nodes of Jean Zay: 30% for the BGK model at $N_g = 4$ and 40% for the MRT model at $N_g = 4$. For this lower computational charge, the overhead associated with increased computing resources outweighs the computational gains, leading to a rapid deterioration in performance.

## 7.3 | Comparison With Existing Scientific LBM Solvers

To contextualize the performance of PyLB relative to other established GPU-based LBM solvers, the present work compares its single-GPU throughput against three scientific frameworks available in the literature: Sailfish [34], Palabos [32], and waLBerla [98]. Since these solvers target different hardware generations, implement distinct collision models, and employ a variety of kernel-optimization and memory-layout strategies, the goal of this comparison is not to establish a strict ranking, but rather

to position PyLB within the broader landscape of GPU LBM implementations.

For the NVIDIA V100 GPU, PyLB achieves a performance of ≈ 0.6 GLUPS for a $256^3$ Taylor-Green vortex test case. On the more recent A100 GPU, PyLB reaches ≈ 1.4 GLUPS and ≈ 1.5 GLUPS for $256^3$ and $512^3$ grids, respectively. On the state-of-the-art NVIDIA H100 device, the performance increases to ≈ 1.9 GLUPS for $256^3$ and ≈ 2.3 GLUPS for $512^3$. These values can be compared, in terms of order of magnitude, with results reported in the literature. Sailfish reports approximately 0.8 GLUPS on an NVIDIA K40 GPU for a $254^3$ lattice [34]. Although this involves an older GPU generation and a different benchmark configuration, it illustrates that PyLB's performance on modern hardware is consistent with expectations for LBM solvers employing similar computational kernels. More recent comparisons are available for the NVIDIA GPUs: waLBerla reports ≈ 1.0 GLUPS on the V100 architecture when using the D3Q19 stencil while Palabos obtains ≈ 3.6 GLUPS for a $256^3$ grid [32] on the A100 accelerator. Although these values are higher than the performance
<figure><img src="image_11.png" /><figcaption>**FIGURE 9** | Weak scalability efficiency, $\eta^{\text{weak}}$, and GLUPS for the flow over a sphere test case with the LBM-MRT model using different computing charges, $C_g = 128^3$, $256^3$, and $512^3$ on the three GPU architectures available on the Jean Zay supercomputer. The violet (—), green (—), and yellow (—) colors represent the V100, A100, and H100 GPUs architectures, respectively. The solid line with the square marker (■—) represents the speedup, while the dashed line with the circle marker (○—) stands for the efficiency. (a) $C_g = 128^3$, (b) $C_g = 256^3$, (c) $C_g = 512^3$.</figcaption></figure>

<table><caption>**TABLE 6** | Number of GPUs, $N_g$, at which the PyLB $\eta^{\text{weak}}$ converges to its minimum value for both BGK and MRT models across the three NVIDIA GPU architectures evaluated in the weak scaling test.</caption><thead><tr><th rowspan="2">BGK</th><th colspan="2">V100</th><th colspan="2">A100</th><th colspan="2">H100</th></tr><tr><th>N<sub>g</sub></th><th>η<sup>weak</sup></th><th>N<sub>g</sub></th><th>η<sup>weak</sup></th><th>N<sub>g</sub></th><th>η<sup>weak</sup></th></tr></thead><tbody><tr><td>C<sub>g</sub></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>128<sup>3</sup></td><td>8</td><td>50%</td><td>32</td><td>55%</td><td>4</td><td>30%</td></tr><tr><td>256<sup>3</sup></td><td>8</td><td>80%</td><td>16</td><td>80%</td><td>4</td><td>70%</td></tr><tr><td>512<sup>3</sup></td><td>—</td><td>—</td><td>2</td><td>85%</td><td>4</td><td>95%</td></tr></tbody></table><table><thead><tr><th rowspan="2">MRT</th><th colspan="2">V100</th><th colspan="2">A100</th><th colspan="2">H100</th></tr><tr><th>N<sub>g</sub></th><th>η<sup>weak</sup></th><th>N<sub>g</sub></th><th>η<sup>weak</sup></th><th>N<sub>g</sub></th><th>η<sup>weak</sup></th></tr></thead><tbody><tr><td>C<sub>g</sub></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>128<sup>3</sup></td><td>8</td><td>60%</td><td>16</td><td>55%</td><td>4</td><td>40%</td></tr><tr><td>256<sup>3</sup></td><td>8</td><td>80%</td><td>16</td><td>80%</td><td>4</td><td>70%</td></tr><tr><td>512<sup>3</sup></td><td>—</td><td>—</td><td>16</td><td>75%</td><td>16</td><td>95%</td></tr></tbody></table>

achieved by PyLB on the V100 and A100 GPUs, they remain within the same order of magnitude. The discrepancy can reasonably be attributed to several factors, including differences in the underlying test cases, algorithmic choices, collision models, data structures, and the varying degree of GPU-specific optimizations applied in each framework.

In particular, Palabos and waLBerla implement highly specialized data layouts, non-trivial fusion strategies, architecture-tailored kernels, and extensive tuning for specific GPU generations. PyLB, in contrast, adopts a portability-oriented design through OpenACC, maintaining a unified codebase for CPU and GPU execution. While this design choice naturally limits the theoretical peak performance on a single GPU, it greatly simplifies multi-architecture development, testing, and maintenance.

Overall, PyLB demonstrates a performance level that is consistent with the expected behavior of state-of-the-art LBM solvers when
<figure><img src="image_12.png" /><figcaption>FIGURE 10 | A CPU/GPU performance comparison for the TGV test case with 512³ grid points.</figcaption></figure>

considering differences in hardware, benchmark setup, and optimization philosophy. These results indicate that PyLB achieves a competitive throughput while preserving a design that favors simplicity, portability, and long-term maintainability.

## 8 | CPU/GPU Power Efficiency Comparison

In the current section, the LBM solver is executed on both the GPU and CPU nodes of the Jean Zay supercomputer to compare its performance, measured in GLUPS, as well as its energy consumption (in kWh) and CO₂ equivalent emissions (kg CO₂e) across different hardware architectures. The carbon emissions are reported based on data from the French Electricity Transmission Network (RTE) [99, 100]. These values account for the varying energy sources contributing to the French national grid and are only used here for comparison. The Taylor-Green vortex (TGV) test case, with 512³ and 1024³ grid resolutions, is used to assess the strong scalability of PyLB-OpenACC on the Jean Zay system. The code is executed on 500, 4000, and 8000 CPU cores across nodes equipped with Intel Cascade Lake 6248 processors. Its performance is compared with calculations using 1, 4, and 8 GPUs from the accelerated architectures available in the Jean Zay supercomputer.

Figure 10 shows the PyLB performance in GLUPS, illustrating scalability as the number of GPUs or CPU cores increases, for the TGV flow configuration using 512³ mesh points. The code demonstrates near-linear scaling across all evaluated configurations. The scalability trends of PyLB-OpenACC on H100 GPUs and CSL CPUs are comparable: 1000 CPU cores and 1 H100 GPU each deliver approximately ≈ 2.5 GLUPS, 4000 CPU cores and 4 H100 GPUs reach ≈ 9 GLUPS, and 8000 CPU cores and 8 H100 GPUs achieve ≈ 16 GLUPS. The A100 and V100 GPU architectures yield ≈ 1 GLUPS per GPU. With 4 GPUs, the A100 configuration reaches ≈ 6 GLUPS and the V100 configuration reaches ≈ 2.5 GLUPS. With 8 GPUs, the A100 and V100 achieve ≈ 11 and ≈ 4 GLUPS, respectively. To reach equivalent performance levels, approximately 1500 CPU cores are required to match the ≈ 4 GLUPS achieved by 8 V100 GPUs, while around 5500 Cascade

Lake cores are necessary to match the ≈ 11 GLUPS delivered by 8 A100 GPUs.

The Jean Zay system employs the CEEMS toolkit [101, 102] to monitor both the energy consumption, in kWh, and the total equivalent carbon dioxide emissions, in kg CO₂e, associated with each SLURM job. These metrics enable a quantitative comparison of the energy efficiency and carbon footprint between GPU- and CPU-based architectures. To this end, the present study performs simulations of the three-dimensional Taylor-Green vortex (TGV) configuration on a 1024³ grid, integrated up to a non-dimensional time of t* = 20 (or Nᵢᵢ = 67456). This flow problem is first executed on 8000 CPU cores distributed across Jean Zay nodes, yielding a performance of ≈ 15 GLUPS. The same test case is then computed using 64 V100 GPUs, 16 A100 GPUs, and 8 H100 GPUs, with each configuration selected to deliver a comparable performance of ≈ 15 GLUPS, as shown in Figure 5d.

The CEEMS data collected from the TGV test cases performed on the Jean Zay supercomputer are normalized by the GLUPS performance metric (see Equation 14) to define energy and carbon emission efficiency coefficients. Since GLUPS accounts for both the number of lattice points and iterations per second, these normalized coefficients, denoted by α, enable a quantitative comparison across architectures. The coefficients are defined as:

$$ 
\begin{aligned}
\alpha_{\text{energy}} \text{ [Wh/LU]} &= \frac{E \text{ [Wh]}}{N_x \times N_y \times N_z \times N_{it}} \\
&= \frac{E \text{ [Wh]}}{\text{GLUPS}_{N_g}^{N_p} \times 10^9 \text{ [LU/s]} \times T \text{ [s]}}, \quad (18)
\end{aligned}
$$

$$ 
\begin{aligned}
\alpha_{\text{emission}} \text{ [gCO}_2\text{e/LU}] &= \frac{e_{\text{CO}_2} \text{ [gCO}_2\text{e]}}{N_x \times N_y \times N_z \times N_{it}} \\
&= \frac{e_{\text{CO}_2} \text{ [gCO}_2\text{e]}}{\text{GLUPS}_{N_g}^{N_p} \times 10^9 \text{ [LU/s]} \times T \text{ [s]}}. \quad (19)
\end{aligned}
$$

Table 7 reports the computational performance, total wall-clock time (T), energy consumption (E), equivalent CO₂ emissions
<table><caption>**TABLE 7** | Computational performance, energy consumption, CO₂ emissions (RTE), and α coefficients, for CPU and GPU jobs running the 1024³ points TGV test case.</caption><thead><tr><th rowspan="2">Architecture</th><th rowspan="2">Resources (N<sub>g</sub>)</th><th rowspan="2">GLUPS<sub>N<sub>p</sub>=1024</sub></th><th rowspan="2">T (s)</th><th colspan="2">E (kWh)</th><th colspan="2">e<sub>CO2</sub> (g CO<sub>2</sub>e)</th><th rowspan="2">α<sub>energy</sub> (Wh/LU)</th><th rowspan="2">α<sub>emission</sub> (g CO<sub>2</sub>/LU)</th></tr><tr><th>Host</th><th>GPU</th><th>Host</th><th>GPU</th></tr></thead><tbody><tr><td>CLS 6248</td><td>8000 cores</td><td>16.840</td><td>4506</td><td>80.4</td><td>—</td><td>1045</td><td>—</td><td>1.06 ×10<sup>-6</sup></td><td>13.772 ×10<sup>-12</sup></td></tr><tr><td>V100</td><td>64 GPUs</td><td>12.880</td><td>5537</td><td>3.56</td><td>12.5</td><td>45.2</td><td>158</td><td>0.225 ×10<sup>-6</sup></td><td>2.850 ×10<sup>-12</sup></td></tr><tr><td>A100</td><td>16 GPUs</td><td>15.453</td><td>4939</td><td>1.85</td><td>6.63</td><td>31.3</td><td>112</td><td>0.111 ×10<sup>-6</sup></td><td>1.878 ×10<sup>-12</sup></td></tr><tr><td>H100</td><td>8 GPUs</td><td>17.670</td><td>4444</td><td>1.21</td><td>5.57</td><td>13.3</td><td>61.3</td><td>0.0863 ×10<sup>-6</sup></td><td>0.950 ×10<sup>-12</sup></td></tr></tbody></table>

(eCO₂, and the corresponding energy and carbon emission coefficients, α energy and α emission, for the CPU and GPU nodes evaluated on Jean Zay.

As shown in Figure 10, Table 7 also demonstrates that GPU-based configurations require fewer computational resources to achieve comparable performance levels relative to CPU-based systems even for different computing charges, 512³ and 1024³ points, respectively. However, direct comparisons of resource counts are not straightforward, as CPU and GPU architectures are fundamentally different and are designed to achieve high computational throughput via distinct architectural paradigms. In contrast, metrics such as energy consumption and CO₂ emissions provide a meaningful basis for evaluating architectural efficiency, as computing centers increasingly prioritize cost-effectiveness and environmental sustainability for a given level of computational performance. In this study, all GPU configurations exhibit significantly lower energy consumption and markedly reduced CO₂ emissions compared to the CPU-based configuration. Furthermore, among the evaluated GPUs, newer architectures demonstrate improved energy and carbon efficiency, highlighting ongoing advancements in accelerator technology. Quantitatively, for the TGV test case with a 1024³ grid, the Cascade Lake architecture requires approximately 5, 10, and 12 times more energy per lattice update than the V100, A100, and H100 architectures, respectively. In terms of CO₂ emissions, the V100, A100, and H100 configurations show approximately 5, 7, and 14 times lower emission levels, respectively, compared to the CSL CPUs.

These observed differences are rooted in the fundamental architectural design of CPUs and GPUs in high-performance computing. CPUs typically feature a smaller number of highly complex cores, optimized for sequential tasks and intricate instruction scheduling, which limits their parallel throughput. In contrast, GPUs consist of many simpler cores designed for data-parallel execution, offering significantly higher throughput for massively parallel workloads, as the LBM algorithm proves to be. This architectural distinction allows GPUs to achieve higher computational performance with greater energy efficiency and lower carbon emissions, particularly for workloads that benefit from extensive parallelism, as reflected in the present results.

9 | Conclusion

Porting of an existing research code to a specific and promising architecture can require significant development time to meet architecture's expectations. The current work presented

the porting and performance evaluation of the PyLB solver, a Lattice Boltzmann Method code originally developed using Python, Mpi4Py, and Fortran, onto GPU-accelerated platforms. The migration leveraged the OpenACC programming paradigm, enabling GPU support in a minimally intrusive manner. The solver architecture was reorganized into Python scripts that manage high-level operations and interface with compiled Fortran-MPI-OpenACC binaries responsible for the computational core. The validity of the GPU-accelerated implementation was assessed through two benchmark cases: the 3-D decaying Taylor-Green vortex and the flow past a 3-D sphere. Results showed good agreement with both reference numerical and experimental data.

To improve performance on GPU architectures, several targeted optimizations were implemented. These included loop restructuring, scalar substitution, memory coalescence, and the elimination of data dependencies to promote compiler vectorization and efficient memory access. A key modification involved reducing the number of distribution function buffers used during the LBM collision and propagation routines from four to two, significantly simplifying data movement and improving memory efficiency. These enhancements translated into notable performance gains, as confirmed through a detailed Roofline analysis on V100, A100, and H100 architectures. The study showed that the optimized implementation substantially increased both arithmetic intensity and computational throughput, particularly on the V100 and A100, and moved the solver closer to the performance ceilings imposed by memory bandwidth and compute capability.

Comprehensive strong and weak scalability tests were carried out on the Jean Zay supercomputer using accelerated architectures. The PyLB solver demonstrated consistent scalability across all generations of accelerators, with particularly favorable results on the H100 platform. Strong scaling efficiency remained above 50% for a wide range of domain sizes, with scalability observed up to 128 GPUs. In weak scaling tests, efficiency converged to stable values as the number of GPUs increased, particularly for larger workloads ($C_g \ge 256^3$), confirming the robustness of the solver. For the TGV test case with a 512³ grid, 8 H100 GPUs achieved approximately 16 GLUPS, the same performance obtained using 8000 CPU cores, demonstrating a 1000:1 resource equivalence ratio.

In addition to raw performance, the study also examined the energy consumption and carbon footprint associated with large-scale simulations. GPU-based configurations significantly outperformed CPU-based setups in both energy efficiency and
CO₂ emissions, underlining the environmental benefits of using modern accelerators for computational fluid dynamics applications based on the LBM formulation.

Overall, this work demonstrates that the OpenACC-enabled PyLB solver is an efficient, portable, and sustainable tool for high-performance LBM simulations on heterogeneous HPC platforms. It also serves as a roadmap for future PyLB users and developers interested in porting HPC applications to accelerator architectures using the OpenACC paradigm.

## Acknowledgments

This work was performed using HPC resources from IDRIS-GENCI (Grant A0172A10636, A0152A10636, A0152A12067, and A0172A12067) and from the IDRIS-NVIDIA Hackathon 2023 (Grant A0132B10947). The authors also acknowledge the ANR-FRANCE (French National Research Agency) for its financial support of the project JENII (ANR-21-DMES-0006).

## Disclosure

The authors have nothing to report.

## Conflicts of Interest

The authors declare no conflicts of interest.

## Data Availability Statement

The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.

## References

1. CEA, "Joliot-Curie, One of the Most Powerful Supercomputer Dedicated to French and European Research," 2025, https://www-hpc.cea.fr/en/Joliot-Curie.html.

2. IDRIS, "Jean Zay: HPE SGI 8600 Supercomputer," 2025, http://www.idris.fr/eng/jean-zay/cpu/jean-zay-cpu-hw-eng.html.

3. CINES, "Adastra Architecture," 2025, https://dci.dci-gitlab.cines.fr/webextranet/architecture/index.html.

4. Top 500, Top 500, The list, 2025, https://www.top500.org/lists/top500/.

5. A. Abdelfattah, V. Barra, N. Beams, et al., "GPU Algorithms for Efficient Exascale Discretizations," *Parallel Computing* 108 (2021): 102841.

6. P. Fischer, S. Kerkemeier, M. Min, et al., "NekRS, a GPU-Accelerated Spectral Element Navier-Stokes Solver," *Parallel Computing* 114 (2022): 102982.

7. R. Farber, *Parallel Programming With OpenACC* (Morgan Kaufmann, 2017).

8. W. Xue and C. J. Roy, "Multi-GPU Performance Optimization of a Computational Fluid Dynamics Code Using OpenACC," *Concurrency and Computation: Practice and Experience* 33, no. 5 (2021): e6036.

9. A. Xu and B.-T. Li, "Multi-GPU Thermal Lattice Boltzmann Simulations Using OpenACC and MPI," *International Journal of Heat and Mass Transfer* 201 (2023): 123649.

10. Y.-H. Qian, D. d'Humières, and P. Lallemand, "Lattice BGK Models for Navier-Stokes Equation," *Europhysics Letters* 17 (1992): 479-484.

11. X. He and L. S. Luo, "A Priori Derivation of the Lattice Boltzmann Equation," *Physical Review E* 55 (1997): R6333.

12. A. Kummerländer, M. Dorn, M. Frank, and M. J. Krause, "Implicit Propagation of Directly Addressed Grids in Lattice Boltzmann Methods," *Concurrency and Computation: Practice and Experience* 35, no. 8 (2023): e7509.

13. X. Shan, X. Yuan, and H. Chen, "Kinetic Theory Representation of Hydrodynamics: A Way Beyond the Navier-Stokes Equation," *Journal of Fluid Mechanics* 550 (2006): 413-441.

14. L. Zhong, S. Feng, P. Dong, and S. Gao, "Lattice Boltzmann Schemes for the Nonlinear Schrödinger Equation," *Physical Review E* 74 (2006): 036704.

15. P. Romatschke, M. Mendoza, and S. Succi, "Fully Relativistic Lattice Boltzmann Algorithm," *Physical Review C* 84 (2011): 034903.

16. P. L. Bhatnagar, E. P. Gross, and M. Krook, "A Model for Collision Process in Gases. I. Small Amplitude Process in Charged and Neutral One-Component Systems," *Physics Review* 94, no. 3 (1954): 511-525.

17. D. d'Humière, I. Ginzburg, Y. Krafczyk, P. Lallemand, and L. S. Luo, "Multiple Relaxation Time Lattice Boltzmann Models in Three Dimensions," *Philosophical Transactions of the Royal Society of London. Series A* 360 (2002): 437.

18. S. Marié, D. Ricot, and P. Sagaut, "Comparison Between Lattice Boltzmann Method and Navier-Stokes High Order Schemes for Computational Aeroacoustics," *Journal of Computational Physics* 228 (2009): 1056-1070.

19. Y. Peng, W. Liao, L.-S. Luo, and L.-P. Wang, "Comparison of the Lattice Boltzmann and Pseudo-Spectral Methods for Decaying Turbulence: Low-Order Statistics," *Computers & Fluids* 39, no. 4 (2010): 568-591.

20. J. D. Sterling and S. Chen, "Stability Analysis of Lattice Boltzmann Methods," *Journal of Computational Physics* 123 (1996): 196-206.

21. C. David and P. Sagaut, "Structural Stability of Lattice Boltzmann Schemes," *Physica A* 444 (2016): 1-8.

22. C. Coreixas, G. Wissocq, B. Chopard, and J. Latt, "Impact of Collision Models on the Physical Properties and the Stability of Lattice Boltzmann Methods," *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 378, no. 2175 (2020): 20190397.

23. M. Geier, M. Sch nherr, A. Pasquali, and M. Krafczyk, "The Cumulant Lattice Boltzmann Equation in Three Dimensions: Theory and Validation," *Computers and Mathematics With Applications* 70, no. 4 (2015): 507-547.

24. J. Latt and B. Chopard, "Lattice Boltzmann Method With Regularized Pre-Collision Distribution Functions," *Mathematics and Computers in Simulation* 72, no. 2 (2006): 165-168.

25. F. Tosi, S. Ubertini, S. Succi, H. Chen, and I. V. Karlin, "Numerical Stability of Entropic Versus Positivity-Enforcing Lattice Boltzmann Schemes," *Mathematics and Computers in Simulation* 72 (2006): 227-231.

26. P. Lallemand and F. Dubois, "Some Results on Energy-Conserving Lattice Boltzmann Models," *Computers and Mathematics With Applications* 65 (2013): 831-844.

27. F. Dubois, B. Graille, and S. V. R. Rao, "A Notion of Non-Negativity Preserving Relaxation for a Mono-Dimensional Three Velocities Scheme With Relative Velocity," *Journal of Computational Science* 47 (2020): 101181.

28. C. Coreixas, B. Chopard, and J. Latt, "Comprehensive Comparison of Collision Models in the Lattice Boltzmann Framework: Theoretical Investigations," *Physical Review E* 100 (2019): 034901.

29. S. Marié and X. Gloerfelt, "Adaptive Filtering for the Lattice Boltzmann Method," *Journal of Computational Physics* 333C (2017): 212-226.

30. L. Vienne, S. Marié, and F. Grasso, "Lattice Boltzmann Method for Miscible Gases: A Forcing-Term Approach," *Physical Review E* 100 (2019): 034902.
31. L. Vienne and S. Marié, "Lattice Boltzmann Study of Miscible Viscous Fingering for Binary and Ternary Mixtures," *Physical Review Fluids* 6 (2021): 053904.

32. J. Latt, O. Malaspinas, D. Kontaxakis, et al., "Palabos: Parallel Lattice Boltzmann Solver," *Computers & Mathematics With Applications* 81 (2021): 334–350.

33. S. Kemmler, C. Rettinger, U. R. de, P. Cuéllar, and H. Köstler, "Efficiency and Scalability of Fully-Resolved Fluid-Particle Simulations on Heterogeneous CPU-GPU Architectures," *International Journal of High Performance Computing Applications* 39, no. 3 (2025): 345–363.

34. M. Januszewski and M. Kostur, "Sailfish: A Flexible Multi-GPU Implementation of the Lattice Boltzmann Method," *Computer Physics Communications* 185, no. 9 (2014): 2350–2368.

35. J. Nickolls, I. Buck, M. Garland, and K. Skadron, "Scalable Parallel Programming With CUDA," in *ACM SIGGRAPH 2008 Classes* (ACM New York, 2008).

36. R. P. Broussard and R. W. Ives, "Parallel Processing for Imaging Applications," in *International Society for Optics and Photonics*, ed. D. Owens John, L. I-Jong, Z. Yu-Jin, and B. Beretta Giordano (SPIE, 2011), 787202.

37. D. B. Kirk and W. H. Wen-Mei, *Programming Massively Parallel Processors: A Hands-On Approach* (Morgan kaufmann, 2016).

38. J. E. Stone, D. Gohara, and G. Shi, "OpenCL: A Parallel Programming Standard for Heterogeneous Computing Systems," *Computing in Science & Engineering* 12, no. 3 (2010): 66–73.

39. S. Chandrasekaran and G. Juckeland, *OpenAC for Programmers: Concepts and Strategies* (Addison-Wesley Professional, 2017).

40. S. Bak, C. Bertoni, S. Boehm, et al., "OpenMP Application Experiences: Porting to Accelerated Nodes," *Parallel Computing* 109 (2022): 102856.

41. J. Huber, M. Cornelius, G. Georgakoudis, et al., "Efficient Execution of OpenMP on GPUs," in *2022 IEEE/ACM International Symposium on Code Generation and Optimization (CGO)* (IEEE, 2022), 41–52.

42. G. Wellein, T. Zeiser, G. Hager, and S. Donath, "On the Single Processor Performance of Simple Lattice Boltzmann Kernels," *Computers & Fluids* 35, no. 8 9 (2006): 910–919.

43. D. Li, C. Xu, Y. Wang, et al., "Parallelizing and Optimizing Large-Scale 3D Multi-Phase Flow Simulations on the Tianhe-2 Supercomputer," *Concurrency and Computation: Practice and Experience* 28, no. 5 (2015): 1678–1692.

44. P. Valera-Lara, "Reducing Memory Requirements for Large Size LBM Simulations on GPUs," *Concurrency and Computation: Practice and Experience* 29, no. 24 (2017): e4221.

45. E. Riegel, T. Indinger, and N. A. Adams, "Implementation of a Lattice-Boltzmann Method for Numerical Fluid Mechanics Using the nVIDIA CUDA Technology," *Computer Science - Research and Development* 23, no. 3 4 (2009): 241–247.

46. F. Kuznik, C. Obrecht, G. Rusaouen, and J.-J. Roux, "LBM Based Flow Simulation Using GPU Computing Processor," *Computers & Mathematics With Applications* 59, no. 7 (2010): 2380–2392.

47. C. Obrecht, F. Kuznik, B. Tourancheau, and J.-J. Roux, "Multi-GPU Implementation of the Lattice Boltzmann Method," *Computers & Mathematics With Applications* 65, no. 2 (2013): 252–261.

48. P. Valero-Lara, J. Vetter, J. Gounley, and A. Randles, "Moment Representation of Regularized Lattice Boltzmann Methods on NVIDIA and AMD GPUs," in *Proceedings of the SC 23 Workshops of the International Conference on High Performance Computing, Network, Storage, and Analysis* (ACM, 2023), 1697–1704.

49. K. Sano, "FPGA-Based Scalable Custom Computing Accelerator for Computational Fluid Dynamics Based on Lattice Boltzmann Method," in *Sustained Simulation Performance 2014*, ed. M. M. Resch, W. Bez, E. Focht, H. Kobayashi, and N. Patel (Springer International Publishing, 2015), 187–201.

50. W. Altoyan and J. J. Alonso, "Accelerating the Lattice Boltzmann Method," in *2023 IEEE Aerospace Conference* (IEEE, 2023), 1–20.

51. H. Schottenhamml, A. A. Sedrakian, F. Blondel, H. Köstler, and U. Rüde, "waLBerla-Wind: A Lattice-Boltzmann-Based High-Performance Flow Solver for Wind Energy Applications," *Concurrency and Computation: Practice and Experience* 36, no. 16 (2024): e8117.

52. P. Lallemand and L. S. Luo, "Theory of the Lattice Boltzmann Method: Dispersion, Dissipation, Isotropy, Galilean Invariance and Stability," *Physical Review E* 61, no. 6 (2000): 6546–6562.

53. T. Février, "PhD Thesis, Univeristé de Paris-Saclay," 2015.

54. O. Malaspinas, "Increasing Stability and Accuracy of the Lattice Boltzmann Scheme: Recursivity and Regularization," *arXiv* e-prints, *arXiv:1505.06900* (2015).

55. Y. Feng, P. Sagaut, and W.-Q. Tao, "A Compressible Lattice Boltzmann Finite Volume Model for High Subsonic and Transonic Flows on Regular Lattices," *Computers & Fluids* 131 (2016): 45–55.

56. A. De Rosis and K. H. Luo, "Role of Higher-Order Hermite Polynomials in the Central-Moments-Based Lattice Boltzmann Framework," *Physical Review E* 99 (2019): 016901.

57. J. Jacob, O. Malaspinas, and P. Sagaut, "A New Hybrid Recursive Regularised Bhatnagar–Gross–Krook Collision Model for Lattice Boltzmann Method-Based Large Eddy Simulation," *Journal of Turbulence* 19, no. 11–12 (2018): 1051–1076.

58. F. Dubois and P. Lallemand, "Towards Higher Order Lattice Boltzmann Schemes," *Journal of Statistical Mechanics* 2009 (2009): P06006.

59. S. Guo, Y. Feng, J. Jacob, F. Renard, and P. Sagaut, "An Efficient Lattice Boltzmann Method for Compressible Aerodynamics on D3Q19 Lattice," *Journal of Computational Physics* 418 (2020): 109570.

60. L. Dalcín, R. Paz, and M. Storti, "MPI for Python," *Journal of Parallel and Distributed Computing* 65, no. 9 (2005): 1108–1115.

61. M. Rogowski, S. Aseeri, D. Keyes, and L. Dalcin, "mpi4py.Futures: MPI-Based Asynchronous Task Execution for Python," *IEEE Transactions on Parallel and Distributed Systems* 34, no. 2 (2023): 611–622.

62. Mpy4Py, "MPI for Python Tutorial," 2025, https://mpi4py.readthedocs.io/en/4.0.3/tutorial.html.

63. M. Haussmann, S. Simonis, H. Nirschl, and M. J. Krause, "Direct Numerical Simulation of Decaying Homogeneous Isotropic Turbulence - Numerical Experiments on Stability, Consistency and Accuracy of Distinct Lattice Boltzmann Methods," *International Journal of Modern Physics C* 30, no. 9 (2019): 1950074.

64. C. Mimeau, S. Marié, and I. Mortazavi, "A Comparison of Semi-Lagrangian Vortex Method and Lattice Boltzmann Method for Incompressible Flows," *Computers and Fluids* 224 (2021): 104946.

65. A. Suss, I. Mary, T. Le Garrec, and S. Marié, "Comprehensive Comparison Between the Lattice Boltzmann and Navier–Stokes Methods for Aerodynamic and Aeroacoustic Applications," *Computers & Fluids* 257 (2023): 105881.

66. P. Skordos, "Initial and Boundary Conditions for the Lattice Boltzmann Method," *Physical Review E* 48, no. 6 (1993): 4823–4842.

67. R. Mei, L.-S. Luo, P. Lallemand, and D. d Humières, "Consistent Initial Conditions for Lattice Boltzmann Simulations," *Computers & Fluids* 35, no. 8 (2006): 855–862.

68. Z. J. Wang, K. Fidkowski, R. Abgrall, et al., "High-Order CFD Methods: Current Status and Perspective," *International Journal for Numerical Methods in Fluids* 72, no. 8 (2013): 811–845.
69. K. Hu, J. Meng, H. Zhang, X.-J. Gu, D. R. Emerson, and Y. Zhang, "Ninth International Conference on Computational Fluid Dynamics (ICCFD9)," *Computers & Fluids* 156 (2017): 1-8.

70. J. Jeong and F. Hussain, "On the Identification of a Vortex," *Journal of Fluid Mechanics* 285 (1995): 69-94.

71. E. Achenbach, "Vortex Shedding From Spheres," *Journal of Fluid Mechanics* 62, no. 2 (1974): 209-221.

72. T. A. Johnson and V. C. Patel, "Flow Past a Sphere up to a Reynolds Number of 300," *Journal of Fluid Mechanics* 378 (1999): 19-70.

73. J. D. Anderson, *Fundamentals of Aerodynamics* (McGraw-Hill, 2011).

74. F. W. Roos and W. W. Willmarth, "Some Experimental Results on Sphere and Disk Drag," *AIAA Journal* 9, no. 2 (1971): 285-291.

75. A. G. Tomboulides and S. A. Orszag, "Numerical Investigation of Transitional and Weak Turbulent Flow Past a Sphere," *Journal of Fluid Mechanics* 416 (2000): 45-73.

76. G. S. Constantinescu and K. Squires, "LES and DES Investigations of Turbulent Flow Over a Sphere at $\bar{\epsilon} = 10,000$," *Flow, Turbulence and Combustion* 70, no. 1-4 (2003): 267-298.

77. D. Kim and H. Choi, "Laminar Flow Past a Sphere Rotating in the Streamwise Direction," *Journal of Fluid Mechanics* 461 (2002): 365-386.

78. P. Ploumhans, G. S. Winckelmans, J. K. Salmon, A. Leonard, and M. S. Warren, "Vortex Methods for Direct Numerical Simulation of Three-Dimensional Bluff Body Flows: Application to the Sphere at $\bar{\epsilon} = 300, 500, \text{ and } 1000$," *Journal of Computational Physics* 178, no. 2 (2002): 427-436.

79. C. Mimeau, I. Mortazavi, and G.-H. Cottet, "Passive Control of the Flow Around a Hemisphere Using Porous Media," *European Journal of Mechanics - B/Fluids* 65 (2017): 213-226.

80. M. S. Birrittella, M. Debbage, R. Huggahalli, et al., "Intel Omni-Path Architecture: Enabling Scalable, High Performance Fabrics," in *Proceedings of the 2015 IEEE 23rd Annual Symposium on High-Performance Interconnects* (IEEE Computer Society, 2015), 1-9.

81. Mellanox, "Introducing 200G HDR InfiniBand Solutions," 2022, https://network.nvidia.com/sites/default/files/pdf/whitepapers/WP_Introducing_200G_HDR_InfiniBand_Solutions.pdf.

82. GENCI, "GENCI and CNRS Choose Eviden to Make the Jean Zay Supercomputer One of the Most Powerful in France," 2024, https://www.genci.fr/en/news/genci-and-cnrs-choose-eviden-make-jean-zay-supercomputer-one-most-powerful-france.

83. NVIDIA, "NVIDIA HPC SDK Documentation," 2025, https://docs.nvidia.com/hpc-sdk/index.html.

84. J. Hursey, E. Mallove, J. M. Squyres, and A. Lumsdaine, "An Extensible Framework for Distributed Testing of MPI Implementations," in *Proceedings, Euro PVM/MPI* (Springer, Berlin, Heidelberg, 2007).

85. NVIDIA, "NVIDIA Nsight Systems," 2025, https://developer.nvidia.com/nsight-systems.

86. S. Williams, A. Waterman, and D. Patterson, "Roofline: An Insightful Visual Performance Model for Multicore Architectures," *Communications of the ACM* 52, no. 4 (2009): 65-76.

87. Berkley Lab, "CS Roofline Toolkit Repository," 2025, https://bitbucket.org/berkeleylab/cs-roofline-toolkit/src/master/.

88. C. Yang, R. Gayatri, T. Kurth, et al., "An Empirical Roofline Methodology for Quantitatively Assessing Performance Portability," in *2018 IEEE/ACM International Workshop on Performance, Portability and Productivity in HPC (P3HPC)* (IEEE, 2018), 14-23.

89. T. Nguyen, S. Williams, M. Siracusa, C. MacLean, D. Doerfler, and N. J. Wright, "The Performance and Energy Efficiency Potential of FPGAs in

Scientific Computing," in *2020 IEEE/ACM Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS)* (IEEE, 2020), 8-19.

90. W. Ertel, "On the Definition of Speedup," in *PARLE'94 Parallel Architectures and Languages Europe*, ed. C. Halatsis, D. Maritsas, G. Philokyprou, and S. Theodoridis (Springer Science & Business Media, 1994), 289-300.

91. R. Winkelmann, J. Häuser, and R. D. Williams, "Strategies for Parallel and Numerical Scalability of CFD Codes," *Computer Methods in Applied Mechanics and Engineering* 174, no. 3 (1999): 433-456.

92. I. Bermejo-Moreno, J. Bodart, J. Larsson, B. M. Barney, J. W. Nichols, and S. Jones, "Solving the Compressible Navier-Stokes Equations on up to 1.97 Million Cores and 4.1 Trillion Grid Points," in *SC 13: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis* (ACM, Association for Computing Machinery, 2013), 1-10.

93. F. Jia, Z. J. Wang, R. Bhaskaran, U. Paliath, and G. M. Laskowski, "Accuracy, Efficiency and Scalability of Explicit and Implicit FR/CPR Schemes in Large Eddy Simulation," *Computers & Fluids* 195 (2019): 104316.

94. C. Junqueira-Junior, J. L. F. Azevedo, J. Panetta, W. R. Wolf, and S. Yamouni, "On the Scalability of CFD Tool for Supersonic Jet Flow Configurations," *Parallel Computing* 93 (2020): 102620.

95. Z. Liu, X. Chu, X. Lv, H. Liu, H. Fu, and G. Yang, "Accelerating Large-Scale CFD Simulations With Lattice Boltzmann Method on a 40-Million-Core Sunway Supercomputer," in *ICPP 23: Proceedings of the 52nd International Conference on Parallel Processing* (ACM, Association for Computing Machinery, 2023), 797-806.

96. W. Verdier and P. Kestener, "Alain Cartalade," *Computer Methods in Applied Mechanics and Engineering* 370, no. 113266 (2020): 113266.

97. G. Jahoski and J. Kupis, "Performance Optimization of Implementation of Lattice Boltzmann Method in ARUZ," in *25th International Conference Mixed Design of Integrated Circuits and System (MIXDES)* (IEEE, 2018), 188-191.

98. M. Bauer, S. Eibl, C. Godenschwager, et al., "Development and Application of Open-Source Software for Problems With Numerical PDEs," *Computers & Mathematics With Applications* 81 (2021): 478-501.

99. France's Transmission System Operator, "eCO_2mix: C CO₂ Emissions per kWh of Electricity Generated in France," in *Réseau de Transport D'Électricité (RTE), Real-Time Data* (Réseau de Transport d'Électricité (RTE) en France, 2025), https://www.rte-france.com/en/eco2mix/co2-emissions.

100. France's Transmission System Operator, "Carbon Intensity of Electricity Production in Metropolitan France in 2023," in *RTE Report: 32 g CO₂ eq/kWh* (Réseau de Transport d'Électricité (RTE) en France, 2024), https://analysesetdonnees.rte-france.com/bilan-electrique-2023/emissions.

101. CEEMS, "Compute Energy & Emissions Monitoring Stack (CEEMS) Git Repository," 2025, https://github.com/mahendrapaipuri/ceems.

102. CEEMS, "Compute Energy & Emissions Monitoring Stack (CEEMS) Documentation," 2025, https://mahendrapaipuri.github.io/ceems/docs/.