Immersed-boundary force reads finite-difference coefficients out of bounds, giving a force that depends on the MPI decomposition

[sbryngelson]

This issue originally attributed the symptom to grid stretching in the force stencil. That diagnosis was wrong and is corrected below; the stretching correlation was a coincidence of which runs happened to put a subdomain edge through the body.

What happens

s_compute_ib_forces (src/simulation/m_ibm.fpp) evaluates the viscous stress tensor at every cell its own stencil reaches:

do l = -fd_number, fd_number
    call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i + l, j, k)

s_compute_viscous_stress_tensor (src/simulation/m_viscous.fpp:1141) then indexes the coefficient array at the cell it was handed:

velocity_gradient_tensor(l, 1) = velocity_gradient_tensor(l, 1) + fd_coeff_x(r, i)*q_prim_vf(...)%sf(i + r, j, k)

so the coefficients are read at i - fd_number through i + fd_number. They are allocated 0:m (src/simulation/m_derived_variables.fpp:43). A body cell within fd_number of a subdomain edge therefore reads fd_coeff off the end of the array, and whatever is adjacent in memory becomes a force.

Nothing reports this. The read is in bounds for a body sitting in the interior of one rank's subdomain, which is why single-rank runs and the golden files never show it.

Reproducer

examples/2D_ibm_force_decomposition — a static cylinder at Re 40, 200×200, 200 steps, fd_order = 4, body at the domain center so an even rank count in a direction puts a subdomain edge through it. Drag from restart_data/ib_state_200.dat:

ranks	Fx
1	1.02352019
4	1.02933438

0.57 % apart, on a uniform grid, same case, same binary. A force that changes when you re-cut the domain.

The size of the discrepancy is set by whatever happens to be adjacent in memory, so it is not bounded by anything. On a 3D sphere at Re 100 on 64 GPU ranks, with the body at the center of the domain and so at the corner of eight subdomains, it appears as a transverse force of 1.08 × the drag on a body that has no transverse force at all: Fx +0.4040, Fy +0.4364, Fz +0.4364. That number is pinned at 0.4364 from the first time step to the fifty-thousandth while Fx relaxes from 1.638 to 0.404 — it is not a flow quantity. Reproducing MFC's own volume sum offline from the restart field gives Fx +0.40459 (0.15 % from MFC's 0.40399) and Fy +0.00002, confirming the local sum itself is right and the spurious part enters through the out-of-bounds read.

Fix

s_compute_finite_difference_coefficients already accepts an offset_s. Passing one and widening the allocation to -fd_number:m + fd_number covers the range that is actually read. x_cc carries at least buff_size = 10 ghost cells whenever ib is set (s_configure_coordinate_bounds), so the extra coefficients are built from real coordinates rather than clamped.

With that change the reproducer gives 1.02352019 at 1 rank and 1.02351629 at 4 — 0.0004 % apart, and the 1-rank answer is bit-identical to before, as it must be.

Scope

Affects ib + viscous on more than one rank whenever a body straddles a subdomain edge. Forces are diagnostic for a static body, so nothing else changes there; for moving_ibm the force drives the motion, so trajectories are affected too.

Separate, smaller observation

The force integral's outer stencil is the centered one, whose coefficients on a non-uniform grid are the uniform-grid weights over a non-uniform denominator. It is first-order there and, more to the point for a force, non-telescoping, so the volume sum of −∇p does not reduce exactly to the surface integral it stands for. Measured on the stretched sphere grid this moves the pressure drag about 2 % (0.20975 vs 0.20541 for a conservative flux-form stencil) and contributes no transverse force. Worth a look, but it is an accuracy question and not the defect above; filing separately if it is wanted.

[sbryngelson]

Correction: the sphere is not this bug.

I rebuilt with the fix and re-ran the sphere from its own restart at step 49780 on the same 64 ranks. The spurious transverse force is unchanged:

step	Fx	Fy	Fz
49780 (original, unfixed)	+0.403987	+0.436413	+0.436399
49781 (fixed build)	+0.404403	+0.436447	+0.436433
49790 (fixed build)	+0.404404	+0.436445	+0.436432

The fix is in the binary that ran — the generated source carries allocate (fd_coeff_x(-fd_number:fd_number, -fd_number:m + fd_number)) and the fd_offset call sites. So the out-of-bounds read is real and worth fixing on its own evidence, but it is not what the sphere is showing, and I should not have put the sphere in this issue as if it were.

What still stands, unchanged, is the 2D reproducer in examples/2D_ibm_force_decomposition: drag 1.02352019 at 1 rank against 1.02933438 at 4, on a uniform grid, going to 1.02351629 after the fix with the 1-rank answer bit-identical. That is a decomposition-dependent force from a read past the end of an array, and it is the whole of what #1862 claims.

The sphere's transverse force is a separate, still-unexplained defect. What is established about it:

• It 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 50,000 steps. Constant to four digits while the flow develops completely — not a flow quantity.
• Fy and Fz agree to five digits (y and z carry identical stretching), and Ty = −0.028511 against Tz = +0.028507, which is what a spatially uniform spurious force density in y and z gives through Σ r × F.
• Reproducing MFC's volume sum offline from the restart field, with the fd_order = 4 stencil the run actually uses, gives Fx +0.40459 against MFC's +0.40399 — 0.15 % — and Fy +0.00002. So the local sum over the field as the restart file holds it carries no transverse force.
• The same case on a uniform grid at 32 ranks gives Fy = 6e-6.

That leaves the field the force loop sees inside the body (which the restart need not hold), or s_communicate_ib_forces, as the remaining candidates. A single-rank run of the same restart is queued to separate them. I will open a separate issue with the answer rather than keep amending this one.