Immersed Boundary force output

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_force_wrt`          | Logical | Write IB force & torque 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_force_wrt` activates the output of data specified by patch_ib(i)%force(:) and patch_ib(i)%torque(:) into a single binary datafile for all IBs at all timesteps. During post_processing, this file is converted into separate force/torque 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`.

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_force_files, &
  s_write_intf_data_file, &
  s_write_energy_data_file, &
  s_close_formatted_database_file, &
@@ -1499,6 +1500,58 @@ contains
 
     end subroutine s_write_lag_variable_to_formatted_database_file
 
+    impure subroutine s_write_ib_force_files()
+
+        character(len=len_trim(case_dir) + 4*name_len) :: in_file, out_file, file_loc
+        integer :: iu_in, ios, i, rec_step, rec_id
+        integer, allocatable, dimension(:) :: iu_out
+        real(wp) :: rec_time
+        real(wp), dimension(3) :: rec_force, rec_torque
+
+        file_loc = trim(case_dir)//'/D'
+
+        in_file = trim(file_loc)//'/ib_force.dat'
+        open (newunit=iu_in, file=trim(in_file), form='unformatted', access='stream', &
+              status='old', action='read', iostat=ios)
+        if (ios /= 0) then
+            print *, 'ERROR: cannot open input file: ', trim(in_file), ' iostat=', ios
+            return
+        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
+                print *, 'ERROR: cannot open output file: ', trim(out_file), ' iostat=', ios
+                close (iu_in)
+                do rec_id = 1, i - 1
+                    close (iu_out(rec_id))
+                end do
+                deallocate (iu_out)
+                return
+            end if
+            write (iu_out(i), '(A)') 'mytime fx fy fz Tau_x Tau_y Tau_z'
+        end do
+
+        do
+            read (iu_in, iostat=ios) rec_step, rec_time, rec_id, rec_force, rec_torque
+            if (ios /= 0) exit  ! EOF (<0) or read error (>0)
+            if (rec_id >= 1 .and. rec_id <= num_ibs) then
+                write (iu_out(rec_id), '(7(ES24.16E3,1X))') rec_time, &
+                    rec_force(1), rec_force(2), rec_force(3), &
+                    rec_torque(1), rec_torque(2), rec_torque(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_force_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)
 

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_force_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_force_wrt = .false.
         lag_txt_wrt = .false.
         lag_header = .true.
         lag_db_wrt = .false.

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_force_wrt']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 

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_force_wrt
 
         ! Inquiring the status of the post_process.inp file
         file_loc = 'post_process.inp'

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_force_wrt) then
+        call s_write_ib_force_files()
+    end if
+
     close (11)
 
     call s_finalize_modules()

src/simulation/m_data_output.fpp

@@ -37,17 +37,20 @@ module m_data_output
               s_open_run_time_information_file, &
               s_open_com_files, &
               s_open_probe_files, &
+              s_open_ib_force_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_force_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_force_file, &
+              s_finalize_data_output_module
 
     real(wp), allocatable, dimension(:, :, :) :: icfl_sf  !< ICFL stability criterion
     real(wp), allocatable, dimension(:, :, :) :: vcfl_sf  !< VCFL stability criterion
@@ -254,6 +257,18 @@ contains
 
     end subroutine s_open_probe_files
 
+    impure subroutine s_open_ib_force_file
+        character(LEN=path_len + 2*name_len) :: file_loc
+
+        write (file_loc, '(A)') 'ib_force.dat'
+        file_loc = trim(case_dir)//'/D/'//trim(file_loc)
+        open (92, FILE=trim(file_loc), &
+              FORM='unformatted', &
+              ACCESS='stream', &
+              STATUS='replace', &
+              POSITION='append')
+    end subroutine s_open_ib_force_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 +1086,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 +1151,22 @@ contains
             call s_write_serial_ib_data(time_step)
         end if
 
-    end subroutine
+    end subroutine s_write_ib_data_file
+
+    !> @brief Writes IB force/torque records to D/ib_force.dat on rank 0.
+    impure subroutine s_write_ib_force_file(time_step)
+
+        integer, intent(in) :: time_step
+        if (proc_rank == 0) then
+            block
+                integer :: i
+                do i = 1, num_ibs
+                    write (92) time_step, mytime, i, patch_ib(i)%force, patch_ib(i)%torque
+                end do
+            end block
+        end if
+
+    end subroutine s_write_ib_force_file
 
     !>  This writes a formatted data file where the root processor
     !!      can write out the CoM information
@@ -1907,6 +1937,12 @@ contains
 
     end subroutine s_close_probe_files
 
+    impure subroutine s_close_ib_force_file
+
+        close (92)
+
+    end subroutine s_close_ib_force_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.

src/simulation/m_global_parameters.fpp

@@ -396,6 +396,7 @@ module m_global_parameters
     !> @{
     logical :: ib
     integer :: num_ibs
+    logical :: ib_force_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_force_wrt = .false.
 
         ! Bubble modeling
         bubbles_euler = .false.

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_force_wrt']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 

src/simulation/m_start_up.fpp

@@ -144,6 +144,7 @@ contains
             x_domain, y_domain, z_domain, &
             hypoelasticity, &
             ib, num_ibs, patch_ib, &
+            ib_force_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, &

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 data file
+        if (proc_rank == 0 .and. ib_force_wrt) then
+            call s_open_ib_force_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_force_wrt .and. s == nstage) then
+                    call s_compute_ib_forces(q_prim_vf, fluid_pp)
+                end if
+
+                ! Write IB forces to file if requested and at the RK final stage
+                if (ib_force_wrt .and. s == nstage) then
+                    call s_write_ib_force_file(t_step)
                 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_force_wrt) then
+            call s_close_ib_force_file()
+        end if
+
     end subroutine s_finalize_time_steppers_module
 
 end module m_time_steppers

toolchain/mfc/params/definitions.py

@@ -214,6 +214,7 @@
     "prim_vars_wrt": "Write primitive variables",
     "cons_vars_wrt": "Write conservative variables",
     "run_time_info": "Print runtime info",
+    "ib_force_wrt": "Write IB force/torque data",
     # Misc
     "case_dir": "Case directory path",
     "cantera_file": "Cantera mechanism file",
@@ -932,7 +933,7 @@ def _load():  # pylint: disable=too-many-locals,too-many-statements
     _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"]:
+              "cons_vars_wrt", "fft_wrt", "ib_force_wrt"]:
         _r(n, LOG, {"output"})
     for n in ["schlieren_wrt", "alpha_rho_wrt", "rho_wrt", "mom_wrt", "vel_wrt",
               "flux_wrt", "E_wrt", "pres_wrt", "alpha_wrt", "kappa_wrt", "gamma_wrt",