ParCFD 2026 37th International Conference on Parallel Computational Fluid Dynamics Nov 03-06 2026, Nantes, FRANCE
SCALING WATERLILY.JL WITH MPI AND AN IMPROVED GEOMETRIC MULTIGRID SOLVER Bernat FONT∗ , Marin LAUBER, Tzu-Yao HUANG, Gabriel D. WEYMOUTH
arXiv:2607.07687v1 [physics.comp-ph] 8 Jul 2026
Delft University of Technology Faculty of Mechanical Engineering Mekelweg 2, 2628 CD Delft, The Netherlands ∗ e-mail: [email protected], web page: https://b-fg.github.io
Key words: distributed-memory parallelism, message-passing interface, geometric multigrid, finite volume, Julia Abstract. We present recent performance-oriented developments in WaterLily.jl, a scale-resolving incompressible flow solver written in pure Julia that runs seamlessly on CPUs and GPUs of any vendor. Supported by the newly added MPI-based parallelism, strong-scalability tests display a near-ideal linear trend, and weak-scaling efficiency is kept above 85% before node memory-concurrency contention dominates parallel performance. Inter-node weak scalability is sustained above 96% with grid size up to 1 billion cells. We further benchmark improvements to the geometric multigrid Poisson solver enabled by an adaptive under-relaxed red-black Gauss–Seidel smoother together with anisotropic coarsening operators. 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
Bernat Font, Marin Lauber, Tzu-Yao Huang, Gabriel D. Weymouth
10 6 1.0 0.8 p
10 4
E = t1/tn
t / Δt [ms]
10 5
10 3 TGV 2563 TGV 5123 TGV 10243 Sphere 256 × 96 × 96 Sphere 512 × 192 × 192 Sphere 1024 × 384 × 384
10 2 10 1 10 0
1
2
4
8
16
0.6
Internode : E = t64/tn
p
1.00
0.4
0.95 0.90
0.2
np
32
64
128 256 512
0.0
3
TGV 64 / rank TGV 1283 / rank
1
2
4
8
0.85 0.80
16
64 128 192 256 320 384 448 512
np
32
64
128 256 512
Figure 1: Strong scaling (left) and weak scaling (right) results on the Rome CPU partition of the Snellius Dutch national supercomputer. The inset shows the inter-node weakscaling efficiency, normalised by the single-node (64-rank) time. 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 ReD = 3700 in a (16 × 6 × 6)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 memoryconcurrency contention discussed below. The trend deteriorates when MPI halo updates dominate the local workload, once the smallest local subdomain dimension approaches ∼16 cells. The weak-scalability results maintain an efficiency (E) above 90% for the small (643 cells) and large (1283 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 109 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,
Bernat Font, Marin Lauber, Tzu-Yao Huang, Gabriel D. Weymouth
Table 1: Solver comparison for the case of a simplified swimming jellyfish. Grid size is 4 × 23p . Cost reported in ns per cell per time-step running on a laptop RTX 4060 GPU. p Solver
# Iter. Final resid.
Cost
p Solver
# Iter. Final resid. Cost
5 PCG GMG-PCG GMG-RBGS
14.90 1.75 1.70
9.45 × 10−5 3.09 × 10−5 5.44 × 10−5
139.84 86.97 38.22
6 PCG GMG-PCG GMG-RBGS
26.09 1.89 2.34
9.43 × 10−5 5.33 × 10−5 5.69 × 10−5
95.93 26.25 16.61
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 lowfrequency 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 immersedboundary flows. In contrast, Red-Black Gauss–Seidel (RBGS) naturally damps highfrequency 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 1) 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. 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 ≥ 4, 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.
Bernat Font, Marin Lauber, Tzu-Yao Huang, Gabriel D. Weymouth
anisotropic
isotropic
N=2
domain: (N:1)
Number of V-cycles
N=4
N=6 256
256
256
domain: (N:1:1)
64
64
16
16
16
4
4
4
1
1 2
3 time step
4
N=10 domain: (N:N:1)
64
1
N=8
1 1
2
3 time step
4
1
2
3
4
time step
Figure 2: Total number of V-cycle (predictor + corrector steps) for the first 4 time steps of canonical impulsive-flow problems: (left) 2D flow around a circular cylinder, (center) 3D flow around a spanwise-periodic cylinder and (right) 3D flow around a sphere. The different domain aspect ratios are given as a small insert in the top-right corner of each subplot. 4
CONCLUSIONS
Recent developments in WaterLily.jl have improved the Poisson solver through anisotropic 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 [1] G. Weymouth and B. Font, “WaterLily.jl: A differentiable and backend-agnostic Julia solver for incompressible viscous flow around dynamic bodies,” Computer Physics Communications, vol. 315, p. 109748, 2025. doi:10.1016/j.cpc.2025.109748, github.com/WaterLily-jl/WaterLily.jl. [2] V. Churavy, “KernelAbstractions.jl.” github.com/JuliaGPU/KernelAbstractions.jl.
doi:10.5281/zenodo.4021259,
[3] S. Omlin, L. Räss, and I. Utkin, “Distributed Parallelization of xPU Stencil Computations in Julia,” JuliaCon Proceedings, vol. 6, p. 137, Nov. 2024. doi:10.21105/jcon.00137, github.com/eth-cscs/ImplicitGlobalGrid.jl. [4] S. Williams, A. Waterman, and D. Patterson, “Roofline: an insightful visual performance model for multicore architectures,” Communications of the ACM, vol. 52, no. 4, pp. 65–76, 2009. doi:10.1145/1498765.1498785.