From 6e4d90c7c5f2b76614846d66cf487a606315517d Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 11 Sep 2026 10:42:11 -0400 Subject: [PATCH 1/3] Record immersed-boundary forces and kinematics every time step `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_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 --- docs/documentation/case.md | 3 +- src/simulation/m_data_output.fpp | 84 +++++++++++++++++++++++++- src/simulation/m_global_parameters.fpp | 1 + src/simulation/m_start_up.fpp | 2 + src/simulation/m_time_steppers.fpp | 2 + toolchain/mfc/case_validator.py | 2 + toolchain/mfc/params/definitions.py | 2 + toolchain/mfc/params/descriptions.py | 1 + 8 files changed, 94 insertions(+), 3 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 0bb752d64..25a81cc3d 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -728,6 +728,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca | `alpha_wrt(i)` | Logical | Add the volume fraction of fluid $i$ to the database | | `gamma_wrt` | Logical | Add the specific heat ratio function to the database | | `heat_ratio_wrt` | Logical | Add the specific heat ratio to the database | +| `ib_force_stride` | Integer | Stride, in time steps, of the per-step immersed-boundary force record (default 1) | | `ib_state_wrt` | Logical | Parameter to handle writing IB state on saves and outputting the state as a point mesh to SILO files. | | `pi_inf_wrt` | Logical | Add the liquid stiffness function to the database | | `pres_inf_wrt` | Logical | Add the liquid stiffness to the formatted database | @@ -795,7 +796,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu - `probe_wrt` activates the output of state variables at coordinates specified by `probe(i)%[x;y,z]`. -- `ib_state_wrt` is used to trigger post-processing of the IB state to be written out as a point mesh in the SILO files. When no IBs are moving, it also triggers force and torque calculation so that those values may be written to the output state files. +- `ib_state_wrt` is used to trigger post-processing of the IB state to be written out as a point mesh in the SILO files. When no IBs are moving, it also triggers force and torque calculation so that those values may be written to the output state files. It also records one line per time step in `D/ib_forces.dat` for each immersed boundary (time step, time, force, torque, velocity, angular velocity, angles, centroid). Records are buffered and written in batches rather than opened per step, and `ib_force_stride` writes only every N-th step for runs long enough that the record itself becomes large. - `output_partial_domain` activates the output of part of the domain specified by `[x,y,z]_output%%beg` and `[x,y,z]_output%%end`. This is useful for large domains where only a portion of the domain is of interest. diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 57f7eaf9f..254cde5bc 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -26,8 +26,9 @@ module m_data_output private public :: s_initialize_data_output_module, s_open_run_time_information_file, s_open_com_files, s_open_probe_files, & & s_write_run_time_information, s_write_data_files, s_write_serial_data_files, s_write_parallel_data_files, & - & s_write_ib_data_file, s_write_com_files, s_write_probe_files, s_write_ib_state_file, s_close_run_time_information_file, & - & s_close_com_files, s_close_probe_files, s_finalize_data_output_module + & s_write_ib_data_file, s_write_com_files, s_write_probe_files, s_write_ib_state_file, s_write_ib_force_files, & + & s_flush_ib_force_files, s_close_run_time_information_file, s_close_com_files, s_close_probe_files, & + & s_finalize_data_output_module real(wp), public, allocatable, dimension(:,:) :: c_mass $:GPU_DECLARE(create='[c_mass]') @@ -42,6 +43,11 @@ module m_data_output type(scalar_field), allocatable, dimension(:) :: q_cons_temp_ds + !> 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 + integer :: ib_force_buf_n = 0 + contains !> Write data files. Dispatch subroutine that replaces procedure pointer. @@ -1100,6 +1106,80 @@ contains end subroutine s_write_serial_ib_state + !> Record the force, torque and kinematic state of every owned immersed boundary for this time step. + !! + !! Rows are accumulated in a rank-local buffer and flushed to D/ib_forces.dat in batches, because opening + !! a file 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 would issue 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: the old + !! owner's pending rows still reach the right file. `ib_force_stride` subsamples very long runs. + impure subroutine s_write_ib_force_files(t_step) + + integer, intent(in) :: t_step + integer :: i, ib_idx, n_write + + if (mod(t_step, max(ib_force_stride, 1)) /= 0) return + + n_write = num_local_ibs + if (num_procs == 1) n_write = num_ibs + if (n_write == 0) return ! ranks holding no body have nothing to record + + $: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) + if (ib_force_buf_n == ib_force_buf_len) call s_flush_ib_force_files() + ib_force_buf_n = ib_force_buf_n + 1 + ib_force_buf(1, ib_force_buf_n) = real(max(patch_ib(ib_idx)%gbl_patch_id, ib_idx), wp) + ib_force_buf(2, ib_force_buf_n) = real(t_step, wp) + ib_force_buf(3, ib_force_buf_n) = mytime + ib_force_buf(4:6,ib_force_buf_n) = patch_ib(ib_idx)%force(1:3) + ib_force_buf(7:9,ib_force_buf_n) = patch_ib(ib_idx)%torque(1:3) + ib_force_buf(10:12,ib_force_buf_n) = patch_ib(ib_idx)%vel(1:3) + ib_force_buf(13:15,ib_force_buf_n) = patch_ib(ib_idx)%angular_vel(1:3) + ib_force_buf(16:18,ib_force_buf_n) = patch_ib(ib_idx)%angles(1:3) + ib_force_buf(19, ib_force_buf_n) = patch_ib(ib_idx)%x_centroid + ib_force_buf(20, ib_force_buf_n) = patch_ib(ib_idx)%y_centroid + ib_force_buf(21, ib_force_buf_n) = patch_ib(ib_idx)%z_centroid + end do + + end subroutine s_write_ib_force_files + + !> Write every buffered immersed-boundary force record and empty the buffer + impure subroutine s_flush_ib_force_files() + + character(LEN=path_len + 2*name_len) :: file_loc + integer :: i, j, ib_id, file_unit + logical :: file_exist + + do i = 1, ib_force_buf_n + ib_id = nint(ib_force_buf(1, i)) + if (ib_id < 0) cycle ! already written as part of an earlier body's pass + + 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' + end if + + do j = i, ib_force_buf_n ! all rows for this body, in time order + if (nint(ib_force_buf(1, j)) /= ib_id) cycle + write (file_unit, '(I10,19ES18.10E3)') nint(ib_force_buf(2, j)), ib_force_buf(3:21,j) + ib_force_buf(1, j) = -1._wp + end do + + close (file_unit) + end do + + ib_force_buf_n = 0 + + end subroutine s_flush_ib_force_files + !> @brief Writes IB state records to restart_data/ib_state.dat. Must be called only on rank 0. impure subroutine s_write_ib_state_file(time_step) diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 447bb05c5..04aaf8e44 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -512,6 +512,7 @@ contains collision_time = dflt_real ib_coefficient_of_friction = dflt_real ib_state_wrt = .false. + ib_force_stride = 1 many_ib_patch_parallelism = .false. ! Bubble modeling (sim-specific) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 485643aff..608f95c4a 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1109,6 +1109,8 @@ contains !> Finalize and deallocate all simulation sub-modules in reverse initialization order impure subroutine s_finalize_modules + if (ib .and. ib_state_wrt) call s_flush_ib_force_files() ! write whatever is still buffered + if (model_eqns == model_eqns_6eq) call s_report_pressure_relaxation() call s_finalize_time_steppers_module() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 95192f7ac..a38018e62 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -477,6 +477,8 @@ contains call s_compute_derived_variables(t_step, q_cons_ts(1)%vf, q_prim_ts1, q_prim_ts2) end if + if (ib_state_wrt) call s_write_ib_force_files(t_step) + if (cfl_dt) then if (mytime >= t_stop) return else diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index b15cbffb9..664d7fb52 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -759,6 +759,8 @@ def check_ibm(self): ) 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 + self.prohibit(ib_force_stride < 1, "ib_force_stride must be >= 1") self.prohibit(many_ib_patch_parallelism and not ib, "many_ib_patch_parallelism requires ib to be enabled") for i in range(1, num_particle_clouds + 1): diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index ca888ab52..63d43694c 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -677,6 +677,7 @@ def _load(): # Output _r("precision", INT, {"output"}) _r("format", INT, {"output"}) + _r("ib_force_stride", INT, {"output", "ib"}) for n in ["parallel_io", "file_per_process", "run_time_info", "prim_vars_wrt", "cons_vars_wrt", "fft_wrt", "ib_state_wrt"]: _r(n, LOG, {"output"}) for n in [ @@ -1343,6 +1344,7 @@ def _decl(targets: set, *names: str) -> None: "prim_vars_wrt", "fd_order", "ib_state_wrt", + "ib_force_stride", "avg_state", "alt_soundspeed", "mixture_err", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index 7d1278e00..f280620d1 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -136,6 +136,7 @@ "cf_wrt": "Write color function field", # Immersed boundaries "ib": "Enable immersed boundary method", + "ib_force_stride": "Write the per-step immersed-boundary force record every N steps (default 1)", "num_ibs": "Number of immersed boundary patches", "num_stl_models": "Number of STL/OBJ model entries in the stl_models array", "num_particle_clouds": "Number of particle bed specifications to generate immersed boundary patches from", From 0d79d0edbc31baa9cf034ea4d68f866d9d2f56ac Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 11 Sep 2026 11:02:50 -0400 Subject: [PATCH 2/3] Separate the columns in the IB force record 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 --- src/simulation/m_data_output.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 254cde5bc..57547d7f9 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1169,7 +1169,7 @@ contains do j = i, ib_force_buf_n ! all rows for this body, in time order if (nint(ib_force_buf(1, j)) /= ib_id) cycle - write (file_unit, '(I10,19ES18.10E3)') nint(ib_force_buf(2, j)), ib_force_buf(3:21,j) + write (file_unit, '(I10,19(1X,ES17.9E3))') nint(ib_force_buf(2, j)), ib_force_buf(3:21,j) ib_force_buf(1, j) = -1._wp end do From d0912a21adb9a61d9015238a2693a7db921c86a5 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 11 Sep 2026 15:01:38 -0400 Subject: [PATCH 3/3] Flush buffered force records before ownership changes hands 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 --- src/simulation/m_data_output.fpp | 5 +++-- src/simulation/m_time_steppers.fpp | 3 +++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 57547d7f9..84130cda4 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1111,8 +1111,9 @@ contains !! Rows are accumulated in a rank-local buffer and flushed to D/ib_forces.dat in batches, because opening !! a file 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 would issue 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: the old - !! owner's pending rows still reach the right file. `ib_force_stride` subsamples very long runs. + !! row carries its own global body id, so a record always reaches the right file. Ownership changes are + !! handled by flushing before the handoff, so a rank that stops owning a body cannot hold stale records and + !! append them after the new owner's newer ones. `ib_force_stride` subsamples very long runs. impure subroutine s_write_ib_force_files(t_step) integer, intent(in) :: t_step diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index a38018e62..d6b7c7a57 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -601,6 +601,9 @@ contains if (ib) then if (moving_immersed_boundary_flag) then + ! Write out what is buffered before ownership can change: a rank that stops owning a body would + ! otherwise hold its records until shutdown and append them after the new owner's newer ones. + if (ib_state_wrt) call s_flush_ib_force_files() call s_wrap_periodic_ibs() ! wraps the positions of IBs to the local proc call s_handoff_ib_ownership() ! recomputes which ranks own which IBs and communicate to neighbors else if (ib_state_wrt) then