Read IB viscous FD coefficients at the nearest interior cell

[sbryngelson]

Closes #1856.

s_compute_ib_forces samples the viscous stress fd_number cells out from each interior cell:

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

and s_compute_viscous_stress_tensor uses that sample index to look up the coefficient itself, fd_coeff_y(r, j). The coefficients exist for interior cells only:

@:ALLOCATE(fd_coeff_x(-fd_number:fd_number, 0:m))
@:ALLOCATE(fd_coeff_y(-fd_number:fd_number, 0:n))
@:ALLOCATE(fd_coeff_z(-fd_number:fd_number, 0:p))

So a body within fd_number cells of a domain boundary reads past the array:

At line 1143 of file src/simulation/m_viscous.fpp
Fortran runtime error: Index '-2' of dimension 2 of array 'fd_coeff_y' below lower bound of 0

Only the reldebug lanes say so, because only they carry bounds checking. Everywhere else this reads whatever sits in front of the array and returns a force that looks plausible, which is the part worth fixing rather than the crash.

This reads the nearest interior cell's coefficients instead. On a uniform grid they are the same coefficients; on a stretched grid they are the stencil the boundary cell itself uses.

i_fd = min(max(i, 0), m)
j_fd = min(max(j, 0), n)
k_fd = min(max(k, 0), p)

Testing

Reproduced and fixed with gfortran under ./mfc.sh test --no-gpu --debug. The case that exercises it is 86A2C59B (3D -> Example -> ibm_flapping_plate), which arrives with #1847 - it is the first viscous IB case in the suite with a body near a domain boundary, which is why this sat unnoticed. Without the change that case dies with the error above; with it, it and both Prescribed Kinematics cases pass. No existing case changes behaviour: every current caller passes an interior index, for which min/max are identities.

The regression coverage therefore lands with #1847 rather than here. I did not add a case for it, since duplicating that example only to cover the same path seemed worse than letting #1847 carry it - say the word if you would rather it were self-contained.

Note on the first commit

4f4987a is unrelated and fixes three doc-reference lint hits that came in with #1855: two single-backtick spans containing a single quote (which Doxygen reads as closing the span) and a link to a #selection-and-execution-pitfalls anchor that was never defined. Master's Lint Toolchain is red on them as of 5953684, so they are in here only to get this branch's own lint gate green. Happy to split them out if you would rather they went separately.

[sbryngelson]

@danieljvickers - flagging this one for you since it sits in the IB force path you have been working in, and it is adjacent to what you hit in #1792.

Both are the same shape: a finite-difference coefficient is asked for outside the range it was computed for. In #1792 it was s_compute_finite_difference_coefficients filling past the end of s_cc; here it is s_compute_viscous_stress_tensor reading fd_coeff_y at a sample index that can sit fd_number cells outside the interior, while the coefficients only span 0:n. Neither shows up outside reldebug, because bounds checking is the only thing that notices - otherwise you get a force built partly from whatever precedes the array.

Two things I would value your read on:

1. Whether clamping to the nearest interior cell is the behaviour you want for the IB force integration, or whether those ghost-adjacent samples should be getting real coefficients. That is what your Debug ibm stability #1792 change was reaching for with the widened range, and if the IB work needs it, the honest fix for both is to compute coefficients over the margin buff_size actually provides rather than to clamp.
2. Whether any other fd_coeff_* consumer indexes at a sample location rather than a cell location. I only audited the viscous path.

Detail on the reproduction and the verification is in #1856 and the description above.

[Copilot]

Warning

Copilot couldn't run its full agentic review because it didn't start before the timeout. Make sure your repository has a runner available, or add a copilot-code-review.yml file specifying one with the runs-on attribute. See the docs for more details.

Pull request overview

Fixes an out-of-bounds access when immersed-boundary (IB) viscous force computation samples finite-difference coefficients near domain boundaries, and updates documentation to resolve Doxygen/lint issues.

Changes:

• Clamp FD coefficient lookup indices to the nearest interior cell to prevent out-of-bounds reads in m_viscous.
• Add a Doxygen anchor in testing docs to fix a broken reference.
• Escape inline-code spans in documentation to avoid Doxygen parsing issues.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 2 comments.

File	Description
src/simulation/m_viscous.fpp	Clamps FD coefficient indices to interior bounds to prevent invalid coefficient reads near boundaries.
docs/documentation/testing.md	Adds the missing Doxygen anchor for the “Selection and Execution Pitfalls” section.
docs/documentation/gpuParallelization.md	Fixes inline-code formatting to avoid Doxygen span termination issues.
docs/documentation/contributing.md	Fixes inline-code formatting for a GPU directive snippet to satisfy doc linting.

---

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

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_viscous.fpp	1025	+4

Directory	Lines	Diff
simulation	28090	+4
total	46344	+4

[sbryngelson]

I found the same bug independently and opened #1862 before seeing this. This PR came first, so I am closing mine and have moved its example case here — examples/2D_ibm_force_decomposition (pushed to this branch, no behavior change).

Measured evidence this fix is right. Static cylinder, Re 40, Ma 0.1, uniform 200×200, fd_order = 4, 200 steps, body at the domain center so an even rank count puts a subdomain edge through it. Drag from restart_data/ib_state_200.dat:

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. Lift stays at 8e-6 in all four configurations, so the drag agreement is not bought with a symmetry error.

One thing worth reconsidering in the comment, and maybe in the fix. The comment says "a body near a domain boundary." It is not only the domain boundary — fd_coeff is allocated over the rank-local interior 0:m, so this fires whenever a body is within fd_number of a subdomain edge. That is the common case, not the rare one: a body at the center of the domain on an even decomposition hits it every time, which is exactly what the table above measures on an interior body with non-reflecting boundaries ten diameters away.

That also bears on clamp-versus-widen. Clamping to the nearest interior cell is right at a true domain edge, where there is no cell beyond. At a subdomain edge the cell at m+1 is a real interior cell of the neighbouring rank with perfectly well-defined coefficients, and the clamp substitutes its neighbour's. On a uniform grid the two agree exactly, so the table above is identical either way; on a stretched grid they differ at first order. s_compute_finite_difference_coefficients already takes an offset_s, so the alternative is three allocation bounds and one local:

@:ALLOCATE(fd_coeff_x(-fd_number:fd_number, -fd_number:m + fd_number))
...
fd_offset%beg = fd_number; fd_offset%end = fd_number
call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, fd_number, fd_order, fd_offset)

x_cc carries at least buff_size = 10 ghost cells whenever ib is set (s_configure_coordinate_bounds), so those coefficients come from real coordinates. That diff is in #1862 if you want to take it; either way the bug is caught and I have no objection to merging this as it stands.

One correction to carry over. I had also attributed a spurious transverse force on a 3D sphere at Re 100 on 64 ranks (Fy = Fz = 1.08 × drag on an axisymmetric wake) to this bug. Rebuilding with the fix and re-running that case from its own restart leaves it unchanged at 0.4364, so that is a separate, still-unexplained defect and should not be claimed for this PR. What this PR fixes is exactly the table above.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 61.33%. Comparing base (dc0aec1) to head (ee1ec6e).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1859      +/-   ##
==========================================
+ Coverage   61.26%   61.33%   +0.06%     
==========================================
  Files          84       84              
  Lines       22330    22333       +3     
  Branches     3265     3265              
==========================================
+ Hits        13680    13697      +17     
+ Misses       6207     6187      -20     
- Partials     2443     2449       +6     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[danieljvickers]

I have already added a fix for this on my debug branch. To ensure rank invariance, you need to sample into the halo region. The fix I added in this vein was minimal and is currently being demonstrated to work on my branch. It is annoying that i have added all of these changes into a monolithic PR, but desperate times... Regardless, I think I have a sufficent answer on my own branch.