Latent out-of-bounds read in mibm central-difference IB drag (s_compute_viscous_stress_tensor) for boundary-adjacent bodies

[sbryngelson]

Summary

s_compute_viscous_stress_tensor (added in #1614 for higher-order central-difference IB drag) performs an out-of-bounds heap read whenever an immersed-boundary body's finite-difference stencil is centered on a cell adjacent to a domain boundary. The bug has been latent in master since #1614 (2026-06-24); it is silent in optimized builds and only manifests as a hard crash under bounds-checking (./mfc.sh ... --debug/reldebug) builds.

cc @danieljvickers (author of #1614 and #1618)

Root cause

In src/simulation/m_viscous.fpp, the velocity-gradient stencil indexes the FD-coefficient arrays by cell:

velocity_gradient_tensor(l, 2) = velocity_gradient_tensor(l, 2) &
    + fd_coeff_y(r, j)*q_prim_vf(...)%sf(i, j + r, k)

But fd_coeff_x/y/z are allocated over the interior only — e.g. fd_coeff_y(-fd_number:fd_number, 0:n) (in m_derived_variables.fpp). The IB drag path (s_compute_ib_forces in m_ibm.fpp) evaluates this stencil at body ghost-point cells, which for a body sitting against a boundary (very common for periodic IB configurations) land at j = -1, j = n+1, etc. Indexing fd_coeff_y(r, -1) reads the word before the interior slab — an out-of-bounds access.

The q_prim_vf%sf field itself is fine to index there (it is allocated with ghost cells and populated by the BC/periodic halo). Only the coefficient array is interior-only.

Why it has been silent

• Optimized builds: the OOB read grabs adjacent heap memory. Garbage, but no crash — the run completes and produces a golden that (being reproducible per-platform) matches. Every normal-build CI lane has been passing while silently reading out of bounds on boundary bodies.
• Bounds-checking (reldebug) builds: the same read becomes a hard Index '-1' ... below lower bound of '0' abort.

Whether the crash triggers additionally depends on whether the (chaotic) 3D_mibm_periodic_collision trajectory happens to place a body against the boundary during the drag evaluation on a given compiler/platform — so it can appear and disappear across unrelated changes. #1618 ("Fix periodic ib issues") shifted that trajectory, and any downstream change to shared code can move it again. This makes it look like an intermittent regression when it is actually a fixed latent bug being surfaced by chaotic sensitivity.

Reproduce

./mfc.sh test --reldebug --no-gpu -o 38B4AED7   # 3D -> Example -> mibm_periodic_collision

on a tree where the chaotic trajectory lands a body on the boundary (currently reproduces on the block-AMR mega branch, PR #1628).

Secondary issue: 3D_mibm_periodic_collision is not a portable regression target

Independent of the OOB: this case is chaotic (stiff, collisional). Even with the OOB fixed, its step-50 fields diverge across compilers/platforms, so its golden file is not portable. It should not be a hard tolerance-gated regression case (it is effectively already treated as marginal on GPU). Recommend marking it non-portable / smoke-only.

Suggested fix

Clamp the coefficient's cell index into the interior (the coefficients are cell-independent on uniform grids, so this is exact there and a reasonable boundary approximation on stretched grids), while keeping the q_prim stencil on the populated ghost/periodic neighbors:

fd_coeff_y(r, min(max(j, 0), n))

A one-sided boundary stencil would be more accurate on stretched grids and is worth considering. I have applied the clamp form in PR #1628 (block-structured AMR) to unblock its bounds-checking CI; happy to split it into a standalone PR against master if you'd prefer to take it directly.