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.