# Performance analysis of the free surface lattice Boltzmann implementation in waLBerla

Jonas Plewinski, Christoph Alt, Harald Köstler, Ulrich Rüde

► To cite this version:

Jonas Plewinski, Christoph Alt, Harald Köstler, Ulrich Rüde. Performance analysis of the free surface lattice Boltzmann implementation in waLBerla. PAMM, 2024, 24, 10.1002/pamm.202400196. hal-04963725

HAL Id: hal-04963725

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

Submitted on 24 Feb 2025

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 Attribution 4.0 International License
# Performance analysis of the free surface lattice Boltzmann implementation in waLBerla

Jonas Plewinski<fnref n="1" /> ID

Christoph Alt<fnref n="1" />,<fnref n="2" /> ID

Harald Köstler<fnref n="1" />,<fnref n="3" /> ID

Ulrich Rüde<fnref n="1" />,<fnref n="4" /> ID

## Correspondence

Jonas Plewinski, Computer Science 10 (System Simulation), Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany.
Email: jonas.plewinski@fau.de

## Funding information

Deutsche Forschungsgemeinschaft,
Grant/Award Number: 434946896

## Abstract

The free-surface lattice Boltzmann method uses a volume of fluid approach to simulate immiscible two-fluid flow problems. It divides the simulation domain into three distinct phases—gas, fluid, and interface—where computation within the gas phase is disregarded. The interface delineates a one-cell-thick layer between the first two phases, validated physically for implementation in the HPC C++ multiphysics framework WALBERLA but lacking an exhaustive performance analysis. This paper aims to shed light on node-level performance on different architectures, employing continuous benchmarking, showing and analyzing weak scaling results on a modern HPC cluster, the Fritz supercomputer, and reporting energy consumption for the current implementation.

# 1 | INTRODUCTION

Free surface flows are a category in fluid dynamics characterized by an interface between a dense liquid and a lighter gas phase, typically air. These flows exhibit a distinct boundary—the free surface—whose position can change over time and under the influence of various forces. Free surface flows can model ocean waves, the flow in rivers and stream, or the behavior of liquids in containers, and the dynamics of droplets and bubbles.

High-performance computing (HPC) is pivotal in advancing our understanding of free surface flows. The complex nature of these three-dimensional flows leads to physical models that often require large-scale numerical simulations to allow researchers to analyze them with high enough resolution. Efficient parallel solvers can help to capture the full complexity of free surface flows, providing otherwise unattainable insights.

This manuscript is focussed on a systematic performance analysis of free surface flows within the waLBerla multi-physics simulation framework. We analyze the current state of the software, highlight the challenges towards reaching maximal performance for free surface simulations, and discuss which additional special HPC's techniques will be needed for additional performance improvements.

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.

© 2024 The Author(s). *Proceedings in Applied Mathematics & Mechanics* published by Wiley-VCH GmbH.

<footnote n="1">¹Computer Science 10 (System Simulation), Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany</footnote>

<footnote n="2">²Paderborn University, Paderborn Center for Parallel Computing, Paderborn, Germany</footnote>

<footnote n="3">³Erlangen National High Performance Computing Center (NHR@FAU), Erlangen, Germany</footnote>

<footnote n="4">⁴CERFACS, Toulouse Cedex 1, France</footnote>
## 2 | NUMERICAL METHODS

### 2.1 | The lattice Boltzmann method

The lattice Boltzmann method (LBM) is a modern computational fluid dynamics approach [1]. It discretizes the Boltzmann equation on a Cartesian grid with spacing $\Delta x \in \mathbb{R}^+$ using DdQq velocity sets (standard sets include D2Q9, D3Q19, or D3Q27 for 2D or 3D problems, respectively). Particle distribution functions (PDFs) $f_i(\mathbf{x}, t) \in \mathbb{R}$ represent the density probability of fluid particles. The discrete LBM equation is given by

$$f_i(\mathbf{x} + c_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t) + \Omega_i(\mathbf{x}, t) + F(\mathbf{x}, t). \quad (1)$$

Here, $\Omega_i(\mathbf{x}, t)$ denotes the collision operator, and $F(\mathbf{x}, t)$ accounts for external forces. To enhance performance, the computation is split into two parts: execution of the right hand side, known as collision and updating the post-collision values to the correct position in neighboring cells, known as streaming.

The collision operator serves as a parameter that can be adjusted to tailor the algorithm to specific problem requirements [2]. The simplest operator is the single relaxation time (SRT) operator, which is defined by

$$\Omega_i(\mathbf{x}, t) = \frac{f_i(\mathbf{x}, t) - f_i^{\text{eq}}(\mathbf{x}, t)}{\tau} \Delta t \quad (2)$$

with relaxation time $\tau \in (\Delta t/2, 2\Delta t)$. Other collision operator includes the two-relaxation time (TRT), multiple-relaxation time (MRT), central moments, cumulants, and entropic collision operators.

The equilibrium

$$f_i^{\text{eq}}(\mathbf{x}, t) = w_i \rho \left( 1 + \frac{\mathbf{u} \cdot \mathbf{c}_i}{c_s^2} + \frac{(\mathbf{u} \cdot \mathbf{c}_i)^2}{2c_s^4} - \frac{\mathbf{u} \cdot \mathbf{u}}{2c_s^2} \right). \quad (3)$$

is derived from the Maxwell-Boltzmann equation.

Macroscopic quantities such as the density

$$\rho(\mathbf{x}, t) = \sum_i f_i(\mathbf{x}, t) \quad (4)$$

and velocity

$$\mathbf{u}(\mathbf{x}, t) = \frac{F(\mathbf{x}, t)\Delta t}{2\rho(\mathbf{x}, t)} + \frac{1}{\rho(\mathbf{x}, t)} \sum_i c_i f_i(\mathbf{x}, t), \quad (5)$$

can be calculated from moments of the particle distribution functions (PDFs).

### 2.2 | The free-surface lattice Boltzmann method

The free-surface lattice Boltzmann method (FSLBM) employed in this work uses the implementation of Schwarzmeier et al. [3–5]. It models a dynamic interface between two immiscible fluids, where the heavier fluid primarily governs the flow dynamics. This reduces the problem to a single-fluid flow with a free boundary, applicable when there is a significant difference in the densities and viscosities of the fluids. The heavier phase is typically considered as the fluid, and lighter one as gas.

The interface is modeled using the volume-of-fluid approach [6], where each lattice cell is assigned a fill level $\varphi(\mathbf{x}, t)$ to indicate its phase. Cells can be categorized as fluid ($\varphi(\mathbf{x}, t) = 1$), gas ($\varphi(\mathbf{x}, t) = 0$), or interface ($\varphi(\mathbf{x}, t) \in (0, 1)$). In the FSLBM, fluid, gas, and interface cells are handled differently. The simulation of the gas phase is neglected entirely; the fluid cells are only executing the pure Lattice Boltzmann Method (LBM), whereas for the free surface, a lot more computations are necessary.
To compute a free surface, we have to compute the liquid mass of a cell, which is determined by

$$m(\mathbf{x}, t) = \varphi(\mathbf{x}, t)\rho(\mathbf{x}, t)\Delta x^3 \quad (6)$$

Then, the mass flux between the interface and other cells during the streaming is computed by

$$\frac{\Delta m_i(\mathbf{x}, t)}{\Delta x^3} = \begin{cases} 0 & \mathbf{x} + \mathbf{c}_i \Delta t \in \text{gas} \\ f_i^*(\mathbf{x} + \mathbf{c}_i \Delta t, t) - f_i^*(\mathbf{x}, t) & \mathbf{x} + \mathbf{c}_i \Delta t \in \text{liquid} \\ \frac{1}{2}(\varphi(\mathbf{x}, t) + \varphi(\mathbf{x} + \mathbf{c}_i \Delta t, t))(f_i^*(\mathbf{x} + \mathbf{c}_i \Delta t, t) - f_i^*(\mathbf{x}, t)) & \mathbf{x} + \mathbf{c}_i \Delta t \in \text{interface}, \end{cases} \quad (7)$$

with $\bar{i}$ defining the reverse lattice direction resulting in $\bar{c}_i = -c_i$.

After determining the mass flux for the interface cells, it is necessary to compute a new mass. Depending on the new fill level, cell conversions are identified with

$$\frac{m_{\text{ex}}(\mathbf{x}, t)}{\rho(\mathbf{x}, t)\Delta x^3} = \begin{cases} \varphi^{\text{conv}}(\mathbf{x}, t) - 1 & \text{if } \mathbf{x} \text{ is converted to liquid} \\ \varphi^{\text{conv}}(\mathbf{x}, t) & \text{if } \mathbf{x} \text{ is converted to gas.} \end{cases} \quad (8)$$

Here $\varepsilon_\varphi = 10^{-2}$ serves as a threshold to prevent oscillatory conversions [7]. The excess mass due to cell conversions has to be accounted for, otherwise mass conservation would be violated. Therefore, it is evenly distributed among neighboring interface cells to conserve total mass. Unnecessary interface cells without gas or liquid neighbors are adjusted in mass flux.

When converting cells, there is no modification of the PDFs. In gas-to-interface conversions, PDFs are initialized based on the equilibrium. The collision and streaming (1) occur in interface and liquid cells.

The macroscopic boundary condition [8, 9] at the free surface is modeled by

$$\begin{align} p(\mathbf{x}, t) - p^G(\mathbf{x}, t) + p^L(\mathbf{x}, t) &= 2\mu\partial_n u_n(\mathbf{x}, t) \nonumber \\ 0 &= \partial_{t_1} u_n(\mathbf{x}, t) + \partial_n u_{t_1}(\mathbf{x}, t) \nonumber \\ 0 &= \partial_{t_2} u_n(\mathbf{x}, t) + \partial_n u_{t_2}(\mathbf{x}, t). \tag{9} \end{align}$$

with gas pressure $p^G(\mathbf{x}, t)$, Laplace pressure $p^L(\mathbf{x}, t)$, tangent vectors $\mathbf{t}_1(\mathbf{x}, t) \in \mathbb{R}^d$, and $\mathbf{t}_2(\mathbf{x}, t) \in \mathbb{R}^d$. At the free surface interface, the LBM anti-bounce-back pressure boundary condition [10]

$$f_i^*(\mathbf{x} - \mathbf{c}_i \Delta t, t) = f_i^{\text{eq}}(\rho^G(\mathbf{x}, t), \mathbf{u}(\mathbf{x}, t)) + f_i^{\text{eq}}(\rho^L(\mathbf{x}, t), \mathbf{u}(\mathbf{x}, t)) - f_i^*(\mathbf{x}, t) \quad (10)$$

is used. Here $\rho^G(\mathbf{x}, t)$ defines the gas density and $\mathbf{u}(\mathbf{x}, t)$ the velocity of the free surface interface. Due to their unavailability, all PDFs streaming from gas cells to interface cells must adhere to the free-surface boundary condition (10).

The gas pressure is expressed as

$$p^G(\mathbf{x}, t) = p^V(t) - p^L(\mathbf{x}, t), \quad (11)$$

where $p^V(t)$ represents volume pressure, and $p^L(\mathbf{x}, t)$ is Laplace pressure. The Laplace pressure is determined by surface tension $\sigma \in \mathbb{R}^+$ and interface curvature $\kappa(\mathbf{x}, t) \in \mathbb{R}$

$$p^L(\mathbf{x}, t) = 2\sigma\kappa(\mathbf{x}, t). \quad (12)$$

In this study, the interface curvature is calculated using finite difference methods [11]

$$\kappa(\mathbf{x}, t) = -\nabla \cdot \hat{\mathbf{n}}(\mathbf{x}, t), \quad (13)$$

where the normalized interface normal is obtained through a weighted central finite difference methods [12]

$$\mathbf{n}(\mathbf{x}, t) = \nabla \varphi(\mathbf{x}, t). \quad (14)$$
<table><caption>ALGORITHM 1 Overview of FSLBM algorithm (based on Schwarzmeier [19])</caption><tr><td>1:</td><td>Initialization of PDFs</td><td>▷ pure LBM</td></tr><tr><td>2:</td><td>fill level, cell-types</td><td></td></tr><tr><td>3:</td><td><b>for each</b> time step <i>t</i> <b>do</b></td><td></td></tr><tr><td>4:</td><td>Smoothing of fill level</td><td></td></tr><tr><td>5:</td><td>Compute normal of free-surface</td><td></td></tr><tr><td>6:</td><td>Compute curvature of free-surface</td><td></td></tr><tr><td>7:</td><td>Apply LBM boundary conditions</td><td>▷ pure LBM</td></tr><tr><td>8:</td><td>Compute force density</td><td></td></tr><tr><td>9:</td><td>Apply free-surface boundary conditions</td><td></td></tr><tr><td>10:</td><td>Reconstruction of missing gas-side PDFs in interface</td><td></td></tr><tr><td>11:</td><td>Perform LBM streaming in interface</td><td></td></tr><tr><td>12:</td><td>Exchange mass in interface cells and update fill level</td><td></td></tr><tr><td>13:</td><td>Execute LBM collision and streaming</td><td>▷ pure LBM</td></tr><tr><td>14:</td><td>Cell conversion and refilling of fluid cells</td><td></td></tr><tr><td>15:</td><td>Distribute excessive mass</td><td></td></tr><tr><td>16:</td><td>Initialization of new free-surface cell-types</td><td></td></tr><tr><td>17:</td><td><b>end for</b></td><td></td></tr></table>

The computation of $n(\mathbf{x}, t) \in \mathbb{R}^d$ is adjusted near solid obstacle cells [13]. The curvature $\kappa(\mathbf{x}, t)$ is effectively computed from the second-order derivative of the fill level $\varphi(\mathbf{x}, t)$. To mitigate errors introduced by the nonsmooth indicator function $\varphi(\mathbf{x}, t)$, a smoothing process using the $K_s$-Kernel [14] with a support radius of 2.0 is applied [11].

## 3 | IMPLEMENTATION IN WALBERLA

WALBERLA [15, 16] is a framework programmed in C++ and Python and specifically designed for simulations on massively parallel systems. Originally focused on simulating fluid flow problems with LBM, WALBERLA has evolved into a versatile multiphysics framework, encompassing rigid particle dynamics, phase-field simulations, and coupled algorithms like the FSLBM. It employs fully distributed data structures on uniform blocks to ensure scalability and supports both pure message passing interface (MPI) and hybrid MPI/OpenMP parallelization. As described in Section 2.1, the simulation domain is discretized in a cartesian grid for the LBM. The cells of this grid are partitioned into uniform blocks that can then be distributed to multiple processes. Communication is realized with MPI.

WALBERLA features a modular architecture and employs advanced code generation techniques *lbmpy* [17, 18] as essential techniques to improve on user productivity, reusability, and maintainability. The Python package *lbmpy* uses metaprogramming techniques to facilitate the efficient implementation and automated generation of various LBM variants, producing highly optimized compute kernels for different hardware. These techniques allow easy switching between different implemented stencil types, streaming patterns, and collision operators.

An overview of the FSLBM introduced in Section 2.2 is now illustrated in Algorithm 1. Depending on the state of a block, that is, whether a block contains only fluid, or only gas, or a mixture of both with an interface in between, the algorithm branches into different operations:

* Only fluid cells present ⇒ pure LBM is executed, namely lines 1, 7, and 11 in Algorithm 1

* Only gas cells present ⇒ nothing is executed

* Interface cells present ⇒ the whole algorithm given in Algorithm 1 is executed

For the fluid blocks, a highly efficient fused streaming and collision kernel generated by *lbmpy* for the LBM is used. The FSLBM parts are currently natively implemented in waLBerla and only support CPUs. In order to get a better understanding of the performance results, which will be shown later on, we take a closer look at how the FSLBM is implemented in waLBerla. Hence, we look at the different parts of the algorithm and how many load/store operations are present to
<table><caption>**TABLE 1** Memory usage of the FSLBM per cell and time step of each part of the algorithm.</caption><thead><tr><td>Sweep</td><td>Memory usage [Bytes]</td></tr></thead><tbody><tr><td>Smoothing of fill level</td><td>28</td></tr><tr><td>Compute normal of free-surface</td><td>60</td></tr><tr><td>Compute curvature of free-surface</td><td>68</td></tr><tr><td>Compute force density</td><td>192</td></tr><tr><td>Apply free-surface boundary condition</td><td>360</td></tr><tr><td>Execute LBM collision</td><td>328</td></tr><tr><td>Execute cell conversion</td><td>160</td></tr><tr><td>Refilling of converted fluid cells</td><td>308</td></tr><tr><td>Distribute excessive mass</td><td>172</td></tr><tr><td>Initialization of new free-surface cell-types</td><td>8</td></tr><tr><td>Total memory usage per cell and time step</td><td>1704</td></tr></tbody></table>

estimate what theoretical performance could be achieved with the current implementation. This information is shown in Table 1. We extracted the information for the table by looking at the code, specifically the data structures, so-called fields in waLBerla, that are used per sweep. Fields in waLBerla are four-dimensional data structures, where the first three dimensions represent the spatial 3D position, and the fourth dimension holds the data for that specific position. That can be either a scalar value (e.g., the macroscopic density), 3D vectors (e.g., the macroscopic velocity), or holding the discrete simulation data of the PDF-fields (e.g., for a D3Q19 stencil, the fourth dimension would hold a vector with 19 values). Additionally, we look at the datatype, loads, and stores. For example, how exactly we calculated the memory usage of a sweep, we look at the most expensive part of the algorithm, the *application of the free-surface boundary conditions*. Hereafter, we list the datatypes used, if these are only loaded, stored, or both, and note the datatype:

* pdfField datatype: double field with 19 values per cell; load and store ⇒ 19 · (8 + 8) Bytes = 304 Bytes

* flagField datatype: scalar int32-field; load and store ⇒ 4 + 4 Bytes = 8 Bytes

* fillField datatype: scalar double field; load and store ⇒ 8 + 8 Bytes = 16 Bytes

* curvatureField datatype: scalar double field; load ⇒ 8 Bytes

* normalField datatype: double field with 3D vector per cell; load ⇒ 3 · 8 Bytes = 24 Bytes

This leads to a memory usage of 360 Bytes per cell for the *application of the free-surface boundary conditions*. The memory usage per cell of the other kernels is shown in Table 1 and is calculated in the same way. This results in 1.704-kB memory transfer per cell and time step for the entire algorithm executed on a block with a free surface.

The benchmark employed in this study emulates a gravity wave, as depicted in Figure 1. Prior investigations have leveraged this and similar test cases to verify the physical accuracy of the FSLBM implementation [3–5]. The fluid is initialized as a wave with an amplitude $a_0$. Additionally, periodic boundary conditions in the x- and z-direction and no-slip boundary conditions in the y-direction are used. Furthermore, gravity is acting on the fluid. This test case is also suitable for a performance analysis because there is an even amount of fluid and gas cells, and the free surface is moving across the entire domain, such that we have no domination of either phase.

# 4 | BENCHMARKING RESULTS

## 4.1 | Node-level results on various architectures

As presented in Alt et al. [20], the gravity wave benchmark described in Section 3 is part of the continuous benchmarking pipeline for WALBERLA. With continuous benchmarking, this and other benchmarks are executed automatically on the Testcluster at NHR@FAU [21] whenever a code change is submitted to the code base. The Testcluster is a special HPC resource and contains compute nodes with different CPU and GPU architectures. All compute nodes are different in this cluster; thus only running single node benchmarks is feasible. Thus, these tests aim at monitoring the single-node
<figure><img src="image_2.png" /><figcaption>**FIGURE 1** Illustrates the initialization of the gravity wave, using the fluid depth *h*, initialized amplitude of the wave *a*₀, wavenumber *k* = 2π/l, wavelength *l*, and gravitational acceleration *g*. In *x*- and *z*-direction, periodic and no slip boundary conditions in *y*-direction were used (based on Schwarzmeier et al. [3]).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**FIGURE 2** Visualization panel, which shows the distribution of the simulation time for the gravity wave benchmark in the continuous benchmarking pipeline executed on different hardware architectures. The red part is for communication, and the green is for computation. The name of the used compute nodes is showed on the *x*-axis. The names correspond to the used CPU architectures: Intel “Cascade Lake,” AMD “Genoa,” Intel “Haswell,” AMD “Bergamo,” two times AMD “Genoa,” and two times Intel “Ice Lake” in that order from left to right.</figcaption></figure>

performance. Though sometimes overlooked, good node level performance is a key prerequisite for achieving good performance also on multinode configuratons. Thus, it is essential that the continuous benchmarking system tracks the impact of all code changes and their impact on the performance for various hardware architectures. Note that this is accomplished automatically, without forcing a developer to run any benchmarks manually. In Figure 2, we see the visualization of some benchmarking results gathered in our continuous benchmarking pipeline for the gravity wave test case. The FSLBM algorithm is complex and consists of multiple substeps (cf. Section 2.2). Figure 2 shows the contribution of the
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Benchmarking results of the FSLBM gravity wave test case on the Fritz supercomputer with up to 384 nodes.</figcaption></figure>

different substeps to the simulation time on different hardware architectures. Communication steps are highlighted in red, and actual computation parts are highlighted in green. This shows that communication takes 35%–45% of the overall simulation time for single-node benchmarks.

## 4.2 | Multinode results on the fritz supercomputer

### 4.2.1 | Performance

We want to show and analyze the performance of the current FSLBM implementation in WALBERLA. Currently, the FSLBM implementation only supports CPUs. Therefore, we chose the Fritz Supercomputer [22] to execute these benchmarks. Fritz uses per node two Intel Xeon Platinum 8360Y ("Ice Lake") CPUs with 36 cores each, resulting in a total of 72 cores per node.

In this article, we use the benchmark described in Section 3, execute it for every node distribution multiple times to get a reproducible and accurate result and simulate 1000 time steps in order to obtain a realistic estimate of the performance of the FSLBM. These are weak scaling benchmarks, meaning the total domain is scaled by the number of cores used for the respective benchmark. Therefore, the workload for every core remains constant. In our case, every core is assigned a block of the simulation domain containing 64³ cells. For this specific type of benchmark, an optimal scaling is reported if the simulation time remains constant for every number of cores, resulting in horizontal graph. The weak scaling performance results are shown in Figure 3; once the *time to solution* values in Figure 3A and once with the standard performance metric for lattice Boltzmann simulations *Mega Lattice Updates per Second (MLUP/s)* in Figure 3B. The latter refers to the rate of lattice updates, that is, for a pure LBM implementation without extensions shows how many streaming and collision operations are performed per second. MLUP/s per CPU-Core are calculated by

$$MLUP/s = \frac{\text{Domain size} \cdot \text{\#time steps}}{\text{Simulation time} \cdot \text{\#cores}} \cdot 1e6 \qquad (15)$$

Clearly, there is an immediate correlation between these two measured values: if the time to solution increases, the MLUP/s drop. Every blue point in the graphs represents one benchmark. The red line represents the median of these
measured values. Additionally, we plotted the maximum performance, including its computation in Figure 3B. In order to evaluate how well the current implementation utilizes the machine's resources, we compare the actual performance with an estimate of what could be maximally achievable. For this, we refer to Section 2.2, where we found that 1704 Bytes must be transferred per lattice update for each interface cell. If we assume having 45% fluid and gas cells and 10% interface cells, this results in an average of $0.45 \cdot 304$ Byte + $0.1 \cdot 1704$ Bytes = 307.2 Bytes needed per lattice update. From this, an optimistic value for the best possible performance can be calculated as the maximal memory bandwidth (per core) divided by the bytes that need to be loaded for every lattice update (cf. Figure 3B). This must be put in context with a realistic estimate of the bandwidth of a compute node. Specifically, on the Fritz supercomputer, we used the copy-benchmark of the likwid tool suite [23] to find a maximal bandwidth of 266.3 GB/s. The optimistic upper bound for the FSLBM performance is then found as $266 \text{ GB/s}/(72 \text{ cores} \cdot 307.2 \text{ Byte/lup}) = 12.15 \text{ MLUP/s per core}$.

Such prepared, we analyze the current performance as follows: First, there are two outliers in the measurement for eight and 64 nodes, and an equal number of values for 256 nodes. This emphasizes the importance of running multiple benchmarks per node distribution since there is always a danger that the results will be affected by some other activity on the machine. Second, the overall trend looks promising since the trend indicates nearly constant performance per core and thus good scalability, except for a jump when progressing from a single node to several nodes. Furthermore, we observe a slight decrease of performance from 16 nodes onwards. Regarding the first drop of performance, this is expected since this is a natural consequence of needing communication between nodes, which causes extra overhead. Since the FSLBM algorithm is compute intensive, we observe here a moderate decrease of performance of about 7%. The reported performance for 384 is likely an outlier, due to technical difficulties of Fritz, it was only possible to run one benchmark. The moderate performance degradation from 16 nodes onwards is a natural consequence of the high volume of communication.

As a third observation, however, we only achieve approximately 7% of the optimistic maximum performance of 12.15 MLUP/s per core. Thus, while the algorithm scales, this should be put in context of the performance obtained for a pure LBM simulation. There 85% of the theoretical peak performance is possible and can be achieved with WALBERLA for CPUs [20] and even with GPUs [24]. Clearly, this indicates a substantial potential for further optimization based on the current FSLBM implementation. Thus, though parallel scalability seems to be good, the node level performance analysis indicates that the method can potentially benefit from further code-tuning.

To investigate this further, we evaluate the percentages of the different parts of the algorithm split into computation (sweeps, background is dotted) and communication. This is now illustrated in Figure 4. Here, we can also make a few deductions: First, from two to 64 nodes, the percentages of the FSLBM algorithm are almost exactly at 50% for the computation and 50% for the communication. From then onward, it increased slightly. Second, we can directly see the main contributors to the computational costs. Applying the free-surface boundary condition (streamReconstructAdvectSweep), fill level smoothing and normal computation are the most expensive parts of computation. Moreover, for communication, the cost is directly proportional to what has to be communicated; an PDF field is much more expensive than a vector or scalar field. Third, approximately, all parts of the algorithm scale, independent of sweep or communication, although the communication seems to increase overall slightly. Nevertheless, returning to the fact that 50% of the simulation time is used for communication: Or rather it seems like that, since we are not using any load balancing here, the blocks and thus processes without a free-surface have significantly less work to do. Since we only make a distinction between communication and computation. Consequently, when these cores have finished their work, they start the communication, while the free-surface blocks still perform computation. Given this, a more intricate analysis will be needed, to distinguish the time spent in synchronization from time spent in communication. These and other possible improvements will be summarized in Section 5.

Here, we conclude by finally presenting an evaluation of the energy consumption for the FSLBM implementation.

## 4.2.2 | Energy consumption

With more and more powerful supercomputers, the energy bill for operating such high-performance machines is financial and environmental becoming a critical factor. Hence, it is crucial to use HPCs as efficiently as possible. Therefore, we also look at the energy consumption of the FSLBM implementation in WALBERLA. To gain further insights, we use Clustercockpit [25]. It is a framework used to monitor job-specific performance and power consumption on distributed HPC clusters. We run the same benchmark as described above on 64 nodes, but this time for half an hour (approximately 5000 time steps); the results regarding the energy consumption are visualized in the screenshot from Clustercockpit in Figure 5.
<figure><img src="image_5.png" /><figcaption>FIGURE 4 Benchmarking results showing the percental distribution of the different parts of the FSLBM algorithm divided into computation (sweeps) (with dots) and communication (no dots).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5 Screenshot of one 30-m benchmark run with 64 nodes on Fritz taken from Clustercockpit regarding the energy consumption per node. On the y-axis, the power consumption per node in Watts is shown, and the x-axis shows the time in minutes.</figcaption></figure>

We see that the power energy consumption per node is oscillating close to 500 W during almost the whole benchmark, except for the beginning and end. This is promising only if the CPU workload is also high; the measurements for that are shown in Figure 6. Here, we see that the workload is also constant (except for the beginning), and we achieve a high workload per node.
<figure><img src="image_7.png" /><figcaption>**FIGURE 6** Screenshot of one 30-m benchmark run with 64 nodes on Fritz taken from Clustercockpit regarding the CPU load per node. On x-axis again the time is noted, and the y-axis shows the workload of per node. Every node consists of 72 core, such that a workload of 72 shows a 100% utilization of every core.</figcaption></figure>

## 5 | CONCLUSION AND OUTLOOK

This paper presents a performance analysis of the FSLBM implementation in WALBERLA. It shows node-level performance results on different architectures using the CB pipeline and weak scaling performance results on the Fritz supercomputer. While the node-level measurements show a significant difference of the results depending on which architecture was chosen, the benchmarks on Fritz give more insights into the performance of different parts of the FSLBM implementation. Even though the different parts of the algorithm are scaling well, there is optimization potential for the entire algorithm. We will further investigate the load imbalance using HPC tracer tools to estimate how bad the imbalance is and how much communication and synchronization we have with the current implementation of using load balancing. It is also planned to partially re-implement the algorithm using meta-programming strategies and include it to the code generation infrastructure in WALBERLA. This enables a vast pool of optimizations, for example, using common subexpression elimination to reduce the operations per computation, using nontemporal stores to reduce the memory traffic, using communication hiding, and advanced streaming patterns for the LBM (currently, it is only possible to use a simple pull-streaming pattern). Furthermore, it allows fusing kernels to reduce the number of communications necessary per time step.

## ACKNOWLEDGMENTS

This project has received funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under the project number 434946896. The authors gratefully acknowledge the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) under the NHR project "HEISSRISSE—Massively Parallel Simulation of the Melt Pool Area during Laser Beam Welding using the Lattice Boltzmann Method." NHR funding is provided by federal and Bavarian state authorities. NHR@FAU hardware is partially funded by the German Research Foundation (DFG)—440719683.

Open access funding enabled and organized by Projekt DEAL.

## ORCID

Jonas Plewinski https://orcid.org/0000-0002-4815-7643

Christoph Alt https://orcid.org/0000-0001-8897-5205

Harald Köstler https://orcid.org/0000-0002-6992-2690

Ulrich Rüde https://orcid.org/0000-0001-8796-8599
## REFERENCES

1. Krüger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Silva, G., & Viggen, E. M. (2017). *The lattice Boltzmann method: Principles and practice*. Springer.

2. Coreixas, C., Chopard, B., & Latt, J. (2019). Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations. *Physical Review E*, *100*, 033305.

3. Schwarzmeier, C., Holzer, M., Mitchell, T., Lehmann, M., Häusl, F., & Rüde, U. (2023). Comparison of free-surface and conservative Allen-Cahn phase-field lattice Boltzmann method. *Journal of Computational Physics*, *473*, 111753.

4. Schwarzmeier, C., & Rüde, U. (2022). Comparison of refilling schemes in the free-surface lattice Boltzmann method. *AIP Advances*, *12*, 115324.

5. Schwarzmeier, C., & Rüde, U. (2023). Analysis and comparison of boundary condition variants in the free-surface lattice Boltzmann method. *International Journal For Numerical Methods in Fluids*, *95*, 820–850.

6. Hirt, C., & Nichols, B. (1981). Volume of fluid (VOF) method for the dynamics of free boundaries. *Journal of Computational Physics*, *39*(1), 201–225.

7. Pohl, T. (2008). *High performance simulation of free surface flows using the lattice Boltzmann method* [PhD thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg].

8. Scardovelli, R., & Zaleski, S. (1999). Direct numerical simulation of free-surface and interfacial flow. *Annual Review of Fluid Mechanics*, *31*(1), 567–603.

9. Bogner, S. (2017). *Direct numerical simulation of liquid-gas-solid flows based on the lattice Boltzmann method* [PhD thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg].

10. Körner, C., Thies, M., Hofmann, T., Thürey, N., & Rüde, U. (2005). Lattice Boltzmann model for free surface flow for modeling foaming. *Journal of Statistical Physics*, *121*(1), 179–196.

11. Bogner, S., Rüde, U., & Harting, J. (2016). Curvature estimation from a volume-of-fluid indicator function for the simulation of surface tension and wetting with a free-surface lattice Boltzmann method. *Physical Review E*, *93*(4), 043302.

12. Parker, B. J., & Youngs, D. L. (1992). *Two and three dimensional Eulerian simulation and fluid flow with material interfaces* (Technical Report 01/92). UK Atomic Weapons Establishment.

13. Donath, S. (2011). *Wetting models for a parallel high-performance free surface lattice Boltzmann method* [PhD thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg].

14. Williams, M. W., Kothe, D. B., & Puckett, E. G. (1999). Accuracy and convergence of continuum surface tension models. In W. Shyy & R. Narayanan (Eds.), *Fluid dynamics at interfaces*. (1st ed., pp. 294–305). Cambridge University Press.

15. Bauer, M., Eibl, S., Godenschwager, C., Kohl, N., Kuron, M., Rettinger, C., Schornbaum, F., Schwarzmeier, C., Thönnes, D., Köstler, H., & Rüde, U. (2021). waLBerla: A block-structured high-performance framework for multiphysics simulations. *Computers & Mathematics with Applications*, *81*, 478–501.

16. Chair for System Simulation, FAU. (2023). *waLBerla* (widely applicable Lattice Boltzmann from Erlangen). (v6.1). Zenodo. https://doi.org/10.5281/zenodo.10054460

17. Bauer, M., Köstler, H., & Rüde, U. (2021). lbmpy: Automatic code generation for efficient parallel lattice Boltzmann methods. *Journal of Computational Science*, *49*, 101269.

18. Hennig, F., Holzer, M., & Rüde, U. (2023). Advanced automatic code generation for multiple relaxation-time lattice Boltzmann methods. *SIAM Journal on Scientific Computing*, *45*(4), C233–C254.

19. Schwarzmeier, C. (2024). *Parallel high-performance simulation of liquid-gas-solid flows with the free-surface lattice Boltzmann method* [PhD thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg].

20. Alt, C., Lanser, M., Plewinski, J., Janki, A., Klawonn, A., Köstler, H., Selzer, M., & Rüde, U. (2024). A continuous benchmarking infrastructure for high-performance computing applications. *International Journal of Parallel, Emergent and Distributed Systems*, *39*(4), 501–523.

21. NHR@FAU Team. (2024). *Testcluster*. https://doc.nhr.fau.de/clusters/testcluster/

22. NHR@FAU Team. (2024). *Fritz — NHR@FAU HPC documentation*. https://doc.nhr.fau.de/clusters/fritz/

23. Gruber, T., Eitzinger, J., Hager, G., & Wellein, G. (2022). LIKWID. Zenodo.

24. Holzer, M., Bauer, M., Köstler, H., & Rüde, U. (2021). Highly efficient lattice Boltzmann multiphase simulations of immiscible fluids at high-density ratios on CPUs and GPUs through code generation. *The International Journal of High Performance Computing Applications*, *35*(4), 413–427.

25. Eitzinger, J., Gruber, T., Afzal, A., Zeiser, T., & Wellein, G. (2019). ClusterCockpit—a web application for job-specific performance monitoring. In *2019 IEEE international conference on cluster computing (CLUSTER)* (pp. 1–7).

**How to cite this article:** Plewinski, J., Alt, C., Köstler, H., & Rüde, U. (2024). Performance analysis of the free surface lattice Boltzmann implementation in waLBerla. *Proceedings in Applied Mathematics and Mechanics*, 24, e202400196. https://doi.org/10.1002/pamm.202400196