[ruru25254649]

Summary

s_compute_moment_of_inertia sets patch_ib%moment, which the two-way-coupled (moving_ibm = 2) update divides the angular momentum by. It has three independent errors.

1. It runs before the body's cells are marked, so non-analytic shapes get NaN.

s_ibm_setup calls it in the first loop, before s_apply_ib_patches fills ib_markers. Shapes without a closed form are:

• airfoil
• STL
• cuboid
• cylinder
• 3D airfoil and STL

For these the routine sums r² over the cells marked with the patch id. At that point it finds none and stores 0/0 = NaN.

• In 2D this is the only time the moment is computed, so a two-way-coupled airfoil or STL body carries a NaN moment for the whole run.
• In 3D, s_propagate_immersed_boundaries multiplies by the stored moment before it recomputes it, so a body started with a nonzero angular velocity goes NaN on the first stage.

The call now runs after the markers and the centre-of-mass offsets (s_compute_centroid_offset) are in place. The rotation-matrix and prescribed-kinematics setup it was grouped with stay where they were, because the markers depend on them.

2. The cell sum divides by count*cell_volume (#1994).

patch%mass is the body's total mass. It is used as force/mass and accel_bf*mass, and the analytic branches use it as 0.4*mass*r². Each of the count cells therefore carries mass/count, and the moment is mass*Σr²/count. The extra 1/cell_volume gives units of mass/length (mass/length² in 2D) and makes the moment scale with the grid.

3. The rectangle uses m(lx² + ly²)/6.

length_x and length_y are full side lengths:

• f_is_inside_cuboid tests |x| ≤ length_x/2.
• The bounding radius is 0.5*sqrt(lx² + ly²).
• The ellipse branch, m(lx² + ly²)/16, already uses that convention.

The polar moment of a rectangle about its centroid is m(lx² + ly²)/12, so the old value was twice too large.

Verification

I printed patch_ib(1)%moment at the end of s_ibm_setup (a temporary print, not part of the PR). The build was GNU 13.4 + Open MPI, run on 1 rank. All three cases use mass = 1e-6 and are derived from examples/2D_mibm_cylinder_in_cross_flow and examples/2D_ibm_airfoil, with moving_ibm = 2:

• 2D rectangle: 0.8 × 0.4 mm on a 201² grid
• 3D cuboid: 0.8 × 0.4 × 0.6 mm on an 81³ grid, rotating about z, angular_vel(3) = 1
• 2D NACA 2415 airfoil: c = 1 mm on a 201² grid

Case	Exact	master	this PR
2D rectangle (analytic branch)	m(lx²+ly²)/12 = 6.667e-14	1.333e-13 (2×)	6.667e-14
2D rectangle, routed through the cell-sum branch	6.667e-14	NaN; the run aborts	6.664e-14 (−0.03%)
3D cuboid (cell sum)	6.667e-14	NaN; the run aborts (icfl = NaN)	6.653e-14 (−0.2%)
2D airfoil (cell sum, 117 cells)	5.652e-14 (polar moment of the NACA 2415 section, integrated separately)	NaN; the run aborts (NaN(s) in timestep output, MPI_ABORT)	5.832e-14 (+3.2%, the 117-cell discretization)

After the fix, the cell sum agrees with the closed forms. Before it, the cell sum could not be evaluated at all.

Test set	Result
All 61 IBM tests (./mfc.sh test -o IBM), GNU 13.4 + Open MPI	61 passed, 0 failed. No golden changes are needed: the circle and sphere formulas are unchanged, and the one two-way-coupled rectangle test (09FDDC81) has no torque
Same, no MPI	61 passed, 0 failed
All 61 IBM tests, NVHPC 25.9 --gpu acc, H100	61 passed, 0 failed
Full regression suite, GNU 13.3 + Open MPI 4.1.6, CPU, master 961cde5 (on a branch merging all candidates)	725 passed, 0 failed (39 convergence cases skipped by default)

Out of scope: the cell sum still covers only the local rank's cells, and periodic markers are still compared undecoded. Both are noted in #1994.

Relation to #1992 (closed, tracked by pr-to-do #2068). That PR decodes periodic markers and replaces 0/0 on a rank holding none of the body with the dummy value 1. It does not fix point 1. With #1992 alone, the setup-time call still finds no marked cells, but now every body that needs the cell sum gets moment = 1 instead of NaN. For the airfoil above, that is 1 kg·m² instead of 5.7e-14, and the body effectively cannot rotate. The two changes are complementary. They overlap on one line, the final moment = ... assignment. If #1992 is reopened first, the rebase keeps its count > 0 guard and drops *cell_volume.

Side effects

• Results change only for two-way-coupled bodies (moving_ibm = 2):
  • rectangles: the moment is halved to the correct value;
  • airfoil, STL, cuboid, cylinder and 3D shapes: the moment is now finite, where before it was NaN and the run aborted.
• Circles and spheres, and every prescribed-motion or static IB, are unchanged.
• No goldens change.
• The setup call moved later in s_ibm_setup, after s_apply_ib_patches and s_compute_centroid_offset. Nothing between the old and new positions reads patch_ib%moment.

Contribution Policy

All contributions are expected to demonstrate:

• A clear understanding of the codebase
• Alignment with product direction
• Thoughtful reasoning behind changes
• Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

---

This PR is from Your Team Name participating in the SC26 Student Cluster Competition.

🤖 Generated with Claude Code

https://claude.ai/code/session_01WPhkY8p1EHhGY7jScizJmC

[github-actions]

Claude Code Review

Head SHA: cae2acb

Files changed:

• 1
• src/simulation/m_ibm.fpp

Findings:

• medium src/simulation/m_ibm.fpp: In s_compute_moment_of_inertia, the sum fallback (geometry not 2/3/6/8) now ends with moment = moment*patch%mass/count. count and the sum run only over this rank's cells (do i = 0, m over ib_markers%sf). Two cases follow.
  • A rank that does not hold any cell of the body has count = 0 and moment = 0, so 0/0 gives NaN. The new loop calls the routine for every moving IB on every rank, and the result is copied to host, so a NaN patch_ib(i)%moment can be present on those ranks. The old count*cell_volume denominator had the same hazard.
  • When a body spans several ranks, each rank divides its partial sum by its own partial count and applies the full mass. The result is a per-rank value, not the body's moment of inertia, and it differs between ranks.
  • Suggested fix: guard count == 0, and reduce the sum and count across ranks, for example in the caller with s_mpi_allreduce_*, before dividing. Confidence: fairly high for the count == 0 case, moderate for the multi-rank case, since I did not check how moment is used downstream.
• medium src/simulation/m_ibm.fpp: The moment of inertia for sum-based shapes is now computed once in the new loop. The axis passed there is patch_ib(i)%angular_vel. If the initial angular velocity is zero in 3D, the routine returns the dummy moment = 1._wp, and it is not recomputed until m_time_steppers line 911 does so. This matches the old behaviour. The time stepper only recomputes it when num_dims == 3, so the 2D sum path keeps this initial value, which is fine. The new code is order-correct: the loop runs after s_apply_ib_patches and the centroid offsets, and the host update follows it. Confidence: low, informational only.
• Positive check: the rectangle moment changes from /6 to /12, which is the correct m(a²+b²)/12. Dropping cell_volume matches the new mass-normalised average of r² over marked cells.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_ibm.fpp	1824	+6

Directory	Lines	Diff
simulation	28569	+6
total	47632	+6