Pith. sign in

REVIEW 2 major objections 5 minor 1 cited by

Scaling WaterLily.jl with MPI and an improved geometric multigrid solver

T0 review · 2 major / 5 minor · reviewed 2026-07-10 · grok-4.5

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read MPI-enabled WaterLily keeps near-ideal scaling to a billion cells while an adaptive anisotropic multigrid solver cuts Poisson cost several-fold.

desk verdict Solid ParCFD engineering note: first multi-node MPI numbers for WaterLily plus measured Poisson cost cuts; claims match the data on the tested CPU setup, with the usual short-paper generalization limits. read the letter →

arxiv 2607.07687 v2 pith:5LGSBEEY submitted 2026-07-08 physics.comp-ph cs.DCphysics.flu-dyn

Bernat Font , Marin Lauber , Tzu-Yao Huang , Gabriel D. Weymouth This is my paper · ORCID
classification physics.comp-phcs.DCphysics.flu-dyn
keywords distributed-memoryparallelismmessage-passinginterfacegeometricmultigridfinitevolumeJuliaincompressibleflowstrongandweakscalingPoissonsolver
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

WaterLily is a compact pure-Julia finite-volume solver for scale-resolving incompressible flow that already runs on any CPU or GPU backend. This paper reports the addition of distributed-memory MPI parallelism and measured strong- and weak-scaling behaviour on a national supercomputer. Strong scaling follows a near-ideal linear drop in time per step; weak-scaling efficiency stays above 85 percent until intra-node memory contention appears, and inter-node efficiency remains above 96 percent out to 512 ranks and roughly one billion cells. In parallel the authors replace the previous Poisson solver options with a geometric multigrid method that uses an adaptive under-relaxed red-black Gauss–Seidel smoother and anisotropic coarsening. On a simplified swimming-jellyfish case the new solver reduces cost per cell by factors of roughly 1.6–5.8 relative to preconditioned conjugate-gradient alternatives, and on high-aspect-ratio impulsive-flow problems it cuts the number of multigrid V-cycles by about an order of magnitude. Together the two advances are presented as making WaterLily ready for large multi-node and multi-GPU scale-resolving simulations.

What carries the argument

Anisotropic geometric multigrid with an adaptive under-relaxed red-black Gauss–Seidel smoother: each grid dimension is coarsened independently so that low-wavenumber errors are fully resolved on the coarsest level, while the relaxation factor is adjusted after each V-cycle according to residual growth or decay, avoiding the global reductions of PCG and preserving monotonic residual decay on heterogeneous coefficients.

What would settle it

Repeat the same strong- and weak-scaling suite and the jellyfish / high-aspect-ratio impulsive-flow residual histories on full-node occupancy, multi-GPU nodes, or a different memory hierarchy; if inter-node efficiency falls well below 96 percent or the multigrid cycle advantage disappears, the central performance claims do not generalise.

Watch

Extended reading notes

Core claim

With newly added MPI parallelism WaterLily exhibits near-ideal strong scaling and sustains inter-node weak-scaling efficiency above 96 percent up to grids of about one billion cells; simultaneously an adaptive under-relaxed red-black Gauss–Seidel geometric multigrid solver with anisotropic coarsening reduces Poisson-solver cost by roughly 1.6–5.8 times versus earlier PCG-based options and markedly lowers V-cycle counts on anisotropic domains.

Load-bearing premise

The reported near-ideal efficiencies and large cycle reductions are assumed to hold beyond the half-node FP32 SIMD CPU tests, the two main flow cases, and the simplified jellyfish and impulsive-flow benchmarks used in the paper.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

2 major / 5 minor

Summary. The manuscript reports two performance-oriented advances in WaterLily.jl, a pure-Julia finite-volume incompressible-flow solver that targets CPUs and GPUs via KernelAbstractions.jl. First, distributed-memory parallelism is added through ImplicitGlobalGrid.jl/MPI.jl; strong-scaling tests on the Snellius Rome CPU partition (SIMD FP32, half-node occupancy) show near-ideal linear reduction in time per step for TGV and sphere cases, while weak scaling remains above ~85% until intra-node memory-concurrency contention appears and inter-node efficiency stays above 96% up to 512 ranks and ~10^9 cells (Fig. 1). Second, the geometric multigrid Poisson solver is improved by an adaptive under-relaxed red-black Gauss–Seidel smoother and anisotropic coarsening; Table 1 reports 1.6–5.8 imes cost reductions versus PCG-based options on a simplified jellyfish case, and Fig. 2 shows roughly an order-of-magnitude fewer V-cycles for impulsive flows on high-aspect-ratio domains. The authors conclude that the solver is now positioned for multi-GPU/exascale work.

Significance. If the reported numbers hold under the stated conditions, the paper supplies concrete, reproducible evidence that a compact Julia CFD code can achieve near-ideal strong scaling and high inter-node weak-scaling efficiency on a national CPU system, together with substantial Poisson-solver cost reductions from relatively simple multigrid changes. The open-source stack (WaterLily, KernelAbstractions, ImplicitGlobalGrid) and the explicit diagnosis of memory-concurrency limits are useful engineering contributions for the Julia HPC and scale-resolving CFD communities. The work is appropriately scoped as a short ParCFD-style performance note rather than a full methods paper; the main value is the measured scaling and cost data on defined benchmarks.

major comments (2)
  1. §2 and Fig. 1: All strong- and weak-scaling results are restricted to the SIMD single-thread CPU backend, FP32, and 50% node occupancy with uniform-stride pinning on dual-socket Rome nodes. The abstract and conclusions claim readiness for many-GPU architectures, yet no multi-GPU or multi-node GPU scaling data are shown despite IGG’s CUDA/ROCm-aware support. At least one GPU weak- or strong-scaling curve (or an explicit statement that GPU results are deferred) is needed so the central “exascale-ready” claim is not overstated relative to the evidence.
  2. §3.1 and Table 1: The adaptive under-relaxed RBGS is a load-bearing algorithmic claim (1.6–2.3 imes vs GMG-PCG, 3.7–5.8 imes vs plain PCG), but the adaptation rule for the relaxation factor is only described qualitatively (“adapts to residual growth or decay after each V-cycle”). Without the update formula, initial value, or bounds, the result cannot be reproduced or assessed for robustness on other residual tolerances or media. A short algorithmic statement or pseudocode is required for the cost-reduction claim to stand.
minor comments (5)
  1. Fig. 1 caption and body: Clarify units and definition of t/Δt (ms per time step) and of the efficiency E = t1/tnp; the inset normalizes by the 64-rank single-node time, which should be stated once in the main text as well.
  2. Table 1: Report the residual tolerance used as the stopping criterion and whether iteration counts are averages over a fixed number of time steps; residual values alone leave the comparison slightly under-specified.
  3. §3.2 / Fig. 2: State the residual tolerance and the number of smoothing steps per level so that the V-cycle counts for anisotropic vs isotropic coarsening can be compared quantitatively by others.
  4. Minor typography: several missing spaces after commas/periods in the abstract and §1–2 (e.g., “distributed-memoryparallelism”, “geometricmulti-grid”); fix for camera-ready.
  5. References: the Williams et al. Roofline citation is appropriate; a brief pointer to a standard anisotropic multigrid reference would help readers place the coarsening choice in the broader literature.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: empirical scaling and multigrid cost measurements on stated benchmarks, not predictions forced by definition or self-citation.

full rationale

WaterLily.jl reports measured wall-clock times, residual histories, V-cycle counts, and parallel efficiencies for MPI strong/weak scaling (Fig. 1) and for an adaptive under-relaxed RBGS geometric multigrid solver with anisotropic coarsening (Table 1, Fig. 2). These are engineering benchmarks on defined cases (TGV, sphere, simplified jellyfish, impulsive cylinder/sphere flows) under stated conditions (SIMD FP32, half-node occupancy, Snellius Rome). Self-citations supply the baseline solver and Julia packages (WaterLily.jl, KernelAbstractions.jl, ImplicitGlobalGrid.jl); they do not force the reported efficiencies or cost ratios by construction. There is no fitted free parameter renamed as a prediction, no uniqueness theorem imported from the authors, and no algebraic identity that makes the scaling or cycle-count claims tautological. The derivation chain is observational measurement, not circular deduction. Score 0 is appropriate.

Assumptions & free parameters 3 free parameters · 4 assumptions · 0 invented entities

Central claims rest on standard incompressible CFD and multigrid practice plus measured runtimes on specific hardware. No new physical entities. Free choices are engineering (occupancy, residual targets, when to adapt under-relaxation) rather than fitted physical constants. Domain assumptions are the usual finite-volume Poisson/multigrid setting and the Cartesian MPI decomposition provided by ImplicitGlobalGrid.

free parameters (3)
  • adaptive under-relaxation factor schedule
    Relaxation factor is adjusted dynamically after each V-cycle based on residual growth/decay; the exact update rule and bounds are not fully specified and act as a tunable control that enables the reported robust convergence.
  • node occupancy / rank pinning (50% occupancy, uniform stride)
    Chosen to mitigate memory-concurrency contention on dual-socket EPYC nodes; scaling efficiencies depend on this operational choice.
  • Poisson residual convergence tolerance (~1e-4 final residual in Table 1)
    Iteration counts and cost comparisons are reported at a chosen residual level; different tolerances would change absolute costs though relative ranking may hold.
assumptions (4)
  • domain assumption Incompressible Navier–Stokes on a structured finite-volume grid with a Poisson pressure solve is an adequate model for the reported scale-resolving cases.
    Background of the entire WaterLily solver; assumed throughout §§1–3.
  • domain assumption Geometric multigrid with smoothers for high-frequency error and coarse-grid correction for low-frequency error converges for the variable-coefficient Poisson problems arising in multiphase/immersed-boundary flows when the smoother is suitably under-relaxed.
    Stated motivation in §3 for preferring RBGS over PCG as smoother.
  • domain assumption Cartesian MPI domain decomposition with halo exchange (via ImplicitGlobalGrid) correctly implements distributed finite-volume operators including multigrid levels.
    §2 integration of IGG; correctness of scaling results depends on correct parallel operators.
  • standard math Standard arithmetic and parallel-timing measurements on Snellius Rome nodes are reliable for efficiency ratios.
    Implicit in Fig. 1 efficiency definitions E = t1/tnp and internode E = t64/tnp.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Scaling WaterLily.jl with MPI and an improved geometric multigrid solver." pith.science (2026). https://pith.science/paper/5LGSBEEY

@misc{pith2026260707687,
  author       = {Pith},
  title        = {Pith review of: Scaling WaterLily.jl with MPI and an improved geometric multigrid solver},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/5LGSBEEY}},
  note         = {Machine review of arXiv:2607.07687}
}
read the original abstract

We present recent performance-oriented developments in WaterLily, 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.

Figures

Figures reproduced from arXiv: 2607.07687 by the authors.

Figure 1
Figure 1. Strong scaling (left) and weak scaling (right) results on the Rome CPU partition [PITH_FULL_IMAGE:figures/full_fig_p002_1.png] view at source ↗
Figure 2
Figure 2. Total number of V-cycle (predictor + corrector steps) for the first 4 time steps [PITH_FULL_IMAGE:figures/full_fig_p004_2.png] view at source ↗

Review history (2 revisions) →

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

  1. Julia for CFD: A Critical Survey of Ecosystem, Performance, and Composability

    cs.CE 2026-08 conditional novelty 3.0 of 10

    Julia is now credible for several research-facing CFD regimes, but its advantage is integration and composability, not unique performance, and it is not yet a complete industrial CFD platform.

Reference graph

Works this paper leans on

4 extracted references · 4 canonical work pages · cited by 1 Pith paper

  1. [1]

    WaterLily.jl: A differentiable and backend-agnostic Julia solver for incompressible viscous flow around dynamic bodies,

    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

    work page doi:10.1016/j.cpc.2025.109748 2025

  2. [2]

    KernelAbstractions.jl

    V. Churavy, “KernelAbstractions.jl.”doi:10.5281/zenodo.4021259, github.com/JuliaGPU/KernelAbstractions.jl

    work page doi:10.5281/zenodo.4021259

  3. [3]

    Distributed Parallelization of xPU Sten- cil Computations in Julia,

    S. Omlin, L. Räss, and I. Utkin, “Distributed Parallelization of xPU Sten- cil Computations in Julia,”JuliaCon Proceedings, vol. 6, p. 137, Nov. 2024. doi:10.21105/jcon.00137,github.com/eth-cscs/ImplicitGlobalGrid.jl

    work page doi:10.21105/jcon.00137 2024

  4. [4]

    Roofline: an insightful visual perfor- mance model for multicore architectures,

    S. Williams, A. Waterman, and D. Patterson, “Roofline: an insightful visual perfor- mance model for multicore architectures,”Communications of the ACM, vol. 52, no. 4, pp. 65–76, 2009.doi:10.1145/1498765.1498785

    work page doi:10.1145/1498765.1498785 2009

Pith tools

Reviewed July 10, 2026 · model on record in the stance chip above.