Skip to content

Commit e9910d9

Browse files
danieljvickersdanieljvickers
authored
IB Particle Mesh Output (#1375)
1 parent 92585b1 commit e9910d9

10 files changed

Lines changed: 305 additions & 100 deletions

File tree

‎examples/2D_mibm_shock_cylinder/case.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,7 @@
8787
"precision": 2,
8888
"prim_vars_wrt": "T",
8989
"E_wrt": "T",
90+
"ib_state_wrt": "T",
9091
"parallel_io": "T",
9192
# Patch: Constant Tube filled with air
9293
# Specify the cylindrical air tube grid geometry
@@ -128,7 +129,7 @@
128129
"patch_ib(1)%angular_vel(1)": 0.0, # x-axis rotational velocity in radians per second
129130
"patch_ib(1)%angular_vel(2)": 0.0, # y-axis rotation
130131
"patch_ib(1)%angular_vel(3)": 0.0, # z-axis rotation
131-
"patch_ib(1)%mass": 0.5, # z-axis rotation
132+
"patch_ib(1)%mass": 0.25, # z-axis rotation
132133
# Fluids Physical Parameters
133134
"fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00), # 2.50(Not 1.40)
134135
"fluid_pp(1)%pi_inf": 0,

‎src/post_process/m_data_output.fpp‎

Lines changed: 133 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,8 @@ module m_data_output
1919
& s_open_intf_data_file, s_open_energy_data_file, s_write_grid_to_formatted_database_file, &
2020
& s_write_variable_to_formatted_database_file, s_write_lag_bubbles_results_to_text, &
2121
& s_write_lag_bubbles_to_formatted_database_file, s_write_ib_state_files, s_write_intf_data_file, &
22-
& s_write_energy_data_file, s_close_formatted_database_file, s_close_intf_data_file, s_close_energy_data_file, &
23-
& s_finalize_data_output_module
22+
& s_write_energy_data_file, s_write_ib_bodies_to_formatted_database_file, s_close_formatted_database_file, &
23+
& s_close_intf_data_file, s_close_energy_data_file, s_finalize_data_output_module
2424

2525
! Include Silo-HDF5 interface library
2626
include 'silo_f9x.inc'
@@ -1341,6 +1341,137 @@ contains
13411341
13421342
end subroutine s_write_energy_data_file
13431343
1344+
!> Read IB state and write a Silo point mesh with per-body scalar fields.
1345+
impure subroutine s_write_ib_bodies_to_formatted_database_file(t_step)
1346+
1347+
integer, intent(in) :: t_step
1348+
character(len=len_trim(case_dir) + 3*name_len) :: file_loc
1349+
1350+
#ifdef MFC_MPI
1351+
integer, parameter :: NFIELDS_PER_IB = 20
1352+
real(wp) :: ib_buf(NFIELDS_PER_IB)
1353+
real(wp), dimension(:,:), allocatable :: ib_data
1354+
logical :: file_exist
1355+
character(LEN=4*name_len), dimension(num_procs) :: meshnames
1356+
integer, dimension(num_procs) :: meshtypes
1357+
integer :: i, ios, file_unit
1358+
integer :: ierr, nBodies
1359+
real(wp), dimension(:), allocatable :: px, py, pz
1360+
real(wp), dimension(:), allocatable :: force_x, force_y, force_z
1361+
real(wp), dimension(:), allocatable :: torque_x, torque_y, torque_z
1362+
real(wp), dimension(:), allocatable :: vel_x, vel_y, vel_z
1363+
real(wp), dimension(:), allocatable :: omega_x, omega_y, omega_z
1364+
real(wp), dimension(:), allocatable :: angle_x, angle_y, angle_z
1365+
real(wp), dimension(:), allocatable :: ib_diameter
1366+
1367+
! Build path to per-timestep IB state file
1368+
write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
1369+
file_loc = trim(case_dir) // trim(file_loc)
1370+
1371+
inquire (FILE=trim(file_loc), EXIST=file_exist)
1372+
if (.not. file_exist) then
1373+
call s_mpi_abort('Restart file ' // trim(file_loc) // ' does not exist!')
1374+
end if
1375+
1376+
nBodies = num_ibs
1377+
1378+
if (nBodies > 0) then
1379+
allocate (ib_data(nBodies, NFIELDS_PER_IB))
1380+
allocate (px(nBodies), py(nBodies), pz(nBodies))
1381+
allocate (force_x(nBodies), force_y(nBodies), force_z(nBodies))
1382+
allocate (torque_x(nBodies), torque_y(nBodies), torque_z(nBodies))
1383+
allocate (vel_x(nBodies), vel_y(nBodies), vel_z(nBodies))
1384+
allocate (omega_x(nBodies), omega_y(nBodies), omega_z(nBodies))
1385+
allocate (angle_x(nBodies), angle_y(nBodies), angle_z(nBodies))
1386+
allocate (ib_diameter(nBodies))
1387+
1388+
if (proc_rank == 0) then
1389+
open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios)
1390+
if (ios /= 0) call s_mpi_abort('Cannot open IB state file: ' // trim(file_loc))
1391+
1392+
do i = 1, nBodies
1393+
read (file_unit, iostat=ios) ib_buf
1394+
if (ios /= 0) call s_mpi_abort('Error reading IB state file')
1395+
ib_data(i,:) = ib_buf(:)
1396+
end do
1397+
1398+
close (file_unit)
1399+
end if
1400+
1401+
call MPI_BCAST(ib_data, nBodies*NFIELDS_PER_IB, mpi_p, 0, MPI_COMM_WORLD, ierr)
1402+
1403+
do i = 1, nBodies
1404+
force_x(i) = ib_data(i, 2); force_y(i) = ib_data(i, 3); force_z(i) = ib_data(i, 4)
1405+
torque_x(i) = ib_data(i, 5); torque_y(i) = ib_data(i, 6); torque_z(i) = ib_data(i, 7)
1406+
vel_x(i) = ib_data(i, 8); vel_y(i) = ib_data(i, 9); vel_z(i) = ib_data(i, 10)
1407+
omega_x(i) = ib_data(i, 11); omega_y(i) = ib_data(i, 12); omega_z(i) = ib_data(i, 13)
1408+
angle_x(i) = ib_data(i, 14); angle_y(i) = ib_data(i, 15); angle_z(i) = ib_data(i, 16)
1409+
px(i) = ib_data(i, 17); py(i) = ib_data(i, 18); pz(i) = ib_data(i, 19)
1410+
ib_diameter(i) = ib_data(i, 20)*2.0_wp
1411+
end do
1412+
1413+
if (proc_rank == 0) then
1414+
do i = 1, num_procs
1415+
write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:ib_bodies'
1416+
meshtypes(i) = DB_POINTMESH
1417+
end do
1418+
err = DBSET2DSTRLEN(len(meshnames(1)))
1419+
err = DBPUTMMESH(dbroot, 'ib_bodies', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, DB_F77NULL, ierr)
1420+
end if
1421+
1422+
err = DBPUTPM(dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
1423+
1424+
call s_write_ib_variable('ib_force_x', t_step, force_x, nBodies)
1425+
call s_write_ib_variable('ib_force_y', t_step, force_y, nBodies)
1426+
call s_write_ib_variable('ib_force_z', t_step, force_z, nBodies)
1427+
call s_write_ib_variable('ib_torque_x', t_step, torque_x, nBodies)
1428+
call s_write_ib_variable('ib_torque_y', t_step, torque_y, nBodies)
1429+
call s_write_ib_variable('ib_torque_z', t_step, torque_z, nBodies)
1430+
call s_write_ib_variable('ib_vel_x', t_step, vel_x, nBodies)
1431+
call s_write_ib_variable('ib_vel_y', t_step, vel_y, nBodies)
1432+
call s_write_ib_variable('ib_vel_z', t_step, vel_z, nBodies)
1433+
call s_write_ib_variable('ib_omega_x', t_step, omega_x, nBodies)
1434+
call s_write_ib_variable('ib_omega_y', t_step, omega_y, nBodies)
1435+
call s_write_ib_variable('ib_omega_z', t_step, omega_z, nBodies)
1436+
call s_write_ib_variable('ib_angle_x', t_step, angle_x, nBodies)
1437+
call s_write_ib_variable('ib_angle_y', t_step, angle_y, nBodies)
1438+
call s_write_ib_variable('ib_angle_z', t_step, angle_z, nBodies)
1439+
call s_write_ib_variable('ib_diameter', t_step, ib_diameter, nBodies)
1440+
1441+
deallocate (ib_data, px, py, pz, force_x, force_y, force_z)
1442+
deallocate (torque_x, torque_y, torque_z, vel_x, vel_y, vel_z)
1443+
deallocate (omega_x, omega_y, omega_z, angle_x, angle_y, angle_z)
1444+
deallocate (ib_diameter)
1445+
end if
1446+
#endif
1447+
1448+
end subroutine s_write_ib_bodies_to_formatted_database_file
1449+
1450+
!> Write a single IB point-variable to the Silo database slave and master files.
1451+
subroutine s_write_ib_variable(varname, t_step, data, nBodies)
1452+
1453+
character(len=*), intent(in) :: varname
1454+
integer, intent(in) :: t_step
1455+
real(wp), dimension(:), intent(in) :: data
1456+
integer, intent(in) :: nBodies
1457+
character(len=4*name_len), dimension(num_procs) :: var_names
1458+
integer, dimension(num_procs) :: var_types
1459+
integer :: ierr, i
1460+
1461+
if (proc_rank == 0) then
1462+
do i = 1, num_procs
1463+
write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:' // trim(varname)
1464+
var_types(i) = DB_POINTVAR
1465+
end do
1466+
err = DBSET2DSTRLEN(len(var_names(1)))
1467+
err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), num_procs, var_names, len_trim(var_names), var_types, &
1468+
& DB_F77NULL, ierr)
1469+
end if
1470+
1471+
err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
1472+
1473+
end subroutine s_write_ib_variable
1474+
13441475
!> Close the formatted database slave file and, for the root process, the master file.
13451476
impure subroutine s_close_formatted_database_file()
13461477

‎src/post_process/m_start_up.fpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -746,6 +746,8 @@ contains
746746
if (lag_db_wrt) call s_write_lag_bubbles_to_formatted_database_file(t_step) ! silo file output
747747
end if
748748
749+
if (ib_state_wrt) call s_write_ib_bodies_to_formatted_database_file(t_step)
750+
749751
if (sim_data .and. proc_rank == 0) then
750752
call s_close_intf_data_file()
751753
call s_close_energy_data_file()

‎src/post_process/p_main.fpp‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -80,10 +80,6 @@ program p_main
8080
end do
8181
! END: Time-Marching Loop
8282

83-
if (proc_rank == 0 .and. ib_state_wrt) then
84-
call s_write_ib_state_files()
85-
end if
86-
8783
close (11)
8884

8985
call s_finalize_modules()

‎src/simulation/m_data_output.fpp‎

Lines changed: 122 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -24,15 +24,13 @@ module m_data_output
2424

2525
private
2626
public :: s_initialize_data_output_module, s_open_run_time_information_file, s_open_com_files, s_open_probe_files, &
27-
& s_open_ib_state_file, s_write_run_time_information, s_write_data_files, s_write_serial_data_files, &
28-
& s_write_parallel_data_files, s_write_ib_data_file, s_write_com_files, s_write_probe_files, s_write_ib_state_file, &
29-
& s_close_run_time_information_file, s_close_com_files, s_close_probe_files, s_close_ib_state_file, &
30-
& s_finalize_data_output_module
31-
32-
integer :: ib_state_unit = -1 !< I/O unit for IB state binary file
33-
real(wp), allocatable, dimension(:,:,:) :: icfl_sf !< ICFL stability criterion
34-
real(wp), allocatable, dimension(:,:,:) :: vcfl_sf !< VCFL stability criterion
35-
real(wp), allocatable, dimension(:,:,:) :: Rc_sf !< Rc stability criterion
27+
& s_write_run_time_information, s_write_data_files, s_write_serial_data_files, s_write_parallel_data_files, &
28+
& s_write_ib_data_file, s_write_com_files, s_write_probe_files, s_write_ib_state_file, s_close_run_time_information_file, &
29+
& s_close_com_files, s_close_probe_files, s_finalize_data_output_module
30+
31+
real(wp), allocatable, dimension(:,:,:) :: icfl_sf !< ICFL stability criterion
32+
real(wp), allocatable, dimension(:,:,:) :: vcfl_sf !< VCFL stability criterion
33+
real(wp), allocatable, dimension(:,:,:) :: Rc_sf !< Rc stability criterion
3634
real(wp), public, allocatable, dimension(:,:) :: c_mass
3735
$:GPU_DECLARE(create='[icfl_sf, vcfl_sf, Rc_sf, c_mass]')
3836

@@ -161,26 +159,6 @@ contains
161159

162160
end subroutine s_open_probe_files
163161

164-
!> Open the immersed boundary state file for binary output
165-
impure subroutine s_open_ib_state_file
166-
167-
character(len=path_len + 2*name_len) :: file_loc
168-
integer :: ios
169-
170-
call s_create_directory(trim(case_dir) // '/restart_data')
171-
write (file_loc, '(A)') 'ib_state.dat'
172-
file_loc = trim(case_dir) // '/restart_data/' // trim(file_loc)
173-
if (t_step_start > 0) then
174-
! On restart, append to existing file to preserve history
175-
open (newunit=ib_state_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', &
176-
& position='append', iostat=ios)
177-
else
178-
open (newunit=ib_state_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios)
179-
end if
180-
if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc))
181-
182-
end subroutine s_open_ib_state_file
183-
184162
!> Write stability criteria extrema to the run-time information file at the given time step
185163
impure subroutine s_write_run_time_information(q_prim_vf, t_step)
186164

@@ -918,16 +896,125 @@ contains
918896

919897
end subroutine s_write_ib_data_file
920898

921-
!> @brief Writes IB state records to restart_data/ib_state.dat. Must be called only on rank 0.
922-
impure subroutine s_write_ib_state_file()
899+
!> Writes the IB state information out to file
900+
subroutine s_write_parallel_ib_state(t_step)
923901

924-
integer :: i
902+
integer, intent(in) :: t_step
903+
904+
#ifdef MFC_MPI
905+
character(LEN=path_len + 2*name_len) :: file_loc
906+
integer(kind=MPI_OFFSET_KIND) :: disp
907+
integer(kind=MPI_OFFSET_KIND) :: WP_MOK
908+
integer :: ifile, ierr
909+
integer, dimension(MPI_STATUS_SIZE) :: status
910+
logical :: file_exist
911+
integer :: i
912+
integer, parameter :: NFIELDS_PER_IB = 20
913+
real(wp) :: ib_buf(NFIELDS_PER_IB)
914+
915+
! Partition IBs across ranks round-robin style
916+
integer :: ib_start, ib_end, nibs_per_rank, remainder
917+
918+
WP_MOK = int(storage_size(0._wp)/8, MPI_OFFSET_KIND)
919+
920+
if (proc_rank == 0) then
921+
call s_create_directory(trim(case_dir) // '/restart_data')
922+
end if
923+
call s_mpi_barrier()
924+
925+
! Divide num_ibs across num_procs
926+
nibs_per_rank = num_ibs/num_procs
927+
remainder = mod(num_ibs, num_procs)
928+
929+
! Ranks < remainder get one extra IB
930+
if (proc_rank < remainder) then
931+
ib_start = proc_rank*(nibs_per_rank + 1) + 1
932+
ib_end = ib_start + nibs_per_rank ! nibs_per_rank + 1 total
933+
else
934+
ib_start = remainder*(nibs_per_rank + 1) + (proc_rank - remainder)*nibs_per_rank + 1
935+
ib_end = ib_start + nibs_per_rank - 1
936+
end if
937+
938+
write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
939+
file_loc = trim(case_dir) // trim(file_loc)
940+
941+
inquire (FILE=trim(file_loc), EXIST=file_exist)
942+
if (file_exist .and. proc_rank == 0) then
943+
call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr)
944+
end if
945+
call s_mpi_barrier()
946+
947+
call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr)
948+
949+
do i = ib_start, ib_end
950+
ib_buf(1) = mytime
951+
ib_buf(2:4) = patch_ib(i)%force(1:3)
952+
ib_buf(5:7) = patch_ib(i)%torque(1:3)
953+
ib_buf(8:10) = patch_ib(i)%vel(1:3)
954+
ib_buf(11:13) = patch_ib(i)%angular_vel(1:3)
955+
ib_buf(14:16) = patch_ib(i)%angles(1:3)
956+
ib_buf(17) = patch_ib(i)%x_centroid
957+
ib_buf(18) = patch_ib(i)%y_centroid
958+
ib_buf(19) = patch_ib(i)%z_centroid
959+
ib_buf(20) = patch_ib(i)%radius
960+
961+
! Global IB index (i) determines position in file
962+
disp = int(i - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK
963+
964+
call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr)
965+
end do
966+
967+
call MPI_FILE_CLOSE(ifile, ierr)
968+
#endif
969+
970+
end subroutine s_write_parallel_ib_state
971+
972+
!> Write IB state data to a per-timestep serial (unformatted) file
973+
subroutine s_write_serial_ib_state(t_step)
974+
975+
integer, intent(in) :: t_step
976+
character(LEN=path_len + 2*name_len) :: file_loc
977+
integer :: i, ios, file_unit
978+
integer, parameter :: NFIELDS_PER_IB = 20
979+
real(wp) :: ib_buf(NFIELDS_PER_IB)
980+
981+
call s_create_directory(trim(case_dir) // '/restart_data')
982+
983+
write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
984+
file_loc = trim(case_dir) // trim(file_loc)
985+
986+
open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios)
987+
if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc))
925988

926989
do i = 1, num_ibs
927-
write (ib_state_unit) mytime, i, patch_ib(i)%force, patch_ib(i)%torque, patch_ib(i)%vel, patch_ib(i)%angular_vel, &
928-
& patch_ib(i)%angles, patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, patch_ib(i)%z_centroid
990+
ib_buf(1) = mytime
991+
ib_buf(2:4) = patch_ib(i)%force(1:3)
992+
ib_buf(5:7) = patch_ib(i)%torque(1:3)
993+
ib_buf(8:10) = patch_ib(i)%vel(1:3)
994+
ib_buf(11:13) = patch_ib(i)%angular_vel(1:3)
995+
ib_buf(14:16) = patch_ib(i)%angles(1:3)
996+
ib_buf(17) = patch_ib(i)%x_centroid
997+
ib_buf(18) = patch_ib(i)%y_centroid
998+
ib_buf(19) = patch_ib(i)%z_centroid
999+
ib_buf(20) = patch_ib(i)%radius
1000+
1001+
write (file_unit) ib_buf
9291002
end do
930-
flush (ib_state_unit)
1003+
1004+
close (file_unit)
1005+
1006+
end subroutine s_write_serial_ib_state
1007+
1008+
!> @brief Writes IB state records to restart_data/ib_state.dat. Must be called only on rank 0.
1009+
impure subroutine s_write_ib_state_file(time_step)
1010+
1011+
integer, intent(in) :: time_step
1012+
1013+
if (parallel_io) then
1014+
call s_write_parallel_ib_state(time_step)
1015+
else
1016+
call s_write_serial_ib_state(time_step)
1017+
end if
9311018

9321019
end subroutine s_write_ib_state_file
9331020

@@ -1543,13 +1630,6 @@ contains
15431630

15441631
end subroutine s_close_probe_files
15451632

1546-
!> Close the immersed boundary state file
1547-
impure subroutine s_close_ib_state_file
1548-
1549-
close (ib_state_unit)
1550-
1551-
end subroutine s_close_ib_state_file
1552-
15531633
!> Initialize the data output module
15541634
impure subroutine s_initialize_data_output_module
15551635

0 commit comments

Comments
 (0)