Fix periodic ib issues

[danieljvickers]

Description

I was made aware of several bugs in the IB MPI communication in periodic cases. This branch resolves all of the issues that I have been made aware of. They include

1. periodic IBs not having the radial vector account for their periodicity
2. Force integration can double-count if the Ib neighbors itself or has the same neighbor periodically
3. Collisions can be missed as neither processor believes they own the collision
4. There is drift in IB information between ranks.

The good news is that the fix for number 2 is also a computational optimization in small cases by limiting the total number of communications to the minimum required to successfully sum. So we have a slight upside.

A new example will be added that should check for this behavior and ensure that it does not break again.

Type of change (delete unused ones)

• Bug fix

[github-actions]

Claude Code Review

Head SHA: 4a2cb4c

Files changed:

• 8
• src/common/m_global_parameters_common.fpp
• src/common/m_mpi_common.fpp
• src/post_process/m_start_up.fpp
• src/simulation/m_collisions.fpp
• src/simulation/m_derived_variables.fpp
• src/simulation/m_ib_patches.fpp
• src/simulation/m_ibm.fpp
• src/simulation/m_time_steppers.fpp

Findings:

1. s_compute_moment_of_inertia: general-case numerical integration is permanently dead (correctness)

src/simulation/m_ibm.fpp — The refactoring renamed the axis dummy argument to moment and made it intent(out), but inside the else (non-analytic geometry) branch the pre-existing line moment = 0._wp now zeros the output variable before the old axis-magnitude check is reached:

else  ! general case
    count = 0
    moment = 0._wp          ! zeros the intent(out) dummy
    ...
    else if (sqrt(sum(moment**2)) < sgm_eps) then   ! 0 < sgm_eps — always true
        moment = 1._wp
        return              ! numerical integration below never runs
    else
        normal_axis = moment/sqrt(sum(moment**2))   ! unreachable
    end if
    ! ... volume-integration loop — never reached for 3D general geometry
    moment = moment*patch%mass/(count*cell_volume)

In the old code sqrt(sum(axis**2)) < sgm_eps checked whether the incoming angular velocity was degenerate. The rename conflated the axis-direction input (formerly intent(in) axis) with the scalar output (formerly a local variable moment). Every 3D IB with non-analytic geometry (not circle/rectangle/ellipse/sphere) now silently returns moment = 1._wp. This is also a Fortran standard violation: reading an intent(out) dummy before defining it.

2. patch_ib(i)%moment struct field never updated after refactoring (correctness)

src/simulation/m_ibm.fpp / src/simulation/m_time_steppers.fpp — Old code stored results via patch_ib(ib_idx)%moment = .... New code returns the result through the moment dummy, whose actual argument at both call sites is patch_ib(i)%angular_vel:

! m_time_steppers.fpp:759
if (num_dims == 3) call s_compute_moment_of_inertia(patch_ib(i), patch_ib(i)%angular_vel)
! m_time_steppers.fpp:761
patch_ib(i)%angular_vel = patch_ib(i)%angular_vel/patch_ib(i)%moment

After the call, angular_vel holds the computed scalar moment broadcast across all 3 components; patch_ib(i)%moment is never written and remains stale from initialisation. The division on line 761 uses the stale denominator, and the angular velocity direction is lost. The same issue applies to the m_ibm.fpp call site.

3. ib_gbl_idx_lookup moved to m_collisions without allocation or deallocation (memory management)

src/simulation/m_collisions.fpp — The variable is declared allocatable with a GPU directive:

integer, dimension(:), allocatable :: ib_gbl_idx_lookup
$:GPU_DECLARE(create='[ib_gbl_idx_lookup]')

s_initialize_collisions_module allocates only collision_lookup and wall_overlap_distances; ib_gbl_idx_lookup has no @:ALLOCATE. s_finalize_collisions_module has no matching @:DEALLOCATE. Any access to ib_gbl_idx_lookup(decoded_pairs(1)) will segfault at runtime. The allocation and GPU-setup code that previously lived in m_ibm must be moved here.

[danieljvickers]

Fixed all of the bugs that I have been made aware of, and now all minimum-reproducers I was given behave as anticipated.

Final list of bugfixes:

1. Updated the n2 collision lookup to take into account periodicity
2. Added an NVHPC/24.5- compiler workaround for the moment of inertia calculation
3. Reduced number of communications to the minimum required to accomplish the correct summation, which prevents periodic cases from wrapping back on themselves and double-counting
4. Uses the correct periodic location when computing the radial vector in various subroutines.

Will add an example to make sure we don't regress on this.

[danieljvickers]

Did a lot of debugging here and determined that this is a machine-precision issue with the levelset. If I compute the levelset on the CPU always and gopy to the GPU, the GPU test passes the test suite. This makes me quite sad. The test is pretty poorly posed. I can get around this for now, but it is another time that this issue has cropped up.

[codecov]

Codecov Report

❌ Patch coverage is 45.16129% with 34 lines in your changes missing coverage. Please review.
✅ Project coverage is 60.43%. Comparing base (9906fa6) to head (1f3438f).
⚠️ Report is 3 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_ibm.fpp	48.64%	13 Missing and 6 partials ⚠️
src/simulation/m_collisions.fpp	38.09%	12 Missing and 1 partial ⚠️
src/simulation/m_ib_patches.fpp	66.66%	0 Missing and 1 partial ⚠️
src/simulation/m_time_steppers.fpp	0.00%	0 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1618      +/-   ##
==========================================
+ Coverage   60.39%   60.43%   +0.04%     
==========================================
  Files          83       83              
  Lines       19854    19868      +14     
  Branches     2955     2956       +1     
==========================================
+ Hits        11990    12007      +17     
+ Misses       5863     5860       -3     
  Partials     2001     2001              

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

@sbryngelson can I get approval for case op tests?

[sbryngelson]

you only need approval for benchmark not case opt. though i think case opt only runs once your other tests pass (automatically regardless).

[sbryngelson]

nice work. bench is clean. code is clean. merging.