Debug ibm stability

[danieljvickers]

There has been a growing amount of technical debt on the immersed boundary code for multi-rank cases since the introduction of the IB neighborhoods. This has led to a host of new potential issues that threaten the stability of simulations being run. As I recently began scaling a relatively-difficult case, in terms of opportunities for instability, I made multiple bugfixes that were latent and untested. Some were relatively innocuous and others were extremely problematic, but all obvious bugs. And explanation of the changes are as follows

1. Added a mask to the interior of immersed boundaries for CFL assessment. Since the values interior to the IBs are meant to be non-physical, using them to check CFL causes violations that otherwise are not problematic. Particularly, the interior can produce odd pressures, which affect the speed of sound computed at the grid cell in question. This has been resolved with a mask during CFL checking.
2. Fixed race condition in IB marker write. The multi-IB parallelism that was introduced also created a race condition for actively-colliding particles. Two IBs could try to claim the same grid cell, but nothing deconflicted this race. In the other IB parallelism, we have left the behavior that the IB with the highest global index claims the grid cell, which I mirrored in the many-ib parallelism state
3. Fixed race condition in ghost-point correction. When two spheres are close together, a ghost point on sphere A may sample an image point from the interior of sphere B, which is also a ghost point. Because the ghost points are written in parallel, this generated a race condition. To deconflict this, I have separated the subroutine into two kernel calls. One kernel interpolates the image points and the following updates the ghost point values. This prevents one ghost point being updated before another image point may be interpolated. The initial pass I had Claude write is somewhat clunky, and I will have to touch it up by hand before the merge. But it works and resolves the race.
4. Fixed non-updated num_ibs. The num_ibs is computed in the ib ownership handoff subroutine, but was never updated to the GPU. This means that if new particles entered the neighborhood, then several subroutines would not iterate over them, including the state integration. This caused cross-rank drift and several issues with integrating the forces themselves. I originally thought that this was an error in the updating of the global IB index lookup, which caused me to make changes there that I should revert before merge.
5. Fixed Latent out-of-bounds read in mibm central-difference IB drag (s_compute_viscous_stress_tensor) for boundary-adjacent bodies #1633 out of bounds memory read when integrating the drag coefficeint of IBs near processor boundaries. Details in that issue. The fix I opted for was to extend the integration rather than use one-sided integrals. This keeps results consistent in multi-rank cases.
6. Fix non-updated num_gps, for similar reasons to num_ibs. If the number of ghost points on a rank goes down, then the update loop will update some grid cell that should no longer be interior to an IB, causing potential instabilities in fast-moving IBs where large discontinuities can occur as grid cells are ejected.
7. Added a density correction that was causing incorrect energy conversion. The reference material that I used to implement the IBM never explicitly declared a density conversion at the ghost points, but naively using a Neumann was accurate previously if the pressure at the ghost points matches that at the image points. However, when this is not the case, a correction based on your EOS is required. I have added this term to the correction to compute the energy, which now takes the density into account when computing energy.
8. Fixed edge case where collisions rarely were not being detected across rank boundaries. This occurred when one rank did not see the collision at all, but the other rank thought that it did not own the collision. The result was that the collision did not get counted until there was more-significant overlap, which caused a violent spring effect and destabilized the collision.
9. Updated unstable image point value calculation when IP is deep inside a neighbor particle. This caused instabilities in high-mach flows. The implementation was due to a misreading of the original source text. I have corrected the implementation to average over neighboring ghost points in this case, matching the reference literature.

The most impactful bug fixes were the update of num_ibs and num_gps, and the checking of interior GP points with a mask during CFL. Current tests show total stability on 8 ranks with 600 IBs that are very light compared to the ambient fluid in a mach 10 shock. Assuming these results hold, then we should have much higher stability, even in non-physically significant regimes. I am currently working on extending this result to higher numbers of IBs and ranks.

---

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

---

PR template credit: junegunn

[github-actions]

Claude Code Review

Head SHA: 88877c1

Files changed:

• 10
• src/common/m_constants.fpp
• src/common/m_derived_types.fpp
• src/simulation/m_data_output.fpp
• src/simulation/m_global_parameters.fpp
• src/simulation/m_ib_patches.fpp
• src/simulation/m_ibm.fpp
• src/simulation/m_mpi_proxy.fpp
• src/simulation/m_particle_cloud.fpp
• src/simulation/m_start_up.fpp
• src/simulation/m_time_steppers.fpp

Findings:

• src/simulation/m_ibm.fpp: The PR wires in "TEMPORARY DEBUG INSTRUMENTATION" that is not gated behind any build/debug flag and runs in every production build. s_debug_log_ib_divergence is called unconditionally from m_time_steppers.fpp whenever moving_immersed_boundary_flag is set, and once it detects 2 distinct timesteps where cross-rank IB state disagrees it calls call s_mpi_abort(...) — any real moving-IB, multi-rank run that hits this condition (which the comments imply is a currently-unresolved, reproducible bug) will hard-abort. This must not ship in a merged PR.
• src/simulation/m_ibm.fpp (s_communicate_ib_forces, both ACCUM and BACKPROP directions): the GPU offload directives ($:GPU_PARALLEL_LOOP / $:END_GPU_PARALLEL_LOOP / $:GPU_UPDATE) around the force/torque pack and unpack loops have been commented out "to test whether it's the source of the ... corruption bug", replaced by per-send/per-recv file open/write/close calls gated on trace_active. This disables GPU acceleration for IB force communication and adds blocking file I/O inside the MPI exchange hot path — a correctness-neutral but severe performance regression, and inconsistent with this same PR raising num_local_ibs_max/num_ib_patches_max_namelist 4x (the "these loops are tiny, no perf concern" rationale no longer holds at that scale).
• src/simulation/m_ibm.fpp: new module-level integer, save :: dbg_t_step, dbg_divergence_count, dbg_last_divergent_t_step variables are global mutable state shared across calls/ranks, which CLAUDE.md's Critical Rules explicitly forbid ("NEVER ... global save variables"). They also make s_compute_ib_forces/s_communicate_ib_forces implicitly stateful across timesteps.

[danieljvickers]

This is not ready for PR, and I will @ maintainers when this is ready.

[sbryngelson]

Merged master to resolve the conflict from #1762. Beyond the marked hunk, three calls in the new stability-report code in m_data_output.fpp still used the old s_compute_enthalpy / ten-argument sound-speed API and would not have compiled; they now use s_compute_cell_state and the six-argument form, and the unused H local is dropped. Verified on the IBM set (56 passed).

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_data_output.fpp	1487	+120
src/simulation/m_ibm.fpp	1309	+56
src/simulation/m_start_up.fpp	1291	+33
src/simulation/m_ib_patches.fpp	559	+19
src/simulation/m_particle_cloud.fpp	430	+19
src/common/m_derived_types.fpp	470	+1
src/simulation/m_global_parameters.fpp	772	+1
src/simulation/m_mpi_proxy.fpp	524	+1

Directory	Lines	Diff
common	10354	+1
simulation	28335	+249
total	46590	+250

[sbryngelson]

The two red reldebug lanes are eight Hypoelasticity tests, all the same bounds violation:

At line 65 of file src/common/m_finite_differences.fpp
Fortran runtime error: Index '-4' of dimension 1 of array 's_cc' below lower bound of -2

Only the reldebug lanes catch it because they compile with bounds checking. The no-debug lanes read past the array and pass.

s_compute_finite_difference_coefficients takes s_cc as dimension(-local_buff_size:q + local_buff_size), and at fd_order = 4 the loop body reads s_cc(i - 2). The new bounds start the loop at lB = -fd_number_in, so with fd_number = 2 the first iteration reads s_cc(-4) while x_cc only carries buff_size = 2 ghost cells. The offset_s branch has the same reach: max(fd_number_in, offset_s%beg) can widen past what the caller allocated.

The window is only available where buff_size exceeds fd_number — true where the IB floor has raised buff_size, not true for these hypoelastic cases where the two are equal. Clamping to that margin keeps the ghost-adjacent coefficients where they exist and collapses to the old 0:q where they do not:

--- a/src/common/m_finite_differences.fpp
+++ b/src/common/m_finite_differences.fpp
@@ -25,6 +25,7 @@
         integer                                                               :: lB, lE  !< loop bounds
+        integer                                                               :: fd_margin  !< usable ghost width
         integer, intent(in)                                                   :: q
@@ -33,16 +34,17 @@
-        ! Coefficients always extend at least fd_number_in beyond the interior on each side, so a stencil centered on a
-        ! ghost-adjacent cell (e.g. an immersed boundary near a domain boundary) has a real coefficient to read instead of
-        ! reading past the caller's allocation. offset_s, when given, widens this further (never narrows it) for callers
-        ! that need more than fd_number_in of margin.
-
-        if (present(offset_s)) then
-            lB = -max(fd_number_in, offset_s%beg)
-            lE = q + max(fd_number_in, offset_s%end)
-        else
-            lB = -fd_number_in
-            lE = q + fd_number_in
-        end if
+        ! A centered stencil at cell i reads s_cc(i - fd_number_in : i + fd_number_in), so coefficients exist only where
+        ! that window stays inside s_cc, which carries local_buff_size ghost cells. That leaves fd_margin cells past the
+        ! interior on each side: zero where buff_size equals fd_number, and positive only where buff_size has been floored
+        ! higher, which is the immersed-boundary case these ghost-adjacent coefficients are for.
+        fd_margin = max(0, local_buff_size - fd_number_in)
+
+        if (present(offset_s)) then
+            lB = -min(offset_s%beg, fd_margin)
+            lE = q + min(offset_s%end, fd_margin)
+        else
+            lB = -fd_margin
+            lE = q + fd_margin
+        end if
--- a/src/simulation/m_hypoelastic.fpp
+++ b/src/simulation/m_hypoelastic.fpp
@@ -43,6 +43,7 @@
         integer :: i
+        integer :: fd_margin_hypo  !< usable ghost width for the coefficient range
@@ -62,13 +63,15 @@
-        ! s_compute_finite_difference_coefficients always extends fd_number beyond the interior on each side
-        @:ALLOCATE(fd_coeff_x_hypo(-fd_number:fd_number,-fd_number:m + fd_number))
+        ! Match the range s_compute_finite_difference_coefficients can actually fill: it stops where the stencil would
+        ! leave x_cc's ghost region, so the margin is what buff_size has over fd_number and is zero when they are equal.
+        fd_margin_hypo = max(0, buff_size - fd_number)
+        @:ALLOCATE(fd_coeff_x_hypo(-fd_number:fd_number,-fd_margin_hypo:m + fd_margin_hypo))
         if (n > 0) then
-            @:ALLOCATE(fd_coeff_y_hypo(-fd_number:fd_number,-fd_number:n + fd_number))
+            @:ALLOCATE(fd_coeff_y_hypo(-fd_number:fd_number,-fd_margin_hypo:n + fd_margin_hypo))
         end if
         if (p > 0) then
-            @:ALLOCATE(fd_coeff_z_hypo(-fd_number:fd_number,-fd_number:p + fd_number))
+            @:ALLOCATE(fd_coeff_z_hypo(-fd_number:fd_number,-fd_margin_hypo:p + fd_margin_hypo))
         end if

The allocation has to move with it, or the margin the routine fills and the margin the caller reserved disagree in the other direction.

Checked on your head c009a9819, gfortran, ./mfc.sh test --debug: without the patch 43CADBB8 reproduces the error above; with it all eight pass (43CADBB8 DC6D7467 5D405BF9 DA44D68D 879C490D E04B6502 C531DC93 CC89283B).

One thing to weigh: nothing reads outside the interior today. Every consumer in m_hypoelastic indexes fd_coeff_x_hypo(r, k) from do k = 0, m, so the widened range is unused either way, and reverting both hunks to the previous 0:m bounds would be a smaller fix. I wrote it as the clamp because it keeps what you were reaching for rather than removing it, but if the ghost-adjacent coefficients are not what the IB work ends up needing, the revert is the cleaner end state.

[sbryngelson]

Following up on the finite-difference diff above - a fuller read of the branch turned up something more serious than that one, plus a few smaller things.

Out-of-bounds writes in s_handoff_ib_ownership whenever ib_neighborhood_radius > 1

The buffers are still sized for a radius-1 stencil (src/simulation/m_ibm.fpp:1478-1482, unchanged by this PR):

! 26 neighbors max in 3D (8 in 2D); each gets its own recv buffer
integer, parameter             :: max_nbrs = 26
character(len=1), allocatable  :: send_buf(:), recv_bufs(:,:)
integer, dimension(2*max_nbrs) :: requests
integer, dimension(max_nbrs)   :: recv_neighbor_list

but the loops this PR changes (:1559-1561 and :1576-1578) now run -ib_neighborhood_radius .. +ib_neighborhood_radius with nbr_idx incremented once per non-zero offset, so the count is (2R+1)**num_dims - 1:

radius	num_dims	neighbours	vs dimension(26)
1	3	26	fits exactly
2	3	124	overflow
3	2	48	overflow
3	3	342	overflow

So recv_neighbor_list(nbr_idx) and requests(nreqs) are written past the end of two stack arrays, and MPI_IRECV(recv_bufs(:, nbr_idx), ...) posts receives past the end of a heap allocation that is (buf_size, 26) (:1528). Stack corruption plus MPI writing into memory it does not own - a crash or a silent wrong answer depending on what the stack happens to hold.

This is reachable rather than hypothetical. ib_neighborhood_radius is a case parameter with {"min": 0} (toolchain/mfc/params/definitions.py:381), and when left at 0 it is auto-computed at src/simulation/m_start_up.fpp:1606:

ib_neighborhood_radius = max(1, ceiling(1.1_wp*max_ib_bound/(min_rank_width)))

which exceeds 1 exactly when a body is large relative to a rank's width - the regime this PR exists to handle. There is also now an example in the tree that sets ib_neighborhood_radius: 3 explicitly.

Everything else on the radius-aware path is sized correctly - ib_neighbor_ranks is allocated (-ax:ax) for ax = ib_neighborhood_radius at m_start_up.fpp:1427 - so it is just these three arrays that were missed. Deriving max_nbrs from the radius (and allocating recv_bufs and requests from it) rather than hardcoding 26 would close it.

Two smaller things in the same hunks:

• The tag formula tag = 200 + (dx+1)*9 + (dy+1)*3 + (dz+1) (:1563, :1580) is a base-3 encoding, injective only for offsets in -1..1. At radius 2 distinct offsets collide - (0,2,-1) and (1,-1,-1) both give 205. Harmless today because every MPI_ISEND sends the same send_buf, so a mismatched pair is indistinguishable, but it becomes a real bug the moment the payload is per-neighbour.
• The unpack bound ((2*ib_neighborhood_radius + 1)**num_dims) - 1 (:1592) disagrees with the enumeration in 1D: the dy loop runs -R..R regardless of dimensionality, so nbr_idx reaches 8 at radius 1 while the unpack loop stops at 2. Harmless now (those entries are MPI_PROC_NULL) but the two should come from one expression.

Worth a look, lower severity

• s_ibm_correct_state now allocates and frees 13 device arrays per RK stage (:231-234, :507-508). @:ALLOCATE expands to allocate + GPU_ENTER_DATA(create=...), and device allocation synchronises on both CUDA and HIP, so at ~3 stages per step this is likely the dominant new cost for a PR aimed at throughput. Hoisting them to module scope sized at max_num_gps, next to ghost_points, would remove it. Relatedly, the phase-1 copy-out and phase-2 copy-in of r_IP/v_IP/pb_IP/mv_IP/nmom_IP/presb_IP/massv_IP/Ys_IP are unconditional while the interpolation that fills them is guarded by bubbles_euler/qbmm/chemistry - about 20 reals per ghost point per stage of uninitialised data copied in every non-bubble non-chemistry case.
• The two-kernel split itself is right: phase 2 reads no cell of q_prim_vf/pb_in/mv_in other than its own (j,k,l), so the read-after-write hazard is genuinely removed, and the race was real (the eta/sum(eta) fallback in s_compute_interpolation_coeffs fires exactly when an image point's stencil lies inside a neighbouring particle). One behavioural note: on CPU the old serial loop let ghost point i see corrections from 1..i-1 in that fallback, and it no longer does, so overlapping-particle cases can shift.
• if (ib) moved inside the CFL kernel (src/simulation/m_data_output.fpp:202-203). Every other ib_markers reference in the tree hoists the flag outside the compute construct (m_rhs.fpp:813). When ib is false, ib_markers%sf is never allocated and @:ACC_SETUP_SFs never ran, so the region references an allocatable component absent from the device present table under default(present). nvfortran will likely tolerate the untaken branch; CCE's lookup is stricter. Also worth noting the mask excludes ghost-point cells as well as deep-interior ones, which is wider than "values interior to the IBs" and means a diverging ghost state no longer trips the stability guard.
• s_report_icfl_violation (m_data_output.fpp:333-457) is called unconditionally from the per-time-step path whenever a rank's local ICFL exceeds 1: ~130 lines of print from every offending rank, a full GPU_UPDATE(host=...) of every q_prim_vf plus ib_markers, and a second global barrier per step. It also reads patch_ib%force/torque/vel on the host with no GPU_UPDATE(host='[patch_ib]'), so for moving IBs the particle state it reports can be a step stale and actively mislead. The PR description says this part is a work in progress - gating it behind MFC_DEBUG or dropping it before merge would be my suggestion.
• s_read_ib_restart_data: nothing checks num_gbl_ibs <= num_ib_patches_max_namelist before patch_ib(i) is written (m_start_up.fpp:1230-1240); the @:PROHIBIT that enforces it lives in s_reduce_ib_patch_array, which runs after. And the broadcast loop issues six MPI_BCAST calls per global IB - 3600 tiny collectives for a 600-particle bed. The underlying fix is right and valuable, though: the old do i = 1, num_ibs read only namelist patches, so restart with particle clouds was silently dropping every particle.
• Pre-existing but in a block this PR rewrote: r is missing from the phase-2 private list (:275-278, used at :486 for QBMM non-polytropic). OpenACC predetermines compute-region scalars private; OpenMP offload does not - the documented trap in .claude/rules/common-pitfalls.md. The list was copied from master, so the PR did not introduce it, but the rewrite is the natural moment.
• ib_gbl_idx_lookup(tmp_patch%gbl_patch_id) = num_ibs at :1604 looks redundant next to s_update_ib_lookup() two lines later, but it is what lets s_get_neighborhood_idx see a patch added earlier in the same unpack loop, de-duplicating one that arrives from two neighbours - which becomes possible as soon as radius > 1 aliases two offsets onto one rank. Worth a comment so it survives a cleanup.
• num_ib_patches_max_namelist 54000 -> 216000 and num_local_ibs_max 2000 -> 8000 are not free: ib_patch_parameters measures 512 bytes, so patch_ib goes from 26 MB to 105 MB in host BSS and as a static device allocation, for every run including a single sphere. Several per-step transfers are sized by the constant rather than by num_ibs - m_ibm.fpp:1496/1521/1608 do whole-array GPU_UPDATEs where other sites already use patch_ib(1:num_ibs), and :1527-1528 allocates ~110 MB of send/recv buffers every step. Sizing those by the actual count is the change that would let the ceiling rise cheaply.

Verified clean, for what it is worth: all four fd_coeff_* allocation sites match the new coefficient bounds including the offset_s path, @:ALLOCATE/@:DEALLOCATE pairing in s_ibm_correct_state is balanced, no wp/stp mixing is introduced, the shell_axis plumbing is complete and self-consistent end to end, s_restart_particle_clouds's id accounting matches both packers, and $:GPU_UPDATE(device='[num_gps]') at :1017 is a genuine fix for a stale device loop bound.

[sbryngelson]

The two failures on this branch look like they are not yours.

Both hypoelasticity cases die the same way, and it is a hard crash rather than a tolerance miss:

At line 65 of file src/common/m_finite_differences.fpp
Fortran runtime error: Index '-4' of dimension 1 of array 's_cc' below lower bound of -2

That is the out-of-bounds finite-difference read tracked in #1856 and #1860: s_compute_ib_forces samples fd_number cells out from each interior cell, and the coefficient lookup follows that index past the end of the array. #1859 fixes it by reading at the nearest interior cell instead.

So the order is #1878 -> #1859 -> here. #1859 is currently red for an unrelated reason, the 23 model_eqns=3 failures every branch inherits from master, which #1878 clears. Once both land, rebasing this branch should take those two failures with them, without any change on your side.

Worth confirming rather than assuming, though: if the crash persists after #1859 lands, it is a second site and I would want to see it.

[sbryngelson]

Correcting my earlier comment: I said #1859 would fix these two failures. #1859 was closed without merging, so that fix is not coming from there, and you should not wait on it.

The diagnosis itself still holds. Both cases die the same way:

At line 65 of file src/common/m_finite_differences.fpp
Fortran runtime error: Index '-4' of dimension 1 of array 's_cc' below lower bound of -2

and the underlying bug is still open as #1856 and #1860. src/simulation/m_ibm.fpp on master still has the do l = -fd_number, fd_number loops that walk the coefficient lookup past the end of the array, so nothing about this has changed on master.

So the two failures here are still not caused by this branch, but they will not clear on their own either. Either #1856 gets a new fix, or this branch needs to work around it.

Sorry for the bad steer — I should have checked that #1859 had actually landed before pointing you at it.

[sbryngelson]

Following up properly, because my last two comments were both partly wrong and this branch is where the fix actually lives.

#1859 was closed in favour of the fix on this branch, so "wait for #1859" was exactly backwards. Your fix is here, and the widening in s_compute_finite_difference_coefficients is the right shape. The crash is a detail inside it.

The loop bounds were widened:

lB = -max(fd_number_in, offset_s%beg)
lE = q + max(fd_number_in, offset_s%end)

but s_cc is still declared over the caller's buffer:

real(wp), dimension(-local_buff_size:q + local_buff_size), intent(in) :: s_cc

and the 4th-order stencil reaches fd_number cells either side of i:

fd_coeff_s(-2, i) = 1._wp/(s_cc(i - 2) - 8._wp*s_cc(i - 1) - s_cc(i + 2) + 8._wp*s_cc(i + 1))

With local_buff_size = 2 and fd_number = 2, the loop now starts at i = -2 and the first thing it does is read s_cc(-4), two below the array. That is the failure verbatim:

At line 65 of file src/common/m_finite_differences.fpp
Fortran runtime error: Index '-4' of dimension 1 of array 's_cc' below lower bound of -2

So the coefficient array got wider but the coordinate array it reads did not, and the loop now runs where the stencil has nothing to stand on. The safe start is i >= -local_buff_size + fd_number_in (0 for this case), with the mirror at the top end.

Two ways out, and the choice is yours because they mean different things physically:

1. Clamp the loop to where the stencil fits, and let coefficients outside that range stay whatever the caller pre-filled. Contained, but a cell that close to the edge still has no real coefficient, so it only moves the problem if the IB force integral genuinely reaches there.
2. Pass a wider s_cc. The callers already hold the coordinate arrays out to the halo, so widening the dummy's declared bounds and the actual passed section gives the stencil real data to read. That is the one that matches "sample into the halo region", which is what you described when the other PR was closed.

Both examples in #1856 are still worth running against whichever you pick: rank-invariance of Fx on 1 vs 4 ranks is the property that actually has to hold, and a clamp that silently leaves an edge coefficient unset can pass a bounds check while still making the force decomposition-dependent.