Give the finite-difference coefficients the buffer range the IB force asks for

[sbryngelson]

Fixes #1860.

The bug

s_compute_ib_forces 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 then indexes the coefficient array at the cell it was handed (m_viscous.fpp:1141):

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

so fd_coeff_x is read at i - fd_number through i + fd_number, while it is allocated 0:m. A body cell within fd_number of a subdomain edge reads off the end of the array and whatever is adjacent in memory becomes a force.

Nothing reports it. The read stays in bounds for a body in the interior of one rank's subdomain, which is why single-rank runs and the goldens never showed it.

The fix

s_compute_finite_difference_coefficients already takes an offset_s. Pass one, widen the allocation to match. x_cc carries at least buff_size = 10 ghost cells whenever ib is set (s_configure_coordinate_bounds), so the extra coefficients come from real coordinates rather than a clamp.

Three lines of allocation, three call sites, one local. Nothing else reads outside 0:m, so no other consumer changes.

How the validation is obtained

examples/2D_ibm_force_decomposition/ — static cylinder, Re 40, Ma 0.1, uniform 200 x 200, fd_order = 4, 200 steps. The body is at the domain center, so an even rank count in a direction puts a subdomain edge through it. The force on a rigid body is a property of the flow, so the same case on a different rank count must give the same answer.

Read restart_data/ib_state_200.dat (20 doubles, [time, Fx, Fy, Fz, Tx, Ty, Tz, ...]):

ranks	Fx before	Fx after
1	1.02352019	1.02352019
4	1.02933438	1.02351629

0.57 % of drag out of adjacent memory, against 0.0004 % after. The single-rank number is bit-identical before and after, as it must be — one rank never reaches outside its own interior. Lift stays at ~8e-6 in all four configurations (the cylinder is symmetric about y = 0), so the drag agreement is not bought with a symmetry error.

The README in that directory carries the same table and the command to reproduce it.

Scope

ib + viscous on more than one rank with a body straddling a subdomain edge. Forces are diagnostic for a static body; for moving_ibm the force drives the motion, so those trajectories change too. Existing IB goldens are single-rank and are unaffected — no regeneration needed.

The discrepancy is not bounded by anything, since it is set by adjacent memory. On a 3D sphere at Re 100 on 64 GPU ranks, body at the center of the domain and so at the corner of eight subdomains, the same read shows up as a transverse force of 1.08 x the drag on a body that has none (Fx +0.4040, Fy +0.4364, Fz +0.4364), pinned at that value from the first time step to the fifty-thousandth while Fx relaxes from 1.638 to 0.404.

https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG

[sbryngelson]

Duplicate of #1859, which was opened first and fixes the same out-of-bounds read by clamping to the nearest interior cell rather than widening the allocation. Closing in its favour.

The example case and the measured before/after table have been moved to that branch, and the clamp-versus-widen tradeoff is written up in my comment there — the diff here stays available if the wider allocation is preferred, since the two differ at first order on a stretched grid where the clamp substitutes a neighbouring cell's stencil.

[Copilot]

🟡 Changes recommended

A critical probe-only out-of-bounds access remains, alongside reproducer inconsistencies and missing automated multi-rank regression coverage.

Get a fresh assessment by requesting another Copilot review.

Pull request overview

This pull request fixes MPI-decomposition-dependent immersed-boundary force errors by extending finite-difference coefficient coverage and adding a 2D reproducer.

Changes:

• Extends coefficient allocation and initialization for IB stencils.
• Adds a static-cylinder decomposition test case.
• Documents reproduction steps and force comparisons.

File summaries

File	Summary
src/simulation/m_derived_variables.fpp	Extends finite-difference coefficient buffers and initialization.
examples/2D_ibm_force_decomposition/README.md	Documents the reproducer and expected results.
examples/2D_ibm_force_decomposition/case.py	Defines the static-cylinder test case.

Review details

Suppressed comments (4)

examples/2D_ibm_force_decomposition/README.md:3

• case.py sets L=5, DX=0.1, then N = 2*L/DX = 100 and m = n = N - 1; because the grid is generated over 0:m, this case has 100×100 cells, not 200×200 as documented here and in the PR description. Please either correct the stated resolution/results provenance or change the default L/DX values so the reproducer matches the documented grid.

A static cylinder in uniform flow at Re 40, Ma 0.1, on a uniform 200 x 200 grid. The body sits at the

examples/2D_ibm_force_decomposition/README.md:14

• Both commands write the same restart_data/ib_state_200.dat, so by the time the extraction command runs only the 4-rank result remains; the documented procedure cannot actually report the 1-rank/4-rank comparison without manually copying the first file. Read or save the result immediately after each run.

./mfc.sh run examples/2D_ibm_force_decomposition/case.py -n 1
./mfc.sh run examples/2D_ibm_force_decomposition/case.py -n 4

examples/2D_ibm_force_decomposition/case.py:18

• With the default values L=5.0 and DX=0.1, this computes N=100 and sets m=n=99, so the simulation has 100 cells per direction (the loops use indices 0:m/0:n), not the 200×200 grid claimed by the README and PR description. Please either adjust the case defaults or correct the stated resolution, since the published force numbers are not for the documented grid.

N = int(2 * L / dx)

src/simulation/m_derived_variables.fpp:45

• This fixes a decomposition-dependent out-of-bounds read, but the regression is not exercised by the automated suite: the new example is not registered in toolchain/mfc/test/cases.py, and the existing IB force case runs with the default single rank while the existing MPI IB case does not enable ib_state_wrt. A future reversion can therefore pass all goldens; please add a small multi-rank regression that checks the force or rank-count agreement.

            @:ALLOCATE(fd_coeff_x(-fd_number:fd_number, -fd_number:m + fd_number))

• Files reviewed: 3/3 changed files
• Comments generated: 1
• Review effort level: Lite

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[github-actions]

Claude Code Review

Head SHA: f0f86b2

Files changed:

• 3
• examples/2D_ibm_force_decomposition/README.md
• examples/2D_ibm_force_decomposition/case.py
• src/simulation/m_derived_variables.fpp

Findings:

• src/simulation/m_derived_variables.fpp:78,82,86: the calls to s_compute_finite_difference_coefficients now pass a 7th actual argument, fd_offset. The changed-files list for this PR contains only these three files, so the file defining s_compute_finite_difference_coefficients (imported via use m_finite_differences) is not modified here. Unless that subroutine's interface already accepted this argument before this PR, all three call sites will fail to compile with an argument-count/type mismatch.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_derived_variables.fpp	443	+2

Directory	Lines	Diff
simulation	28088	+2
total	46342	+2