Record immersed-boundary forces and kinematics every time step

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<id>_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.

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,81 @@ 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<id>_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 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
+        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,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
+
+            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)
 

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)

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()

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
@@ -599,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

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):

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",

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",