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

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:

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.