Output Immersed Boundary load and state data (#1302)

mrvandenboom · claude · sbryngelson · web-flow · commit 9b7a503784f0 · 2026-03-16T09:11:08.000-04:00
Co-authored-by: Claude Opus 4.6 (1M context) &lt;noreply@anthropic.com&gt;
Co-authored-by: Spencer Bryngelson &lt;sbryngelson@gmail.com&gt;
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -615,6 +615,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_state_wrt`          | Logical | Write IB state and loads to a datafile at each time step |
 | `pi_inf_wrt`            | Logical | Add the liquid stiffness function to the database |
 | `pres_inf_wrt`          | Logical | Add the liquid stiffness to the formatted database	 |
 | `c_wrt`                 | Logical | Add the sound speed to the database	 |
@@ -675,6 +676,8 @@ 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` activates the output of data specified by patch_ib(i)%force(:) (and torque, vel, angular_vel, angles, [x,y,z]_centroid) into a single binary datafile for all IBs at all timesteps. During post_processing, this file is converted into separate time histories for each IB.
+
 - `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.
 It is not supported when `precision = 1` and `format = 1`.
diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90
@@ -164,7 +164,7 @@ impure subroutine s_read_ib_data_files(file_loc_base, t_step)
         if (parallel_io) then
             write (file_loc, '(A)') trim(file_loc_base)//'ib.dat'
         else
-            write (file_loc, '(A)') trim(file_loc_base)//'/ib.dat'
+            write (file_loc, '(A)') trim(file_loc_base)//'/ib_data.dat'
         end if
         inquire (FILE=trim(file_loc), EXIST=file_exist)
 
diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp
@@ -30,6 +30,7 @@ module m_data_output
  s_write_variable_to_formatted_database_file, &
  s_write_lag_bubbles_results_to_text, &
  s_write_lag_bubbles_to_formatted_database_file, &
+ s_write_ib_state_files, &
  s_write_intf_data_file, &
  s_write_energy_data_file, &
  s_close_formatted_database_file, &
@@ -1501,6 +1502,61 @@ contains
 
     end subroutine s_write_lag_variable_to_formatted_database_file
 
+    impure subroutine s_write_ib_state_files()
+
+        character(len=len_trim(case_dir) + 4*name_len) :: in_file, out_file, file_loc
+        integer :: iu_in, ios, i, rec_id
+        integer, allocatable, dimension(:) :: iu_out
+        real(wp) :: rec_time
+        real(wp), dimension(3) :: rec_force, rec_torque
+        real(wp), dimension(3) :: rec_vel, rec_angular_vel
+        real(wp), dimension(3) :: rec_angles, rec_centroid
+
+        file_loc = trim(case_dir)//'/D'
+
+        in_file = trim(file_loc)//'/ib_state.dat'
+        open (newunit=iu_in, file=trim(in_file), form='unformatted', access='stream', &
+              status='old', action='read', iostat=ios)
+        if (ios /= 0) then
+            call s_mpi_abort('Cannot open IB state input file: '//trim(in_file))
+        end if
+
+        allocate (iu_out(num_ibs))
+        do i = 1, num_ibs
+            write (out_file, '(A,I0,A)') trim(file_loc)//'/ib_', i, '.txt'
+            open (newunit=iu_out(i), file=trim(out_file), form='formatted', status='replace', action='write', iostat=ios)
+            if (ios /= 0) then
+                call s_mpi_abort('Cannot open IB state output file: '//trim(out_file))
+            end if
+            write (iu_out(i), '(A)') &
+                'mytime fx fy fz Tau_x Tau_y Tau_z vx vy vz omega_x omega_y omega_z angle_x angle_y angle_z x_c y_c z_c'
+        end do
+
+        do
+            read (iu_in, iostat=ios) rec_time, rec_id, &
+                rec_force, rec_torque, rec_vel, rec_angular_vel, rec_angles, &
+                rec_centroid(1), rec_centroid(2), rec_centroid(3)
+            if (ios /= 0) exit
+
+            if (rec_id >= 1 .and. rec_id <= num_ibs) then
+                write (iu_out(rec_id), '(19(ES24.16E3,1X))') rec_time, &
+                    rec_force(1), rec_force(2), rec_force(3), &
+                    rec_torque(1), rec_torque(2), rec_torque(3), &
+                    rec_vel(1), rec_vel(2), rec_vel(3), &
+                    rec_angular_vel(1), rec_angular_vel(2), rec_angular_vel(3), &
+                    rec_angles(1), rec_angles(2), rec_angles(3), &
+                    rec_centroid(1), rec_centroid(2), rec_centroid(3)
+            end if
+        end do
+
+        close (iu_in)
+        do i = 1, num_ibs
+            close (iu_out(i))
+        end do
+        deallocate (iu_out)
+
+    end subroutine s_write_ib_state_files
+
     !> @brief Extract the volume-fraction interface contour from primitive fields and write the coordinates to the interface data file.
     impure subroutine s_write_intf_data_file(q_prim_vf)
 
diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp
@@ -261,6 +261,7 @@ module m_global_parameters
     logical :: schlieren_wrt
     logical :: cf_wrt
     logical :: ib
+    logical :: ib_state_wrt
     logical :: chem_wrt_Y(1:num_species)
     logical :: chem_wrt_T
     logical :: lag_header
@@ -494,6 +495,7 @@ contains
         sim_data = .false.
         cf_wrt = .false.
         ib = .false.
+        ib_state_wrt = .false.
         lag_txt_wrt = .false.
         lag_header = .true.
         lag_db_wrt = .false.
diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp
@@ -101,7 +101,7 @@ contains
             & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt',          &
             & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange',        &
             & 'output_partial_domain', 'relativity', 'cont_damage', 'bc_io',   &
-            & 'down_sample','fft_wrt', 'hyper_cleaning' ]
+            & 'down_sample','fft_wrt', 'hyper_cleaning', 'ib_state_wrt']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp
@@ -116,7 +116,7 @@ contains
             lag_rad_wrt, lag_rvel_wrt, lag_r0_wrt, lag_rmax_wrt, &
             lag_rmin_wrt, lag_dphidt_wrt, lag_pres_wrt, lag_mv_wrt, &
             lag_mg_wrt, lag_betaT_wrt, lag_betaC_wrt, &
-            alpha_rho_e_wrt
+            alpha_rho_e_wrt, ib_state_wrt
 
         ! Inquiring the status of the post_process.inp file
         file_loc = 'post_process.inp'
diff --git a/src/post_process/p_main.fpp b/src/post_process/p_main.fpp
@@ -93,6 +93,10 @@ program p_main
     end do
     ! END: Time-Marching Loop
 
+    if (proc_rank == 0 .and. ib_state_wrt) then
+        call s_write_ib_state_files()
+    end if
+
     close (11)
 
     call s_finalize_modules()
diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
@@ -40,6 +40,8 @@ contains
 
         call s_check_inputs_time_stepping
 
+        @:PROHIBIT(ib_state_wrt .and. .not. ib, "ib_state_wrt requires ib to be enabled")
+
     end subroutine s_check_inputs
 
     !> Checks constraints on compiler options
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
@@ -37,17 +37,22 @@ module m_data_output
               s_open_run_time_information_file, &
               s_open_com_files, &
               s_open_probe_files, &
+              s_open_ib_state_file, &
               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_close_ib_state_file, &
+              s_finalize_data_output_module
+
+    integer :: ib_state_unit = -1 !< I/O unit for IB state binary file
 
     real(wp), allocatable, dimension(:, :, :) :: icfl_sf  !< ICFL stability criterion
     real(wp), allocatable, dimension(:, :, :) :: vcfl_sf  !< VCFL stability criterion
@@ -254,6 +259,20 @@ contains
 
     end subroutine s_open_probe_files
 
+    impure subroutine s_open_ib_state_file
+        character(len=path_len + 2*name_len) :: file_loc
+        integer :: ios
+
+        write (file_loc, '(A)') 'ib_state.dat'
+        file_loc = trim(case_dir)//'/D/'//trim(file_loc)
+        open (newunit=ib_state_unit, file=trim(file_loc), &
+              form='unformatted', &
+              access='stream', &
+              status='replace', &
+              iostat=ios)
+        if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: '//trim(file_loc))
+    end subroutine s_open_ib_state_file
+
     !>  The goal of the procedure is to output to the run-time
         !!      information file the stability criteria extrema in the
         !!      entire computational domain and at the given time-step.
@@ -1071,7 +1090,7 @@ contains
         write (t_step_dir, '(A,I0,A,I0)') trim(case_dir)//'/p_all'
         write (t_step_dir, '(a,i0,a,i0)') trim(case_dir)//'/p_all/p', &
             proc_rank, '/', time_step
-        write (file_path, '(A,I0,A)') trim(t_step_dir)//'/ib.dat'
+        write (file_path, '(A,I0,A)') trim(t_step_dir)//'/ib_data.dat'
 
         open (2, FILE=trim(file_path), &
               FORM='unformatted', &
@@ -1136,7 +1155,26 @@ contains
             call s_write_serial_ib_data(time_step)
         end if
 
-    end subroutine
+    end subroutine s_write_ib_data_file
+
+    !> @brief Writes IB state records to D/ib_state.dat. Must be called only on rank 0.
+    impure subroutine s_write_ib_state_file()
+
+        integer :: i
+
+        do i = 1, num_ibs
+            write (ib_state_unit) mytime, i, &
+                patch_ib(i)%force, &
+                patch_ib(i)%torque, &
+                patch_ib(i)%vel, &
+                patch_ib(i)%angular_vel, &
+                patch_ib(i)%angles, &
+                patch_ib(i)%x_centroid, &
+                patch_ib(i)%y_centroid, &
+                patch_ib(i)%z_centroid
+        end do
+
+    end subroutine s_write_ib_state_file
 
     !>  This writes a formatted data file where the root processor
     !!      can write out the CoM information
@@ -1907,6 +1945,12 @@ contains
 
     end subroutine s_close_probe_files
 
+    impure subroutine s_close_ib_state_file
+
+        close (ib_state_unit)
+
+    end subroutine s_close_ib_state_file
+
     !>  The computation of parameters, the allocation of memory,
         !!      the association of pointers and/or the execution of any
         !!      other procedures that are necessary to setup the module.
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
@@ -396,6 +396,7 @@ module m_global_parameters
     !> @{
     logical :: ib
     integer :: num_ibs
+    logical :: ib_state_wrt
 
     type(ib_patch_parameters), dimension(num_patches_max) :: patch_ib
     type(vec3_dt), allocatable, dimension(:) :: airfoil_grid_u, airfoil_grid_l
@@ -715,6 +716,7 @@ contains
         ! Immersed Boundaries
         ib = .false.
         num_ibs = dflt_int
+        ib_state_wrt = .false.
 
         ! Bubble modeling
         bubbles_euler = .false.
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
@@ -114,7 +114,7 @@ contains
             & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'surface_tension',       &
             & 'shear_stress', 'bulk_stress', 'bubbles_lagrange',                &
             & 'hyperelasticity', 'down_sample', 'int_comp','fft_wrt', &
-            & 'hyper_cleaning' ]
+            & 'hyper_cleaning', 'ib_state_wrt']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
@@ -144,6 +144,7 @@ contains
             x_domain, y_domain, z_domain, &
             hypoelasticity, &
             ib, num_ibs, patch_ib, &
+            ib_state_wrt, &
             fluid_pp, bub_pp, probe_wrt, prim_vars_wrt, &
             fd_order, probe, num_probes, t_step_old, &
             alt_soundspeed, mixture_err, weno_Re_flux, &
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
@@ -463,6 +463,11 @@ contains
             call s_open_run_time_information_file()
         end if
 
+        ! Opening and writing the header of the ib state data file
+        if (proc_rank == 0 .and. ib_state_wrt) then
+            call s_open_ib_state_file()
+        end if
+
         if (cfl_dt) then
             @:ALLOCATE(max_dt(0:m, 0:n, 0:p))
         end if
@@ -630,6 +635,14 @@ contains
                 ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms
                 if (moving_immersed_boundary_flag) then
                     call s_propagate_immersed_boundaries(s)
+                    ! compute ib forces for fixed immersed boundaries if requested for output
+                else if (ib_state_wrt .and. s == nstage) then
+                    call s_compute_ib_forces(q_prim_vf, fluid_pp)
+                end if
+
+                ! Write IB state to file if requested and at the RK final stage
+                if (proc_rank == 0 .and. ib_state_wrt .and. s == nstage) then
+                    call s_write_ib_state_file()
                 end if
 
                 ! update the ghost fluid properties point values based on IB state
@@ -1056,6 +1069,11 @@ contains
             call s_close_run_time_information_file()
         end if
 
+        ! Writing the footer of and closing the IB data file
+        if (proc_rank == 0 .and. ib_state_wrt) then
+            call s_close_ib_state_file()
+        end if
+
     end subroutine s_finalize_time_steppers_module
 
 end module m_time_steppers
diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py
@@ -558,9 +558,12 @@ def check_ibm(self):
         n = self.get("n", 0)
         num_ibs = self.get("num_ibs", 0)
 
+        ib_state_wrt = self.get("ib_state_wrt", "F") == "T"
+
         self.prohibit(ib and n <= 0, "Immersed Boundaries do not work in 1D (requires n > 0)")
         self.prohibit(ib and (num_ibs <= 0 or num_ibs > 1000), "num_ibs must be between 1 and num_patches_max (1000)")
         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")
 
     def check_stiffened_eos(self):
         """Checks constraints on stiffened equation of state fluids parameters"""
diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py
@@ -215,6 +215,7 @@
     "prim_vars_wrt": "Write primitive variables",
     "cons_vars_wrt": "Write conservative variables",
     "run_time_info": "Print runtime info",
+    "ib_state_wrt": "Write IB state and load data",
     # Misc
     "case_dir": "Case directory path",
     "cantera_file": "Cantera mechanism file",
@@ -905,7 +906,7 @@ def _load():
     _r("precision", INT, {"output"})
     _r("format", INT, {"output"})
     _r("schlieren_alpha", REAL, {"output"})
-    for n in ["parallel_io", "file_per_process", "run_time_info", "prim_vars_wrt", "cons_vars_wrt", "fft_wrt"]:
+    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 [
         "schlieren_wrt",
