Immersed-boundary force gains a large spurious transverse component on multiple ranks

[sbryngelson]

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.

[sbryngelson]

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:

	C_D	C_L
8 ranks	1.8982	+0.080406
1 rank	1.9108	−0.001278

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:

t     5    10    15    20    25    30    35    40    45    50    55    60
C_L  .0804 .0759 .0727 .0728 .0742 .0752 .0772 .0751 .0740 .0787 .0753 .0776
C_D  1.898 1.630 1.633 1.613 1.600 1.594 1.590 1.587 1.586 1.591 1.583 1.585

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]

Narrowed to a precise trigger. Same restart field, same binary, same case — only the rank count changes, which changes only the decomposition:

ranks	topology	F_x	F_y	F_z
8	2×2×2	+0.404403	+0.000023	+0.000010
16	4×2×2	+0.404403	+0.000023	+0.000010
64	4×4×4	+0.404403	+0.436446	+0.436433

F_x is identical to six digits in all three. The defect needs three or more ranks in a direction, i.e. more than one hop in s_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

x, 4 ranks:  -1.069, +2.029, +7.926     (body spans [-0.5, +0.5] -- no split inside it)
y, 4 ranks:  -2.754, +0.030, +2.831     (+0.030 is inside the body)
z, 4 ranks:  -2.754, +0.030, +2.831     (+0.030 is inside the body)

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:

• The hop cap is not it. min(2*ib_neighborhood_radius, num_procs_<dir> - 1) truncates when a direction carries more than 2*radius+1 ranks. 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. Raising ib_neighborhood_radius from 3 to 16 changes nothing (F_y = 0.436446), as it must, since min(32, 3) = 3 either way.
• It is not grid stretching. The single-rank run on the identical stretched grid gives F_y = +0.000024.

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_snap double-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 is recv_forces_snap/recv_torques_snap, which are written on the device at hop k and read there at hop k+1 while appearing in no copy/copyin clause, 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]

Two more hypotheses tested and disproved, recorded so nobody spends time on them again.

Not recv_forces_snap device 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 no copy/copyin clause — 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:

$:GPU_PARALLEL_LOOP(private='[i, j]', copy='[forces, torques, recv_forces_snap, recv_torques_snap]')

	F_x	F_y	F_z
4×4×4 as-is	+0.404403	+0.436447	+0.436433
4×4×4, snap forced through copy	+0.403987	+0.436413	+0.436399

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:

ranks	topology	C_D	C_L, true value zero
1	1×1	2.047040	0.000017
4	2×2	2.047033	0.000016
16	4×4	2.047033	0.000016

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]

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:

ranks	topology	F_x	F_y, true value 0	F_z, true value 0
8	2×2×2	1.601324	0.000000	−0.000002
64	4×4×4	1.601324	0.000000	−0.000002

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:

grid	ranks	F_y
stretched	1	+0.000024
uniform	64, 4×4×4	0.000000
uniform	32, 8×2×2	+0.000029
stretched	64, 4×4×4	+0.436446

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_bounds builds neighbor_domain_{x,y,z} by propagating each rank's physical edges outward ib_neighborhood_radius rank hops, and s_handoff_ib_ownership then deletes any patch whose centroid falls outside that box. A deleted patch makes s_get_neighborhood_idx return −1, and the reduction's if (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_radius from 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]

ib_neighborhood_radius has no effect at all. Same restart, same binary, 64 ranks (4×4×4), only the radius changed:

ib_neighborhood_radius	F_x	F_y	F_z
1	0.404403	0.436446	0.436433
3	0.404403	0.436446	0.436433
16	0.404403	0.436446	0.436433

Bit-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_bounds grows neighbor_domain_* by ib_neighborhood_radius rank hops and s_handoff_ib_ownership deletes any patch whose centroid falls outside, which would make s_get_neighborhood_idx return −1 and the reduction's if (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%force does 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:

candidate	how it was ruled out
the hop cap truncating the sweep	the clean 8×2×2 case is the truncated one; radius 1 vs 3 vs 16 bit-identical
recv_forces_snap device residency	forced through the loop's copy clause; F_y unchanged at 0.436413
grid stretching alone	single rank on the same stretched grid gives F_y = +0.000024
multiple ranks alone	uniform 128³ at 4×4×4 gives F_y = 0.000000
the neighbourhood box / patch deletion	radius 1 vs 16, bit-identical

What 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]

The trigger stated earlier survives a deliberate attempt to falsify it. Four topologies now, same restart field, same binary, only the rank count changing:

ranks	topology	F_x	F_y	F_z
8	2×2×2	0.404403	0.000023	0.000010
16	4×2×2	0.404403	0.000023	0.000010
32	8×2×2	0.404403	0.000023	0.000010
64	4×4×4	0.404403	0.436446	0.436433

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:

three or more ranks in a direction whose interior splits fall inside the body, corrupting the force components along those directions, with F_x exact to six digits throughout.

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_radius bit-identical at 1, 3 and 16, snap buffers forced through a copy clause unchanged, and an offline reproduction of the local volume sum matching MFC's F_x to 0.15 % while giving F_y = +0.00002.