Euler–Lagrange two-way coupling with moving bubbles gives rank-count-dependent results

[wilfonba]

Euler–Lagrange two-way coupling gives rank-count-dependent results — investigated: round-off seed amplified by chaos, not a decomposition bug

Summary

Moving Lagrangian bubbles with two-way coupling (lag_params%vel_model = 2, lag_params%solver_approach = 2) produce visibly different results when the same case is run with different MPI rank counts. A systematic bisection campaign (attached tarball) traced this to a one-ulp (2.2e-16) summation-order difference in the two-way smearing deposition wherever a bubble's kernel straddles a subdomain boundary — the halo-accumulate reduce in s_mpi_reduce_beta_variables_buffers sums the same contributions in a decomposition-dependent order. The strongly chaotic coupled-bubble dynamics then amplify that seed exponentially (measured e-folding time ≈ 0.38 nondim time units) to O(1) trajectory differences within a few hundred steps.

Verdict: there is no decomposition bug in the moving-bubble / two-way-coupling path. Every mechanism tested in isolation (Eulerian solver, one-way coupling, kernel deposition away from boundaries, particle migration across ranks) is bitwise decomposition-invariant. The rank-count discrepancy is round-off seeded and chaos-amplified — numerically benign, but it means trajectory-wise cross-rank validation of these cases is fundamentally impossible.

Background

Observed symptom: examples/3D_moving_lag_particles-style cases (moving Keller–Miksis bubbles, two-way coupled) diverge between e.g. 1-rank and 8-rank runs. Because a cloud of coupled oscillating bubbles is chaotic, any machine-precision difference grows exponentially — so divergence alone is not evidence of a bug. The campaign was designed to distinguish:

• a real decomposition bug (first-step differences ≫ round-off, or divergence growing faster than the system's own chaos rate, or O(1) jumps at migration events), from
• benign round-off seeding (first differences at 1 ulp, growing exactly at the independently measured chaos rate).

Method

Self-contained campaign (attached tarball; run 2026-07-05 on wingtip-gpu3, nvfortran 25.11, CPU+MPI and OpenACC builds). All code the campaign exercises (lag_params%vel_model, solver_approach, kahan_summation) is in upstream master as of PR #1290 (Moving EL Bubbles with MPI Decomposition); the campaign was originally run on the pre-merge development branch. Diagnostic-scale configuration mirroring the failing case: 50³ cells, 64 Keller–Miksis bubbles, tri-periodic box, central 10× pressure sphere, thermal = 3, polytropic = F.

Two methodology upgrades over eyeballing ParaView output:

• Bitwise field comparison. compare_fields.py reads the parallel-IO restart files (lustre_<step>.dat) directly — raw doubles in a global layout identical across rank counts — so runs are compared bit-for-bit and the worst-differing cell is located. compare_particles.py parses the native lag_bubble_evol_<rank>.dat time series and fits divergence growth rates.
• A control for chaos. Phase 2 measures the system's intrinsic sensitivity by perturbing bubble ICs by 1e-14 at fixed rank count, giving a baseline exponential growth rate to compare the 1-rank-vs-8-rank divergence against.

Phases:

phase	question
0	Is the same config bitwise reproducible run-to-run? (CPU 1r/8r, GPU 1r/4r)
1	Feature bisection M0 (no bubbles) / M1 (static + two-way) / M2 (moving + one-way) / M3 (moving + two-way), each at 1/2/8 ranks, per-step saves
2	Chaos baseline (1e-14 IC perturbation at fixed ranks) vs rank divergence (1r vs 8r), 500 steps
3	Single-bubble forensics: kernel interior vs straddling a rank boundary vs migrating across ranks mid-run
—	Kahan probe: does lag_params%kahan_summation remove the seed?

Evidence

test	result
P0: same config run twice (CPU 1r, CPU 8r, GPU 1r, GPU 4r)	all bitwise identical — comparisons well-posed, no GPU-atomics noise (deposition is a cell-centric gather)
P1 M0: no bubbles, 1 vs 2 & 1 vs 8 ranks, 20 steps	bitwise identical — Eulerian solver is decomposition-exact
P1 M2: moving + one-way	bitwise identical, fields and particle trajectories
P1 M1: static + two-way	16–65 cells differ by 1 ulp at step 1; ~4e-14 by step 20
P1 M3: moving + two-way	same seed as M1 (1 ulp); ~4e-13 by step 20 (coupling feedback)
P3 interior: single bubble, kernel away from boundary, two-way, 100 steps	bitwise identical — the deposition itself is exact
P3 straddle: single bubble on the 2-rank boundary	1-ulp seed the moment the kernel touches the boundary; ~1e-13 after 100 steps
P3 migrate: single bubble crossing ranks mid-run	first diff (1e-16) when the kernel reaches the boundary, no jump at the crossing — migration handoff exact to round-off
P2 chaos baseline: 8r vs 8r with 1e-14 IC perturbation, 500 steps	smooth exponential growth, e-folding time ≈ 0.38; punctuated burst once amplitude reaches ~1e-8
P2 rank divergence: 8r vs 1r, identical ICs	starts at 1 ulp, same growth character and rate, bursts at the same ~1e-7 amplitude, saturates at O(1) by t ≈ 7
Kahan probe: M1 with lag_params%kahan_summation = T	identical 1-ulp seed — Kahan compensation cannot equalize different summation orders

(Plots: p2_chaos_baseline.png and p2_rank_divergence.png in the tarball — the two divergence curves are visually indistinguishable in growth rate and burst character.)

The Phase 2 "bursts" occur at a divergence amplitude threshold (~1e-8–1e-7), not at a fixed time: collapse/rebound timing shifts by one step once the perturbation is large enough. The single-bubble migrate test rules out a migration kick as the cause.

Practical implications

• Cross-rank-count validation of two-way EL cases must be statistical (integrals, spectra, ensemble envelopes), not trajectory-wise: trajectories decorrelate at the chaos rate (e-folding ≈ 0.4 nondim time here) from an irreducible 1-ulp seed.
• Bitwise decomposition invariance would require a canonical-order reduction in s_mpi_reduce_beta_variables_buffers (accumulate straddling-kernel contributions in a fixed global order) — a performance/complexity trade, not a correctness fix.
• One-way coupling (solver_approach = 1) is bitwise decomposition-invariant and remains a good regression anchor.

Reproducing

The attached tarball is the complete self-contained campaign; it is not part of the upstream tree, so unpack it as examples/3D_lagrange_coupling_matrix/ in an MFC checkout (master at or after PR #1290). From that directory:

./run_all.sh          # builds + all phases + Kahan probe (~1.5–2 h, ~20 GB)
GPU=0 ./run_all.sh    # CPU only: skips the GPU build and GPU determinism pairs

run_all.sh writes phase{0,1,2,3}.log and kahan.log and ends with a verdict summary to check against FINDINGS.md. Phases run individually (./run_phase1_matrix.sh, …), and any single matrix cell by hand:

./run_case.sh out_test 8 --vel-model 2 --solver-approach 2 --steps 50 --save-every 5
python3 compare_fields.py out_p1_M3_1rank out_test      # per-step field diff
python3 compare_particles.py out_p1_M3_1rank out_test   # trajectory divergence + growth fit

Note: the scripts source a machine-local compiler environment (/fastscratch/bwilfong3/NVEnv.sh) — edit run_case.sh and run_all.sh for another machine. Runs are sequential by design (they share the case directory).

Tarball contents

el_rank_divergence_campaign.tar.gz

file	purpose
README.md	campaign design, decision tree, per-phase runtimes
FINDINGS.md	results and verdict (reproduced above)
case_matrix.py	one parameterized case; every physics axis is a CLI flag
run_all.sh	full campaign: builds, phases 0–3, Kahan probe, verdict summary
run_case.sh	run one cell (<outdir> <nranks> [flags]), snapshot outputs
run_phase0_determinism.sh	same config twice → must be bitwise identical
run_phase1_matrix.sh	M0–M3 feature bisection × {1,2,8} ranks
run_phase2_chaos.sh	1e-14 IC perturbation rate vs 1-vs-8-rank divergence rate
run_phase3_boundary.sh	single bubble on/off the rank boundary + migration
compare_fields.py	restart-file diff (bitwise-capable, locates worst cell)
compare_particles.py	per-particle divergence over time + growth-rate fit
p2_chaos_baseline.png, p2_rank_divergence.png	Phase 2 divergence plots

---

AI/LLM Disclosure: This report and the attached code and scripts were generated with the assistance of Claude Code