Force reshape

[joshgillis]

Description

Refactored IB force and torque storage to avoid reshaping during communication.

Forces and torques are now stored in ib_ft as (6, num_ibs), where entries 1:3 are for the force components, and 4:6 are for the torque components. Additional related arrays such as recv_forces/torques and recv_forces/torques_snap are now also combined to match the formatting, labeled as recv_ft and recv_ft_snap.

Additional changes were made to m_collisions to reflect these updates, as well as an additional GPU_UPDATE in m_time_steppers for multi-rank GPU runs.

Closes #1452 .

Type of change (delete unused ones)

• Refactor

Testing

Tested with test suite, and compared results from a non-mpi and mpi on GPU to ensure results matched.

Checklist

Check these like this [x] to indicate which of the below applies.

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

See the developer guide for full coding standards.

GPU changes (expand if you modified src/simulation/)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

AI code reviews

Reviews are not retriggered automatically. To request a review, comment on the PR:

• @claude full review — Claude full review (also triggers on PR open/reopen/ready)
• Or add label claude-full-review — Claude full review via label

[sbryngelson]

The core refactor is sound — I verified it end to end. Collapsing forces(num_ibs,3)/torques(num_ibs,3) into ib_ft(6,num_ibs) preserves the index mapping everywhere (forces(i,l)→ib_ft(l,i), torques(i,l)→ib_ft(l+3,i)), the MPI path is equivalent (the old send_ft was already (6,num_ibs) column-major, so packing ib_ft directly is bit-identical), the snap-buffer merge into recv_ft_snap(:,j) matches the old per-component update, and the alloc/dealloc and host/device data movement all line up. The added m_time_steppers host update before the periodic-wrap/ownership handoff looks correct for multi-rank GPU.

One issue: there's an unrelated, broken change bundled in — the raw_pairs declaration in m_collisions.fpp was transposed but its accesses weren't. Details inline. I'd treat that as a blocker; the rest is good to go once it's resolved.

[joshgillis]

Sounds good, definitely an accidental change. Went ahead and reverted it. Will submit for review once I verify my CPU results.

[danieljvickers]

This looks good. I see you were not able to delete one of the array assignments on the GPU, but it is very clear why you had to leave it. When my main review comments are about naming conventions, it is a good day. I do think the name should be considered being changed. But that is about it. Otherwise, this is pretty much how it should be added.

[danieljvickers]

This is a point of some disagreement in coding philosophy, but I very-much like when variables are named extremely explicity. Physicists/mathemeticians/engineers are notorious for naming variiables things like mu_0 instead of permeability_of_free_space. The former saves much space when it is used in several equations, but the latter is extremely explicity about what the value is, even to non-technical users.

it is clear that ib_ft is "immersed boundary force and torque", but it could just as easily be something like forces_and_torques without causing too much clutter. But a name like this makes the code read more like a book. No one needs to go up a layer of abstraction to see what is being passed into this subroutine and understand what the array has in it. It is obvious from the name. This is a good thought exercises to participate in when naming variables. How to maximize clarity without making your equations bloated. And if you have to choose one or the other, which is more important?

[joshgillis]

I see, I definitely fall into the pattern of saving space, but I also agree that sometimes gets messy when it comes to legibility. What about ib_forces_torques for the name? I wanted to preserve the "ib" part of the name just so it doesn't get confused with other force/torque arrays in m_collisions if they are added later on, unless if you don't see that being necessary.

Edit: Do you also want this to extend to recv_ft and recv_ft_snap?

[danieljvickers]

I don't know if this was intentional, but you did a good job of respecting fortran's column-major ordering for memory layout. Something I forgot to consider in my original implementation.

[joshgillis]

Funny, I was looking into array shaping and saw the column-major preference Fortran has and thought that is what you had been implying with "no particular reason for their original shape" in your issue. Glad I could make an improvement.

[danieljvickers]

This would be a good place to put a comment about the contents of this array. Like I said, I already think it would be best renamed. But it is not clear immediately to a reader how memory is written here. A comment briefly describing that it is f_x, f_y, f_x, tau_x, tau_y, tau_z contiguously would make it easier for anyone reading the code to understand how they are supposed to fill in the array. Although the array bounds help clarify this a little for you.

[codecov]

Codecov Report

❌ Patch coverage is 42.42424% with 19 lines in your changes missing coverage. Please review.
✅ Project coverage is 60.52%. Comparing base (be974fc) to head (8eccf24).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_ibm.fpp	40.90%	13 Missing ⚠️
src/simulation/m_collisions.fpp	45.45%	6 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1616      +/-   ##
==========================================
+ Coverage   60.49%   60.52%   +0.02%     
==========================================
  Files          83       83              
  Lines       19902    19893       -9     
  Branches     2951     2951              
==========================================
  Hits        12040    12040              
+ Misses       5868     5859       -9     
  Partials     1994     1994              

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

[danieljvickers]

@joshgillis I saw you had a single failure on Phoenix that looked spurious to me. I reran just that test. If it fails in the same manner, it is something in your branch. If it passes, ping me here and I will restart the rest of the tests

Edit: Rerun still has the failing. So this is a real failure somehow. But out of bounds errors are some of the easiest to track down. Phoenix uses NVHPC/24.5.