Skip to content

Ib collisions#1348

Open
danieljvickers wants to merge 63 commits intoMFlowCode:masterfrom
danieljvickers:ib-collisions
Open

Ib collisions#1348
danieljvickers wants to merge 63 commits intoMFlowCode:masterfrom
danieljvickers:ib-collisions

Conversation

@danieljvickers
Copy link
Copy Markdown
Member

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

@github-actions
Copy link
Copy Markdown

github-actions bot commented Apr 1, 2026

Claude Code Review

Head SHA: a7874f1

Files changed:

  • 19
  • src/common/m_constants.fpp
  • src/simulation/m_collisions.fpp
  • src/simulation/m_ibm.fpp
  • src/simulation/m_global_parameters.fpp
  • src/simulation/m_start_up.fpp
  • src/simulation/m_data_output.fpp
  • src/simulation/m_mpi_proxy.fpp
  • src/simulation/m_time_steppers.fpp
  • src/pre_process/m_global_parameters.fpp
  • src/pre_process/m_mpi_proxy.fpp

Findings

1. num_patches_max reduced from 1000 to 10 — breaks existing IC patch cases

File: src/common/m_constants.fpp

-    integer, parameter  :: num_patches_max = 1000              !< Maximum number of IC patches
+    integer, parameter  :: num_patches_max = 10                !< Maximum number of IC patches

num_patches_max governs the patch_icpp(num_patches_max) array used for initial condition patches. Reducing it from 1000 to 10 is a regression: any existing case file referencing patch_icpp(11) or above will crash with an array-bounds violation at runtime (or silent memory corruption if bounds checking is off). The intent of separating the IB patch limit (num_ib_patches_max = 50000) from the IC patch limit is correct in principle, but the new IC limit of 10 is far too small — the old value of 1000 should be retained.


2. Debug print * statements left in s_read_ib_restart_data

File: src/simulation/m_start_up.fpp

print *, "IB Restart File Exists: ", file_exist
...
print *, "iostat ", ios

These two print * statements are diagnostic artifacts that will fire on rank 0 for every IB restart. They bypass the @:LOG macro (which is gated behind MFC_DEBUG) and will pollute production run output.


3. Torque sign error for the second IB in IB-IB collisions

File: src/simulation/m_collisions.fpps_apply_ib_collision_forces_soft_sphere

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

torque is computed as cross(normal_vector * pid1%radius, tangental_force) — the lever arm is pid1's radius. For pid2, both the contact normal (-normal_vector) and the friction force (-tangental_force) reverse sign, giving cross(-normal_vector * pid2%radius, -tangental_force) = +cross(normal_vector * pid2%radius, tangental_force), which equals +torque * (pid2%radius / pid1%radius). The current expression negates this, resulting in a torque that has the wrong sign on the second IB.


4. wall_overlap_distances not reset each time step — stale overlap values

File: src/simulation/m_collisions.fpps_detect_wall_collisions

The function writes to wall_overlap_distances(patch_id, idx) only when the corresponding wall BC matches BC_SLIP_WALL or BC_NO_SLIP_WALL. Boundaries with other BC types (e.g. outflow) never zero out their entry. Since the array is initialized once at module startup and never zeroed on subsequent calls, an IB that overlapped a wall in one time step but moves away will retain the stale non-zero distance, causing a phantom collision force to be applied every subsequent step. The function should zero the relevant entries at the top of the loop before recomputing.

@danieljvickers danieljvickers marked this pull request as ready for review April 6, 2026 18:18
@qodo-code-review
Copy link
Copy Markdown
Contributor

Review Summary by Qodo

Add soft-sphere collision model and restart capability for immersed boundaries

✨ Enhancement

Grey Divider

Walkthroughs

Description
• Added soft-sphere collision model for immersed boundary particles
  - Detects collisions between IBs and walls using overlap distances
  - Detects collisions between IBs using ghost-point neighbor search
  - Applies normal and tangential forces with friction
• Implemented IB restart capability from previous simulation states
  - Reads IB kinematic state from restart_data/ib_state.dat
  - Preserves collision history by appending to state file
• Separated immersed boundary patch array from initial condition patches
  - Increased num_ib_patches_max to 50000 for scalability
  - Reduced num_patches_max to 10 for initial conditions only
• Added collision parameters and configuration options
  - collision_model, coefficient_of_restitution, collision_time, ib_coefficient_of_friction
Diagram
flowchart LR
  A["IB Collision Detection"] -->|wall overlaps| B["Wall Collision Forces"]
  A -->|IB-IB overlaps| C["IB Collision Forces"]
  B --> D["Apply Forces & Torques"]
  C --> D
  D --> E["Update IB Motion"]
  E --> F["Write IB State File"]
  F --> G["Restart from State"]
Loading

Grey Divider

File Changes

1. src/simulation/m_collisions.fpp ✨ Enhancement +445/-0

New collision detection and force application module

src/simulation/m_collisions.fpp


2. src/simulation/m_ibm.fpp ✨ Enhancement +15/-24

Integrate collision forces into IBM workflow

src/simulation/m_ibm.fpp


3. src/simulation/m_start_up.fpp ✨ Enhancement +64/-1

Add IB restart and collision initialization

src/simulation/m_start_up.fpp


View more (18)
4. src/simulation/m_global_parameters.fpp ⚙️ Configuration changes +21/-7

Add collision parameters and preserve original BCs

src/simulation/m_global_parameters.fpp


5. src/simulation/m_time_steppers.fpp ✨ Enhancement +11/-12

Update IB state writing for restart capability

src/simulation/m_time_steppers.fpp


6. src/simulation/m_mpi_proxy.fpp ✨ Enhancement +3/-2

Broadcast collision parameters across MPI ranks

src/simulation/m_mpi_proxy.fpp


7. src/simulation/m_ib_patches.fpp 🐞 Bug fix +9/-3

Use preserved BCs and skip out-of-domain particles

src/simulation/m_ib_patches.fpp


8. src/simulation/m_data_output.fpp ✨ Enhancement +11/-3

Write IB state to restart_data directory with append mode

src/simulation/m_data_output.fpp


9. src/common/m_helper.fpp ✨ Enhancement +17/-2

Add cross product subroutine for collision calculations

src/common/m_helper.fpp


10. src/common/m_constants.fpp ⚙️ Configuration changes +2/-1

Separate IC and IB patch array bounds

src/common/m_constants.fpp


11. src/pre_process/m_global_parameters.fpp ⚙️ Configuration changes +6/-6

Update patch_ib array bound to use num_ib_patches_max

src/pre_process/m_global_parameters.fpp


12. src/pre_process/m_mpi_proxy.fpp 🐞 Bug fix +3/-2

Fix IB parameter broadcast loop structure

src/pre_process/m_mpi_proxy.fpp


13. src/post_process/m_data_input.f90 🐞 Bug fix +3/-2

Fix MPI offset type declarations for IB data reading

src/post_process/m_data_input.f90


14. toolchain/mfc/params/definitions.py ⚙️ Configuration changes +17/-3

Register collision parameters and separate IB patch limit

toolchain/mfc/params/definitions.py


15. toolchain/mfc/case_validator.py ⚙️ Configuration changes +3/-3

Validate collision parameters and use num_ib_patches_max

toolchain/mfc/case_validator.py


16. toolchain/mfc/params/descriptions.py 📝 Documentation +4/-0

Add descriptions for collision parameters

toolchain/mfc/params/descriptions.py


17. docs/documentation/case.md 📝 Documentation +12/-0

Document collision parameters and IB configuration

docs/documentation/case.md


18. docs/module_categories.json 📝 Documentation +2/-1

Add m_collisions to module categories

docs/module_categories.json


19. examples/3D_mibm_sphere_head_on_collision/case.py 🧪 Tests +129/-0

New example case for sphere collision verification

examples/3D_mibm_sphere_head_on_collision/case.py


20. tests/AA49A8BC/golden-metadata.txt 🧪 Tests +120/-0

Test metadata for collision example case

tests/AA49A8BC/golden-metadata.txt


21. tests/AA49A8BC/golden.txt Additional files +66/-0

...

tests/AA49A8BC/golden.txt


Grey Divider

Qodo Logo

@qodo-code-review
Copy link
Copy Markdown
Contributor

qodo-code-review bot commented Apr 6, 2026

Code Review by Qodo

🐞 Bugs (7)   📘 Rule violations (2)   📎 Requirement gaps (0)   🎨 UX Issues (0)
🐞\ ≡ Correctness (2) ☼ Reliability (4) ➹ Performance (1)
📘\ ≡ Correctness (1) ☼ Reliability (1)

Grey Divider


Action required

1. 1.0e-10 missing _wp 📘
Description
The new collision code assigns unkinded numeric literals (e.g., 1.0e-10, 0) into real(wp)
variables, breaking the project precision discipline and risking inconsistent precision across
builds. This violates the requirement to use wp/stp-kinded literals and generic intrinsics
consistently.
Code

src/simulation/m_collisions.fpp[R416-421]

+                            ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it
+                            if (collision_location(${ID}$) < ${X}$_domain%beg) then
+                                projected_location(${ID}$) = ${X}$_domain%beg
+                            else if (${X}$_domain%end < collision_location(${ID}$)) then
+                                projected_location(${ID}$) = ${X}$_domain%end - 1.0e-10
+                            end if
Evidence
PR Compliance ID 3 requires wp/stp kinds and *_wp-kinded literals; the added code uses unkinded
literals in real(wp) computations/assignments (e.g., subtracting 1.0e-10 from a real(wp)
domain bound and setting a real(wp) array to 0).

CLAUDE.md
src/simulation/m_collisions.fpp[416-421]
src/simulation/m_collisions.fpp[31-44]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
`real(wp)` code in `m_collisions` uses unkinded numeric literals (e.g., `1.0e-10`, `0`), which violates the project precision discipline.
## Issue Context
The code operates on `real(wp)` state (domain bounds, overlaps) and should use kinded literals like `0._wp` / `1.0e-10_wp` to avoid precision inconsistencies.
## Fix Focus Areas
- src/simulation/m_collisions.fpp[31-44]
- src/simulation/m_collisions.fpp[416-421]

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


2. collision_model lacks validation 📘
Description
New collision parameters are introduced but no corresponding physics/validity constraints are
enforced in toolchain/mfc/case_validator.py, allowing invalid values (e.g.,
coefficient_of_restitution <= 0 leading to log(e) issues). This violates the requirement to add
new parameters to required validation locations when constraints exist.
Code

toolchain/mfc/params/definitions.py[R915-918]

+    _r("collision_model", INT, {"ib"})
+    _r("coefficient_of_restitution", REAL, {"ib"})
+    _r("collision_time", REAL, {"ib"})
+    _r("ib_coefficient_of_friction", REAL, {"ib"})
Evidence
PR Compliance ID 7 requires updating case_validator.py when new parameters have physics
constraints; the documentation specifies ranges/allowed values for the new collision parameters, but
case_validator.py contains no checks for them while the params are registered in the toolchain.

CLAUDE.md
docs/documentation/case.md[368-374]
toolchain/mfc/params/definitions.py[915-918]
toolchain/mfc/case_validator.py[554-571]
src/simulation/m_collisions.fpp[33-38]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
Collision parameters (`collision_model`, `coefficient_of_restitution`, `collision_time`, `ib_coefficient_of_friction`) are added, but the Python `case_validator.py` does not validate their allowed ranges/values.
## Issue Context
`m_collisions` computes `log(e)` where `e = coefficient_of_restitution`, so invalid `e` (<= 0, > 1) can cause NaNs or undefined behavior. Documentation states allowed ranges/values, implying physics constraints that should be enforced.
## Fix Focus Areas
- toolchain/mfc/params/definitions.py[915-918]
- toolchain/mfc/case_validator.py[554-571]
- docs/documentation/case.md[368-374]
- src/simulation/m_collisions.fpp[33-38]

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


3. IC patches capped at 10 🐞
Description
num_patches_max was reduced to 10, capping patch_icpp and other IC-sized arrays to 10 and breaking
any case with >10 initial-condition patches (validation and/or runtime array bounds). This also
desynchronizes the toolchain’s fallback constants, which still assume 1000 when Fortran sources are
unavailable.
Code

src/common/m_constants.fpp[R26-27]

+    integer, parameter  :: num_patches_max = 10                !< Maximum number of IC patches
+    integer, parameter  :: num_ib_patches_max = 50000          !< Maximum number of immersed boundary patches (patch_ib)
Evidence
The Fortran compile-time cap is now 10, and IC patch data structures directly depend on
num_patches_max; the toolchain fallback constants are not updated accordingly.

src/common/m_constants.fpp[24-28]
src/common/m_derived_types.fpp[187-210]
toolchain/mfc/params/namelist_parser.py[506-513]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
`num_patches_max` was reduced to 10, which caps **initial condition** patch arrays (`patch_icpp`, `alter_patch`, etc.) to 10 and breaks cases needing more IC patches; toolchain fallback constants also become inconsistent.
### Issue Context
This PR adds a separate `num_ib_patches_max` for IB patches, but `num_patches_max` still governs IC patches.
### Fix Focus Areas
- src/common/m_constants.fpp[24-28]
- toolchain/mfc/params/namelist_parser.py[506-513]
- (If intentional) add explicit release-note/docs update; otherwise restore prior value (e.g., 1000).

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


View more (5)
4. IB checks stop at 10 🐞
Description
The pre_process IB patch validator still loops to num_patches_max, so with the new IC/IB separation
and num_patches_max=10 it only validates the first 10 IB patches and can miss invalid IB definitions
beyond that. This can allow out-of-range/invalid IB parameters through pre_process checks and fail
later at runtime.
Code

src/pre_process/m_global_parameters.fpp[R174-175]

+    type(ib_patch_parameters), dimension(num_ib_patches_max) :: patch_ib  !< Immersed boundary patch parameters
+    type(vec3_dt), allocatable, dimension(:)                 :: airfoil_grid_u, airfoil_grid_l
Evidence
patch_ib is now sized by num_ib_patches_max, but the validator still iterates using num_patches_max;
with num_patches_max reduced, the validator no longer covers all possible IB indices.

src/pre_process/m_global_parameters.fpp[169-176]
src/pre_process/m_check_ib_patches.fpp[32-71]
src/common/m_constants.fpp[24-28]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
`m_check_ib_patches` validates IB patches using `do i = 1, num_patches_max`, which now (a) is an IC patch limit and (b) was reduced to 10. As a result, IB patches beyond 10 are not validated.
### Issue Context
`patch_ib` is now dimensioned by `num_ib_patches_max` in pre_process, so validation must follow that bound (or at least validate all `1..num_ibs` and assert the rest are default).
### Fix Focus Areas
- src/pre_process/m_check_ib_patches.fpp[32-71]
- src/common/m_constants.fpp[24-28]
- src/pre_process/m_global_parameters.fpp[169-176]

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


5. Fallback constants missing NIB 🐞
Description
The toolchain now requires num_ib_patches_max, but the fallback constants used when
src/common/m_constants.fpp isn’t shipped do not define it, causing the parameter registry
initialization to throw at runtime. This breaks installed/Homebrew-style toolchain usage.
Code

toolchain/mfc/params/definitions.py[R30-31]

+NUM_PATCHES_MAX = _fc("num_patches_max")  # patch_icpp (Fortran array bound)
+NIB = _fc("num_ib_patches_max")  # patch_ib (Fortran array bound)
Evidence
definitions.py now hard-requires num_ib_patches_max via _fc(), but namelist_parser’s fallback
constants dict does not include it, so toolchain boot fails when it can’t parse Fortran sources.

toolchain/mfc/params/definitions.py[15-32]
toolchain/mfc/params/namelist_parser.py[506-513]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
Toolchain startup can fail in environments where `src/common/m_constants.fpp` is not shipped, because fallback constants omit `num_ib_patches_max` while `definitions.py` now requires it.
### Issue Context
`_fc("num_ib_patches_max")` raises when the constant is absent.
### Fix Focus Areas
- toolchain/mfc/params/namelist_parser.py[506-513]
- toolchain/mfc/params/definitions.py[15-32]
### Expected change
Add `"num_ib_patches_max": <Fortran value>` to `_FALLBACK_CONSTANTS` and keep it in sync with `src/common/m_constants.fpp`.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


6. raw_pairs stack allocation 🐞
Description
s_detect_ib_collisions declares raw_pairs(num_gps,2) as an automatic local array, which can overflow
stack or cause large per-step overhead for large num_gps. This is called during force evaluation, so
it risks hard crashes on large runs.
Code

src/simulation/m_collisions.fpp[R239-241]

+        ! Temporary array to hold all detected pairs (with potential duplicates)
+        integer, dimension(num_gps, 2) :: raw_pairs
+        integer                        :: num_raw, local_num_raw
Evidence
raw_pairs is sized at runtime by num_gps and created as a local automatic array, which is
stack-backed on typical Fortran compilers and not safe for large sizes in a hot timestep loop.

src/simulation/m_collisions.fpp[227-277]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
`raw_pairs(num_gps,2)` is an automatic local array and can blow the stack or incur heavy per-timestep allocation costs for large `num_gps`.
### Issue Context
This routine runs in the timestep loop (via force computation).
### Fix Focus Areas
- src/simulation/m_collisions.fpp[227-277]
### Expected change
Move `raw_pairs` to an allocatable module-scope buffer sized to a safe upper bound (or allocate/deallocate on the heap inside the routine), and reuse it across calls.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


7. Zero-distance collision NaNs 🐞
Description
When two IB centroids coincide, norm2(normal_vector) is zero and the normal normalization divides by
zero, producing NaNs that propagate into forces/torques. This can occur from bad initial conditions
or restart data and is not guarded against.
Code

src/simulation/m_collisions.fpp[R110-114]

+            normal_vector = centroid_2 - centroid_1
+            overlap_distance = patch_ib(pid1)%radius + patch_ib(pid2)%radius - norm2(normal_vector)
+            if (overlap_distance > 0._wp) then  ! if the two patches are close enough to collide
+                normal_vector = normal_vector/norm2(normal_vector)
+                if (f_local_rank_owns_collision(centroid_1)) then
Evidence
The code normalizes the centroid-to-centroid vector whenever overlap_distance>0, but
overlap_distance can be positive even when the vector magnitude is zero (coincident centroids).

src/simulation/m_collisions.fpp[109-117]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
Normal vector normalization can divide by zero when centroids are identical.
### Issue Context
`overlap_distance > 0` does not imply `norm2(normal_vector) > 0`.
### Fix Focus Areas
- src/simulation/m_collisions.fpp[109-117]
### Expected change
Compute `dist = norm2(normal_vector)` once and:
- if `dist <= eps`, skip collision or choose a stable fallback normal
- else normalize by `dist`.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


8. Restart reader out-of-bounds 🐞
Description
s_read_ib_restart_data does not check open() iostat for failure and does not bounds-check ib_id
before writing patch_ib(ib_id), risking memory corruption on malformed files or when num_ibs differs
between runs. It also prints debug output unconditionally on restart.
Code

src/simulation/m_start_up.fpp[R1147-1169]

+        if (proc_rank == 0) then
+            file_loc = trim(case_dir) // '/restart_data/ib_state.dat'
+            inquire (FILE=trim(file_loc), EXIST=file_exist)
+            print *, "IB Restart File Exists: ", file_exist
+            if (file_exist) then
+                open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
+                print *, "iostat ", ios
+            else
+                call s_mpi_abort('Cannot open IB state file for restart: ' // trim(file_loc))
+            end if
+
+            ! Read all records; the last num_ibs records are the final state
+            do
+                read (file_unit, iostat=ios) time_read, ib_id, force_read, torque_read, vel_read, angular_vel_read, angles_read, &
+                      & xc_read, yc_read, zc_read
+                if (ios /= 0) exit
+                patch_ib(ib_id)%vel = vel_read
+                patch_ib(ib_id)%angular_vel = angular_vel_read
+                patch_ib(ib_id)%angles = angles_read
+                patch_ib(ib_id)%x_centroid = xc_read
+                patch_ib(ib_id)%y_centroid = yc_read
+                patch_ib(ib_id)%z_centroid = zc_read
+            end do
Evidence
The reader prints iostat but does not abort on ios!=0 and directly indexes patch_ib(ib_id) for every
record read. The writer emits one record per i=1..num_ibs, so changing num_ibs or partial/corrupt
files can lead to invalid ib_id values and out-of-bounds writes.

src/simulation/m_start_up.fpp[1147-1172]
src/simulation/m_data_output.fpp[919-930]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
Restart reader can corrupt memory by writing `patch_ib(ib_id)` without validating `ib_id`, and ignores open() failures.
### Issue Context
Restart files can be incomplete/corrupt; users may also change `num_ibs` between runs.
### Fix Focus Areas
- src/simulation/m_start_up.fpp[1135-1185]
- src/simulation/m_data_output.fpp[919-930]
### Expected change
- After `open(..., iostat=ios)`, abort if `ios /= 0`.
- Before assigning `patch_ib(ib_id)%...`, verify `1 <= ib_id <= num_ibs`; otherwise skip/abort with a clear message.
- Remove unconditional `print *` debug lines or gate behind a verbosity flag.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools



Remediation recommended

9. Collision lookup O(n^2) alloc 🐞
Description
collision_lookup is allocated as num_ibs*(num_ibs-1)/2, which becomes extremely large for the
intended big-particle cases and can cause per-rank OOM or heavy memory pressure. This allocation is
also disproportionate to the algorithm, which only produces O(num_gps) candidate pairs per step.
Code

src/simulation/m_collisions.fpp[R40-41]

+        @:ALLOCATE(collision_lookup(num_ibs * (num_ibs-1) / 2, 4))
+        @:ALLOCATE(wall_overlap_distances(num_ibs, 6))
Evidence
The module allocates a dense all-pairs capacity array irrespective of the actual number of detected
collision candidates, which is bounded by ghost points / local neighborhood detection.

src/simulation/m_collisions.fpp[31-44]
src/simulation/m_collisions.fpp[225-321]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
`collision_lookup` is allocated with worst-case all-pairs size `num_ibs*(num_ibs-1)/2`, which becomes very large for 10k+ IBs and can OOM.
### Issue Context
`s_detect_ib_collisions` can only generate up to O(num_gps) raw pairs per step; allocating all-pairs capacity is unnecessary.
### Fix Focus Areas
- src/simulation/m_collisions.fpp[31-44]
- src/simulation/m_collisions.fpp[225-321]
### Expected change
Allocate collision storage proportional to expected collisions (e.g., `num_gps` or a tunable cap with graceful fallback) rather than O(num_ibs^2).

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


Grey Divider

ⓘ The new review experience is currently in Beta. Learn more

Grey Divider

Qodo Logo

Comment on lines +416 to +421
! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it
if (collision_location(${ID}$) < ${X}$_domain%beg) then
projected_location(${ID}$) = ${X}$_domain%beg
else if (${X}$_domain%end < collision_location(${ID}$)) then
projected_location(${ID}$) = ${X}$_domain%end - 1.0e-10
end if
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

1. 1.0e-10 missing _wp 📘 Rule violation ≡ Correctness

The new collision code assigns unkinded numeric literals (e.g., 1.0e-10, 0) into real(wp)
variables, breaking the project precision discipline and risking inconsistent precision across
builds. This violates the requirement to use wp/stp-kinded literals and generic intrinsics
consistently.
Agent Prompt
## Issue description
`real(wp)` code in `m_collisions` uses unkinded numeric literals (e.g., `1.0e-10`, `0`), which violates the project precision discipline.

## Issue Context
The code operates on `real(wp)` state (domain bounds, overlaps) and should use kinded literals like `0._wp` / `1.0e-10_wp` to avoid precision inconsistencies.

## Fix Focus Areas
- src/simulation/m_collisions.fpp[31-44]
- src/simulation/m_collisions.fpp[416-421]

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +915 to +918
_r("collision_model", INT, {"ib"})
_r("coefficient_of_restitution", REAL, {"ib"})
_r("collision_time", REAL, {"ib"})
_r("ib_coefficient_of_friction", REAL, {"ib"})
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

2. collision_model lacks validation 📘 Rule violation ☼ Reliability

New collision parameters are introduced but no corresponding physics/validity constraints are
enforced in toolchain/mfc/case_validator.py, allowing invalid values (e.g.,
coefficient_of_restitution <= 0 leading to log(e) issues). This violates the requirement to add
new parameters to required validation locations when constraints exist.
Agent Prompt
## Issue description
Collision parameters (`collision_model`, `coefficient_of_restitution`, `collision_time`, `ib_coefficient_of_friction`) are added, but the Python `case_validator.py` does not validate their allowed ranges/values.

## Issue Context
`m_collisions` computes `log(e)` where `e = coefficient_of_restitution`, so invalid `e` (<= 0, > 1) can cause NaNs or undefined behavior. Documentation states allowed ranges/values, implying physics constraints that should be enforced.

## Fix Focus Areas
- toolchain/mfc/params/definitions.py[915-918]
- toolchain/mfc/case_validator.py[554-571]
- docs/documentation/case.md[368-374]
- src/simulation/m_collisions.fpp[33-38]

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +174 to +175
type(ib_patch_parameters), dimension(num_ib_patches_max) :: patch_ib !< Immersed boundary patch parameters
type(vec3_dt), allocatable, dimension(:) :: airfoil_grid_u, airfoil_grid_l
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

4. Ib checks stop at 10 🐞 Bug ≡ Correctness

The pre_process IB patch validator still loops to num_patches_max, so with the new IC/IB separation
and num_patches_max=10 it only validates the first 10 IB patches and can miss invalid IB definitions
beyond that. This can allow out-of-range/invalid IB parameters through pre_process checks and fail
later at runtime.
Agent Prompt
### Issue description
`m_check_ib_patches` validates IB patches using `do i = 1, num_patches_max`, which now (a) is an IC patch limit and (b) was reduced to 10. As a result, IB patches beyond 10 are not validated.

### Issue Context
`patch_ib` is now dimensioned by `num_ib_patches_max` in pre_process, so validation must follow that bound (or at least validate all `1..num_ibs` and assert the rest are default).

### Fix Focus Areas
- src/pre_process/m_check_ib_patches.fpp[32-71]
- src/common/m_constants.fpp[24-28]
- src/pre_process/m_global_parameters.fpp[169-176]

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +30 to +31
NUM_PATCHES_MAX = _fc("num_patches_max") # patch_icpp (Fortran array bound)
NIB = _fc("num_ib_patches_max") # patch_ib (Fortran array bound)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

5. Fallback constants missing nib 🐞 Bug ☼ Reliability

The toolchain now requires num_ib_patches_max, but the fallback constants used when
src/common/m_constants.fpp isn’t shipped do not define it, causing the parameter registry
initialization to throw at runtime. This breaks installed/Homebrew-style toolchain usage.
Agent Prompt
### Issue description
Toolchain startup can fail in environments where `src/common/m_constants.fpp` is not shipped, because fallback constants omit `num_ib_patches_max` while `definitions.py` now requires it.

### Issue Context
`_fc("num_ib_patches_max")` raises when the constant is absent.

### Fix Focus Areas
- toolchain/mfc/params/namelist_parser.py[506-513]
- toolchain/mfc/params/definitions.py[15-32]

### Expected change
Add `"num_ib_patches_max": <Fortran value>` to `_FALLBACK_CONSTANTS` and keep it in sync with `src/common/m_constants.fpp`.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +239 to +241
! Temporary array to hold all detected pairs (with potential duplicates)
integer, dimension(num_gps, 2) :: raw_pairs
integer :: num_raw, local_num_raw
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

6. Raw_pairs stack allocation 🐞 Bug ☼ Reliability

s_detect_ib_collisions declares raw_pairs(num_gps,2) as an automatic local array, which can overflow
stack or cause large per-step overhead for large num_gps. This is called during force evaluation, so
it risks hard crashes on large runs.
Agent Prompt
### Issue description
`raw_pairs(num_gps,2)` is an automatic local array and can blow the stack or incur heavy per-timestep allocation costs for large `num_gps`.

### Issue Context
This routine runs in the timestep loop (via force computation).

### Fix Focus Areas
- src/simulation/m_collisions.fpp[227-277]

### Expected change
Move `raw_pairs` to an allocatable module-scope buffer sized to a safe upper bound (or allocate/deallocate on the heap inside the routine), and reuse it across calls.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +110 to +114
normal_vector = centroid_2 - centroid_1
overlap_distance = patch_ib(pid1)%radius + patch_ib(pid2)%radius - norm2(normal_vector)
if (overlap_distance > 0._wp) then ! if the two patches are close enough to collide
normal_vector = normal_vector/norm2(normal_vector)
if (f_local_rank_owns_collision(centroid_1)) then
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

7. Zero-distance collision nans 🐞 Bug ☼ Reliability

When two IB centroids coincide, norm2(normal_vector) is zero and the normal normalization divides by
zero, producing NaNs that propagate into forces/torques. This can occur from bad initial conditions
or restart data and is not guarded against.
Agent Prompt
### Issue description
Normal vector normalization can divide by zero when centroids are identical.

### Issue Context
`overlap_distance > 0` does not imply `norm2(normal_vector) > 0`.

### Fix Focus Areas
- src/simulation/m_collisions.fpp[109-117]

### Expected change
Compute `dist = norm2(normal_vector)` once and:
- if `dist <= eps`, skip collision or choose a stable fallback normal
- else normalize by `dist`.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

Comment on lines +1147 to +1169
if (proc_rank == 0) then
file_loc = trim(case_dir) // '/restart_data/ib_state.dat'
inquire (FILE=trim(file_loc), EXIST=file_exist)
print *, "IB Restart File Exists: ", file_exist
if (file_exist) then
open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
print *, "iostat ", ios
else
call s_mpi_abort('Cannot open IB state file for restart: ' // trim(file_loc))
end if

! Read all records; the last num_ibs records are the final state
do
read (file_unit, iostat=ios) time_read, ib_id, force_read, torque_read, vel_read, angular_vel_read, angles_read, &
& xc_read, yc_read, zc_read
if (ios /= 0) exit
patch_ib(ib_id)%vel = vel_read
patch_ib(ib_id)%angular_vel = angular_vel_read
patch_ib(ib_id)%angles = angles_read
patch_ib(ib_id)%x_centroid = xc_read
patch_ib(ib_id)%y_centroid = yc_read
patch_ib(ib_id)%z_centroid = zc_read
end do
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Action required

8. Restart reader out-of-bounds 🐞 Bug ☼ Reliability

s_read_ib_restart_data does not check open() iostat for failure and does not bounds-check ib_id
before writing patch_ib(ib_id), risking memory corruption on malformed files or when num_ibs differs
between runs. It also prints debug output unconditionally on restart.
Agent Prompt
### Issue description
Restart reader can corrupt memory by writing `patch_ib(ib_id)` without validating `ib_id`, and ignores open() failures.

### Issue Context
Restart files can be incomplete/corrupt; users may also change `num_ibs` between runs.

### Fix Focus Areas
- src/simulation/m_start_up.fpp[1135-1185]
- src/simulation/m_data_output.fpp[919-930]

### Expected change
- After `open(..., iostat=ios)`, abort if `ios /= 0`.
- Before assigning `patch_ib(ib_id)%...`, verify `1 <= ib_id <= num_ibs`; otherwise skip/abort with a clear message.
- Remove unconditional `print *` debug lines or gate behind a verbosity flag.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools

@coderabbitai
Copy link
Copy Markdown
Contributor

coderabbitai bot commented Apr 6, 2026

📝 Walkthrough

Walkthrough

This pull request introduces collision physics modeling for immersed boundary simulations. A new module m_collisions implements collision detection and force application for inter-body and wall collisions using a soft-sphere model. Four new immersed boundary parameters—coefficient_of_restitution, collision_model, collision_time, and ib_coefficient_of_friction—are added to case file configuration and global parameters. The maximum immersed boundary patch capacity is increased from 1000 to 50000 via a new constant num_ib_patches_max. Infrastructure updates include a cross-product subroutine in the helper module, immersed boundary restart file persistence, periodic boundary condition storage, and MPI broadcasting of collision parameters. Documentation is updated to specify the new parameters, and a 3D sphere collision example case is provided.

🚥 Pre-merge checks | ✅ 5
✅ Passed checks (5 passed)
Check name Status Explanation
Title check ✅ Passed The title 'Ib collisions' is directly related to the main changes; it succinctly describes the core feature being added.
Description check ✅ Passed The PR description covers motivation (soft-sphere collision model), testing approach (10k particle case and example), and acknowledges incomplete checklist items and documentation needs.
Linked Issues check ✅ Passed The PR implements both collision model functionality [m_collisions.fpp] and restart-data support [m_start_up.fpp, m_data_output.fpp] as required by issue #1108.
Out of Scope Changes check ✅ Passed All changes align with objectives: collision implementation, parameter additions, restart support, and documentation updates directly support soft-sphere collisions and restart capability.
Docstring Coverage ✅ Passed Docstring coverage is 100.00% which is sufficient. The required threshold is 80.00%.

✏️ Tip: You can configure your own custom pre-merge checks in the settings.


Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands and usage tips.

Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actionable comments posted: 19

Caution

Some comments are outside the diff and can’t be posted inline due to platform limitations.

⚠️ Outside diff range comments (1)
src/simulation/m_ibm.fpp (1)

570-577: ⚠️ Potential issue | 🔴 Critical

Don't leave encoded periodic IDs in ib_markers.

ghost_points_in(local_idx)%ib_patch_id is decoded here, but ib_markers%sf now keeps the encoded periodic-image value. The rest of this module still does raw lookups like patch_ib(ib_markers%sf(...)) in s_ibm_correct_state and s_compute_ib_forces, so periodic ghost cells will pick up the wrong patch state as soon as an IB crosses a periodic boundary.

🧹 Nitpick comments (2)
tests/AA49A8BC/golden-metadata.txt (1)

13-119: Reduce machine-specific metadata in committed golden artifacts.

This file records host/local-environment identifiers (e.g., Line 13 hostname, Line 29 absolute user path, Line 119 CPU model). Since metadata is informational and not used for golden comparison, committing these details adds noisy diffs and leaks local context. Prefer excluding this metadata file from versioned test artifacts (or sanitize generation in a follow-up) and keep only validation-relevant golden data.

Based on learnings: Golden files are generated with ./mfc.sh test --generate and should not be manually edited.

src/common/m_helper.fpp (1)

300-311: Avoid duplicated cross-product logic between f_cross and s_cross_product.

This is functionally correct, but keeping two separate implementations invites divergence. Reuse f_cross in s_cross_product.

Refactor suggestion
     subroutine s_cross_product(a, b, c)

         $:GPU_ROUTINE(parallelism='[seq]')
         real(wp), intent(in)  :: a(3), b(3)
         real(wp), intent(out) :: c(3)

-        c(1) = a(2)*b(3) - a(3)*b(2)
-        c(2) = a(3)*b(1) - a(1)*b(3)
-        c(3) = a(1)*b(2) - a(2)*b(1)
+        c = f_cross(a, b)

     end subroutine s_cross_product

ℹ️ Review info
⚙️ Run configuration

Configuration used: Path: .coderabbit.yaml

Review profile: CHILL

Plan: Pro

Run ID: 701a3269-2add-4226-9ccf-dc2b22957260

📥 Commits

Reviewing files that changed from the base of the PR and between 23bbd4f and a7874f1.

📒 Files selected for processing (21)
  • docs/documentation/case.md
  • docs/module_categories.json
  • examples/3D_mibm_sphere_head_on_collision/case.py
  • src/common/m_constants.fpp
  • src/common/m_helper.fpp
  • src/post_process/m_data_input.f90
  • src/pre_process/m_global_parameters.fpp
  • src/pre_process/m_mpi_proxy.fpp
  • src/simulation/m_collisions.fpp
  • src/simulation/m_data_output.fpp
  • src/simulation/m_global_parameters.fpp
  • src/simulation/m_ib_patches.fpp
  • src/simulation/m_ibm.fpp
  • src/simulation/m_mpi_proxy.fpp
  • src/simulation/m_start_up.fpp
  • src/simulation/m_time_steppers.fpp
  • tests/AA49A8BC/golden-metadata.txt
  • tests/AA49A8BC/golden.txt
  • toolchain/mfc/case_validator.py
  • toolchain/mfc/params/definitions.py
  • toolchain/mfc/params/descriptions.py

Comment on lines +334 to +337
| `coefficient_of_restitution` | Real | A number 0 to 1 describing how elastic IB collisions are |
| `collision_model` | Integer | Integer to select the collision model being used for IB collisions. |
| `collision_time` | Real | Amount of simulation time used to resolve collisions |
| `ib_coefficient_of_friction` | Real | Coefficient of friction used in IB collisions |
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Clean up IB collision parameter wording (typos + consistency).

There are several spelling errors and slightly inconsistent range wording in these new parameter docs.

📝 Suggested edit
-- `coefficient_of_restitution` is a number from 0 (exclusive) to 1 (inclusive) describing how elastic IB collisions are. 0 is for perfectly inellastic collisions while 1 is for perfectly ellastic collisions.
+- `coefficient_of_restitution` is a number from 0 (exclusive) to 1 (inclusive) describing how elastic IB collisions are. 0 is perfectly inelastic, while 1 is perfectly elastic.

-- `collision_model` is an integer to select the collision model being used for IB collisions. Using 0 disables collisions and collisiono checking. 1 enables the soft-sphere collision model, where all IBs must be circles or sphere and those IBs can collide with each other as well as walls.
+- `collision_model` is an integer selecting the IB collision model. Using 0 disables collisions and collision checking. Using 1 enables the soft-sphere model, where IBs represented as circles/spheres can collide with each other and walls.

-- `collision_time` is approximately the amount of simulation time used to resolve collisions. This is handled by modifying the spring gonstant used to apply collision forces.
+- `collision_time` is approximately the amount of simulation time used to resolve collisions. This is handled by modifying the spring constant used to apply collision forces.

Also applies to: 368-374

Comment on lines +110 to +117
"patch_ib(1)%x_centroid": -1.0 * lead_distance * math.sin(collision_angle_radians), # get a lead up distance to the collision
"patch_ib(1)%y_centroid": radius + lead_distance * math.sin(collision_angle_radians),
"patch_ib(1)%z_centroid": 0.0,
"patch_ib(1)%radius": radius,
"patch_ib(1)%slip": "F",
"patch_ib(1)%mass": 1.0e6, # arbitrarily high mass to ignore fluid
"patch_ib(1)%vel(1)": velocity * math.sin(collision_angle_radians),
"patch_ib(1)%vel(2)": -velocity * math.cos(collision_angle_radians),
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Use cos(theta) for the initial wall-normal offset.

vel(1)=velocity*sin(theta) and vel(2)=-velocity*cos(theta) define the approach direction. Backing the sphere away by lead_distance from first contact therefore needs a wall-normal offset of lead_distance*cos(theta), not lead_distance*sin(theta). As written, the example reaches the wall sooner than the simulation_time = 2 * lead_distance / velocity + collision_time derivation assumes, so the rebound-angle validation is measuring the wrong setup.

🛠️ Minimal fix
-            "patch_ib(1)%y_centroid": radius + lead_distance * math.sin(collision_angle_radians),
+            "patch_ib(1)%y_centroid": radius + lead_distance * math.cos(collision_angle_radians),

Comment on lines +26 to +27
integer, parameter :: num_patches_max = 10 !< Maximum number of IC patches
integer, parameter :: num_ib_patches_max = 50000 !< Maximum number of immersed boundary patches (patch_ib)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Restore num_patches_max to avoid breaking existing multi-patch cases.

Line 26 drops the IC patch cap from 1000 to 10, which will reject or break many valid cases unrelated to IB collisions. Keep the legacy IC cap and use the new IB-specific cap only for patch_ib.

Proposed fix
-    integer, parameter  :: num_patches_max = 10                !< Maximum number of IC patches
+    integer, parameter  :: num_patches_max = 1000              !< Maximum number of IC patches
     integer, parameter  :: num_ib_patches_max = 50000          !< Maximum number of immersed boundary patches (patch_ib)
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
integer, parameter :: num_patches_max = 10 !< Maximum number of IC patches
integer, parameter :: num_ib_patches_max = 50000 !< Maximum number of immersed boundary patches (patch_ib)
integer, parameter :: num_patches_max = 1000 !< Maximum number of IC patches
integer, parameter :: num_ib_patches_max = 50000 !< Maximum number of immersed boundary patches (patch_ib)

Comment on lines +24 to +27
integer, allocatable, dimension(:,:) :: collision_lookup
real(wp), allocatable, dimension(:,:) :: wall_overlap_distances
real(wp) :: spring_stiffness, damping_parameter
$:GPU_DECLARE(create='[spring_stiffness, damping_parameter]')
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Add the collision arrays to GPU_DECLARE.

collision_lookup and wall_overlap_distances are both pushed to the device and then read inside GPU kernels, but only the scalar coefficients are in the module device-data declaration. On accelerator builds that leaves these allocatables outside the expected device lifetime.

Suggested change
-    $:GPU_DECLARE(create='[spring_stiffness, damping_parameter]')
+    $:GPU_DECLARE(create='[collision_lookup, wall_overlap_distances, spring_stiffness, damping_parameter]')

Based on learnings, GPU array lifecycle: declare with $:GPU_DECLARE(create='[...]'), allocate with @:ALLOCATE(), setup pointers with @:ACC_SETUP_SFs() or @:ACC_SETUP_VFs(), deallocate with @:DEALLOCATE().

Comment on lines +40 to +41
@:ALLOCATE(collision_lookup(num_ibs * (num_ibs-1) / 2, 4))
@:ALLOCATE(wall_overlap_distances(num_ibs, 6))
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🔴 Critical

🧩 Analysis chain

🏁 Script executed:

# Check the overall structure and key sections of m_collisions.fpp
head -60 src/simulation/m_collisions.fpp | cat -n

Repository: MFlowCode/MFC

Length of output: 2683


🏁 Script executed:

# Look at lines 279-280 and surrounding context
sed -n '270,290p' src/simulation/m_collisions.fpp | cat -n

Repository: MFlowCode/MFC

Length of output: 1089


🏁 Script executed:

# Check integer declarations and num_ibs constraints
rg -n "integer.*collision_lookup|integer.*wall_overlap|integer.*num_ibs" src/simulation/m_collisions.fpp | head -20

Repository: MFlowCode/MFC

Length of output: 121


🏁 Script executed:

# Check the integer kind used in the file
rg -n "implicit none|integer\s*::" src/simulation/m_collisions.fpp | head -30

Repository: MFlowCode/MFC

Length of output: 823


🏁 Script executed:

# Find num_ibs definition and constraints
rg -n "num_ibs" src/simulation/m_collisions.fpp | head -20

Repository: MFlowCode/MFC

Length of output: 571


🏁 Script executed:

# Search for num_ibs in global parameters
rg -n "num_ibs" src/common/m_global_parameters.fpp | head -20

Repository: MFlowCode/MFC

Length of output: 131


🏁 Script executed:

# Look at how collision_lookup is actually used/indexed
rg -n "collision_lookup" src/simulation/m_collisions.fpp -B 1 -A 1

Repository: MFlowCode/MFC

Length of output: 1808


🏁 Script executed:

# Check the context around where collision_lookup is zeroed and populated
sed -n '275,330p' src/simulation/m_collisions.fpp | cat -n

Repository: MFlowCode/MFC

Length of output: 3167


🏁 Script executed:

# Check if there's any validation on num_ibs
rg -n "num_ibs.*>" src/ | grep -E "(constraint|max|limit|check)" | head -10

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Look for num_ibs constraints and documentation
fd -t f "case_validator|parameters" src/common docs/ 2>/dev/null | head -20

Repository: MFlowCode/MFC

Length of output: 79


🏁 Script executed:

# Search for num_ibs limits/constraints in documentation or validation
rg -i "num_ibs|maximum.*ib|ib.*limit" docs/ | head -20

Repository: MFlowCode/MFC

Length of output: 156


🏁 Script executed:

# Verify the integer overflow calculation more carefully
python3 << 'EOF'
import math

# 32-bit signed integer max
INT_MAX_32 = 2**31 - 1

# Find the threshold where n*(n-1)/2 exceeds INT_MAX_32
# Solving: n^2 - n > 2*INT_MAX_32
# n = (1 + sqrt(1 + 8*INT_MAX_32)) / 2

n_threshold = (1 + math.sqrt(1 + 8*INT_MAX_32)) / 2
print(f"32-bit integer overflow threshold: {n_threshold:.0f}")

# Check specific values
for test_n in [46341, 46342, 65574, 65575]:
    result = test_n * (test_n - 1) // 2
    fits = "FITS" if result <= INT_MAX_32 else "OVERFLOWS"
    print(f"  n={test_n}: {test_n}*({test_n}-1)/2 = {result:,} [{fits}]")
EOF

Repository: MFlowCode/MFC

Length of output: 315


Correct the integer kind for collision_lookup and reconsider the O(n²) allocation strategy.

The default 32-bit integer overflows at num_ibs ≈ 65,536, not 46,341. However, the core issue remains: allocating num_ibs * (num_ibs-1) / 2 rows (line 40) and zeroing the entire array every step (line 280) creates O(n²) memory allocation and access. Since only num_considered_collisions rows are actually populated, this table should either use a 1D array with dynamic resizing or be bounded by a reasonable upper limit on detected pairs rather than the maximum possible pair count.

Comment on lines +572 to +577
!
if (ib) then
if (moving_immersed_boundary_flag) call s_wrap_periodic_ibs()
if (ib_state_wrt .and. (.not. moving_immersed_boundary_flag)) then
call s_compute_ib_forces(q_prim_vf, fluid_pp)
end if
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Refresh q_prim_vf before the post-RK IB force sample.

By the time this block runs, q_cons_ts(1)%vf has already been advanced through the last RK stage, but q_prim_vf still reflects the state from the last s_compute_rhs(...) call. s_compute_ib_forces() reads q_prim_vf directly, so the force/torque snapshot computed here is one stage behind unless you reconvert the final conservative state first.

Comment on lines +564 to +567
num_ib_patches_max = get_fortran_constants().get("num_ib_patches_max", 100000)
self.prohibit(
ib and num_ibs > num_patches_max,
f"num_ibs must be <= {num_patches_max} (num_patches_max in m_constants.fpp)",
ib and num_ibs > num_ib_patches_max,
f"num_ibs must be <= {num_ib_patches_max} (num_ib_patches_max in m_constants.fpp)",
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Align num_ib_patches_max fallback with Fortran constant.

The fallback is 100000, but this PR defines num_ib_patches_max as 50000. When constants are unavailable at runtime, validation can incorrectly allow cases that exceed allocated IB capacity.

Proposed fix
-        num_ib_patches_max = get_fortran_constants().get("num_ib_patches_max", 100000)
+        num_ib_patches_max = get_fortran_constants().get("num_ib_patches_max", 50000)
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
num_ib_patches_max = get_fortran_constants().get("num_ib_patches_max", 100000)
self.prohibit(
ib and num_ibs > num_patches_max,
f"num_ibs must be <= {num_patches_max} (num_patches_max in m_constants.fpp)",
ib and num_ibs > num_ib_patches_max,
f"num_ibs must be <= {num_ib_patches_max} (num_ib_patches_max in m_constants.fpp)",
num_ib_patches_max = get_fortran_constants().get("num_ib_patches_max", 50000)
self.prohibit(
ib and num_ibs > num_ib_patches_max,
f"num_ibs must be <= {num_ib_patches_max} (num_ib_patches_max in m_constants.fpp)",

Comment on lines +716 to +720
"collision_model": {
"when_set": {
"requires": ["ib", "coefficient_of_restitution", "collision_time"],
}
},
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Use when_value for collision_model dependencies.

when_set also fires for collision_model = 0, so an explicit “collisions off” input still has to provide ib, restitution, and collision time. Gate these requirements on the enabled model IDs instead, and this is also the right place to hang the numeric validator checks for the new collision knobs.

As per coding guidelines, "New parameters in toolchain/mfc/params/definitions.py must have correct types, constraints, and tags; validation in case_validator.py must cover interdependencies."

attrs=_ib_attrs,
tags=_ib_tags,
max_index=NUM_PATCHES_MAX,
max_index=NIB,
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Propagate NIB into the num_ibs constraint too.

patch_ib is now bounded by num_ib_patches_max, but CONSTRAINTS["num_ibs"] is still unbounded. That lets the toolchain accept num_ibs > NIB, and the Fortran loops over patch_ib(1:num_ibs) will walk past the array bound at runtime.

As per coding guidelines, "New parameters in toolchain/mfc/params/definitions.py must have correct types, constraints, and tags; validation in case_validator.py must cover interdependencies."

Comment on lines +253 to +256
"coefficient_of_restitution": "Real number describing the elasticity of collisions from 0 (perfectly ineleastic) to 1 (perfectly elastic)",
"collision_model": "Integer selecting the collision model being used. 0 for no collision. 1 for soft-sphere collisions",
"collision_time": "Amount of simulation time each collision will take to resolve",
"ib_coefficient_of_friction": "coefficient of friction used in IB collisions",
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Fix typos and wording in new parameter descriptions.

The new descriptions include misspellings and inconsistent phrasing in user-facing text.

✏️ Proposed text cleanup
-    "coefficient_of_restitution": "Real number describing the elasticity of collisions from 0 (perfectly ineleastic) to 1 (perfectly elastic)",
-    "collision_model": "Integer selecting the collision model being used. 0 for no collision. 1 for soft-sphere collisions",
-    "collision_time": "Amount of simulation time each collision will take to resolve",
-    "ib_coefficient_of_friction": "coefficient of friction used in IB collisions",
+    "coefficient_of_restitution": "Real number describing collision elasticity from 0 (perfectly inelastic) to 1 (perfectly elastic)",
+    "collision_model": "Integer selecting the collision model: 0 disables collisions, 1 enables soft-sphere collisions",
+    "collision_time": "Approximate simulation time used to resolve each collision",
+    "ib_coefficient_of_friction": "Coefficient of friction used in IB collisions",
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
"coefficient_of_restitution": "Real number describing the elasticity of collisions from 0 (perfectly ineleastic) to 1 (perfectly elastic)",
"collision_model": "Integer selecting the collision model being used. 0 for no collision. 1 for soft-sphere collisions",
"collision_time": "Amount of simulation time each collision will take to resolve",
"ib_coefficient_of_friction": "coefficient of friction used in IB collisions",
"coefficient_of_restitution": "Real number describing collision elasticity from 0 (perfectly inelastic) to 1 (perfectly elastic)",
"collision_model": "Integer selecting the collision model: 0 disables collisions, 1 enables soft-sphere collisions",
"collision_time": "Approximate simulation time used to resolve each collision",
"ib_coefficient_of_friction": "Coefficient of friction used in IB collisions",

@codecov
Copy link
Copy Markdown

codecov bot commented Apr 6, 2026

Codecov Report

❌ Patch coverage is 42.20183% with 126 lines in your changes missing coverage. Please review.
✅ Project coverage is 64.76%. Comparing base (23bbd4f) to head (caa2274).

Files with missing lines Patch % Lines
src/simulation/m_collisions.fpp 36.80% 84 Missing and 7 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 71.42% 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.88%   64.76%   -0.12%     
==========================================
  Files          70       71       +1     
  Lines       18251    18442     +191     
  Branches     1508     1524      +16     
==========================================
+ Hits        11842    11944     +102     
- Misses       5446     5531      +85     
- Partials      963      967       +4     

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

MIBM Restart Data Support

1 participant