1 INTRODUCTION↩︎

High-performance computing (HPC) is a cornerstone of computational fluid dynamics, enabling an ever-growing resolution in scale-resolving simulations. Driven by hardware advances such as general-purpose processor (CPU) and accelerator-based (GPU) nodes, CFD software is rapidly evolving to harvest the increase in computational power. The trend comes at an added software-complexity cost, and CFD developers now heavily rely on HPC libraries to offload compute-expensive kernels and parallelisation.

In this direction, we present performance-oriented developments and the first parallel scalability results of WaterLily.jl [1]. Written in pure Julia, WaterLily is designed to be a compact and minimalist finite-volume incompressible-flow solver that runs on any CPU or GPU computing backend. This feature is enabled by separating computing workloads (typically for loops) from the algorithmic space-time discretization. By using Julia’s metaprogramming support, the solver ultimately specialises each kernel to the available computing backend using KernelAbstractions.jl[2].

2 MPI PARALLELISATION AND SCALING↩︎

Further leveraging Julia’s scientific computing ecosystem, distributed-memory parallelism has been recently integrated in WaterLily via the ImplicitGlobalGrid.jl (IGG) package [3]. Backed by the MPI.jl wrapper, IGG uses a Cartesian MPI communicator naturally mapping to regular grids, and offering CUDA-aware or ROCm-aware MPI support for GPU-based applications. Additionally, IGG works with different grid sizes, thus enabling halo updates within geometric-multigrid (GMG) levels without added complexity from the IGG-user side.

Strong and weak scalability results are presented in Figure 1. The tests are conducted with the SIMD single-thread CPU backend with single precision (FP32) and half-node occupancy (64 of the 128 cores of the dual-socket AMD EPYC 7H12 nodes). The test cases comprise a Taylor–Green vortex (TGV) at \(Re=1600\) in a triple-periodic cubic domain, and flow past a sphere at \(Re_D=3700\) in a \((16\times6\times6)D\) domain. The MPI domain decomposition always minimizes the halo surface for a given number of processes. The strong-scalability results closely follow the ideal linear reduction in time per time-step, even super-linearly for the largest grids as shrinking per-rank load relaxes the memory-concurrency contention discussed below. The trend deteriorates when MPI halo updates dominate the local workload, once the smallest local subdomain dimension approaches \({\sim}16\) cells. The weak-scalability results maintain an efficiency \((E)\) above 90% for the small (\(64^3\) cells) and large (\(128^3\) cells) loads up to 16 ranks. Beyond, efficiency drops due to intra-node memory-concurrency contention, motivating the 50% node occupancy. This reflects the memory-bound nature of the solver’s stencil kernels [4]. While memory bandwidth is not saturated (checked with hardware counters), ranks compete to access the per-chiplet link to the memory controllers. Ranks are hence pinned with a uniform stride across the node. Beyond the node, weak scaling is nearly ideal, remaining above 96% up to 512 ranks and \(10^9\) cells (Figure 1, inset).

3 ANISOTROPIC GEOMETRIC MULTIGRID SOLVER WITH ADAPTIVE UNDER-RELAXED RED-BLACK GAUSS–SEIDEL SMOOTHER↩︎

WaterLily’s general-purpose Poisson solver is Preconditioned Conjugate Gradient (PCG), which applies to arbitrary grid sizes and extends naturally to parallel backends. However, its global inner products require costly reductions in MPI simulations, and the number of iterations increases with problem size, resulting in poor weak scaling. For this reason, the default WaterLily solver uses GMG, which accelerates convergence by combining smoothers to remove high-frequency errors with coarse-grid correction (CGC) for low-frequency errors. The GMG solver has been further optimized with the development of an adaptive smoother and anisotropic coarsening.

3.1 Adaptive under-relaxed Red-Black Gauss–Seidel↩︎

While PCG can be used as a smoother, its non-local behaviour conflicts with CGC, particularly for the variable-coefficient Poisson equations arising in multiphase and immersed-boundary flows. In contrast, Red-Black Gauss–Seidel (RBGS) naturally damps high-frequency errors and avoids the global reductions required by PCG, making it well suited to parallel multigrid. To ensure monotonic residual decay in highly heterogeneous media, we introduce a dynamically under-relaxed RBGS smoother whose relaxation factor adapts to residual growth or decay after each V-cycle. Numerical tests (Table ¿tbl:tab:side95by95side95comparison?) show that this adaptive GMG-RBGS solver maintains robust convergence while reducing solver cost by 1.6–2.3× compared with a PCG smoother, and by 3.7–5.8× compared with the generic PCG solver.

4pt

3.2 Anisotropic multigrid coarsening↩︎

With isotropic coarsening, the smallest domain dimension limits the number of multigrid levels, bounding the decay of low-frequency error. For highly anisotropic domains, this can severely penalize the convergence of the Poisson solver. This is especially important for impulsive flow cases, where the low wavenumber axial pressure wave remains at an intermediate wavenumber on the coarsest level, hindering the smoother’s efficacy. Instead, anisotropic coarsening downsamples every dimension independently and reaches the coarsest domain size where the low wavenumber content is easily removed by the smoother. Implemented in WaterLily, the anisotropic multigrid coarsening leverages the Cartesian structure of the grid with minimal overhead compared to the classical isotropic approach. In contrast to isotropic coarsening, we find that this approach allows us to effectively reduce the number of multigrid cycles for impulsive flow problems on high aspect-ratio rectangular grids, as displayed in Figure 2. For domains with aspect ratio \(N\ge4\), approximately an order of magnitude fewer V-cycles are required to converge to the same residual. This trend is maintained for all cases and aspect ratios tested.

4 CONCLUSIONS↩︎

Recent developments in WaterLily.jl have improved the Poisson solver through aniso-
tropic GMG coarsening and an RBGS smoother, targeting both faster convergence and total solve wall-time reduction. Together with the newly implemented distributed-memory parallelism, tested in strong and weak-scalability setups, WaterLily is now ready to leverage many-GPU architectures and embrace the exascale computing challenge. This opens new research avenues in scale-resolving simulations including vortex-dominated flows, fluid-structure interaction problems, and optimization.

References↩︎