Ib collisions

[danieljvickers]

Description

Added a soft-sphere collision model for sphere/circle IB patches to collide with each other and walls. Collisions with friction were validated against references to obtain reasonable results. This code was also run on 64 MPI ranks to simulate a 10k particle case in a shock tube. Adding this feature also required the ability to restart from previous data, which was added and integrated with existing file reading. That required some refactoring of the existing IB state writing, and will require more modifications to write in parallel.

Fixes #1108

Type of change

• New feature
• Refactor
• Documentation
• Other: describe

Testing

Ran a 10k particle case as well as the new example to compare rebound angles and collision energies.

Checklist

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

See the developer guide for full coding standards.

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

[codecov]

Codecov Report

❌ Patch coverage is 34.25926% with 142 lines in your changes missing coverage. Please review.
✅ Project coverage is 64.62%. Comparing base (ec94064) to head (4dd02a3).
⚠️ Report is 4 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_collisions.fpp	25.17%	102 Missing and 5 partials ⚠️
src/simulation/m_start_up.fpp	17.64%	28 Missing ⚠️
src/simulation/m_time_steppers.fpp	63.63%	2 Missing and 2 partials ⚠️
src/simulation/m_ibm.fpp	66.66%	1 Missing and 1 partial ⚠️
src/simulation/m_ib_patches.fpp	75.00%	0 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1348      +/-   ##
==========================================
- Coverage   64.90%   64.62%   -0.28%     
==========================================
  Files          70       71       +1     
  Lines       18254    18407     +153     
  Branches     1508     1516       +8     
==========================================
+ Hits        11847    11895      +48     
- Misses       5444     5555     +111     
+ Partials      963      957       -6     

☔ View full report in Codecov by Sentry.
📢 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.

[sbryngelson]

sync with master (hopefully gracefully)

[danieljvickers]

@sbryngelson this is ready for review

[sbryngelson]

Code review

Thorough review of the soft-sphere collision model. Found 17 issues across 4 severity levels.

---

Critical

1. #:if defined('MFC_MPI') instead of #ifdef MFC_MPI — collision forces multiplied by num_procs on multi-rank runs

f_local_rank_owns_collision uses Fypp syntax #:if defined('MFC_MPI') to guard the multi-rank ownership logic. But MFC_MPI is a C preprocessor define, not a Fypp variable — the entire codebase (96 occurrences across 23 files) uses #ifdef MFC_MPI. The Fypp conditional always evaluates to false, so the function unconditionally compiles to owns_collision = .true. on every rank. After s_mpi_allreduce_vectors_sum, collision forces are multiplied by num_procs. Single-rank runs are unaffected, which masks this during testing.

Fix: replace #:if defined('MFC_MPI') / #:else / #:endif with #ifdef MFC_MPI / #else / #endif.

https://github.com/MFlowCode/MFC/blob/957fcd2adb9b807fce99821739cc1ea737fd2d3a/src/simulation/m_collisions.fpp#L710-L712

2. Wrong sign on pid2 torque — Newton's third law violated for rotational DOF

In s_apply_ib_collision_forces_soft_sphere:

torques(pid2, l) = torques(pid2, l) - torque(l)*patch_ib(pid2)%radius/patch_ib(pid1)%radius

Sphere 2 receives force -F_t at contact arm (-r2 * n_hat). The torque is (-r2*n_hat) x (-F_t) = r2*(n_hat x F_t) — the same sign as sphere 1's torque, not opposite. The minus sign causes sphere 2 to spin the wrong direction during oblique collisions.

Fix: change - to +.

https://github.com/MFlowCode/MFC/blob/957fcd2adb9b807fce99821739cc1ea737fd2d3a/src/simulation/m_collisions.fpp#L453

3. Division by zero when particle centroids coincide

normal_vector = centroid_2 - centroid_1
overlap_distance = R1 + R2 - norm2(normal_vector)
if (overlap_distance > 0._wp) then
    normal_vector = normal_vector/norm2(normal_vector)

If two spheres are initialized at the same location (or drift there numerically), norm2 = 0 and this produces NaN/Inf that propagates through all forces and velocities. No minimum-distance guard exists.

MFC/src/simulation/m_collisions.fpp

Lines 414 to 417 in 957fcd2

	! catch the edge case where th collision lies just outside the computational domain
	#:for X, ID in [('x', 1), ('y', 2), ('z', 3)]
	if (num_dims >= ${ID}$) then
	if (ib_bc_${X}$%beg /= BC_PERIODIC) then

4. ib_bc_x/y/z%end not GPU-declared for OpenACC

The GPU_DECLARE block only places %beg on device:

$:GPU_DECLARE(create='[ib_bc_x%beg, ib_bc_y%beg, ib_bc_z%beg]')

But s_detect_wall_collisions reads ib_bc_x%end, ib_bc_y%end, ib_bc_z%end inside a GPU_PARALLEL_LOOP. On OpenACC GPU builds, wall collision detection reads uninitialized device memory, silently producing wrong results.

Fix: add ib_bc_x%end, ib_bc_y%end, ib_bc_z%end to both the GPU_DECLARE and the corresponding GPU_UPDATE calls.

MFC/src/simulation/m_global_parameters.fpp

Lines 807 to 808 in 957fcd2

	! Variables to handle moving immersed boundaries, defaulting to no movement
	patch_ib(i)%moving_ibm = 0

---

Major

5. num_patches_max reduced from 1000 to 10 — silent breaking regression

The PR correctly introduces num_ib_patches_max = 50000 for IB patches, but also reduces num_patches_max from 1000 to 10. This constant bounds the patch_icpp IC patch array, which is unrelated to IB patches. Any existing case with more than 10 initial condition patches will silently fail or crash at namelist read time.

Fix: leave num_patches_max = 1000 unchanged; num_ib_patches_max already handles the IB case.

MFC/src/common/m_constants.fpp

Lines 28 to 29 in 957fcd2

	integer, parameter :: num_bc_patches_max = 10 !< Maximum number of boundary condition patches
	integer, parameter :: max_2d_fourier_modes = 10 !< Max Fourier mode index for 2D modal patch (geometry 13)

6. collision_lookup sized num_ibs * 8 — unguarded buffer overflow

@:ALLOCATE(collision_lookup(num_ibs * 8, 4))

The maximum number of unique IB-IB pairs is num_ibs*(num_ibs-1)/2. For num_ibs >= 18, this exceeds 8*num_ibs. The PR description mentions a 10k particle case — with 10,000 IBs, the maximum pairs is ~50 million while the allocated size is only 80,000. There is no bounds check before writing to collision_lookup(num_considered_collisions, ...), so exceeding the bound silently corrupts memory.

Fix: either size analytically (e.g., num_ibs*(num_ibs-1)/2) with a practical cap, or add @:ASSERT(num_considered_collisions <= size(collision_lookup, 1)) before each write.

MFC/src/simulation/m_collisions.fpp

Line 343 in 957fcd2

distance_vec = centroid_2 - centroid_1

7. raw_pairs(num_gps, 2) — stack-allocated VLA risks overflow

integer, dimension(num_gps, 2) :: raw_pairs

num_gps can reach hundreds of thousands for dense 3D simulations. This is a local stack allocation that will overflow on systems with limited stack size, and is problematic as a local array inside a GPU parallel region. Should be heap-allocated via @:ALLOCATE / @:DEALLOCATE.

https://github.com/MFlowCode/MFC/blob/957fcd2adb9b807fce99821739cc1ea737fd2d3a/src/simulation/m_collisions.fpp#L544

8. Restart reader always loads most recent state, not the requested checkpoint

s_read_ib_restart_data scans ib_state.dat from beginning to EOF, applying every record (last wins). Restarting from an intermediate checkpoint loads the most recent IB state rather than the one matching the fluid checkpoint at t_step_start, making fluid and IB state inconsistent. This defeats the stated purpose of the restart feature.

MFC/src/simulation/m_start_up.fpp

Lines 1133 to 1143 in 957fcd2

	end subroutine s_finalize_modules
	!> @brief Reads IB kinematic state from restart_data/ib_state.dat on restart. Rank 0 reads the last num_ibs records and
	!! broadcasts to all ranks. Overwrites patch_ib vel, angular_vel, angles, and centroid.
	impure subroutine s_read_ib_restart_data()
	character(len=path_len + 2*name_len) :: file_loc
	integer :: i, ios, file_unit, ib_id, ierr
	real(wp) :: time_read
	real(wp), dimension(3) :: force_read, torque_read

9. Commented-out ib_markers overwrite breaks periodic IB cases

The line ib_markers%sf(i, j, k) = patch_id was commented out in m_ibm.fpp without explanation. For periodic IB cases (xp or yp != 0), ghost-point cells retain encoded patch IDs (which include periodicity offsets) instead of plain patch IDs. s_compute_ib_forces indexes forces(ib_idx, l) where ib_idx comes from ib_markers — for periodic cases the encoded ID exceeds num_ibs, causing out-of-bounds array access in forces and patch_ib.

MFC/src/simulation/m_ibm.fpp

Line 955 in 957fcd2

end if

---

Moderate

10. Missing validation for coefficient_of_restitution and collision_time

damping_parameter = -2._wp*log(e)/collision_time
spring_stiffness = (pi**2 + log(e)**2)/(collision_time**2)

If coefficient_of_restitution <= 0, log(0) = -Inf, producing NaN forces. If collision_time <= 0, division produces Inf. Neither case_validator.py nor Fortran-side @:PROHIBIT checks enforce e in (0, 1] or collision_time > 0. The CLAUDE.md parameter rules require Fortran-side runtime validation for physics constraints.

11. Rotation matrix removal — unexplained physics change

The PR removes matmul(rotation_matrix, angular_vel) and matmul(rotation_matrix_inverse, torques) from m_ibm.fpp at two call sites without any explanation in the code or PR description. If angular velocity tracking changed from body-frame to inertial-frame, this is a physics change affecting all existing moving-IB cases, not just the new collision model.

MFC/src/simulation/m_ibm.fpp

Line 934 in 957fcd2

end if

12. Torque arm inconsistency

Velocity at the contact point uses the corrected arm R - 0.5*overlap, but the torque computation uses the full radius R. These should be consistent — either both use the corrected arm or both use the full radius.

13. Coulomb friction cap uses total normal force instead of spring-only component

tangental_force = -ib_coefficient_of_friction*norm2(normal_force)*tangental_vector

Standard DEM applies the Coulomb cap using only the elastic (spring) component: mu * k * delta. Using norm2(normal_force) (spring + dashpot) inflates the friction bound during approach. Not catastrophically wrong, but deviates from the Thornton reference and is non-standard.

14. O(N^2) collision deduplication on CPU

After the GPU-parallel detection phase, the coalescing loop does a linear scan over collision_lookup for each raw pair to check for duplicates — O(num_raw^2). For the 10k particle case mentioned in the PR description, this will be a significant bottleneck. A hash set or sort-then-unique approach would be standard.

---

Minor / Style

15. Dead code: s_detect_ib_collisions_n2 never called

This subroutine is defined but never referenced. It also doesn't populate collision_lookup columns 3/4 (encoded IDs), so calling it would produce uninitialized data in the force routine. Should be removed before merge.

https://github.com/MFlowCode/MFC/blob/957fcd2adb9b807fce99821739cc1ea737fd2d3a/src/simulation/m_collisions.fpp#L627-L660

16. Debug print * left in production code

s_read_ib_restart_data has two bare print * statements that will spam stdout on every restart run. The codebase convention is @:LOG(...) gated behind MFC_DEBUG.

MFC/src/simulation/m_start_up.fpp

Lines 1124 to 1127 in 957fcd2

	call s_finalize_viscous_module()
	end if
	call s_finalize_mpi_proxy_module()

17. Module docstring is wrong

m_collisions.fpp has @brief Ghost-node immersed boundary method: locates ghost/image points... which is m_ibm's description, not a collisions module description. Also, a hardcoded 1.0e-10_wp epsilon in f_local_rank_owns_collision for domain boundary projection assumes domain scales where 1e-10 is meaningful; micro-scale simulations would need a relative tolerance.

---

Generated with Claude Code

_{- If this code review was useful, please react with 👍. Otherwise, react with 👎.}

[danieljvickers]

I looked at the AI review, and none of it looks particularly plausible. It sugests things like the type of macro for MFC_MPI mattering, even though that was not demonstrated in any test, and chenging signs to be wrong a violate conservation of momentum. It is also simply wrong about the upper limit of array sizes. I don't see any comments here that make me think I need to modify the code now.

[sbryngelson]

just delete...