Record immersed-boundary force history to one shared file

docs/documentation/case.md

@@ -741,6 +741,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	 |
@@ -808,7 +809,31 @@ 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 the force, torque and kinematics of every immersed boundary in a single shared text file, `D/ib_forces.dat`, described below. `ib_force_stride` writes only every N-th step, for runs long enough that the history itself becomes large.
+
+#### Immersed-boundary force history {#sec-ib-force-history}
+
+`D/ib_forces.dat` holds one fixed-width record per body per written step. Its twenty columns are
+
+| Columns | Quantity |
+| ---:    | :---     |
+| 1       | body id (the global `patch_ib` index) |
+| 2       | time |
+| 3–5     | force, x/y/z |
+| 6–8     | torque, x/y/z |
+| 9–11    | velocity, x/y/z |
+| 12–14   | angular velocity, x/y/z |
+| 15–17   | angles about x/y/z |
+| 18–20   | centroid, x/y/z |
+
+The file carries no header line, because every record sits at a computed byte offset and a header would shift them all. Each record is exactly 353 bytes including its newline (`I10` followed by nineteen `1X,ES17.9E3` fields), so the whole file loads with `numpy.loadtxt` and a single body or step can be read without scanning it:
+
+```
+row    = t_step / ib_force_stride
+offset = (row * num_ibs + ib_id - 1) * 353
+```
+
+Rows are written in global body-id order, so the file is byte-identical however the domain is decomposed, and no merge step is needed after a parallel run.
 
 - `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_history, &
+        & s_close_ib_force_history, 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,18 @@ module m_data_output
 
     type(scalar_field), allocatable, dimension(:) :: q_cons_temp_ds
 
+    !> One fixed-width text record per immersed body per recorded step, in D/ib_forces.dat.
+    !!
+    !! IB_REC_FMT is fixed width by construction: every ES descriptor right-justifies in its field,
+    !! including for negatives, three-digit exponents, NaN and Inf, so a record is always
+    !! IB_REC_BODY characters. That is what lets a rank compute a byte offset for (step, body) and
+    !! write there directly, giving one shared text file with no gather and no per-rank shards.
+    !! The two must be edited together: widening the format without IB_REC_LEN shears the file.
+    character(len=*), parameter :: IB_REC_FMT = '(I10,19(1X,ES17.9E3))'
+    integer, parameter          :: IB_REC_BODY = 10 + 19*18      !< characters the format emits
+    integer, parameter          :: IB_REC_LEN = IB_REC_BODY + 1  !< plus the newline
+    integer                     :: ib_hist_file = -1             !< held open for the run; -1 until first write
+
 contains
 
     !> Write data files. Dispatch subroutine that replaces procedure pointer.
@@ -1100,6 +1113,127 @@ contains
 
     end subroutine s_write_serial_ib_state
 
+    !> Record every immersed body's force, torque and kinematics for this step.
+    !!
+    !! One shared text file, D/ib_forces.dat, opened once for the run. Each rank writes only the
+    !! bodies it owns, at a byte offset computed from the step and the global body id, so the file
+    !! is byte-identical however the domain is decomposed and needs no merge step. Writing a file
+    !! per body instead costs an inquire, open and close per body per rank per step, which is
+    !! 3e5 filesystem metadata operations per step at 1000 ranks holding 100 bodies each.
+    !!
+    !! Layout: row r = t_step/ib_force_stride holds all num_gbl_ibs bodies in global id order, so
+    !! body g occupies bytes ((r*num_gbl_ibs) + g - 1)*IB_REC_LEN. The file carries no header line,
+    !! which would shift every offset after it; the columns are listed in docs/documentation/case.md.
+    impure subroutine s_write_ib_force_history(t_step)
+
+        integer, intent(in)                  :: t_step
+        character(LEN=IB_REC_LEN)            :: rec
+        character(LEN=path_len + 2*name_len) :: file_loc
+        real(wp)                             :: fields(19)
+        integer                              :: i, ib_idx, n_write, row
+
+#ifdef MFC_MPI
+        integer(kind=MPI_OFFSET_KIND) :: disp
+        integer                       :: ierr, status(MPI_STATUS_SIZE)
+#endif
+
+        if (.not. ib_state_wrt) return
+        if (mod(t_step, max(ib_force_stride, 1)) /= 0) return
+
+        n_write = num_local_ibs
+        if (num_procs == 1) n_write = num_ibs
+        row = t_step/max(ib_force_stride, 1)
+
+        $:GPU_UPDATE(host='[patch_ib(1:num_ibs)]')
+
+        call s_open_ib_force_history()
+
+        do i = 1, n_write
+            ib_idx = i
+            if (num_procs > 1) ib_idx = local_ib_patch_ids(i)
+
+            fields(1) = mytime
+            fields(2:4) = patch_ib(ib_idx)%force(1:3)
+            fields(5:7) = patch_ib(ib_idx)%torque(1:3)
+            fields(8:10) = patch_ib(ib_idx)%vel(1:3)
+            fields(11:13) = patch_ib(ib_idx)%angular_vel(1:3)
+            fields(14:16) = patch_ib(ib_idx)%angles(1:3)
+            fields(17) = patch_ib(ib_idx)%x_centroid
+            fields(18) = patch_ib(ib_idx)%y_centroid
+            fields(19) = patch_ib(ib_idx)%z_centroid
+
+            write (rec, IB_REC_FMT) patch_ib(ib_idx)%gbl_patch_id, fields
+            rec(IB_REC_LEN:IB_REC_LEN) = new_line('a')
+
+#ifdef MFC_MPI
+            disp = (int(row, MPI_OFFSET_KIND)*int(num_gbl_ibs, MPI_OFFSET_KIND) + int(patch_ib(ib_idx)%gbl_patch_id - 1, &
+                    & MPI_OFFSET_KIND))*int(IB_REC_LEN, MPI_OFFSET_KIND)
+            call MPI_FILE_WRITE_AT(ib_hist_file, disp, rec, IB_REC_LEN, MPI_CHARACTER, status, ierr)
+#else
+            write (ib_hist_file, rec=row*num_gbl_ibs + patch_ib(ib_idx)%gbl_patch_id) rec
+#endif
+        end do
+
+    end subroutine s_write_ib_force_history
+
+    !> Open the shared history file. Done once for the run.
+    impure subroutine s_open_ib_force_history
+
+        character(LEN=path_len + 2*name_len) :: file_loc
+        character(LEN=IB_REC_LEN)            :: probe
+        integer                              :: i
+
+#ifdef MFC_MPI
+        integer :: ierr
+        logical :: file_exist
+#endif
+
+        if (ib_hist_file /= -1) return
+
+        ! Every offset below assumes the format emits exactly IB_REC_BODY characters. Measure it once
+        ! rather than trusting that the format and the constant were edited together: a format one
+        ! character wider would shear every record past the first without any other symptom.
+        write (probe, IB_REC_FMT) 0, [(0._wp, i=1, 19)]
+        @:PROHIBIT(len_trim(probe) /= IB_REC_BODY, &
+                   & "IB force record width disagrees with IB_REC_BODY;  IB_REC_FMT and IB_REC_BODY must be changed together")
+
+        file_loc = trim(case_dir) // '/D/ib_forces.dat'
+#ifdef MFC_MPI
+        ! MPI_MODE_CREATE does not truncate, so a shorter run following a longer one in the same
+        ! directory would keep the old tail past its last record. Delete first, as the ib_state
+        ! writer does, then barrier so no rank opens before the delete lands.
+        inquire (FILE=trim(file_loc), EXIST=file_exist)
+        if (file_exist .and. proc_rank == 0) call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr)
+
+        ! Collective: every rank opens, including one holding no body this step.
+        call s_mpi_barrier()
+        call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ib_hist_file, ierr)
+#else
+        ! Unformatted: the record is already a formatted string, so this writes its bytes verbatim and
+        ! produces the same file the MPI branch does. A formatted direct-access write would need a
+        ! format and would pad rather than emit the string as-is.
+        open (newunit=ib_hist_file, file=trim(file_loc), form='unformatted', access='direct', recl=IB_REC_LEN, status='replace')
+#endif
+
+    end subroutine s_open_ib_force_history
+
+    !> Close the history file. Nothing is buffered, so there is nothing to flush first.
+    impure subroutine s_close_ib_force_history
+
+#ifdef MFC_MPI
+        integer :: ierr
+#endif
+
+        if (ib_hist_file == -1) return
+#ifdef MFC_MPI
+        call MPI_FILE_CLOSE(ib_hist_file, ierr)
+#else
+        close (ib_hist_file)
+#endif
+        ib_hist_file = -1
+
+    end subroutine s_close_ib_force_history
+
     !> @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

@@ -514,6 +514,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

@@ -1121,6 +1121,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_close_ib_force_history()
+
         if (model_eqns == model_eqns_6eq) call s_report_pressure_relaxation()
 
         call s_finalize_time_steppers_module()

src/simulation/m_time_steppers.fpp

@@ -478,6 +478,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_history(t_step)
+
                 if (cfl_dt) then
                     if (mytime >= t_stop) return
                 else