Porous-bed parallel scaling

Stokes permeability through DEM-grown sphere packings — cut-cell and ghost-cell IBM, workstation to 32 H100

Weak-scaling study of the full IBM pressure path on a physical porous medium: spheres grown to a target solid fraction with peclet.dem (XPBD growth mode, periodic box), then creeping flow through the bed with peclet.flow, cut-cell and ghost-cell IBM. Companion to the all-fluid TGV study; where TGV measures the solver floor, this measures what a real geometry costs — and adds a physics observable, the bed permeability $k/R^2$, whose convergence the ladders expose.

Pipeline

pack_bed.py — grow $N$ spheres ($R{=}1$) in a periodic box to solid fraction $\phi = 0.50$; the artifact is the sphere list (centers + scales, a small npz), not an SDF.
spheres_bench.py — every MPI rank samples the analytic union-of-spheres SDF over its own block (periodic images included), so one packing serves any grid resolution and nothing global is gathered. Phase A measures per-step cost (warmup-excluded, per-phase breakdown, pressure iterations, allreduce time); phase B marches to steady state and records $k/R^2$.

Ladders

Upscale (snellius/spheres_weak_gpu.sh): 256³ cells/GPU held fixed, the domain — and the packing in it — doubles with the GPU count (1→32 GPUs, 16.8M→537M cells, 489→15.6k spheres, sphere radius fixed at 16 cells, per-rung seeds). Blocks tile exactly, so the decomposition is perfectly balanced at every rung. Permeability should be domain-size independent here (an RVE check); iterations/step should be flat — the agglomerated coarse solve (set_pressure_bottom auto) handles the growing global bottom exactly.

Refine (snellius/spheres_refine_gpu.sh): one fixed physical bed (16³ R-units, seed 100), grid refined 256³→1024³ with the GPU count; the sphere radius spans 16→64 cells. Cells/GPU varies (3-D refinement can’t hold it), so weak efficiency reads from per-GPU throughput. The physics payoff: $k(N) \to k_\infty$ per IBM with its observed order — and the cut-cell and ghost-cell estimates of $k_\infty$ must agree.

Levers (at the largest allocation): legacy smoothed bottom vs the agglomerated default, GPU-aware vs host-staged MPI.

Results (Snellius gpu_h100, 2026-08)

Upscale. Cut-cell pressure iterations are flat across the ladder (21–33/step, 1→32 GPUs, 16.8M→537M cells) — the agglomerated coarse solve + exactly-tiling decomposition deliver the grid-independent solve on a real porous geometry. Per-GPU throughput drops 36→13 Mcell/s from 1→8 GPUs (halo + allreduce cost appearing), then holds flat 8→32: the scaling loss is a fixed communication toll, not a growing one. Permeability is domain-size independent as an RVE should be: every converged rung sits in $k/R^2 = 0.0108$–$0.0112$ (per-rung random beds; the scatter is bed statistics, not solver noise). Ghost-cell tracks cut-cell physics on the rungs that converge, at roughly 4× the iterations and ~⅓ the throughput.

Halo–compute overlap + communication-avoiding smoothing (solid vs dashed above; the communication-scaling work of suite/docs/COMMUNICATION_SCALING.md, flow 1e9c5db + b96dd74, 2026-08). The baseline’s fixed toll was diagnosed as latency-bound halo events — ~700 per step (one per smoother colour per MG level per iteration), not bandwidth and not allreduces. Two bit-exact changes attack the event count: every distributed RB-GS sweep now overlaps its exchange with the interior sweep, and communication-avoiding smoothing (PECLET_FLOW_CA) exchanges a 2-deep ghost layer once per red–black pair, re-smoothing the 1-deep ghost ring so the second colour needs no exchange at all. The result at np=32 cut-cell: 1323 → 759 ms/step (12.7 → 22.1 Mcell/s/GPU), lifting weak efficiency 35% → 62% — with pressure iterations and $k/R^2$ identical to every printed digit at every rung, the numerics are provably untouched. The ablation point (×) shows overlap alone at 13.8 Mcell/s/GPU: the win is almost entirely the CA halving of event count, confirming the latency diagnosis. The ghost IBM gains the same way where its march is stable (np=16: 6966→4281 ms/step).

One anomaly, traced by a np=8 attribution matrix (PECLET_FLOW_CA=mom|mg, host-staged, and UCX-threshold runs; all numerically bit-identical): the np=8 rungs’ projection phase is slower under CA (cut-cell 342→465 ms) even as momentum improves. It is reproducible, bisects entirely to the MG coarse-level CA (CA=mom leaves projection at baseline while keeping the momentum win), and is transport-independent (host-staged MPI shows it too; forcing a large UCX_RNDV_THRESH makes it 2× worse — not a protocol-threshold effect). Reading: MG-CA carries a small per-V-cycle overhead that every rung pays; np=8 — two nodes, with baseline halo events half the price of np=16’s — has nothing for the event-halving to save, so only the overhead shows, while np=16/32 gain 40–47 %. Practical note: cut-cell full CA is still net faster at np=8 (988 vs 1239 ms/step); for the projection-dominated ghost IBM at ≤2 nodes, PECLET_FLOW_CA=mom is the better setting (4917 vs 5312 baseline vs 6337 full-CA ms/step).

Refine. The physics headline: on one fixed bed refined 256³→1024³ (sphere radius 16→64 cells), cut-cell $k/R^2$ is already converged at the coarsest grid — flat at $0.01085$ to four digits — while ghost-cell approaches the same value from above at observed order ≈1.9. The two IBMs’ independent $k_\infty$ estimates agree to 0.1% ($0.01086$ vs $0.01085$) — a strong cross-validation of both discretizations. Open markers are marches that hit the 400-step cap on the plateau (1024³); they sit on the fit.

Levers (32 GPUs, 537M cells): host-staged MPI is ~7% faster per step than GPU-aware UCX (1236 vs 1323 ms) at this size — the GPU-aware path buys nothing here; and the smoothed vs agglomerated bottom is a wash (24.2 iters/step both) at this box/depth — the agglomerated bottom’s win lives in deeper hierarchies and long boxes (see the channel-scaling study), while costing nothing here.

Collocated schemes on the fixed beds (colcmp / colcmp060, 2026-08)

The same two fixed beds carry the collocated-scheme A/B (snellius/collocated_ab_gpu.sh -> analyze_collocated_ab.py results/snellius-h100 [colcmp060]), R = 5..32 cells/radius, staggered cut-cell as the same-rung reference. Error in $k/R^2$ against the staggered reference and the finest-rung cost:

scheme	φ=0.50: R=8 → 32 err%	φ=0.60: R=8 → 32 err%	p.iters (R32?)	ms/step (R32?)
staggered cut-cell (ref)	—	—	37	1185
collocated `plain` (mode 0)	−4.52 → −1.23 (order ~1)	−7.19 → −1.72 (order ~1)	45	1335
collocated `gauge-exact` (mode 9)	−0.77 → +0.29	−2.25 → +0.39	43	1285
collocated `ghost`	−0.34 → +0.14	−1.43 → +0.22	92	4260

Notes:

Interpretation changed 2026-08 (read flow/doc/collocated_invisible_subspace.md): the gauge-exact “+0.3–0.4 % asymptote” is not a truncation ceiling — it is the bias of an attractor family of steady states (a support-inconsistency of the gradient/constraint pair) plus, at these march settings, a rotational-update instability. The ghost scheme is family-free (unique, unconditionally stable, protocol-independent), and its asymptote here (+0.14/+0.22 %) was reproduced to three digits by the independent clean-protocol campaign — it is the ghost’s own small discretization bias, now the production-candidate trade: robustness + ~0.2 % bias at 2.3–3.6× the cost. set_collocated_scheme("ghost") selects it since flow ec693a1.
Face-flux estimator (*_flux reruns): the flux-based $k$ agrees with the cell-mean $k$ to −0.11 % (φ=0.50) / −0.22 % (φ=0.60) at R=32, converging ~O(h²) — the two estimators measure the same limit, closing the estimator item of the refutation catalogue.
Embed variants (modes 6/7) diverged at every bed rung (k = nan): root-caused post-hoc (flow 4732b17) — a degenerate-sliver wall-drag gain at grid/surface-grazing cells, fixed by a resolution-scale floor; the repaired modes remain unstable on beds for a deeper reason (non-adjoint pressure pair) and are retired as scheme candidates. The archived nan runs are the incident record.

Incidents & open issues (found by this study)

dem packing corruption on H100 (RESOLVED 2026-08-18, dem 17288f8): a silent out-of-bounds ghost-slot write — ensureCapacity never resized materialId, so the periodic ghost generator wrote 1-byte strays into the neighbouring allocation. Harmless on the workstation allocator’s layout, catastrophic on the H100s’ (contacts silently unresolved, φ_voxel 0.40 instead of 0.50, while the solver’s own overlap metric — which only sees tracked contacts — stayed quiet). The seven-round probe ladder (snellius/probe_dem*.sh) eliminated budget, sleeping, optimizer, toolkit, and every execution mode before compute-sanitizer named the kernel; the post-fix probe-2 rerun passes all 12 configurations on H100. The independent voxel-fraction gates added at both pipeline ends contained the incident and remain as permanent defenses.
Ghost-projection march instability at scale (open, under active de-confounding 2026-08-23): the ghost IBM steady march diverges at the np=16/32 upscale rungs with per-solve convergence intact. Progress: the recorded rungs ran the PURE (2,2) closures (not the quarantined (1,2) mixed mode — verified from the JSONs’ gporder); the same 7823-sphere np=16 bed marched SINGLE-RANK at R=8 is stable (s116probe), R≥12 exceeds one H100 (the ghost overlay’s ~1.6 KB/cell); the divergence needs ≥7823 spheres at R=16 (np=8/3911 converges, refine np=32/489 converges). The remaining bed-size-vs-rank-count confounder is split by the submitted np=8-on-the-np=16-bed rung (s116np8). Tracking: flow doc/collocated_paper_plan.md row 46.
The perf-phase numbers of the diverged ghost runs are valid (divergence builds over hundreds of march steps) and are included in the scaling panels; their $k$ values are excluded.

Reproducing

# workstation shakedown (any MPI flow build; dem build for the packing)
GNX=256 GNY=256 GNZ=256 RCELLS=16 SEED=100 OUT=pack.npz python pack_bed.py
PACK=pack.npz mpirun -np 2 python spheres_bench.py

# Snellius: one job per rung (queue-parallel safe), argument = GPU count; the optional second
# argument is a result TAG (required after any solver change -- rungs skip existing JSONs)
sbatch --nodes=1 snellius/spheres_weak_gpu.sh 1
sbatch --nodes=8 snellius/spheres_weak_gpu.sh 32 overlap
sbatch --nodes=8 snellius/spheres_weak_gpu.sh levers
sbatch --nodes=1 snellius/spheres_refine_gpu.sh 1

Pre-flight any grid/rank-count change with flow/scripts/check_decomposition.py before spending queue time.

--- title: "Porous-bed parallel scaling" subtitle: "Stokes permeability through DEM-grown sphere packings — cut-cell and ghost-cell IBM, workstation to 32 H100" format: html: toc: true --- Weak-scaling study of the full IBM pressure path on a physical porous medium: spheres grown to a target solid fraction with `peclet.dem` (XPBD growth mode, periodic box), then creeping flow through the bed with `peclet.flow`, cut-cell **and** ghost-cell IBM. Companion to the [all-fluid TGV study](../parallel-scaling/index.qmd); where TGV measures the solver floor, this measures what a real geometry costs — and adds a physics observable, the bed permeability $k/R^2$, whose convergence the ladders expose. ## Pipeline 1. **`pack_bed.py`** — grow $N$ spheres ($R{=}1$) in a periodic box to solid fraction $\phi = 0.50$; the artifact is the *sphere list* (centers + scales, a small npz), not an SDF. 2. **`spheres_bench.py`** — every MPI rank samples the analytic union-of-spheres SDF over its own block (periodic images included), so one packing serves any grid resolution and nothing global is gathered. Phase A measures per-step cost (warmup-excluded, per-phase breakdown, pressure iterations, allreduce time); phase B marches to steady state and records $k/R^2$. ## Ladders **Upscale** (`snellius/spheres_weak_gpu.sh`): 256³ cells/GPU held fixed, the domain — and the packing in it — doubles with the GPU count (1→32 GPUs, 16.8M→537M cells, 489→15.6k spheres, sphere radius fixed at 16 cells, per-rung seeds). Blocks tile exactly, so the decomposition is perfectly balanced at every rung. Permeability should be domain-size independent here (an RVE check); iterations/step should be flat — the agglomerated coarse solve (`set_pressure_bottom auto`) handles the growing global bottom exactly. **Refine** (`snellius/spheres_refine_gpu.sh`): one fixed physical bed (16³ R-units, seed 100), grid refined 256³→1024³ with the GPU count; the sphere radius spans 16→64 cells. Cells/GPU varies (3-D refinement can't hold it), so weak efficiency reads from per-GPU throughput. The physics payoff: $k(N) \to k_\infty$ per IBM with its observed order — and the cut-cell and ghost-cell estimates of $k_\infty$ must agree. **Levers** (at the largest allocation): legacy smoothed bottom vs the agglomerated default, GPU-aware vs host-staged MPI. ## Results (Snellius gpu_h100, 2026-08) ![Upscale weak scaling](results/snellius-h100/porous_weak.png) **Upscale.** Cut-cell pressure iterations are *flat* across the ladder (21–33/step, 1→32 GPUs, 16.8M→537M cells) — the agglomerated coarse solve + exactly-tiling decomposition deliver the grid-independent solve on a real porous geometry. Per-GPU throughput drops 36→13 Mcell/s from 1→8 GPUs (halo + allreduce cost appearing), then holds flat 8→32: the scaling loss is a fixed communication toll, not a growing one. Permeability is domain-size independent as an RVE should be: every converged rung sits in $k/R^2 = 0.0108$–$0.0112$ (per-rung random beds; the scatter is bed statistics, not solver noise). Ghost-cell tracks cut-cell physics on the rungs that converge, at roughly 4× the iterations and ~⅓ the throughput. **Halo–compute overlap + communication-avoiding smoothing** (solid vs dashed above; the communication-scaling work of `suite/docs/COMMUNICATION_SCALING.md`, flow `1e9c5db` + `b96dd74`, 2026-08). The baseline's fixed toll was diagnosed as *latency-bound halo events* — ~700 per step (one per smoother colour per MG level per iteration), not bandwidth and not allreduces. Two bit-exact changes attack the event count: every distributed RB-GS sweep now overlaps its exchange with the interior sweep, and communication-avoiding smoothing (`PECLET_FLOW_CA`) exchanges a 2-deep ghost layer once per red–black *pair*, re-smoothing the 1-deep ghost ring so the second colour needs no exchange at all. The result at np=32 cut-cell: **1323 → 759 ms/step** (12.7 → 22.1 Mcell/s/GPU), lifting weak efficiency **35% → 62%** — with pressure iterations and $k/R^2$ identical *to every printed digit* at every rung, the numerics are provably untouched. The ablation point (×) shows overlap alone at 13.8 Mcell/s/GPU: the win is almost entirely the CA halving of event count, confirming the latency diagnosis. The ghost IBM gains the same way where its march is stable (np=16: 6966→4281 ms/step). One anomaly, **traced by a np=8 attribution matrix** (`PECLET_FLOW_CA=mom|mg`, host-staged, and UCX-threshold runs; all numerically bit-identical): the np=8 rungs' *projection* phase is slower under CA (cut-cell 342→465 ms) even as momentum improves. It is reproducible, bisects entirely to the **MG coarse-level CA** (`CA=mom` leaves projection at baseline while keeping the momentum win), and is transport-independent (host-staged MPI shows it too; forcing a large `UCX_RNDV_THRESH` makes it 2× worse — not a protocol-threshold effect). Reading: MG-CA carries a small per-V-cycle overhead that every rung pays; np=8 — two nodes, with baseline halo events half the price of np=16's — has nothing for the event-halving to save, so only the overhead shows, while np=16/32 gain 40–47 %. Practical note: cut-cell full CA is still net faster at np=8 (988 vs 1239 ms/step); for the projection-dominated ghost IBM at ≤2 nodes, `PECLET_FLOW_CA=mom` is the better setting (4917 vs 5312 baseline vs 6337 full-CA ms/step). ![Refine ladder](results/snellius-h100/porous_refine.png) **Refine.** The physics headline: on one fixed bed refined 256³→1024³ (sphere radius 16→64 cells), cut-cell $k/R^2$ is already converged at the coarsest grid — flat at $0.01085$ to four digits — while ghost-cell approaches the same value from above at observed order ≈1.9. The two IBMs' independent $k_\infty$ estimates agree to **0.1%** ($0.01086$ vs $0.01085$) — a strong cross-validation of both discretizations. Open markers are marches that hit the 400-step cap on the plateau (1024³); they sit on the fit. **Levers** (32 GPUs, 537M cells): host-staged MPI is ~7% *faster* per step than GPU-aware UCX (1236 vs 1323 ms) at this size — the GPU-aware path buys nothing here; and the smoothed vs agglomerated bottom is a wash (24.2 iters/step both) at this box/depth — the agglomerated bottom's win lives in deeper hierarchies and long boxes (see the channel-scaling study), while costing nothing here. ## Collocated schemes on the fixed beds (colcmp / colcmp060, 2026-08) The same two fixed beds carry the collocated-scheme A/B (`snellius/collocated_ab_gpu.sh` -> `analyze_collocated_ab.py results/snellius-h100 [colcmp060]`), R = 5..32 cells/radius, staggered cut-cell as the same-rung reference. Error in $k/R^2$ against the staggered reference and the finest-rung cost: | scheme | φ=0.50: R=8 → 32 err% | φ=0.60: R=8 → 32 err% | p.iters @R32 | ms/step @R32 | |---|---|---|---|---| | staggered cut-cell (ref) | — | — | 37 | 1185 | | collocated `plain` (mode 0) | −4.52 → −1.23 (order ~1) | −7.19 → −1.72 (order ~1) | 45 | 1335 | | collocated `gauge-exact` (mode 9) | −0.77 → **+0.29** | −2.25 → **+0.39** | 43 | 1285 | | collocated `ghost` | −0.34 → **+0.14** | −1.43 → **+0.22** | 92 | 4260 | Notes: - **Interpretation changed 2026-08 (read `flow/doc/collocated_invisible_subspace.md`):** the gauge-exact "+0.3–0.4 % asymptote" is not a truncation ceiling — it is the bias of an *attractor family* of steady states (a support-inconsistency of the gradient/constraint pair) plus, at these march settings, a rotational-update instability. The ghost scheme is family-free (unique, unconditionally stable, protocol-independent), and its asymptote here (+0.14/+0.22 %) was *reproduced to three digits* by the independent clean-protocol campaign — it is the ghost's own small discretization bias, now the production-candidate trade: robustness + ~0.2 % bias at 2.3–3.6× the cost. `set_collocated_scheme("ghost")` selects it since flow `ec693a1`. - **Face-flux estimator** (`*_flux` reruns): the flux-based $k$ agrees with the cell-mean $k$ to −0.11 % (φ=0.50) / −0.22 % (φ=0.60) at R=32, converging ~O(h²) — the two estimators measure the same limit, closing the estimator item of the refutation catalogue. - **Embed variants (modes 6/7) diverged at every bed rung** (`k = nan`): root-caused post-hoc (flow `4732b17`) — a degenerate-sliver wall-drag gain at grid/surface-grazing cells, fixed by a resolution-scale floor; the repaired modes remain unstable on beds for a deeper reason (non-adjoint pressure pair) and are retired as scheme candidates. The archived nan runs are the incident record. ## Incidents & open issues (found by this study) - **dem packing corruption on H100** (**RESOLVED** 2026-08-18, dem `17288f8`): a silent out-of-bounds ghost-slot write — `ensureCapacity` never resized `materialId`, so the periodic ghost generator wrote 1-byte strays into the neighbouring allocation. Harmless on the workstation allocator's layout, catastrophic on the H100s' (contacts silently unresolved, φ_voxel 0.40 instead of 0.50, while the solver's own overlap metric — which only sees tracked contacts — stayed quiet). The seven-round probe ladder (`snellius/probe_dem*.sh`) eliminated budget, sleeping, optimizer, toolkit, and every execution mode before `compute-sanitizer` named the kernel; the post-fix probe-2 rerun passes all 12 configurations on H100. The independent voxel-fraction gates added at both pipeline ends contained the incident and remain as permanent defenses. - **Ghost-projection march instability at scale** (open, under active de-confounding 2026-08-23): the ghost IBM steady march diverges at the np=16/32 upscale rungs with per-solve convergence intact. Progress: the recorded rungs ran the PURE (2,2) closures (not the quarantined (1,2) mixed mode — verified from the JSONs' `gporder`); the same 7823-sphere np=16 bed marched SINGLE-RANK at R=8 is stable (`s116probe`), R≥12 exceeds one H100 (the ghost overlay's ~1.6 KB/cell); the divergence needs ≥7823 spheres *at R=16* (np=8/3911 converges, refine np=32/489 converges). The remaining bed-size-vs-rank-count confounder is split by the submitted np=8-on-the-np=16-bed rung (`s116np8`). Tracking: flow `doc/collocated_paper_plan.md` row 46. - The perf-phase numbers of the diverged ghost runs are valid (divergence builds over hundreds of march steps) and are included in the scaling panels; their $k$ values are excluded. ## Reproducing ```bash # workstation shakedown (any MPI flow build; dem build for the packing) GNX=256 GNY=256 GNZ=256 RCELLS=16 SEED=100 OUT=pack.npz python pack_bed.py PACK=pack.npz mpirun -np 2 python spheres_bench.py # Snellius: one job per rung (queue-parallel safe), argument = GPU count; the optional second # argument is a result TAG (required after any solver change -- rungs skip existing JSONs) sbatch --nodes=1 snellius/spheres_weak_gpu.sh 1 sbatch --nodes=8 snellius/spheres_weak_gpu.sh 32 overlap sbatch --nodes=8 snellius/spheres_weak_gpu.sh levers sbatch --nodes=1 snellius/spheres_refine_gpu.sh 1 ``` Pre-flight any grid/rank-count change with `flow/scripts/check_decomposition.py` before spending queue time.