Repository navigation
Immersed-boundary force gains a large spurious transverse component on multiple ranks #1863
Description
Investigated with Claude Code.
What happens
A static sphere at Re 100, Ma 0.1, in uniform flow along x, returns a transverse force larger than its drag on a flow whose wake is axisymmetric. Grid 405×270×270, 64 GPU ranks, body at the center of the domain:
Fx +0.4040 Fy +0.4364 Fz +0.4364
Fy and Fz agree to five digits, and the torques mirror them (Ty −0.028511, Tz +0.028507), which is what a spatially uniform spurious force density in y and z produces through Σ r × F.
It is not a flow quantity
Fy is 0.4367 at time step 1, on a field that is still essentially uniform, and 0.4364 at step 49,793, while Fx relaxes from 1.638 to 0.4040 over those same 50,000 steps. Constant to four significant figures while the flow develops completely.
Where it comes from
Restarting the same field, with the same binary, changing only the rank count:
| ranks | Fx | Fy | Fz |
|---|---|---|---|
| 64 | +0.404403 | +0.436447 | +0.436433 |
| 1 | +0.404403 | +0.000024 | +0.000010 |
Fx is identical to six digits. The transverse force collapses by four orders of magnitude. s_communicate_ib_forces returns immediately when num_procs == 1, so the single-rank path skips the reduction entirely — as does every multi-rank effect on the ghost-cell fill.
Independently, reproducing s_compute_ib_forces' volume sum offline from the restart file — same fd_order = 4 stencil, same body mask, pressure and viscous terms — gives Fx +0.40459 against MFC's +0.40399 (0.15 %) and Fy +0.00002. So the per-cell sum itself is right and carries no transverse force; the spurious part enters somewhere between that sum and what lands in patch_ib%force.
That leaves s_communicate_ib_forces and the multi-rank ghost-cell fill as the candidates. The reduction is a hand-rolled dimension-by-dimension prefix sum with a hop count capped at min(2*ib_neighborhood_radius, num_procs_<dir> - 1), which truncates whenever a direction carries more than 2*ib_neighborhood_radius + 1 ranks.
Why it has not been caught
Every immersed-boundary golden file is single-rank, where s_communicate_ib_forces returns on its first line. The force is diagnostic for a static body, so a wrong value is silent; for moving_ibm it drives the motion.
Reproducing
Any static, symmetric body on more than one rank, positioned so it spans a subdomain boundary, with the force written per step (ib_state_wrt). The symmetric body is what makes it visible: the true transverse force is zero, so anything non-zero is the defect. examples/2D_ibm_force_decomposition (on #1859) is the 2D version of that setup and runs in seconds; it shows the smaller, related decomposition dependence in drag that #1859 fixes, which is a different bug — rebuilding with that fix leaves the sphere's 0.4364 unchanged.
Not the same as #1860
#1860 is an out-of-bounds fd_coeff read that makes the force depend on the decomposition at the 0.5 % level. It is real and fixed by #1859, and it is not this. Verified by rebuilding with that fix and re-running the sphere from its own restart: Fy stays at 0.4364.
Activity
sbryngelson commentedon Sep 12, 2026
A second, much cheaper case shows the same thing, and it is 2D — so this is not a 3D or a stretched-grid effect.
A static cylinder at Re 40, Ma 0.1, uniform 1000×1000, body at the domain centre. Same case, same binary, only the rank count differs. At the matched time t = 5.001:
The cylinder is symmetric about y = 0, so the true lift is zero. It falls by a factor of 63 when the reduction is skipped.
And on 8 ranks it is a pedestal rather than a transient. C_L over the full run:
Flat to ±0.003 across the whole run while the drag relaxes by 17 %. Same signature as the sphere: a constant, layout-dependent transverse force that ignores the flow entirely.
This one costs a single GPU-hour to reproduce, so it is a better regression case than the sphere. The practical consequence for anyone using IB forces: a spurious lift of order 5 % of drag on a body that has none, silently, on any multi-rank run, and it does not average away.
sbryngelson commentedon Sep 12, 2026
Narrowed to a precise trigger. Same restart field, same binary, same case — only the rank count changes, which changes only the decomposition:
F_xis identical to six digits in all three. The defect needs three or more ranks in a direction, i.e. more than one hop ins_communicate_ib_forces' sweep for that direction, and it corrupts the force component along that direction.The 4×2×2 row is what makes this specific. It runs three hops in x and one in y and z — and it is clean in every component, including x. So multi-hop alone is not sufficient. The distinguishing fact is where the body sits: with this grid the subdomain splits land at
So the trigger is a direction that runs more than one hop and has the body straddling one of its interior splits. In 4×2×2 the body straddles only y and z, which run a single hop each; in 4×4×4 those same directions run three.
Two hypotheses I tested and can rule out:
min(2*ib_neighborhood_radius, num_procs_<dir> - 1)truncates when a direction carries more than2*radius+1ranks. The clean 8×2×2 case (the uniform-grid sphere at 32 ranks, F_y = 2.9e-5) is the truncated one — 7 hops needed in x, 6 taken. The broken 4×4×4 case is complete: 3 needed, 3 taken. Raisingib_neighborhood_radiusfrom 3 to 16 changes nothing (F_y = 0.436446), as it must, sincemin(32, 3) = 3either way.The prefix-sum-then-scatter algorithm is correct in exact arithmetic — I traced it through for 4 ranks with body cells on ranks 1 and 2, and both the
recv_forces_snapdouble-count subtraction and the backward replace give the right answer. So the fault looks like an implementation detail of the multi-hop path rather than the scheme: the most likely candidate isrecv_forces_snap/recv_torques_snap, which are written on the device at hop k and read there at hop k+1 while appearing in nocopy/copyinclause, relying entirely on residency. With a single hop the snap is only ever read as the zero it was initialised to, which is exactly why one-hop directions never show this.I have not confirmed that last step, so I am reporting the trigger rather than the cause.
sbryngelson commentedon Sep 12, 2026
Two more hypotheses tested and disproved, recorded so nobody spends time on them again.
Not
recv_forces_snapdevice residency. I suspected the snap buffers, since they are written on the device at hop k and read there at hop k+1 while appearing in nocopy/copyinclause — and a single-hop direction never reads them as anything but the zero they were initialised to, which fit the trigger exactly. Adding them to the accumulate loop's copy clause, so residency cannot matter:Unchanged. (The small F_x difference is unrelated — that build does not carry the out-of-bounds fix from #1859.)
Not reproducible in 2D. The obvious cheap regression case does not show it. A 2D cylinder, uniform 200×200, body at the domain centre so the mid-domain split lands inside it, at three rank counts:
4×4 runs three hops in y with the body straddling a y split — the same conditions that break the 3D case — and it is clean. So the trigger needs the third dimension, not merely a multi-hop sweep over a body-spanning direction.
Where that leaves it: the trigger is pinned (3D, ≥3 ranks in the directions the body spans, corrupting the components along those directions, F_x exact throughout), the sweep algorithm is correct in exact arithmetic, and three plausible mechanisms are eliminated. I have not found the cause and am not going to guess at a fourth. The sphere restart reproduces it in about twenty seconds on 8 nodes if someone wants to bisect it.
sbryngelson commentedon Sep 12, 2026
Correction to my earlier claim that stretching is not involved. It is — I drew the wrong conclusion from the right measurement.
I built a clean 3D control: a sphere at Re 100, Ma 0.1 on a uniform 128³ grid, body at the origin so the mid-domain split lands inside it, 50 steps. Same geometry, same physics, same topologies as the failing case:
Identical to six digits. A uniform grid at 4×4×4 is clean, where the stretched grid at the same 4×4×4 gives F_y = +0.436446.
So the correct statement is that neither stretching nor multiple ranks is sufficient alone — it needs both:
My earlier "not grid stretching" was based only on the single-rank row, which does not license that conclusion. Apologies for the noise.
This makes a rank-width-dependent quantity the obvious suspect, and there is one:
s_get_neighbor_boundsbuildsneighbor_domain_{x,y,z}by propagating each rank's physical edges outwardib_neighborhood_radiusrank hops, ands_handoff_ib_ownershipthen deletes any patch whose centroid falls outside that box. A deleted patch makess_get_neighborhood_idxreturn −1, and the reduction'sif (j > 0)then silently discards a relayed contribution passing through that rank. On a stretched grid the ranks near the body are physically narrow, so a fixed number of hops covers far less distance than on a uniform grid — which is exactly the dependence observed.I could not confirm it: raising
ib_neighborhood_radiusfrom 3 to 16 (enough hops to saturate a 4-rank direction and make every rank keep the patch) leaves F_y at 0.436446. So either that is not the path, or something else re-narrows the neighbourhood. I am leaving it there rather than guessing further.Running tally of eliminated mechanisms: the hop cap, snap-buffer device residency, stretching alone, multi-rank alone. The reproducer remains the stretched sphere restart, ~20 s on 8 nodes.
sbryngelson commentedon Sep 12, 2026
ib_neighborhood_radiushas no effect at all. Same restart, same binary, 64 ranks (4×4×4), only the radius changed:ib_neighborhood_radiusBit-identical to six digits in every component. That rules out two more things at once, and one of them is a surprise:
Not the neighbourhood.
s_get_neighbor_boundsgrowsneighbor_domain_*byib_neighborhood_radiusrank hops ands_handoff_ib_ownershipdeletes any patch whose centroid falls outside, which would makes_get_neighborhood_idxreturn −1 and the reduction'sif (j > 0)silently drop a relayed contribution. Radius 1 versus 16 changes that box by a factor of sixteen and changes nothing.Not the hop count either, and this is the informative part. The accumulate and scatter loops run
min(2*ib_neighborhood_radius, num_procs_<dir> - 1)hops. At radius 1 that is 2 hops over a 4-rank direction; at radius 3 it is 3. A 4-rank direction needs 3 for the sweep to complete, so radius 1 runs a provably incomplete reduction — and produces a bit-identical answer.If the reported force were the output of that sweep, an incomplete sweep would have to change it. It does not change in the last digit. So whatever reaches
patch_ib%forcedoes not depend on the sweep's hop count at all, which is hard to reconcile with the reduction being where the spurious value enters.Running tally of what is eliminated, each by measurement:
recv_forces_snapdevice residencycopyclause; F_y unchanged at 0.436413What still holds: it needs both stretching and ≥3 ranks in the directions the body spans; F_x is exact throughout; and the offline reproduction of the local volume sum from the restart field matches MFC's F_x to 0.15 % while giving F_y = +0.00002.
I am out of hypotheses I can test cheaply and am leaving it here rather than guessing at a sixth. The reproducer is the stretched sphere restart, about twenty seconds on eight nodes, and the three-row table at the top of this comment is a one-command bisection target.
sbryngelson commentedon Sep 12, 2026
The trigger stated earlier survives a deliberate attempt to falsify it. Four topologies now, same restart field, same binary, only the rank count changing:
The 8×2×2 row was run specifically because the rule predicts it clean and a different answer would have falsified it. With this grid the body spans only y and z, and at 32 ranks those carry two ranks each — one hop — while x carries eight. It comes back clean.
That row also closes the hop-count question from the other side. 8×2×2 runs seven hops in x and is clean; 4×4×4 runs three in y and z and is broken. Hop count alone predicts nothing; what predicts it is whether a multi-hop direction is one the body straddles.
So the trigger is, across four topologies:
Stated so it can be bisected against. Everything else remains as reported: single rank clean, uniform grid clean at the same 4×4×4,
ib_neighborhood_radiusbit-identical at 1, 3 and 16, snap buffers forced through acopyclause unchanged, and an offline reproduction of the local volume sum matching MFC's F_x to 0.15 % while giving F_y = +0.00002.sbryngelson commentedon Sep 12, 2026
Correction: my inference from the
ib_neighborhood_radiustest was wrong, and it wrongly steered away from the reduction.I wrote that radius 1 runs a "provably incomplete" reduction and that a bit-identical answer therefore meant the reported force does not depend on the sweep at all — "hard to reconcile with the reduction being where the spurious value enters". That does not follow.
With 4 ranks in a direction,
min(2*ib_neighborhood_radius, num_procs-1)gives 2 hops at radius 1 and 3 at radius 3. The body straddles the split at y = +0.030, so it occupies ranks 1 and 2 — adjacent, in the middle. Working the sweep through for both hop counts:Two hops already carries both body ranks up to rank 3 and back down to rank 0. Both configurations run a complete reduction, so bit-identical output is exactly what a correct reduction would produce. The test says nothing about the reduction either way.
What the radius test does still rule out is the neighbourhood box — that changes by a factor of sixteen between radius 1 and 16, which would change which ranks keep the patch, and nothing moves.
So the corrected tally:
s_communicate_ib_forcesgenerallyrecv_forces_snapdevice residencycopyclause, unchanged)A test that would actually discriminate needs the hop count to be genuinely insufficient, which means a body-spanning direction with more ranks between the body and the far end than the sweep covers — eight or more ranks in that direction with the body in the middle, at radius 1. Every topology I ran had the body on adjacent ranks in a 2- or 4-rank direction, where two hops is always enough. The 8-rank direction I did test (8×2×2) was x, which this body does not span.
Apologies for the noise. The empirical trigger from the four-topology table stands unchanged; it is the mechanistic inference drawn from the radius rows that was unsound.
sbryngelson commentedon Sep 12, 2026
s_communicate_ib_forcesis not the site, and this is a proof rather than an absence of evidenceThe issue text above names the hop cap
min(2*ib_neighborhood_radius, num_procs_<dir> - 1)as a truncation "whenever a direction carries more than2*ib_neighborhood_radius + 1ranks." That is wrong and I withdraw it. The cap is exactly right, for a reason that is only visible when the routine is read together withs_reduce_ib_patch_array.The relay is gated by neighbourhood membership
On a multi-rank run
num_ibsis not the global patch count.m_start_up.fpp:1306-1323compactspatch_ibdown to the patches whose centroid lies in this rank's subdomain grown byib_neighborhood_radiushops, andib_gbl_idx_lookupis-1for every patch that failed that test. Both phases of the sweep guard their absorb ons_get_neighborhood_idx(...) > 0, and a rank that does not track a patch never puts that patch into its own send list either. So a contribution cannot propagate past the edge of the neighbourhood: the relay chain is confined to ranks[c-R, c+R]per direction, wherecis the rank holding the centroid.That is precisely the hop budget. With
K = min(2R, num_procs_d - 1)over a neighbourhood[n_lo, n_hi] = [c-R, c+R], accumulate leavesA[i] = sum over [max(n_lo, i-K), i]and scatter leavesS[i] = A[min(i+K, n_hi)]. For anyiin the neighbourhood,so every tracking rank ends with the sum over every tracking rank, for any R. The cap is not a truncation;
2Ris the width of the region that can contribute.A replay of the routine with the gate in place agrees exhaustively — rank counts 2–14, every centroid position, body spans 0–4, radii 1–7, periodic ring and non-periodic: 0 failures.
Measurements, on a body straddling a split in all three directions
A static sphere, Re 100, uniform 128³, body at the domain centre so every even split lands inside it,
ib_neighborhood_radius = 3, built with the #1859 out-of-bounds fix:Nothing moves. 64 CPU ranks is the same 4×4×4 topology as the failing sphere. The GPU legs matter because
s_communicate_ib_forcesreturns on its first line whennum_procs == 1, so a single rank is the one configuration in which none of its kernels run — that was the obvious remaining suspect, and it is clean at 8 GPU ranks.The 2D cylinder of
examples/2D_ibm_force_decompositionsays the same at 1, 4 and 16 CPU ranks after #1859: Fx 1.02352019 / 1.02351629 / 1.02351629.(The 1-vs-8 GPU gap in Fx, 5.9e-4 relative, is run-to-run atomic ordering in the force sum at
m_ibm.fpp:1164-1166, not a decomposition effect — the two CPU legs agree exactly and the two GPU legs do not.)What is left
Rank count, the reduction, and the GPU build are each individually excluded on a body straddling splits in every direction. Every one of those runs used a uniform grid, and the failing sphere does not. Combined with the earlier result that stretching alone on one rank is clean, the trigger is the conjunction: grid stretching together with a rank split through the body.
The coordinate halo is the obvious mechanism and it is not it —
s_populate_grid_bc_direction(m_boundary_common.fpp:419) routes every interior rank interface throughs_mpi_sendrecv_grid_variable_buffer, sox_ccanddxin the buffer hold the neighbour's real stretched values; uniform extrapolation from the edge cell happens only at a physical boundary. So the remaining question is what else the force integral reads across a rank interface that is only wrong when cells on the two sides differ in width.Investigated with Claude Code.
https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG
sbryngelson commentedon Sep 13, 2026
Every mechanism I proposed is falsified, and the minimal case does not reproduce this at all
Posting this because I have published several conjectures on this issue and they are now all dead by measurement. The record should not be left with them standing.
The test case is a static sphere, Re 100, D = 1, 128³, body at the domain centre so every even split lands inside it,
ib_neighborhood_radius = 3, built with the #1859 fix. One binary across every leg, so each row differs from its neighbours in exactly one variable.u8andu64are bit-identical, as ares8ands64. Going from 8 to 64 ranks changes nothing to the last digit on either grid.And the restart path, which is how this issue's own numbers were taken — 1 rank against 64 on the same restart file:
What that retires
s_communicate_ib_forcesThe only defect actually found and fixed along the way is the out-of-bounds
fd_coeffread (#1859), which is real but is not this — as the issue text already says.(
Fxmoving 1.60 → 1.33 between uniform and stretched is resolution, not a defect: the stretched core puts 14 cells across the sphere against 16 uniform. The 5.9e-4 and 2.1e-3 spreads within each triple are run-to-run atomic ordering in the force sum atm_ibm.fpp:1164-1166, which is a separate known issue.)What I would do next
I cannot reproduce this with a 128³ sphere under any combination of decomposition, build, stretching or restart. Continuing to guess at mechanisms from a case that does not trigger the defect is not productive.
The remaining difference is the failing case itself: 405×270×270, a non-cubic grid, run 50,000 steps. The useful next step is to bisect that case down — halve the grid, shorten the run, drop to a cube — and find the smallest configuration that still shows Fy ≈ Fx. Without a reproducer the mechanism hunt is unbounded.
Investigated with Claude Code.
https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG
sbryngelson commentedon Sep 18, 2026
Bisection of the real case: the transverse offset is rank-invariant, and it is
F_xthat breaksMy last comment ended by saying the 405×270×270 case needed bisecting down, and that guessing at mechanisms from a 128³ case that does not trigger the defect was not productive. That bisection has now happened, incidentally, while building a grid-convergence ladder for this same deck.
Four grid spacings, several rank counts, one binary throughout — master plus #1887 and the #1859 clamp (
i_fd = min(max(i, 0), m)), so the out-of-boundsfd_coeffread is not in play. The body is static, so none of the moving-IB changes on that branch are either. Static sphere, Re 100, Ma 0.1, D = 1, body at the domain centre; trueF_yandF_zare zero.Table 1 — the same deck, matched times
F_xF_yF_zThree things in there are new, and two of them contradict claims I left standing above.
1. DX 0.05 is clean at the rank count that breaks DX 0.0333
Same deck, same binary, same requested topology — only the grid spacing differs. At 64 ranks the 0.05 grid gives
F_y= 2.7e−05 and the 0.0333 grid gives +0.4365. The 0.05 grid is also bit-identical at 16 and 64 ranks in all three components.That is a two-node run against an eight-node run, under an hour each, and it is the smallest pair I have that brackets the defect.
2. "A uniform grid is clean" does not generalise — it is rank-independent, not transverse-free
I reported a uniform 128³ sphere with
F_y= 0.000000 at 4×4×4 and concluded the trigger needs stretching. The same deck made uniform (361×241×241, same DX 0.0333 core spacing) givesF_xbit-identical across 8, 64 and 256 ranks — a factor of 32 in rank count, not one digit moves — andF_y=F_z= +0.578 at every one of them, including the smallest.So there look to be two separable effects on this geometry: a transverse offset that does not depend on the decomposition at all, and something that additionally corrupts the force once ranks are added on a stretched grid. The 128³ control was clean for a reason that does not carry over.
3.
F_xis not exact — and it is the transverse components that stay putEvery earlier comment on this issue says
F_xsurvives to six digits. On the 0.0222 grid, at one instant, three rank counts:F_xF_yF_yagrees to five digits across all three whileF_xchanges sign. Run out to steady state, the 256-rank case givesC_D= −1.136, constant from the first step to t = 37. The 0.0333 grid does the same thing more mildly:F_x0.5520 at 8 and 64 ranks, 0.1327 at 256.The issue title says "spurious transverse component on multiple ranks". On this evidence the transverse component is the rank-invariant part, and the rank dependence is in
F_x.4. It does not reach a comparable production deck
A static AR-4 flat plate at 45°, Re 300, stretched grid 904×677×567 (347 M cells, boundaries 8 c laterally and 36 c downstream), 216 steps:
F_xF_y(zero by symmetry)F_zIdentical to every digit printed, with the residual
F_yat 0.02 % ofF_x— discretisation asymmetry, not this. So it is not "stretched grid plus many ranks" in general; whatever selects it is narrower than that.Reproducing
The deck is the sphere case this issue has used throughout. I have added an
ASTRenvironment knob to it so the uniform control is the same deck with one parameter changed rather than a separate case:ASTR=0builds the same domain and the same core spacing on a uniform grid and leaves everything else untouched. The cheapest bracketing pair is stretched DX 0.05 at 64 ranks (clean) against stretched DX 0.0333 at 64 ranks (broken).No mechanism proposed — this thread has had enough of those from me. These are measurements, including two that retract earlier ones.
Measured with Claude Code.