IB Particle Mesh Output

examples/2D_mibm_shock_cylinder/case.py

@@ -87,6 +87,7 @@
             "precision": 2,
             "prim_vars_wrt": "T",
             "E_wrt": "T",
+            "ib_state_wrt": "T",
             "parallel_io": "T",
             # Patch: Constant Tube filled with air
             # Specify the cylindrical air tube grid geometry
@@ -128,7 +129,7 @@
             "patch_ib(1)%angular_vel(1)": 0.0,  # x-axis rotational velocity in radians per second
             "patch_ib(1)%angular_vel(2)": 0.0,  # y-axis rotation
             "patch_ib(1)%angular_vel(3)": 0.0,  # z-axis rotation
-            "patch_ib(1)%mass": 0.5,  # z-axis rotation
+            "patch_ib(1)%mass": 0.25,  # z-axis rotation
             # Fluids Physical Parameters
             "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50(Not 1.40)
             "fluid_pp(1)%pi_inf": 0,

src/post_process/m_data_output.fpp

@@ -19,8 +19,8 @@ module m_data_output
         & s_open_intf_data_file, s_open_energy_data_file, s_write_grid_to_formatted_database_file, &
         & 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, s_close_intf_data_file, s_close_energy_data_file, &
-        & s_finalize_data_output_module
+        & s_write_energy_data_file, s_write_ib_bodies_to_formatted_database_file, s_close_formatted_database_file, &
+        & s_close_intf_data_file, s_close_energy_data_file, s_finalize_data_output_module
 
     ! Include Silo-HDF5 interface library
     include 'silo_f9x.inc'
@@ -1341,6 +1341,137 @@ contains
 
     end subroutine s_write_energy_data_file
 
+    !> Read IB state and write a Silo point mesh with per-body scalar fields.
+    impure subroutine s_write_ib_bodies_to_formatted_database_file(t_step)
+
+        integer, intent(in)                            :: t_step
+        character(len=len_trim(case_dir) + 3*name_len) :: file_loc
+
+#ifdef MFC_MPI
+        integer, parameter                              :: NFIELDS_PER_IB = 20
+        real(wp)                                        :: ib_buf(NFIELDS_PER_IB)
+        real(wp), dimension(:,:), allocatable           :: ib_data
+        logical                                         :: file_exist
+        character(LEN=4*name_len), dimension(num_procs) :: meshnames
+        integer, dimension(num_procs)                   :: meshtypes
+        integer                                         :: i, ios, file_unit
+        integer                                         :: ierr, nBodies
+        real(wp), dimension(:), allocatable             :: px, py, pz
+        real(wp), dimension(:), allocatable             :: force_x, force_y, force_z
+        real(wp), dimension(:), allocatable             :: torque_x, torque_y, torque_z
+        real(wp), dimension(:), allocatable             :: vel_x, vel_y, vel_z
+        real(wp), dimension(:), allocatable             :: omega_x, omega_y, omega_z
+        real(wp), dimension(:), allocatable             :: angle_x, angle_y, angle_z
+        real(wp), dimension(:), allocatable             :: ib_diameter
+
+        ! Build path to per-timestep IB state file
+        write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
+        file_loc = trim(case_dir) // trim(file_loc)
+
+        inquire (FILE=trim(file_loc), EXIST=file_exist)
+        if (.not. file_exist) then
+            call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
+        end if
+
+        nBodies = num_ibs
+
+        if (nBodies > 0) then
+            allocate (ib_data(nBodies, NFIELDS_PER_IB))
+            allocate (px(nBodies), py(nBodies), pz(nBodies))
+            allocate (force_x(nBodies), force_y(nBodies), force_z(nBodies))
+            allocate (torque_x(nBodies), torque_y(nBodies), torque_z(nBodies))
+            allocate (vel_x(nBodies), vel_y(nBodies), vel_z(nBodies))
+            allocate (omega_x(nBodies), omega_y(nBodies), omega_z(nBodies))
+            allocate (angle_x(nBodies), angle_y(nBodies), angle_z(nBodies))
+            allocate (ib_diameter(nBodies))
+
+            if (proc_rank == 0) then
+                open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
+                if (ios /= 0) call s_mpi_abort('Cannot open IB state file: ' // trim(file_loc))
+
+                do i = 1, nBodies
+                    read (file_unit, iostat=ios) ib_buf
+                    if (ios /= 0) call s_mpi_abort('Error reading IB state file')
+                    ib_data(i,:) = ib_buf(:)
+                end do
+
+                close (file_unit)
+            end if
+
+            call MPI_BCAST(ib_data, nBodies*NFIELDS_PER_IB, mpi_p, 0, MPI_COMM_WORLD, ierr)
+
+            do i = 1, nBodies
+                force_x(i) = ib_data(i, 2); force_y(i) = ib_data(i, 3); force_z(i) = ib_data(i, 4)
+                torque_x(i) = ib_data(i, 5); torque_y(i) = ib_data(i, 6); torque_z(i) = ib_data(i, 7)
+                vel_x(i) = ib_data(i, 8); vel_y(i) = ib_data(i, 9); vel_z(i) = ib_data(i, 10)
+                omega_x(i) = ib_data(i, 11); omega_y(i) = ib_data(i, 12); omega_z(i) = ib_data(i, 13)
+                angle_x(i) = ib_data(i, 14); angle_y(i) = ib_data(i, 15); angle_z(i) = ib_data(i, 16)
+                px(i) = ib_data(i, 17); py(i) = ib_data(i, 18); pz(i) = ib_data(i, 19)
+                ib_diameter(i) = ib_data(i, 20)*2.0_wp
+            end do
+
+            if (proc_rank == 0) then
+                do i = 1, num_procs
+                    write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:ib_bodies'
+                    meshtypes(i) = DB_POINTMESH
+                end do
+                err = DBSET2DSTRLEN(len(meshnames(1)))
+                err = DBPUTMMESH(dbroot, 'ib_bodies', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, DB_F77NULL, ierr)
+            end if
+
+            err = DBPUTPM(dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
+
+            call s_write_ib_variable('ib_force_x', t_step, force_x, nBodies)
+            call s_write_ib_variable('ib_force_y', t_step, force_y, nBodies)
+            call s_write_ib_variable('ib_force_z', t_step, force_z, nBodies)
+            call s_write_ib_variable('ib_torque_x', t_step, torque_x, nBodies)
+            call s_write_ib_variable('ib_torque_y', t_step, torque_y, nBodies)
+            call s_write_ib_variable('ib_torque_z', t_step, torque_z, nBodies)
+            call s_write_ib_variable('ib_vel_x', t_step, vel_x, nBodies)
+            call s_write_ib_variable('ib_vel_y', t_step, vel_y, nBodies)
+            call s_write_ib_variable('ib_vel_z', t_step, vel_z, nBodies)
+            call s_write_ib_variable('ib_omega_x', t_step, omega_x, nBodies)
+            call s_write_ib_variable('ib_omega_y', t_step, omega_y, nBodies)
+            call s_write_ib_variable('ib_omega_z', t_step, omega_z, nBodies)
+            call s_write_ib_variable('ib_angle_x', t_step, angle_x, nBodies)
+            call s_write_ib_variable('ib_angle_y', t_step, angle_y, nBodies)
+            call s_write_ib_variable('ib_angle_z', t_step, angle_z, nBodies)
+            call s_write_ib_variable('ib_diameter', t_step, ib_diameter, nBodies)
+
+            deallocate (ib_data, px, py, pz, force_x, force_y, force_z)
+            deallocate (torque_x, torque_y, torque_z, vel_x, vel_y, vel_z)
+            deallocate (omega_x, omega_y, omega_z, angle_x, angle_y, angle_z)
+            deallocate (ib_diameter)
+        end if
+#endif
+
+    end subroutine s_write_ib_bodies_to_formatted_database_file
+
+    !> Write a single IB point-variable to the Silo database slave and master files.
+    subroutine s_write_ib_variable(varname, t_step, data, nBodies)
+
+        character(len=*), intent(in)                    :: varname
+        integer, intent(in)                             :: t_step
+        real(wp), dimension(:), intent(in)              :: data
+        integer, intent(in)                             :: nBodies
+        character(len=4*name_len), dimension(num_procs) :: var_names
+        integer, dimension(num_procs)                   :: var_types
+        integer                                         :: ierr, i
+
+        if (proc_rank == 0) then
+            do i = 1, num_procs
+                write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:' // trim(varname)
+                var_types(i) = DB_POINTVAR
+            end do
+            err = DBSET2DSTRLEN(len(var_names(1)))
+            err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), num_procs, var_names, len_trim(var_names), var_types, &
+                            & DB_F77NULL, ierr)
+        end if
+
+        err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
+
+    end subroutine s_write_ib_variable
+
     !> Close the formatted database slave file and, for the root process, the master file.
     impure subroutine s_close_formatted_database_file()
 

src/post_process/m_start_up.fpp

@@ -746,6 +746,8 @@ contains
             if (lag_db_wrt) call s_write_lag_bubbles_to_formatted_database_file(t_step)  ! silo file output
         end if
 
+        if (ib_state_wrt) call s_write_ib_bodies_to_formatted_database_file(t_step)
+
         if (sim_data .and. proc_rank == 0) then
             call s_close_intf_data_file()
             call s_close_energy_data_file()

src/post_process/p_main.fpp

@@ -80,10 +80,6 @@ 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()

src/simulation/m_data_output.fpp

@@ -24,15 +24,13 @@ 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_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_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
-    real(wp), allocatable, dimension(:,:,:)       :: Rc_sf               !< Rc stability criterion
+        & 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
+
+    real(wp), allocatable, dimension(:,:,:)       :: icfl_sf  !< ICFL stability criterion
+    real(wp), allocatable, dimension(:,:,:)       :: vcfl_sf  !< VCFL stability criterion
+    real(wp), allocatable, dimension(:,:,:)       :: Rc_sf    !< Rc stability criterion
     real(wp), public, allocatable, dimension(:,:) :: c_mass
     $:GPU_DECLARE(create='[icfl_sf, vcfl_sf, Rc_sf, c_mass]')
 
@@ -161,26 +159,6 @@ contains
 
     end subroutine s_open_probe_files
 
-    !> Open the immersed boundary state file for binary output
-    impure subroutine s_open_ib_state_file
-
-        character(len=path_len + 2*name_len) :: file_loc
-        integer                              :: ios
-
-        call s_create_directory(trim(case_dir) // '/restart_data')
-        write (file_loc, '(A)') 'ib_state.dat'
-        file_loc = trim(case_dir) // '/restart_data/' // trim(file_loc)
-        if (t_step_start > 0) then
-            ! On restart, append to existing file to preserve history
-            open (newunit=ib_state_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', &
-                  & position='append', iostat=ios)
-        else
-            open (newunit=ib_state_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios)
-        end if
-        if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc))
-
-    end subroutine s_open_ib_state_file
-
     !> Write stability criteria extrema to the run-time information file at the given time step
     impure subroutine s_write_run_time_information(q_prim_vf, t_step)
 
@@ -918,16 +896,125 @@ contains
 
     end subroutine s_write_ib_data_file
 
-    !> @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()
+    !> Writes the IB state information out to file
+    subroutine s_write_parallel_ib_state(t_step)
 
-        integer :: i
+        integer, intent(in) :: t_step
+
+#ifdef MFC_MPI
+        character(LEN=path_len + 2*name_len) :: file_loc
+        integer(kind=MPI_OFFSET_KIND)        :: disp
+        integer(kind=MPI_OFFSET_KIND)        :: WP_MOK
+        integer                              :: ifile, ierr
+        integer, dimension(MPI_STATUS_SIZE)  :: status
+        logical                              :: file_exist
+        integer                              :: i
+        integer, parameter                   :: NFIELDS_PER_IB = 20
+        real(wp)                             :: ib_buf(NFIELDS_PER_IB)
+
+        ! Partition IBs across ranks round-robin style
+        integer :: ib_start, ib_end, nibs_per_rank, remainder
+
+        WP_MOK = int(storage_size(0._wp)/8, MPI_OFFSET_KIND)
+
+        if (proc_rank == 0) then
+            call s_create_directory(trim(case_dir) // '/restart_data')
+        end if
+        call s_mpi_barrier()
+
+        ! Divide num_ibs across num_procs
+        nibs_per_rank = num_ibs/num_procs
+        remainder = mod(num_ibs, num_procs)
+
+        ! Ranks < remainder get one extra IB
+        if (proc_rank < remainder) then
+            ib_start = proc_rank*(nibs_per_rank + 1) + 1
+            ib_end = ib_start + nibs_per_rank  ! nibs_per_rank + 1 total
+        else
+            ib_start = remainder*(nibs_per_rank + 1) + (proc_rank - remainder)*nibs_per_rank + 1
+            ib_end = ib_start + nibs_per_rank - 1
+        end if
+
+        write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
+        file_loc = trim(case_dir) // trim(file_loc)
+
+        inquire (FILE=trim(file_loc), EXIST=file_exist)
+        if (file_exist .and. proc_rank == 0) then
+            call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr)
+        end if
+        call s_mpi_barrier()
+
+        call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr)
+
+        do i = ib_start, ib_end
+            ib_buf(1) = mytime
+            ib_buf(2:4) = patch_ib(i)%force(1:3)
+            ib_buf(5:7) = patch_ib(i)%torque(1:3)
+            ib_buf(8:10) = patch_ib(i)%vel(1:3)
+            ib_buf(11:13) = patch_ib(i)%angular_vel(1:3)
+            ib_buf(14:16) = patch_ib(i)%angles(1:3)
+            ib_buf(17) = patch_ib(i)%x_centroid
+            ib_buf(18) = patch_ib(i)%y_centroid
+            ib_buf(19) = patch_ib(i)%z_centroid
+            ib_buf(20) = patch_ib(i)%radius
+
+            ! Global IB index (i) determines position in file
+            disp = int(i - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK
+
+            call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr)
+        end do
+
+        call MPI_FILE_CLOSE(ifile, ierr)
+#endif
+
+    end subroutine s_write_parallel_ib_state
+
+    !> Write IB state data to a per-timestep serial (unformatted) file
+    subroutine s_write_serial_ib_state(t_step)
+
+        integer, intent(in)                  :: t_step
+        character(LEN=path_len + 2*name_len) :: file_loc
+        integer                              :: i, ios, file_unit
+        integer, parameter                   :: NFIELDS_PER_IB = 20
+        real(wp)                             :: ib_buf(NFIELDS_PER_IB)
+
+        call s_create_directory(trim(case_dir) // '/restart_data')
+
+        write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
+        file_loc = trim(case_dir) // trim(file_loc)
+
+        open (newunit=file_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))
 
         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
+            ib_buf(1) = mytime
+            ib_buf(2:4) = patch_ib(i)%force(1:3)
+            ib_buf(5:7) = patch_ib(i)%torque(1:3)
+            ib_buf(8:10) = patch_ib(i)%vel(1:3)
+            ib_buf(11:13) = patch_ib(i)%angular_vel(1:3)
+            ib_buf(14:16) = patch_ib(i)%angles(1:3)
+            ib_buf(17) = patch_ib(i)%x_centroid
+            ib_buf(18) = patch_ib(i)%y_centroid
+            ib_buf(19) = patch_ib(i)%z_centroid
+            ib_buf(20) = patch_ib(i)%radius
+
+            write (file_unit) ib_buf
         end do
-        flush (ib_state_unit)
+
+        close (file_unit)
+
+    end subroutine s_write_serial_ib_state
+
+    !> @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)
+
+        integer, intent(in) :: time_step
+
+        if (parallel_io) then
+            call s_write_parallel_ib_state(time_step)
+        else
+            call s_write_serial_ib_state(time_step)
+        end if
 
     end subroutine s_write_ib_state_file
 
@@ -1543,13 +1630,6 @@ contains
 
     end subroutine s_close_probe_files
 
-    !> Close the immersed boundary state file
-    impure subroutine s_close_ib_state_file
-
-        close (ib_state_unit)
-
-    end subroutine s_close_ib_state_file
-
     !> Initialize the data output module
     impure subroutine s_initialize_data_output_module