[platelett]

Fixes #1994.

Two-way-coupled immersed boundaries (moving_ibm = 2) get a wrong moment of inertia:

• airfoil, STL, cuboid and cylinder bodies get NaN, and the run aborts;
• rectangles get twice the correct value;
• a slowly spinning 3D cuboid, cylinder or STL body diverges in its first step.

This PR computes the moment after the body's cells are marked, fixes the two formula errors, and fixes the zero-rotation check. Circles, spheres, and prescribed-motion or static IBs are unchanged. No goldens change.

Changes (all in src/simulation/m_ibm.fpp)

1. Compute it after the cells are marked.
   • s_ibm_setup called s_compute_moment_of_inertia in its first loop, before s_apply_ib_patches fills ib_markers. Shapes without a closed form sum r² over the body's marked cells. At that point there are none, so they got 0/0 = NaN.
   • In 2D this is the only time the moment is computed. In 3D, s_propagate_immersed_boundaries multiplies by the stored value before recomputing it.
   • The call now runs in its own loop, after the markers and the centroid offsets. The rotation-matrix and prescribed-kinematics setup stay where they were, because the markers depend on them.
2. Drop cell_volume from the cell sum.
   • patch%mass is the body's total mass: it is used as force/mass, and the analytic branches use 0.4*mass*r².
   • Each of the count cells therefore carries mass/count, so the moment is mass*Σr²/count.
   • Dividing by count*cell_volume gave the wrong units and made the moment scale with the grid.
3. Rectangle: /12, not /6.
   • length_x and length_y are full side lengths: f_is_inside_cuboid tests |x| ≤ length_x/2, and the ellipse branch uses the same convention.
   • The polar moment about the centroid is m(lx² + ly²)/12.
4. No 0/0 on a rank holding no cell of the body. It returns the dummy value 1, the same guard and value as Decode periodic IB markers in the moment of inertia sum #1992 (the non-rotating case already returns 1).
5. Only a zero axis counts as "not rotating".
   • In 3D the time stepper passes the angular momentum as the axis. Its norm was compared with sgm_eps = 1e-16, a threshold with units, so a small, slowly spinning body got the dummy moment 1.
   • The next RK stage multiplied by that dummy: a cuboid with I = 6e-14 kg·m² at 1e-3 rad/s went to ω = 8.3e9 rad/s in one step, and the run stopped on ICFL.
   • Now only the exact zero vector returns the dummy value. The axis is scaled by its largest component before normalising, so a small one does not underflow.

Side effects

• Results change only for two-way-coupled bodies:
  • rectangles: the moment is halved to the correct value;
  • shapes that use the cell sum: the moment is finite, where the run used to abort.
• The setup call moved later in s_ibm_setup. Nothing between the old and new positions reads patch_ib%moment: its only readers are the time stepper and the output.

Scope, and relation to #1992

#1992 (closed, tracked in pr-to-do #2068) decodes periodic markers in this sum and adds the count > 0 guard. This PR uses the same guard. The two PRs still overlap in this routine on the cell_volume lines (declaration, computation, final division), and a rebase drops them.

#1992 alone does not fix point 1. The setup-time call would still find no marked cells, so every cell-sum body would get the dummy moment of 1. For the airfoil below that is 1 kg·m² instead of 5.7e-14, so the body could barely rotate.

Not changed here, as noted in #1994 and #2068:

• The sum covers only the local rank's cells, so a body split across ranks gets a rank-dependent moment.
• Fixing that needs the sum moved to host code with an allreduce, because the routine is a GPU_ROUTINE called inside device loops.
• The sum gives every marked cell equal mass, which matches a uniform-density body only on a uniform grid. On a stretched grid both this sum and the centroid offset (s_compute_centroid_offset) would need volume weighting. That fits with the host-side rewrite above, rather than in this PR.

Verification

I printed patch_ib(1)%moment at the end of s_ibm_setup with a temporary print, not part of the PR. Setup:

• GNU 13.4 + Open MPI, 1 rank, mass = 1e-6, moving_ibm = 2.
• Cases derived from examples/2D_mibm_cylinder_in_cross_flow and examples/2D_ibm_airfoil.

Case	Exact	master	this PR
2D rectangle 0.8 × 0.4 mm, 201² (analytic branch)	m(lx²+ly²)/12 = 6.667e-14	1.333e-13 (2×)	6.667e-14
Same rectangle, routed through the cell sum	6.667e-14	NaN, run aborts	6.664e-14 (−0.03%)
3D cuboid 0.8 × 0.4 × 0.6 mm, 81³, rotating about z (cell sum)	6.667e-14	NaN, run aborts (icfl = NaN)	6.653e-14 (−0.2%)
2D NACA 2415 airfoil, c = 1 mm, 201² (cell sum, 117 cells)	5.652e-14 (section polar moment, integrated separately)	NaN, run aborts (MPI_ABORT)	5.832e-14 (+3.2%, 117-cell discretization)

Test set	Result
All 61 IBM tests (./mfc.sh test -o IBM), final branch, GNU 13.3 + Open MPI 4.1.6	61 passed. The one two-way-coupled rectangle test (09FDDC81) has no torque, so its golden does not move
Same, no MPI (before the count > 0 guard, which only matters with MPI)	61 passed
All 61 IBM tests, NVHPC 25.9 --gpu acc, H100 (before the guard)	61 passed
Full regression suite, GNU 13.3 + Open MPI 4.1.6, CPU (on a branch that also includes our other fixes; before the guard)	725 passed, 0 failed (39 convergence cases skipped by default)

The axis check (change 5).

• Case: examples/3D_rotating_sphere with the sphere replaced by a 0.8 × 0.4 × 0.6 mm cuboid, moving_ibm = 2, mass = 1e-6, angular_vel(3) = 1e-3, still fluid, inviscid, 50³ grid, 1 rank.
• Without the fix, ω_z goes to 8.3e9 rad/s after one step, and the run stops on ICFL.
• With the fix, ω_z = 1.000e-3 for all 20 steps. The torque stays below 2e-17.

The count > 0 guard on 4 ranks.

• Case: the airfoil above, moved to (1, 1) mm so that it lies on one rank, with a temporary per-rank print.
• Without the guard, the three ranks that hold none of the body carry moment = NaN. With it, they carry 1.
• The body's own rank has 5.832e-14 either way, the same as on 1 rank.
• The flow solution is byte-identical over 20 steps.

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 Diablo (Tsinghua University) participating in the SC26 Student Cluster Competition.

[github-actions]

Lines of Code

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

Directory	Lines	Diff
simulation	28569	+6
total	47638	+6

[github-actions]

Claude Code Review

Head SHA: c0a2c4e

Files changed:

• 1
• src/simulation/m_ibm.fpp

Findings:

• medium src/simulation/m_ibm.fpp: The moment-of-inertia sum in s_compute_moment_of_inertia runs over ib_markers%sf(0:m,0:n,0:p), which only holds the current rank's cells. count and moment are therefore per-rank partial values. A body that spans several ranks gets a different moment on each rank, scaled as sum_local(r^2)*mass/count_local. A rank that holds none of the body gets the new dummy 1._wp. The new count > 0 guard avoids the 0/0, but it does not make the value consistent across ranks. The new call site, after s_apply_ib_patches, is the first point where the markers are populated, so this now matters. No s_mpi_allreduce of count and the moment sum appears in the hunk. If the torque or angular acceleration is applied with this value on each rank, a multi-rank run would differ from a single-rank run. Consider reducing the sum of r^2 and count across ranks before dividing, as is done for n_need just above. Confidence is moderate, because I did not trace where patch%moment is consumed. The same limitation existed before, but the old call site ran before the markers were set.
• medium src/simulation/m_ibm.fpp: The rectangle moment changes from /6 to /12, and the cell-volume division is removed. Both are physically correct: m(a^2+b^2)/12 is the right formula, and the old 1/(count*cell_volume) had inconsistent units. They do change the results of every moving-IB case that uses a rectangle or a numerically summed shape. The PR needs the verification and golden-file updates that CLAUDE.md requires for changes to CFD results. Confidence is high that the change is correct, and I could not check from the diff whether the goldens were regenerated.

[codecov]

Codecov Report

❌ Patch coverage is 75.00000% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.98%. Comparing base (ee99151) to head (c0a2c4e).

Files with missing lines	Patch %	Lines
src/simulation/m_ibm.fpp	75.00%	1 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2109      +/-   ##
==========================================
+ Coverage   61.95%   61.98%   +0.02%     
==========================================
  Files          86       86              
  Lines       22906    22908       +2     
  Branches     3378     3378              
==========================================
+ Hits        14192    14200       +8     
+ Misses       6206     6200       -6     
  Partials     2508     2508              

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

[sbryngelson]

claude code review isn't happy with this one