Local Aware IBM

[danieljvickers]

Description

The current code has an issue scaling much past 10k particles due to limitations in the MPIAllReduceSum in the IB force computation. This PR attempts to alleviate this by limiting the number of IBs any given rank can be aware of to its neighbors. This turns the AllReduce compute to a MPI neighbor computation, removing the communication bottlneck. To support this, a massive overhaul of IB ownership between ranks was required.

Type of change

• Refactor

Testing

TBD

Checklist

• I added or updated tests for new behavior

AI code reviews

Reviews are not triggered automatically. To request a review, comment on the PR:

• @coderabbitai review — incremental review (new changes only)
• @coderabbitai full review — full review from scratch
• /review — Qodo review
• /improve — Qodo code suggestions
• @claude full review — Claude full review (also triggers on PR open/reopen/ready)
• Add label claude-full-review — Claude full review via label

[github-actions]

Claude Code Review

Head SHA: 9879df3

Files changed:

• 21
• src/simulation/m_ibm.fpp
• src/simulation/m_start_up.fpp
• src/simulation/m_collisions.fpp
• src/simulation/m_global_parameters.fpp
• src/common/m_constants.fpp
• src/common/m_derived_types.fpp
• src/common/m_mpi_common.fpp
• src/post_process/m_data_output.fpp
• src/simulation/m_ib_patches.fpp
• src/simulation/m_time_steppers.fpp

Findings:

Critical: Stack overflow in s_reduce_ib_patch_array (src/simulation/m_start_up.fpp)

subroutine s_reduce_ib_patch_array()
    type(ib_patch_parameters), dimension(num_ib_patches_max) :: patch_ib_gbl

patch_ib_gbl is a local (stack-allocated) variable dimensioned by num_ib_patches_max, which was just raised from 50,000 to 2,050,000 in src/common/m_constants.fpp. ib_patch_parameters contains a character(LEN=pathlen_max) field (pathlen_max=400) plus ~500 bytes of numeric fields — roughly 900 bytes per element. At 2,050,000 elements this is ~1.8 GB on the call stack, which will segfault on every platform (typical stack limits are 8–64 MB). This variable must be heap-allocated (allocate/deallocate).

---

High: Debug print left in production code (src/simulation/m_ibm.fpp)

print *, proc_rank, " New Owner ", patch_ib(k)%gbl_patch_id  ! TODO :: REMOVE THIS DEBUG PRINT

This fires on every ownership transfer for every locally-tracked IB, generating unbounded output at scale. The ! TODO :: REMOVE THIS DEBUG PRINT marker confirms it is unintentional.

---

Medium: Commented-out code in m_collisions.fpp leaves pid2 lookup absent (src/simulation/m_collisions.fpp)

! call s_get_neighborhood_idx(pid1, pid1) ! global patch ID -> local index call s_get_neighborhood_idx(pid2, pid2)
if (pid1 <= 0 .or. pid2 <= 0) cycle

The comment text contains call s_get_neighborhood_idx(pid2, pid2) — what appears to be a second intended lookup that was accidentally folded into the comment instead of being left as executable code. Neither lookup actually executes: pid1 and pid2 are the raw decoded global IDs from s_decode_patch_periodicity, not local indices. The guard if (pid1 <= 0 .or. pid2 <= 0) cycle would never trigger for valid global IDs (which are ≥ 1), meaning the subsequent patch_ib(pid1) and patch_ib(pid2) accesses use global IDs as local array indices, which is out-of-bounds when num_ibs < num_gbl_ibs. If this is intentional (global IDs equal local indices in the no-MPI / single-rank path), that should be documented; otherwise both lookups need to be uncommented.

---

Low: GPU_UPDATE(device=[patch_ib(ib_idx)%moment]) removed without replacement (src/simulation/m_ibm.fpp)

-            patch_ib(ib_marker)%moment = moment*patch_ib(ib_marker)%mass/(count*cell_volume)
-            $:GPU_UPDATE(device='[patch_ib(ib_marker)%moment]')
+            patch_ib(ib_idx)%moment = moment*patch_ib(ib_idx)%mass/(count*cell_volume)

The per-field device update was removed. If patch_ib(i)%moment is consumed on the GPU before the next full GPU_UPDATE(device=[patch_ib(1:num_ibs)]) (e.g., in a subsequent call to s_ibm_correct_state within the same time step), the device will use a stale value. Verify that a covering full-struct update always precedes any GPU use of %moment after s_compute_moment_of_inertia is called from s_propagate_immersed_boundaries.

[sbryngelson]

I generated and pushed the missing golden file for BA186DDF (2D -> Example -> mibm_particle_bed) directly to your branch — that should fix the CI failure.

While reviewing the diff I found several other issues worth addressing:

---

Critical (wrong results)

f_neighborhood_ranks_own_location guard is num_procs > 2 instead of num_procs > 1 (m_collisions.fpp)
On a 2-rank run the condition is false, so both ranks think they own every collision and forces get double-applied. Should match f_local_rank_owns_location's guard.

recv_forces_snap not zeroed before y/z accumulation passes (m_ibm.fpp, s_communicate_ib_forces)
recv_forces_snap = 0._wp lives only inside the x-dimension block. For y and z it accumulates stale values from the previous dimension, making the double-count subtraction remove too much force.

s_handoff_ib_ownership can silently drop newly received patches (m_ibm.fpp)
After the compact-and-delete loop, ib_gbl_idx_lookup still reflects the old layout (it is rebuilt after unpacking). s_get_neighborhood_idx is therefore testing the pre-handoff table — a just-evicted patch may still appear as resident, causing the newly received authoritative state to be discarded.

---

High (MPI / GPU correctness)

recv_forces_snap, recv_torques_snap, recv_ids, recv_ft allocated with plain allocate, not @:ALLOCATE (m_ibm.fpp)
These arrays are used inside GPU_PARALLEL_LOOP via copyin clauses, but plain allocate skips GPU_ENTER_DATA. On OpenACC builds the GPU sees uninitialized device memory.

patch_ib re-allocated with plain allocate/deallocate while GPU_DECLARE is active (m_global_parameters.fpp + m_start_up.fpp)
patch_ib was changed from a static array to allocatable, but the device pointer registered by GPU_DECLARE at startup goes stale every time s_reduce_ib_patch_array or the MPI proxy resizes it. Each reallocation should use @:ALLOCATE/@:DEALLOCATE.

neighbor_domain_x/y/z computed on host but never GPU_UPDATE'd to device (m_start_up.fpp)
These are GPU_DECLARE'd and used inside f_neighborhood_ranks_own_location (a GPU_ROUTINE). After get_neighbor_bounds fills them on the CPU, a GPU_UPDATE(device=[neighbor_domain_x, neighbor_domain_y, neighbor_domain_z]) is needed.

s_compute_ib_neighbor_ranks: nreqs not reset between outer loop iterations (m_start_up.fpp)
After MPI_WAITALL, nreqs retains its end value going into the next k iteration. With ib_neighborhood_radius > 2 this overflows the dimension(52) requests array.

Post-process SILO writes called from all ranks (post_process/m_data_output.fpp)
DBPUTMVAR to dbroot is now called from all ranks, not just rank 0. In file_per_process=F mode this produces duplicate or corrupted multi-mesh entries.

---

Medium (robustness / performance)

num_local_ibs_max = 2000 hard cap with no overflow check — local_ib_patch_ids(num_local_ibs) = num_ibs writes out of bounds silently if more than 2000 particles land on one rank. Needs an @:PROHIBIT.

ib_force_send_buf/recv_buf heap-allocated on every RK stage call (m_ibm.fpp) — these should be module-level buffers sized once in s_ibm_setup rather than allocated and freed on every timestep.

moving_immersed_boundary_flag not recomputed after s_handoff_ib_ownership — set once at startup, then stale after each ownership handoff cycle when moving IBs migrate between ranks.

---

Low (code quality)

• s_get_neighborhood_idx has an unused local variable i (compiler warning)
• f_xorshift can return exactly 1.0 (violates the [0, 1) contract in the docstring — particle placed at xmax)
• particle_bed_parameters type defined in simulation-only m_global_parameters.fpp; convention is to put derived types in m_derived_types.fpp

---

The two highest-priority fixes are the num_procs > 2 off-by-one (one-character fix, guaranteed wrong on 2-rank jobs) and recv_forces_snap not being zeroed per-dimension (wrong forces in any 2D/3D multi-rank run).

[danieljvickers]

Addressing AI comments:

All of the critical issues are incorrect for one reason or another. I have gone through and verified that the relevant chunks of code are in place to address the concerns it has.

As for the neighbor boundaries, they are never called in a GPU region, so this comment is mute. I have removed the GPU routine call to help with compilation in response.

The file per process F claim is just false, and I tested that the code runs find without it. But it is true that we do not need other ranks to be doing any of this calculation regardless. I extended the rank 0 check in response.

nreqs is in fact reset between loop iterations.

Snapshot arrays are in fact reset between loop iterations

It is confused about the need to GPU allocate arrays for patch_ib. Because the array is used with a declare, but never actually allocated on the GPU, we just have the GPU pointer. When we perform the first copy to the GPU to allocate the array size, we have already done the reshape. This means that the GPU array is properly sized. This is important because the size of the patch_ib array is significant in the new code and could causes issues with limiting problem size because the GPU must be able to hold the global patch_ib array, even though it would then immediately be deallocated and reallocated to a smaller array. If the number of ib patches grows to the extent that it costs many for GB of memory, this could cause memory limitations. The current implementation should be the correct way to go.

the moving IBM flag should not be recomputed. You want it to be the same on all ranks. There is argument that it should be moved to the case level for case optimizaiton later, but that is a future PR.

Force arrays getting moved is going to be handled in a future PR that already has a github issue I created yesterday.

Prohibit for number of local IBs resolved.

In summary, two of these I made minor adjustments for. One I reject for reasons that I believe are valid and better for long-term code support. All others are incorrect.

[danieljvickers]

General Note: Mac errors are on sections of the code that were not touched by this PR (bubbles), and it is a random dyl linker error. I get the same error on reldebug on a chemsitry case. I get the same error locally, but on a different set of tests. This seems like a MacOS instability issue, and not related to the particular PR. Going to attempt to rerun. If ti fails again, but on different tests, then it looks like the MacOS test suite checks may be broken.

Edit: Rerunning the job changed the jobs which failed. This appears to be a spurious failure for MacOS specifically. It appears to be specific to the branch, as other branch are not experiencing it. Will investigate.