IBM force calculation simplification (Dr Bala recommendations)

[mrvandenboom]

User description

(cherry picked from commit 2a034dc) (cherry picked from commit 1d8b109)

Description

This change takes advantage of the symmetric and asymmetric attributes of the tensors to simplify the calculation of force and torque on immersed boundaries.

Type of change

• Bug fix
• New feature
• [ X] Refactor
• Documentation
• Other: describe

Testing

ran ./mfc.sh lint
ran ./mfc.sh format
ran ./mfc.sh test -a on Hipergator (CPU). All tests passed.

---

CodeAnt-AI Description

Compute torque from combined pressure and viscous forces and reduce atomic updates

What Changed

• Torque is computed once from the final local force (pressure + viscous) instead of from many intermediate cross-product terms
• Viscous contributions are added directly to the force divergence rows before torque is computed, ensuring torque reflects combined pressure and viscous effects
• Global force and torque arrays are updated atomically only after each point's force and torque are fully assembled, reducing repeated atomic writes

Impact

✅ Correct torque from combined pressure and viscous contributions
✅ Fewer atomic updates during force/torque assembly
✅ Lower contention during per-point force assembly

💡 Usage Guide

Checking Your Pull Request

Every time you make a pull request, our system automatically looks through it. We check for security issues, mistakes in how you're setting up your infrastructure, and common code problems. We do this to make sure your changes are solid and won't cause any trouble later.

Talking to CodeAnt AI

Got a question or need a hand with something in your pull request? You can easily get in touch with CodeAnt AI right here. Just type the following in a comment on your pull request, and replace "Your question here" with whatever you want to ask:

@codeant-ai ask: Your question here

This lets you have a chat with CodeAnt AI about your pull request, making it easier to understand and improve your code.

Example

@codeant-ai ask: Can you suggest a safer alternative to storing this secret?

Preserve Org Learnings with CodeAnt

You can record team preferences so CodeAnt AI applies them in future reviews. Reply directly to the specific CodeAnt AI suggestion (in the same thread) and replace "Your feedback here" with your input:

@codeant-ai: Your feedback here

This helps CodeAnt AI learn and adapt to your team's coding style and standards.

Example

@codeant-ai: Do not flag unused imports.

Retrigger review

Ask CodeAnt AI to review the PR again, by typing:

@codeant-ai: review

Check Your Repository Health

To analyze the health of your code repository, visit our dashboard at https://app.codeant.ai. This tool helps you identify potential issues and areas for improvement in your codebase, ensuring your repository maintains high standards of code health.

[codeant-ai]

Nitpicks 🔍

🔒 No security issues identified

⚡ Recommended areas for review Torque correctness The PR changes torque assembly: it no longer computes derivatives of the cross-product of the viscous stress (previous approach) and instead computes the cross-product once from the final, combined local force contribution. This is mathematically valid (torque = r x F) but can produce subtle numerical differences compared to the old per-component cross-derivative method. Verify equivalence across test cases (especially near boundaries and with non-uniform viscosity) and check that s_compute_viscous_stress_tensor returns the tensor rows in the ordering assumed here. Atomic update / precision check Forces and torques are updated using GPU atomics. Ensure the target GPU atomic implementation supports the required precision (real(wp)) and that the $:GPU_ATOMIC pragmas correctly protect concurrent updates for both forces and torques. Also validate that local accumulators are private in the parallel region on all backends.

[Copilot]

Pull request overview

Refactors immersed-boundary (IBM) force/torque accumulation in simulation to simplify the force/torque calculation by leveraging tensor symmetry/asymmetry characteristics, per Dr. Bala’s recommendations.

Changes:

• Simplifies viscous-force contribution assembly by removing per-component cross-product derivatives.
• Computes torque via a single r × F cross product after all force contributions are accumulated.
• Reorders/clarifies some inline comments around viscous contribution and atomic updates.

[Copilot]

Inline comment is incorrect: this statement computes the z-derivative of the third row (index 3) of the viscous stress tensor, not the second row. Please update the comment to avoid confusion for future maintenance.

Suggested change

	viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the second-row of viscous stress tensor
	viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the third-row of viscous stress tensor

[Copilot]

This comment says only the force is updated atomically, but the loop performs atomic updates for both forces and torques. Consider updating the comment to reflect both arrays to prevent misleading documentation.

Suggested change

	! Update the force values atomically to prevent race conditions
	! Update the force and torque values atomically to prevent race conditions

[cubic-dev-ai]

No issues found across 1 file

Confidence score: 5/5

• Automated review surfaced no issues in the provided summaries.
• No files require special attention.

[coderabbitai]

Actionable comments posted: 1

🤖 Prompt for all review comments with AI agents

Verify each finding against the current code and only fix it if needed.

Inline comments:
In `@src/simulation/m_ibm.fpp`:
- Line 1130: Remove the trailing whitespace on the blank line flagged by CI (a
line that contains only spaces) in the m_ibm.fpp source so the formatter passes;
locate the empty/blank line in the file, delete the trailing spaces (convert it
to a truly empty line), run ./mfc.sh format to verify no formatting diffs
remain, and commit the change.

[codecov]

Codecov Report

❌ Patch coverage is 71.42857% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 44.95%. Comparing base (7c806be) to head (47afbfe).

Files with missing lines	Patch %	Lines
src/simulation/m_ibm.fpp	71.42%	2 Missing ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1234   +/-   ##
=======================================
  Coverage   44.95%   44.95%           
=======================================
  Files          70       70           
  Lines       20503    20488   -15     
  Branches     1946     1946           
=======================================
- Hits         9217     9211    -6     
+ Misses      10164    10155    -9     
  Partials     1122     1122           

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

[github-actions]

Claude Code Review

1 file changed: src/simulation/m_ibm.fpp

Summary

• Refactors s_compute_ib_forces in m_ibm.fpp to simplify the torque calculation from a per-direction cross-product accumulation to a single cross product applied to the fully-assembled force vector.
• The mathematical equivalence holds by linearity of the cross product: r × F_p + r × F_v = r × (F_p + F_v), so computing torque from the combined force after adding viscous contributions is equivalent to the old per-direction approach.
• Eliminates ~6 cross-product calls per IB grid point per time step in 3D viscous cases, plus removes the explicit loop over rows l — a meaningful GPU compute savings on dense IBM meshes.
• The new code narrows each viscous_stress_div update from a full 3×3 matrix write to a single-row write per direction, which is a cleaner representation of ∇·σ.
• PR states CPU tests passed on HiPerGator; no GPU test run is mentioned.

---

Findings

Issue 1 — Dead variables left in declaration and GPU private clause (m_ibm.fpp)

viscous_cross_1 and viscous_cross_2 are no longer assigned or read anywhere in the new code, but they remain in two places:

Line ~1053 (variable declaration):

real(wp), dimension(1:3, 1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2

Line ~1074 (GPU_PARALLEL_LOOP private clause):

$:GPU_PARALLEL_LOOP(private='[... viscous_cross_1, viscous_cross_2, ...]', ...)

These should be removed. On GPU builds the unused private allocations are harmless but waste per-thread stack space, and they create misleading code. The l loop variable can also be examined — it is still used in the atomic update loop below, so it stays.

Recommendation: Remove viscous_cross_1 and viscous_cross_2 from both the declaration and the private list.

---

Issue 2 — Incorrect comment on z-direction block (m_ibm.fpp)

The new z-direction line says "second-row" but indexes row 3:

viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the second-row of viscous stress tensor

This is a copy-paste error from the y-direction comment. It should read "third-row".

---

Issue 3 — GPU test coverage not confirmed

This PR modifies code in src/simulation/ inside a GPU parallel loop (GPU_PARALLEL_LOOP). The PR description notes only CPU tests on HiPerGator. For a change to GPU-offloaded IBM force/torque assembly, a GPU regression run (e.g., on Phoenix or HiPerGator GPU partition) would close the loop, since CPU and GPU path divergences are possible (e.g., private variable scoping, atomics behavior).

---

Improvement Opportunities

1. Remove the ! TODO :: This is really bad code comment if this refactor resolves the concern it referred to, or update it to reflect what remains suboptimal about the viscous stress stencil approach.

2. Unify the y-block comment style — the x-block says "first-row", y says "second-row", z (after fix) would say "third-row". Consider a single parametric comment pattern or a brief note above the block explaining the ∇·σ assembly pattern to make the index-to-direction mapping self-evident.

3. Consider extracting the viscous force divergence assembly into a helper subroutine (e.g., s_compute_viscous_force_divergence) to reduce the per-dimension repetition and make s_compute_ib_forces easier to read — the current TODO comment suggests this was already considered.

[qodo-code-review]

Review Summary by Qodo

Simplify IBM force and torque calculations with combined contributions

✨ Enhancement

Walkthroughs

Description

• Simplify force and torque calculation by computing torque once from combined pressure and viscous
  forces
• Remove redundant cross-product calculations for intermediate viscous stress terms
• Reduce atomic update contention by assembling complete force before updating global arrays
• Move torque computation after all force contributions are accumulated

Diagram

flowchart LR
  A["Pressure gradient\ncontribution"] --> B["Accumulate all\nforce components"]
  C["Viscous stress\nderivatives"] --> B
  B --> D["Compute torque\nonce from final force"]
  D --> E["Atomic update\nglobal arrays"]

Loading

File Changes

1. src/simulation/m_ibm.fpp ✨ Enhancement +11/-31

Refactor force and torque assembly logic

• Removed intermediate cross-product calculations for viscous stress tensor rows
• Consolidated viscous force contributions directly into local_force_contribution before torque
 computation
• Moved s_cross_product call to after all force components are accumulated, ensuring torque
 reflects combined pressure and viscous effects
• Simplified viscous stress divergence calculations by computing only the required row derivatives
 instead of full tensor cross-products
• Reduced atomic update operations by assembling complete force and torque before updating global
 arrays

src/simulation/m_ibm.fpp

---

[qodo-code-review]

Code Review by Qodo

🐞 Bugs (1) 📘 Rule violations (0) 📎 Requirement gaps (0)

1. Dead GPU privates left 🐞 Bug ➹ Performance

Description

viscous_cross_1/viscous_cross_2 are no longer used after the refactor but are still declared and
listed as GPU-private variables. On GPUs, this can increase register/private-memory pressure and
reduce kernel occupancy/performance for the immersed-boundary force kernel.

Code

src/simulation/m_ibm.fpp[R1063-1078]

+                            ! get the linear force components first
                            call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, j, k)
                            call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, j, k)
-                            viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dx) ! get the x derivative of the viscous stress tensor
-                            local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1, 1:3) ! add te x components of the derivative to the force
-                            do l = 1, 3
-                                ! take the cross products for the torque component
-                                call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
-                                call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
-                            end do
-
-                            viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dx) ! get the x derivative of the cross product
-                            local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(1, 1:3) ! apply the cross product derivative to the torque
+                            viscous_stress_div(1, 1:3) = (viscous_stress_div_2(1, 1:3) - viscous_stress_div_1(1, 1:3))/(2._wp*dx) ! get x derivative of the first-row of viscous stress tensor
+                            local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1, 1:3) ! add the x components of the divergence to the force

                            call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j - 1, k)
                            call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j + 1, k)
-                            viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dy)
-                            local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2, 1:3)
-                            do l = 1, 3
-                                call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
-                                call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
-                            end do
-
-                            viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dy)
-                            local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(2, 1:3)
+                            viscous_stress_div(2, 1:3) = (viscous_stress_div_2(2, 1:3) - viscous_stress_div_1(2, 1:3))/(2._wp*dy) ! get y derivative of the second-row of viscous stress tensor
+                            local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2, 1:3) ! add the y components of the divergence to the force

                            if (num_dims == 3) then
                                call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j, k - 1)
                                call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j, k + 1)
-                                viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dz)
-                                local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3, 1:3)
-                                do l = 1, 3
-                                    call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
-                                    call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
-                                end do
-                                viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dz)
-                                local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(3, 1:3)
+                                viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the second-row of viscous stress tensor
+                                local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3, 1:3) ! add the z components of the divergence to the force

Evidence

The PR removed the viscous cross-product torque path (no references to viscous_cross_* remain in
the viscous divergence computation), but the variables are still declared and explicitly privatized
in the GPU parallel loop, so they can still contribute to per-thread storage requirements.

src/simulation/m_ibm.fpp[995-1035]
src/simulation/m_ibm.fpp[1058-1091]

Agent prompt

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

### Issue description
`viscous_cross_1` and `viscous_cross_2` are no longer used after the torque simplification refactor, but they are still declared and included in the GPU loop's `private` list. This can increase per-thread register/private-memory usage in a hot kernel.

### Issue Context
The refactor removed the viscous cross-product derivative approach and now computes torque from the assembled force. The cross-product temporaries are therefore dead.

### Fix Focus Areas
- src/simulation/m_ibm.fpp[1000-1025]
- src/simulation/m_ibm.fpp[1058-1079]

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

---

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

[github-actions]

Claude Code Review

Head SHA: 47afbfe5ee673bd32847b2a257d2f8c8c9d40c4d

Files changed: 1

• src/simulation/m_ibm.fpp

Summary:

• Refactors IBM force/torque assembly in s_compute_ibm_forces (or equivalent) to collect all force contributions (pressure + viscous) first, then compute torque as one cross product at the end instead of accumulating partial torques.
• Removes 3×N redundant cross-product calls and associated FD loops over viscous stress stencil points, reducing GPU atomic contention.
• Narrows viscous stress FD to only the required row per direction (row 1 for x, row 2 for y, row 3 for z) instead of the full 3×3 FD.
• Mathematical equivalence holds because radial_vector is fixed at the marker point, so r × (f_p + f_v) = r × f_p + r × f_v.

---

Findings

1. Comment typo — z-direction row index (minor)
src/simulation/m_ibm.fpp, z-direction block (added line ~+26 in the diff):

viscous_stress_div(3, 1:3) = ... ! get z derivative of the second-row of viscous stress tensor

The comment says "second-row" but the code accesses index 3 (third row). Should read "third-row".

2. Potentially dead local arrays — viscous_cross_1 / viscous_cross_2
The only consumers of viscous_cross_1 and viscous_cross_2 were the removed loops:

! OLD (removed)
call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
...

If these arrays are no longer referenced anywhere else in the subroutine they should be removed from the local variable declarations and any associated @:ALLOCATE / @:DEALLOCATE pairs. Leaving them allocated wastes device memory and may confuse future readers.

3. Correctness — verified correct
The original code computed local_torque_contribution = r × f_pressure (cross product called before viscous contributions were added), then accumulated viscous torque terms as r × ∂σ_i/∂x_i via the FD loops. The new code accumulates all force contributions first, then calls s_cross_product(radial_vector, local_force_contribution, local_torque_contribution) once. Since radial_vector is a constant at each marker point, the two formulations are algebraically identical. No correctness regression.

4. "TODO :: This is really bad code" comment preserved
This pre-existing comment (not introduced by this PR) is in the viscous block. Not a blocker for this PR, but worth tracking for a follow-up cleanup.

---

Overall: Approve with the comment typo fix and a check that viscous_cross_1/viscous_cross_2 are cleaned up if unused. The simplification is sound and the reported test coverage (CPU full suite on HiPerGator) is appropriate for this refactor.

[coderabbitai]

📝 Walkthrough

Walkthrough

The pull request refactors the force and torque computation logic in the IBM simulation module. Previously, atomic updates occurred inside the per-IBPoint loop. The change consolidates this by accumulating viscous stress derivatives into local force and torque contribution variables, then performing a single cross-product calculation and atomic update block after all local contributions are computed. Comments are adjusted to reflect the updated computation sequence focusing on force divergence contributions rather than per-step operations.

🚥 Pre-merge checks | ✅ 3

✅ Passed checks (3 passed)

Check name	Status	Explanation
Title check	✅ Passed	The title clearly and specifically summarizes the main refactor: simplifying IBM force calculations based on recommendations, which aligns with the primary change in the changeset.
Description check	✅ Passed	The description is mostly complete with a clear summary of changes, motivation (using tensor symmetry/asymmetry), proper type selection (Refactor), and testing information showing lint, format, and test execution on CPU.
Docstring Coverage	✅ Passed	No functions found in the changed files to evaluate docstring coverage. Skipping docstring coverage check.

_{✏️ 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

• X
• Mastodon
• Reddit
• LinkedIn

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

[coderabbitai]

🧹 Nitpick comments (1)

src/simulation/m_ibm.fpp (1)

1077-1077: Clarify two inline comments to match behavior.

Line 1077 says “second-row” but this is the third-row derivative, and Line 1085 mentions only force while both force and torque are updated atomically.

✏️ Suggested comment-only diff

-                                viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the second-row of viscous stress tensor
+                                viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the third-row of viscous stress tensor
@@
-                        ! Update the force values atomically to prevent race conditions
+                        ! Update the force and torque values atomically to prevent race conditions

Also applies to: 1085-1085

---

ℹ️ Review info

Configuration used: Path: .coderabbit.yaml

Review profile: CHILL

Plan: Pro

📥 Commits

Reviewing files that changed from the base of the PR and between 7c806be and 47afbfe.

📒 Files selected for processing (1)

• src/simulation/m_ibm.fpp