Parallel scaling: peclet.flow vs the field

Weak scaling of the incompressible Navier–Stokes solver on CPU (MPI × OpenMP) and multi-GPU, from a 24-core workstation to 32 H100s on Snellius — benchmarked head-to-head against CaNS, incflo (AMReX) and OpenFOAM on the identical triply periodic case.

flow

HPC

MPI

multi-GPU

CPU

benchmark

scaling

Author

Peclet

Published

August 10, 2026

The questions

Three things a prospective user of a parallel CFD solver actually wants to know, each answered here by measurement:

Does weak scaling survive the interconnect? On-node parallelism is cheap; the test of a distributed incompressible solver is what happens to its global pressure solve when halo exchanges and reductions start crossing nodes.
What does the same code deliver on CPU and on GPU, from a workstation to a national supercomputer (Snellius: Genoa nodes with 192 cores; H100 nodes with 4 GPUs)?
How does it compare with leading open-source codes on the identical case? We run three references, chosen to bracket the design space:
- CaNS — the state of the art in raw speed for trivial geometry: explicit RK3 + an exact FFT Poisson solve, possible only on boxes.
- incflo — the DOE-exascale (AMReX) solver in peclet’s own method family: finite-volume incompressible NS, geometric-multigrid projections, embedded-boundary capable.
- OpenFOAM (ESI v2412) — the ubiquitous general-purpose unstructured code, icoFoam + GAMG.

The instrument: a boundary-free Taylor–Green box

Scaling studies want a case that grows without complications. We tile a triply periodic box with Taylor–Green vortices (one vortex per $64^3$ tile, $\mathrm{Re}_\text{tile}=100$): no walls, every cell does identical work, and the domain extends in any direction by adding tiles. The TGV initial condition is also a Stokes eigenfunction, giving the driver a built-in analytic correctness check (matched to $2\times10^{-3}$, identical across rank counts and backends).

tgv_bench.py runs the case distributed and reports warmup-excluded per-step wall time with a per-phase breakdown (predictor / momentum solve / pressure projection), the pressure iteration count, and the accumulated time and count of every global reduction inside the pressure solve — so algorithmic effects and communication effects are never conflated.

Fairness rules, applied symmetrically. Every code runs the identical grids at matched per-step tolerances (pressure ~1e-5/1e-6 class; CaNS’s FFT solve is exact); timings are steady-state per-step wall times, excluding each code’s setup (meshing, FFT plans, grid initialization). References run expert-tuned, not default: incflo gets max_grid_size=64 and MLMG tolerances matched to the study (its AMReX defaults — $32^3$ box shredding and ~1e-11 tolerances — cost it a measured 1.53×, which we refuse to charge it). OpenFOAM uses GAMG with a 64-bit-label build (required above ~2×10⁸ cells). CaNS runs its recommended adaptive CFL 0.95 — its RK3 step advances ~4× more physical time per step than the fixed-CFL codes, which per-step throughput comparisons deliberately do not penalize. On the cluster, single measurements are noise: production points are pinned (--cpu-bind=cores), and repeated draws then reproduce to within 5–10 % — against up to 2× node-set variability measured for unpinned runs, so pinning is part of the method, and differences below ~10 % between single draws should not be over-read.

Workstation (24-core Threadripper PRO 5965WX, RTX 5080)

Hybrid MPI × OpenMP mix

Left: throughput by MPI ranks × OpenMP threads on 24 cores (fixed $192^3$ TGV). Right: the per-phase step decomposition.

The mix barely matters: every configuration lands within 14 % (16.4 Mcell/s at 4 ranks × 6 threads, down to 14.4 for pure MPI or pure OpenMP), a 12.6× speedup over serial on 24 cores. The phase split shows why: the momentum and pressure sweeps are memory-bandwidth bound, so once all memory channels are fed it makes little difference whether the parallelism comes from ranks or threads. Mid-hybrid mixes (4×6, 12×2) edge ahead by avoiding both extremes — 24 ranks’ worth of halo surface on one side, a single rank’s thread-synchronisation on the other. The same flat landscape reappears on Genoa (below).

Weak scaling within one node: four codes, one wall

Weak scaling at 2.1 M cells/rank, pure MPI, up to 24 cores. Left: efficiency vs. own 1-rank baseline. Right: absolute throughput (log scale); the dotted line is the single RTX 5080 running the same case.

With per-rank work fixed, all four codes decay together — and this is not communication (peclet’s pressure iteration count is pinned at 4.0 at every rank count). A single node’s weak scaling measures shared memory-bandwidth saturation, nothing else: the efficiency ordering (incflo 43 %, peclet 42 %, CaNS 31 %, OpenFOAM 29 % at 24 ranks) is the inverse of per-core bandwidth demand — the codes that stream least per core-second saturate the shared memory system last. Absolute throughput on 24 cores: CaNS 41.4 Mcell/s, peclet 13.9, incflo 4.4, OpenFOAM 3.3 — peclet is 4.2× OpenFOAM and 3.2× incflo at identical settings, while CaNS’s remaining 3× lead is the price of generality made visible: an exact FFT solve that exists only because the box is trivial, against machinery that handles arbitrary cut-cell geometry at large d$t$. And the single RTX 5080 at 118 Mcell/s outruns CaNS’s entire 24-core node by 2.9× — the design point of the suite: the CPU path is a capable oracle, the GPU path is the product.

Single-GPU ladder (RTX 5080, 16 GB): three codes on the same device

Both GPU-capable references run here from the same sources as their CPU curves: incflo via its AMReX CUDA backend, CaNS via its OpenACC + cuDecomp backend (NVIDIA HPC SDK toolchain, cuFFT). OpenFOAM has no native GPU solver path. All three at their GPU-fair settings.

grid	cells	peclet	incflo	CaNS
$128^3$	2.1 M	99.7	16.1	131.6
$160^3$	4.1 M	104.7	17.9	136.8
$192^3$	7.1 M	118.1	18.6	138.3
$208^3$	9.0 M	100.1	18.7	129.7

(Mcell/s, higher is better.) Two clean results:

peclet is 6.3× incflo on the same GPU — a far wider margin than on CPU, the GPU-first design paying off where it is meant to.
CaNS keeps the single-device crown by just 17 % (138 vs 118 at $192^3$) — with no communication to pay, the exact FFT solve sets the throughput ceiling, and a general-geometry cut-cell multigrid comes within a sixth of it. What happens to that lead once the FFT’s all-to-all meets an interconnect is the multi-GPU question — measured below.

peclet holds ~100–118 Mcell/s across the whole ladder until 16 GB fills, iteration counts flat at 4–5 regardless of grid divisibility.

Does the domain decomposition care about awkward rank counts?

Sweeping prime and composite rank counts on a fixed $192^3$ box (pure MPI, 5965WX):

ranks	ms/step	Mcell/s	pressure iters	max-block imbalance
5	1765	4.0	4.0	1.04×
7 (prime)	1331	5.3	4.0	1.02×
8	1160	6.1	4.0	1.00×
13 (prime)	795	8.9	4.0	1.08×
16	665	10.6	4.0	1.00×
23 (prime)	544	13.0	4.0	1.20×

No. The pressure iteration count is exactly 4.0 at every rank count — prime or not — because the multigrid’s convergence is decomposition-independent by construction. What remains is the ORB’s geometric imbalance, at worst 20 % of the max block at np = 23, with throughput on the same trend line as the neighbouring powers of two. Users can pick rank counts by hardware, not by numerology.

Cross-node CPU weak scaling on Snellius Genoa

The central experiment: 188.7 M cells per node (768×640×384 each), grown with the node count on Genoa (192 cores/node), every code on the identical case. A single Genoa node peaks at 75 Mcell/s (48 ranks × 4 threads; every hybrid mix from 12×16 to 192×1 lands within 20 %, iterations pinned at 4.0 throughout), 5× the whole 24-core workstation. peclet’s cross-node curves run in two hybrid flavors — the classic 96×2 and the fat-rank 12×16, which uses 8× fewer MPI ranks and therefore suffers less straggler skew at the collectives; the references run their natural pure-MPI mode.

Cross-node weak scaling on Snellius Genoa, 188.7 Mcells/node fixed. Left: efficiency vs own 1-node baseline. Right: aggregate throughput; the dotted line is ideal scaling at peclet’s slope.

nodes	peclet 96×2	peclet 12×16	CaNS	incflo	OpenFOAM
1	63 Mcell/s	67	109	16	9.4
2	134 (106 %)	109 (82 %)	20 (9 %)	32 (99 %)	—
4	259 (103 %)	225 (84 %)	153 (35 %)	67 (103 %)	—
8	437 (87 %)	459 (86 %)	45 (5 %)	117 (90 %)	—

Three structural results:

The geometric-multigrid codes scale; the FFT code does not. peclet and incflo both hold near-flat cross-node curves (single-draw wiggles above 100 % are node-set noise) — halo exchanges plus a handful of latency-bound reductions per iteration survive the interconnect. At 8 nodes the fat-rank 12×16 flavor wins outright (459 vs 437 Mcell/s) with visibly less collective-wait time (249 vs 382 ms/step) — fewer ranks, less straggler skew, exactly as the skew analysis predicts.
CaNS’s single-node crown does not survive the interconnect. 1.6× faster than peclet’s best flavor on one node — and by two nodes it has collapsed to 9 % efficiency (reproduced in independent runs), at eight nodes to 5 %: eight CaNS nodes deliver well under one peclet node’s throughput. Its FFT solve requires global all-to-all transposes every substep — nearly free in shared memory, ruinous across nodes. (Caveat: CaNS ran its automatic pencil-grid choice; expert dims tuning may soften, but cannot remove, the all-to-all.) The crossover lands before two nodes.
Within the method family, peclet leads on speed at equal scaling: ~4× incflo per node with the same flat curve, and 7× OpenFOAM on its single-node point. (OpenFOAM’s multi-node points are limited by its serial mesh generation at these sizes — an operational cost of the unstructured design that the structured codes don’t have.)

Multi-GPU weak scaling: 1–32 H100, three codes

The same case at 47.2 M cells per GPU (4 H100 per node), 1 → 32 GPUs — from one device to eight nodes, with both GPU-capable references on the identical grids:

Left: GPU weak scaling, 47.2 Mcells/GPU fixed — peclet at 90 % efficiency and 7.7 Gcell/s aggregate on 32 H100s. Right: the pressure-solver ablation at the inter-node points (8 and 16 GPUs).

H100 GPUs	peclet	incflo	CaNS
1	269 Mcell/s	61	609
4	998 (92 %)	225 (92 %)	856 (35 %)
8	1968 (92 %)	429 (88 %)	297 (7.6 %)
16	3936 (92 %)	848 (85 %)	977 (10 %)
32	7748 (90 %)	1665 (85 %)	1032 (5.9 %)

The GPU fabric replays the CPU story, at higher stakes. One H100 runs peclet at 269 Mcell/s — four Genoa nodes’ worth — and 32 H100s sustain 7.7 Gcell/s at 90 % weak efficiency: a 1.5-billion-cell grid advancing five time steps per second, with iterations pinned at 4.0 and the entire global-reduction budget under 2 % of the step. incflo scales just as flatly (85 %) at 4.4× less per-GPU speed — the method family’s signature. And CaNS repeats its arc precisely: the single-device crown (609 Mcell/s, 2.3× peclet — cuFFT with nothing to communicate), already sagging within one NVLink node (35 % at 4 GPUs), collapsing at the first node boundary (7.6 %), and flatlining near 1 Gcell/s while peclet pulls away 7.5×. The crossover sits between 4 and 8 GPUs. (CaNS-GPU caveat: cuDecomp ran its autotuned MPI transpose backends — the NVHPC-bundled NCCL errors on this fabric — which particularly understates its intra-node point.)

The ablation on the right is how the solver’s defaults are chosen — measured, at scale, and re-validated after every solver change:

MG-PCG with fine-scoped nullspace removal (the default) wins at both 8 and 16 GPUs (191.8 ms/step).
Removing the mean on every multigrid level instead costs 4 % — pure reduction latency.
The host-staged halo costs 16 % vs GPU-aware MPI (device pointers straight to the interconnect).
Chebyshev comes within 15 % (5 iterations vs PCG’s 4) — a respectable reduction-free fallback, but PCG’s convergence still beats its communication savings at every scale measured.
The agglomerated algebraic bottom solve (GraphAMG) pays far more than it saves on GPU.

What the study says, in one paragraph

On the identical boundary-free case, peclet.flow is the fastest general-geometry incompressible solver in the comparison on every platform tested — 4.2× OpenFOAM and 3.2× incflo per CPU node on the workstation, ~4× incflo per Genoa node, 6.3× incflo on the same GPU — and, together with incflo, the only code family whose weak scaling survives the interconnect: ~87 % across 8 CPU nodes and 90 % across 32 H100s, where the FFT state of the art collapses below 10 % on both fabrics. CaNS keeps the single-device crowns (by 3× on a CPU node, 17 % on a consumer GPU, 2.3× on one H100) — the honest price of general geometry against an exact solve that only exists on trivial boxes — and loses them the moment the problem spans two of anything: the crossover lands before two nodes on CPU and between 4 and 8 GPUs on H100. The same source runs 118 Mcell/s on a consumer RTX 5080 (2.9× CaNS’s entire 24-core node), 269 Mcell/s on one H100, and 7.7 Gcell/s at 90 % weak efficiency on 32 H100s — a 1.5-billion-cell incompressible flow advancing five steps per second.

Reproducing

All drivers, case definitions, runner scripts for the three reference codes, the Snellius job pack, raw result JSONs and the figure scripts live in benchmarks/parallel-scaling/ — tgv_bench.py (peclet), run_cans.sh, run_incflo.sh, openfoam-tgv/, snellius/README.md (exact submission commands), plot_workstation.py / plot_snellius.py.

--- title: "Parallel scaling: peclet.flow vs the field" subtitle: "Weak scaling of the incompressible Navier–Stokes solver on CPU (MPI × OpenMP) and multi-GPU, from a 24-core workstation to 32 H100s on Snellius — benchmarked head-to-head against CaNS, incflo (AMReX) and OpenFOAM on the identical triply periodic case." author: "Peclet" date: "2026-08-10" categories: [flow, HPC, MPI, multi-GPU, CPU, benchmark, scaling] --- ## The questions Three things a prospective user of a parallel CFD solver actually wants to know, each answered here by measurement: 1. **Does weak scaling survive the interconnect?** On-node parallelism is cheap; the test of a distributed incompressible solver is what happens to its global pressure solve when halo exchanges and reductions start crossing nodes. 2. **What does the same code deliver on CPU and on GPU**, from a workstation to a national supercomputer (Snellius: Genoa nodes with 192 cores; H100 nodes with 4 GPUs)? 3. **How does it compare with leading open-source codes** on the identical case? We run three references, chosen to bracket the design space: - [**CaNS**](https://github.com/CaNS-World/CaNS) — the state of the art in raw speed for *trivial geometry*: explicit RK3 + an exact FFT Poisson solve, possible only on boxes. - [**incflo**](https://github.com/AMReX-Fluids/incflo) — the DOE-exascale (AMReX) solver in peclet's *own method family*: finite-volume incompressible NS, geometric-multigrid projections, embedded-boundary capable. - **OpenFOAM** (ESI v2412) — the ubiquitous general-purpose unstructured code, `icoFoam` + GAMG. ## The instrument: a boundary-free Taylor–Green box Scaling studies want a case that grows without complications. We tile a triply periodic box with **Taylor–Green vortices** (one vortex per $64^3$ tile, $\mathrm{Re}_\text{tile}=100$): no walls, every cell does identical work, and the domain extends in any direction by adding tiles. The TGV initial condition is also a Stokes eigenfunction, giving the driver a built-in analytic correctness check (matched to $2\times10^{-3}$, identical across rank counts and backends). [`tgv_bench.py`](tgv_bench.py) runs the case distributed and reports warmup-excluded per-step wall time with a per-phase breakdown (predictor / momentum solve / pressure projection), the pressure iteration count, and the accumulated time and count of every global reduction inside the pressure solve — so algorithmic effects and communication effects are never conflated. **Fairness rules, applied symmetrically.** Every code runs the identical grids at matched per-step tolerances (pressure ~1e-5/1e-6 class; CaNS's FFT solve is exact); timings are steady-state per-step wall times, excluding each code's setup (meshing, FFT plans, grid initialization). References run *expert-tuned*, not default: incflo gets `max_grid_size=64` and MLMG tolerances matched to the study (its AMReX defaults — $32^3$ box shredding and ~1e-11 tolerances — cost it a measured 1.53×, which we refuse to charge it). OpenFOAM uses GAMG with a 64-bit-label build (required above ~2×10⁸ cells). CaNS runs its recommended adaptive CFL 0.95 — its RK3 step advances ~4× more physical time per step than the fixed-CFL codes, which per-step throughput comparisons deliberately do not penalize. On the cluster, single measurements are noise: production points are pinned (`--cpu-bind=cores`), and repeated draws then reproduce to within 5–10 % — against up to 2× node-set variability measured for *unpinned* runs, so pinning is part of the method, and differences below ~10 % between single draws should not be over-read. ## Workstation (24-core Threadripper PRO 5965WX, RTX 5080) ### Hybrid MPI × OpenMP mix ![Left: throughput by MPI ranks × OpenMP threads on 24 cores (fixed $192^3$ TGV). Right: the per-phase step decomposition.](workstation_mix.png) **The mix barely matters: every configuration lands within 14 %** (16.4 Mcell/s at 4 ranks × 6 threads, down to 14.4 for pure MPI or pure OpenMP), a 12.6× speedup over serial on 24 cores. The phase split shows why: the momentum and pressure sweeps are memory-bandwidth bound, so once all memory channels are fed it makes little difference whether the parallelism comes from ranks or threads. Mid-hybrid mixes (4×6, 12×2) edge ahead by avoiding both extremes — 24 ranks' worth of halo surface on one side, a single rank's thread-synchronisation on the other. The same flat landscape reappears on Genoa (below). ### Weak scaling within one node: four codes, one wall ![Weak scaling at 2.1 M cells/rank, pure MPI, up to 24 cores. Left: efficiency vs. own 1-rank baseline. Right: absolute throughput (log scale); the dotted line is the single RTX 5080 running the same case.](workstation_weak.png) With per-rank work fixed, all four codes decay together — and this is not communication (peclet's pressure iteration count is pinned at 4.0 at every rank count). A single node's weak scaling measures **shared memory-bandwidth saturation**, nothing else: the efficiency ordering (incflo 43 %, peclet 42 %, CaNS 31 %, OpenFOAM 29 % at 24 ranks) is the *inverse* of per-core bandwidth demand — the codes that stream least per core-second saturate the shared memory system last. Absolute throughput on 24 cores: **CaNS 41.4 Mcell/s, peclet 13.9, incflo 4.4, OpenFOAM 3.3** — peclet is 4.2× OpenFOAM and 3.2× incflo at identical settings, while CaNS's remaining 3× lead is the price of generality made visible: an exact FFT solve that exists *only because* the box is trivial, against machinery that handles arbitrary cut-cell geometry at large d$t$. **And the single RTX 5080 at 118 Mcell/s outruns CaNS's entire 24-core node by 2.9×** — the design point of the suite: the CPU path is a capable oracle, the GPU path is the product. ### Single-GPU ladder (RTX 5080, 16 GB): three codes on the same device Both GPU-capable references run here from the same sources as their CPU curves: incflo via its AMReX CUDA backend, CaNS via its OpenACC + cuDecomp backend (NVIDIA HPC SDK toolchain, cuFFT). OpenFOAM has no native GPU solver path. All three at their GPU-fair settings. | grid | cells | peclet | incflo | CaNS | |---|---|---|---|---| | $128^3$ | 2.1 M | 99.7 | 16.1 | **131.6** | | $160^3$ | 4.1 M | 104.7 | 17.9 | **136.8** | | $192^3$ | 7.1 M | 118.1 | 18.6 | **138.3** | | $208^3$ | 9.0 M | 100.1 | 18.7 | **129.7** | *(Mcell/s, higher is better.)* Two clean results: - **peclet is 6.3× incflo on the same GPU** — a far wider margin than on CPU, the GPU-first design paying off where it is meant to. - **CaNS keeps the single-device crown by just 17 %** (138 vs 118 at $192^3$) — with no communication to pay, the exact FFT solve sets the throughput ceiling, and a general-geometry cut-cell multigrid comes within a sixth of it. What happens to that lead once the FFT's all-to-all meets an interconnect is the multi-GPU question — measured below. peclet holds ~100–118 Mcell/s across the whole ladder until 16 GB fills, iteration counts flat at 4–5 regardless of grid divisibility. ## Does the domain decomposition care about awkward rank counts? Sweeping prime and composite rank counts on a fixed $192^3$ box (pure MPI, 5965WX): | ranks | ms/step | Mcell/s | pressure iters | max-block imbalance | |---|---|---|---|---| | 5 | 1765 | 4.0 | 4.0 | 1.04× | | 7 (prime) | 1331 | 5.3 | 4.0 | 1.02× | | **8** | 1160 | 6.1 | 4.0 | 1.00× | | 13 (prime) | 795 | 8.9 | 4.0 | 1.08× | | **16** | 665 | 10.6 | 4.0 | 1.00× | | 23 (prime) | 544 | 13.0 | 4.0 | 1.20× | **No.** The pressure iteration count is exactly 4.0 at every rank count — prime or not — because the multigrid's convergence is decomposition-independent by construction. What remains is the ORB's geometric imbalance, at worst 20 % of the max block at np = 23, with throughput on the same trend line as the neighbouring powers of two. Users can pick rank counts by hardware, not by numerology. ## Cross-node CPU weak scaling on Snellius Genoa The central experiment: 188.7 M cells *per node* (768×640×384 each), grown with the node count on Genoa (192 cores/node), every code on the identical case. A single Genoa node peaks at **75 Mcell/s** (48 ranks × 4 threads; every hybrid mix from 12×16 to 192×1 lands within 20 %, iterations pinned at 4.0 throughout), 5× the whole 24-core workstation. peclet's cross-node curves run in two hybrid flavors — the classic 96×2 and the fat-rank 12×16, which uses **8× fewer MPI ranks** and therefore suffers less straggler skew at the collectives; the references run their natural pure-MPI mode. ![Cross-node weak scaling on Snellius Genoa, 188.7 Mcells/node fixed. Left: efficiency vs own 1-node baseline. Right: aggregate throughput; the dotted line is ideal scaling at peclet's slope.](snellius_weak.png) | nodes | peclet 96×2 | peclet 12×16 | CaNS | incflo | OpenFOAM | |---|---|---|---|---|---| | 1 | 63 Mcell/s | 67 | **109** | 16 | 9.4 | | 2 | 134 (106 %) | 109 (82 %) | 20 (9 %) | 32 (99 %) | — | | 4 | 259 (103 %) | 225 (84 %) | 153 (35 %) | 67 (103 %) | — | | 8 | 437 (87 %) | **459 (86 %)** | 45 (**5 %**) | 117 (90 %) | — | Three structural results: 1. **The geometric-multigrid codes scale; the FFT code does not.** peclet and incflo both hold near-flat cross-node curves (single-draw wiggles above 100 % are node-set noise) — halo exchanges plus a handful of latency-bound reductions per iteration survive the interconnect. At 8 nodes the **fat-rank 12×16 flavor wins outright** (459 vs 437 Mcell/s) with visibly less collective-wait time (249 vs 382 ms/step) — fewer ranks, less straggler skew, exactly as the skew analysis predicts. 2. **CaNS's single-node crown does not survive the interconnect.** 1.6× faster than peclet's best flavor on one node — and by two nodes it has collapsed to 9 % efficiency (reproduced in independent runs), at eight nodes to 5 %: **eight CaNS nodes deliver well under one peclet node's throughput.** Its FFT solve requires global all-to-all transposes every substep — nearly free in shared memory, ruinous across nodes. (Caveat: CaNS ran its automatic pencil-grid choice; expert `dims` tuning may soften, but cannot remove, the all-to-all.) **The crossover lands before two nodes.** 3. **Within the method family, peclet leads on speed at equal scaling**: ~4× incflo per node with the same flat curve, and 7× OpenFOAM on its single-node point. (OpenFOAM's multi-node points are limited by its serial mesh generation at these sizes — an operational cost of the unstructured design that the structured codes don't have.) ## Multi-GPU weak scaling: 1–32 H100, three codes The same case at 47.2 M cells per GPU (4 H100 per node), 1 → 32 GPUs — from one device to eight nodes, with both GPU-capable references on the identical grids: ![Left: GPU weak scaling, 47.2 Mcells/GPU fixed — peclet at 90 % efficiency and 7.7 Gcell/s aggregate on 32 H100s. Right: the pressure-solver ablation at the inter-node points (8 and 16 GPUs).](snellius_gpu.png) | H100 GPUs | peclet | incflo | CaNS | |---|---|---|---| | 1 | 269 Mcell/s | 61 | **609** | | 4 | 998 (92 %) | 225 (92 %) | 856 (35 %) | | 8 | 1968 (92 %) | 429 (88 %) | 297 (7.6 %) | | 16 | 3936 (92 %) | 848 (85 %) | 977 (10 %) | | 32 | **7748 (90 %)** | 1665 (85 %) | 1032 (5.9 %) | The GPU fabric replays the CPU story, at higher stakes. **One H100 runs peclet at 269 Mcell/s — four Genoa nodes' worth — and 32 H100s sustain 7.7 Gcell/s at 90 % weak efficiency**: a 1.5-billion-cell grid advancing five time steps per second, with iterations pinned at 4.0 and the entire global-reduction budget under 2 % of the step. incflo scales just as flatly (85 %) at 4.4× less per-GPU speed — the method family's signature. And CaNS repeats its arc precisely: the single-device crown (609 Mcell/s, 2.3× peclet — cuFFT with nothing to communicate), already sagging *within* one NVLink node (35 % at 4 GPUs), collapsing at the first node boundary (7.6 %), and flatlining near 1 Gcell/s while peclet pulls away 7.5×. **The crossover sits between 4 and 8 GPUs.** (CaNS-GPU caveat: cuDecomp ran its autotuned MPI transpose backends — the NVHPC-bundled NCCL errors on this fabric — which particularly understates its intra-node point.) The ablation on the right is how the solver's defaults are chosen — measured, at scale, and re-validated after every solver change: - **MG-PCG with fine-scoped nullspace removal** (the default) wins at both 8 and 16 GPUs (191.8 ms/step). - Removing the mean on every multigrid level instead costs 4 % — pure reduction latency. - The host-staged halo costs 16 % vs GPU-aware MPI (device pointers straight to the interconnect). - **Chebyshev** comes within 15 % (5 iterations vs PCG's 4) — a respectable reduction-free fallback, but PCG's convergence still beats its communication savings at every scale measured. - The agglomerated algebraic bottom solve (GraphAMG) pays far more than it saves on GPU. ## What the study says, in one paragraph On the identical boundary-free case, `peclet.flow` is the fastest general-geometry incompressible solver in the comparison on every platform tested — 4.2× OpenFOAM and 3.2× incflo per CPU node on the workstation, ~4× incflo per Genoa node, 6.3× incflo on the same GPU — and, together with incflo, the only code family whose weak scaling survives the interconnect: ~87 % across 8 CPU nodes and 90 % across 32 H100s, where the FFT state of the art collapses below 10 % on both fabrics. CaNS keeps the single-device crowns (by 3× on a CPU node, 17 % on a consumer GPU, 2.3× on one H100) — the honest price of general geometry against an exact solve that only exists on trivial boxes — and loses them the moment the problem spans two of anything: **the crossover lands before two nodes on CPU and between 4 and 8 GPUs on H100**. The same source runs 118 Mcell/s on a consumer RTX 5080 (2.9× CaNS's entire 24-core node), 269 Mcell/s on one H100, and **7.7 Gcell/s at 90 % weak efficiency on 32 H100s** — a 1.5-billion-cell incompressible flow advancing five steps per second. ## Reproducing All drivers, case definitions, runner scripts for the three reference codes, the Snellius job pack, raw result JSONs and the figure scripts live in [`benchmarks/parallel-scaling/`](https://github.com/computational-chemical-engineering/peclet-examples/tree/main/benchmarks/parallel-scaling) — `tgv_bench.py` (peclet), `run_cans.sh`, `run_incflo.sh`, `openfoam-tgv/`, `snellius/README.md` (exact submission commands), `plot_workstation.py` / `plot_snellius.py`.