Skip to content

Record immersed-boundary forces and kinematics every time step - #1846

Open
sbryngelson wants to merge 6 commits into
masterfrom
feat/ib-force-history
Open

Record immersed-boundary forces and kinematics every time step#1846
sbryngelson wants to merge 6 commits into
masterfrom
feat/ib-force-history

Conversation

@sbryngelson

Copy link
Copy Markdown
Member

Why

ib_state_wrt records the force, torque and kinematic state of each immersed boundary only at snapshot intervals, so the force history is sampled at the field-output cadence — typically every few hundred steps. That is too coarse to compare a transient load against an experiment or a reference computation, or to drive a reduced-order model from a force signal. Every force comparison in the verification work behind this PR needed a per-step record.

What

One record per time step in D/ib<id>_forces.dat: time step, time, force, torque, velocity, angular velocity, angles, centroid. Gated by the existing ib_state_wrt, so no new switch to enable it.

Scaling

The obvious implementation — open, append, close, per body, per step — is a metadata operation per step on a parallel filesystem and does not scale. A particle bed of a thousand bodies over a long run would issue on the order of a hundred million of them. Instead records accumulate in a small fixed rank-local buffer (1024 rows) and are flushed in batches, with a final flush at shutdown. Each buffered row carries its own global body id, so a body changing owner mid-run needs no special handling: the previous owner's pending rows still reach the right file.

ib_force_stride (default 1) subsamples runs long enough for the record itself to become large. At the default the record is about 350 bytes per step per body — 25 MB for a 50k-step case, which is nothing next to the field output, but a thousand-body bed would want a stride.

Testing

Exercised on Frontier across 2D and 3D immersed-boundary cases, static and moving, on 4, 8 and 64 ranks. The per-step forces were checked against an independent control-volume momentum balance and, for a static cylinder at Re 40, against the published drag coefficient. ./mfc.sh precheck passes.

https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG

`ib_state_wrt` writes the force, torque and kinematic state of each immersed boundary only
at snapshot intervals, so the force history is sampled at the field-output cadence. That is
far too coarse to compare a transient load against an experiment or a reference
computation, or to drive a reduced-order model from a force signal: typical cases here
write a field every few hundred steps.

Write one record per time step to `D/ib<id>_forces.dat` (time step, time, force, torque,
velocity, angular velocity, angles, centroid), under the existing `ib_state_wrt` flag.

Records are buffered per rank and flushed in batches rather than opened per body per step:
opening a file per body per step is a metadata operation per step on a parallel filesystem
and does not scale, and a particle bed of a thousand bodies would issue on the order of a
hundred million of them over a long run. Each buffered row carries its own global body id,
so a body changing owner mid-run needs no special handling, and `ib_force_stride`
subsamples runs long enough for the record itself to become large.

Claude-Session: https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG
Copilot AI lite review requested due to automatic review settings September 11, 2026 14:42

Copilot AI left a comment

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.

Warning

Copilot couldn't run its full agentic review because it didn't start before the timeout. Make sure your repository has a runner available, or add a copilot-code-review.yml file specifying one with the runs-on attribute. See the docs for more details.

Pull request overview

Adds per-time-step immersed-boundary force/kinematics logging (buffered + flushed in batches) so IB force histories are recorded at high temporal resolution, with optional subsampling via ib_force_stride.

Changes:

  • Introduces ib_force_stride parameter (docs + toolchain registration + validation).
  • Writes buffered per-step IB force/torque/kinematics rows and flushes at shutdown.
  • Hooks per-step recording into the time-step loop when ib_state_wrt is enabled.

Reviewed changes

Copilot reviewed 8 out of 8 changed files in this pull request and generated 6 comments.

Show a summary per file
File Description
toolchain/mfc/params/descriptions.py Adds description string for ib_force_stride.
toolchain/mfc/params/definitions.py Registers ib_force_stride in parameter definitions.
toolchain/mfc/case_validator.py Validates ib_force_stride value in IBM checks.
src/simulation/m_time_steppers.fpp Calls per-step IB force logging routine.
src/simulation/m_start_up.fpp Flushes buffered IB force records at shutdown.
src/simulation/m_global_parameters.fpp Sets default ib_force_stride = 1.
src/simulation/m_data_output.fpp Implements buffered IB force record writing + flush to D/ib<id>_forces.dat.
docs/documentation/case.md Documents ib_force_stride and new per-step IB force record behavior.

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

Comment on lines +1160 to +1167
write (file_loc, '(A,I0,A)') '/D/ib', ib_id, '_forces.dat'
file_loc = trim(case_dir) // trim(file_loc)
inquire (file=trim(file_loc), exist=file_exist)
if (file_exist) then
open (newunit=file_unit, file=trim(file_loc), form='formatted', status='old', position='append')
else
open (newunit=file_unit, file=trim(file_loc), form='formatted', status='new')
write (file_unit, '(A)') '# t_step time Fx Fy Fz Tx Ty Tz vx vy vz wx wy wz ax ay az xc yc zc'
Comment on lines +46 to +48
!> Buffered immersed-boundary force records: (id, t_step, time, force, torque, vel, angular_vel, angles, centroid)
integer, parameter :: ib_force_buf_len = 1024
real(wp), dimension(21, ib_force_buf_len) :: ib_force_buf
Comment on lines +1127 to +1131
$:GPU_UPDATE(host='[patch_ib(1:num_ibs)]')

do i = 1, n_write
ib_idx = i
if (num_procs > 1) ib_idx = local_ib_patch_ids(i)
Comment thread toolchain/mfc/case_validator.py Outdated
)
self.prohibit(not ib and num_ibs > 0, "num_ibs is set, but ib is not enabled")
self.prohibit(ib_state_wrt and not ib, "ib_state_wrt requires ib to be enabled")
ib_force_stride = self.get("ib_force_stride", 1) or 1
integer :: i, j, ib_id, file_unit
logical :: file_exist

do i = 1, ib_force_buf_n
write (file_unit, '(A)') '# t_step time Fx Fy Fz Tx Ty Tz vx vy vz wx wy wz ax ay az xc yc zc'
end if

do j = i, ib_force_buf_n ! all rows for this body, in time order
ES18.10E3 is exactly wide enough for a negative value with a three-digit exponent
(-1.2345678901E+003), so adjacent values were written with no space between them and the
file could not be read as whitespace-separated columns. Use an explicit 1X separator.

Claude-Session: https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG
@github-actions

Copy link
Copy Markdown

Claude Code Review

Head SHA: 0d79d0e

Files changed:

  • 8
  • docs/documentation/case.md
  • src/simulation/m_data_output.fpp
  • src/simulation/m_global_parameters.fpp
  • src/simulation/m_start_up.fpp
  • src/simulation/m_time_steppers.fpp
  • toolchain/mfc/case_validator.py
  • toolchain/mfc/params/definitions.py
  • toolchain/mfc/params/descriptions.py

Findings:

  • toolchain/mfc/case_validator.py (check_ibm): ib_force_stride = self.get("ib_force_stride", 1) or 1 silently coerces a user-supplied 0 into 1 before the very next line checks ib_force_stride < 1. Since 0 is falsy in Python, the or 1 fallback fires for the exact invalid value (0) this check exists to catch, so ib_force_stride = 0 passes validation instead of being prohibited. Negative values are still caught (they're truthy), so only 0 slips through silently. Use self.get("ib_force_stride", 1) without the or 1 fallback (or check for None explicitly) so 0 still hits the < 1 prohibit.
  • src/simulation/m_time_steppers.fpp: s_write_ib_force_files(t_step) is called unconditionally whenever ib_state_wrt is true (line ~480), but per the adjacent existing code (if (moving_immersed_boundary_flag) ... else if (ib_state_wrt) call s_compute_ib_forces(...), lines ~602-609) and the doc text this PR itself edits in case.md ("When no IBs are moving, it also triggers force and torque calculation..."), patch_ib(...)%force/%torque are only refreshed when IBs are not moving. When moving_immersed_boundary_flag is true, the new D/ib<id>_forces.dat writer will keep recording stale/unrefreshed force and torque values every step instead of the actual current-step forces, with no indication in the output that the data is invalid for that run configuration.

A rank that stops owning a body kept whatever it had buffered and wrote it at shutdown,
after the new owner's newer records, so the file came out non-monotonic in time: on a
64-rank flapping case the time jumped from 3.24 back to 0 partway through.

The claim in the original comment -- that a body changing owner needs no special handling
because each row carries its own body id -- is true for which file a record lands in, but
not for its order. Flush before the handoff so a rank cannot hold stale records, and say so
correctly in the comment.

Found on Frontier, 8 nodes, flapping plate whose centroid crosses a rank boundary during
the stroke. Built with --gpu mp.

Claude-Session: https://claude.ai/code/session_01HMJ7cycfo7kTFSFq5yhHLG
@sbryngelson

Copy link
Copy Markdown
Member Author

Two updates pushed after running this on a real case.

A correctness fix. My original note claimed that because each buffered row carries its own body id, a body changing owner mid-run needed no special handling. That is true for which file a record lands in, but not for its order: a rank that stops owning a body kept its buffered records and wrote them at shutdown, after the new owner's newer ones. On a 64-rank flapping case the recorded time jumped from 3.24 back to 0 partway through the file. Fixed by flushing before the ownership handoff, and the comment now says the right thing.

The column separator. ES18.10E3 is exactly wide enough for a negative value with a three-digit exponent, so adjacent columns ran together and the file could not be read as whitespace-separated. Now uses an explicit separator.

Both were found by using the output rather than by review, which is the argument for landing this with a case that exercises it. The branch builds clean on Frontier with --gpu mp; I have now made a habit of building each of these before pushing, after learning that precheck lints but does not compile.

@codecov

codecov Bot commented Sep 12, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 76.19048% with 10 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.29%. Comparing base (dc0aec1) to head (17854d0).

Files with missing lines Patch % Lines
src/simulation/m_data_output.fpp 78.94% 1 Missing and 7 partials ⚠️
src/simulation/m_start_up.fpp 0.00% 0 Missing and 1 partial ⚠️
src/simulation/m_time_steppers.fpp 50.00% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1846      +/-   ##
==========================================
+ Coverage   61.26%   61.29%   +0.02%     
==========================================
  Files          84       84              
  Lines       22330    22372      +42     
  Branches     3265     3276      +11     
==========================================
+ Hits        13680    13712      +32     
- Misses       6207     6209       +2     
- Partials     2443     2451       +8     

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

@github-actions

Copy link
Copy Markdown

Lines of Code

File Lines Diff
src/simulation/m_data_output.fpp 1422 +55
src/simulation/m_time_steppers.fpp 863 +2
src/simulation/m_global_parameters.fpp 772 +1
src/simulation/m_start_up.fpp 1259 +1
Directory Lines Diff
simulation 28145 +59
total 46399 +59

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.

2 participants