1454-projection-optimization

[danieljvickers]

Description

In many-particle cases, the limiting factor in the IBM compute is by far the time it takes to project the immersed boundaries onto the grid. This is fundamentally rooted in how we parallelize the work. The current code parallelizes of x, y, and z, but sequentially iterates through the IB patches. In cases where there are many IBs that are small, we are launching several (thousands) of GPU kernels each time step, but each kernel only works on hundreds of grid cells. Adding an option to parallelize over the thousands of particles should be significantly larger parallelism and thus optimize the projection.

In order to maintain the functionality of both parallelism methods, I need a separate set of geometry bounding checks. Since we perform a check in icpp patches and now 2 in IB patches, this is a significant amount of redundant code that must be maintained. To be somewhat forward-looking, I opted to merge all geometry checking into a single module that can be called from both forms of IB parallelism and the icpp pre_processing code. This should clean up the code nicely and significantly reduce code maintenance going forward. Since we can now change cylinder orientation with angles, I also deprecate the unneeded cylinder length checks.

The end result is the creation of a new module in common, the deletion of duplicate code, and a new parallelism path for IB patches when utilizing GPU compute.

Closes #1454, #1532, #1543

Type of change (delete unused ones)

• New feature
• Refactor

Testing

All tests pass on GNU compiler

Checklist

Check these like this [x] to indicate which of the below applies.

• I added or updated tests for new behavior
• I updated documentation if user-facing behavior changed

See the developer guide for full coding standards.

GPU changes (expand if you modified src/simulation/)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

AI code reviews

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

• @claude full review — Claude full review (also triggers on PR open/reopen/ready)
• Or add label claude-full-review — Claude full review via label

[github-actions]

Claude Code Review

Head SHA: aa507f7

Files changed:

• 15
• src/simulation/m_ib_patches.fpp
• src/common/m_patch_geometries.fpp
• src/common/m_model.fpp
• src/pre_process/m_icpp_patches.fpp
• src/simulation/m_ibm.fpp
• src/simulation/m_time_steppers.fpp
• src/simulation/m_checker.fpp
• src/simulation/m_global_parameters.fpp
• src/pre_process/m_global_parameters.fpp
• src/post_process/m_global_parameters.fpp

Findings

1. Data race in s_apply_ib_patches_ib_parallelism — non-deterministic IB marker writes when patches overlap

File: src/simulation/m_ib_patches.fpp — s_apply_ib_patches_ib_parallelism

Both the 2D and 3D variants GPU-parallelize over patch_id (outermost loop), then serially iterate over i, j, k inside each GPU thread. If two patches contain the same cell, two distinct GPU threads will write to the same ib_markers%sf(i, j, k) location concurrently without synchronization — a data race that produces non-deterministic results.

$:GPU_PARALLEL_LOOP(private='[i, il, ir, j, jl, jr, k, kl, kr, ...]', copyin='[xp, yp, zp]')
do patch_id = 1, num_ibs          ! parallelized — each thread owns a patch_id
    ...
    do k = kl, kr
        do j = jl, jr
            do i = il, ir
                ...
                if (...) ib_markers%sf(i, j, k) = encoded_patch_id  ! races with other threads

The old s_apply_ib_patches_grid_cell_parallelism serialized the patch loop (outer serial, inner GPU), so later patches deterministically overwrote earlier ones. The new IB-parallelism path breaks that invariant. The checker in m_checker.fpp only validates that ib is enabled — it does not require non-overlapping patches. Any case where patches touch (including particle-cloud spheres that may be placed adjacently) will silently produce wrong IB markers.

---

2. Dropped k <= Np_local upper-bound guard in f_is_inside_airfoil — latent out-of-bounds while-loop

File: src/common/m_patch_geometries.fpp — f_is_inside_airfoil

The deleted s_ib_airfoil guarded its upper-surface while-loop with k <= Np_local:

! old s_ib_airfoil — guarded:
do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < xy_local(1) .and. k <= Np_local)

The replacement function omits the guard entirely:

! new f_is_inside_airfoil — unguarded:
do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < x)
    k = k + 1
end do

The early-return check x <= ib_airfoil(airfoil_id)%c provides protection in exact arithmetic, but if a grid cell centre lands at c only by floating-point accumulation (rounding slightly above the stored chord endpoint), the loop will walk off the end of the array. Because f_is_inside_airfoil is tagged GPU_ROUTINE(parallelism='[seq]') and called inside GPU_PARALLEL_LOOP, an out-of-bounds access here silently corrupts GPU memory.

[danieljvickers]

All AI review comments are now irrelevant because of significant changes that have occurred between now and that review.

I doubt that this will pass the test suite on the first try, but on the off chance that it does we should not yet merge it. A data product of computational optimization performance should be presented before this PR is merged. Otherwise, it is unnecessary refactoring of the code.

[codecov]

Codecov Report

❌ Patch coverage is 59.55056% with 108 lines in your changes missing coverage. Please review.
✅ Project coverage is 60.65%. Comparing base (825adb2) to head (d1764ed).

Files with missing lines	Patch %	Lines
src/simulation/m_ib_patches.fpp	63.08%	58 Missing and 21 partials ⚠️
src/common/m_patch_geometries.fpp	52.77%	9 Missing and 8 partials ⚠️
src/simulation/m_ibm.fpp	0.00%	6 Missing ⚠️
src/pre_process/m_icpp_patches.fpp	50.00%	2 Missing and 2 partials ⚠️
src/common/m_model.fpp	0.00%	1 Missing ⚠️
src/simulation/m_time_steppers.fpp	0.00%	0 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1549      +/-   ##
==========================================
- Coverage   60.95%   60.65%   -0.30%     
==========================================
  Files          82       83       +1     
  Lines       19926    19856      -70     
  Branches     2924     2945      +21     
==========================================
- Hits        12145    12043     -102     
- Misses       5805     5820      +15     
- Partials     1976     1993      +17     

☔ 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]

Successful runtime comparison on Frontier. This plot shows run times for the s_apply_ib_patches subroutines.

In this test, I ran with 8k particles being projected onto the grid. Master is the current master branch (before this is merged) as a baseline of performance. The other two plots show the two types of parallelism now available from this branch.

the grid parallelism is the current Master branch functionality. This demonstrates that the performance is equivalent. The IB parallelism is the new parallelism introduced by this feature. In this particular case of 8k sphere with only about 300 grid cells per sphere, there was a 150x performance improvement over the previous parallelism. This this subroutine is called 3x per time step in the test, this is a 327 ms time reduction PER TIME STEP.

This plot should demonstrate that this feature is useful and should be merged. I will now move on to resolving any remaining issues with the test suite and hopefully getting this on master.

[sbryngelson]

I've pushed a merge of master into this branch (ac1681a, a true merge commit — no history rewritten) to resolve the conflicts from the recently-merged global-parameters/registry refactors (#1550–#1556). What changed and why:

src/post_process/m_global_parameters.fpp — master moved shear_*, proc_coords, start_idx into the new shared m_global_parameters_common; kept your ib_airfoil/ib_airfoil_grids stubs that the new common m_patch_geometries module needs to compile for post_process.

src/simulation/m_global_parameters.fpp — the GPU_DECLARE lines for ib, num_ibs, and ib_coefficient_of_friction now live in m_global_parameters_common, so the manual list reduces to ib_airfoil_grids. Your addition of stl_models to the GPU declare list was relocated to the new parameter pipeline: derived-type declarations (and their sim-side GPU_DECLARE) are auto-generated from the TYPED_DECLS table in toolchain/mfc/params/definitions.py — I flipped the stl_models gpu flag there to True, which emits the same declare in the generated decls.

Heads-up for future commits on this branch: Fortran parameter declarations, namelist bindings, and MPI broadcasts are now generated from toolchain/mfc/params/definitions.py at build time (see docs/documentation/contributing.md). Your many_ib_patch_parallelism registration came through cleanly and its broadcast is auto-generated — no action needed.

Verification on my end: ./mfc.sh format, ./mfc.sh precheck, the toolchain pytest suite (342 passed), and a full CPU build of all three targets all pass on the merged tree. Nothing remains on your side; CI should run on the updated head.

[danieljvickers]

I am very glad you have merged this for me because I screwed up the previous one significantly, lol.

[danieljvickers]

@sbryngelson I assume you are going to want to wait merging this until after we resolve the master branch AMD issue with the 2D viscous shock tube case. So I am going to just table it here. But otherwise it looks like we pass the full test suite and it is ready* to merge.

[sbryngelson]

@sbryngelson I assume you are going to want to wait merging this until after we resolve the master branch AMD issue with the 2D viscous shock tube case. So I am going to just table it here. But otherwise it looks like we pass the full test suite and it is ready* to merge.

@danieljvickers noted!