!> 读写粒子二进制文件 module spx_pif_h5part use spx_kinds, only: rk use spx_command_line, only: cli_obj use spx_particle, only: particle_type use spx_configuration, only: ioc_obj use spx_pif_namelist, only: pnl_obj use spx_env, only: stat, float_obj, rgn_obj use spx_logger, only: lgr_obj, tmr_obj use fffc_module, only: operator(.join.), terminal_obj, nowtime use h5part_module, only: h5pt_writefileattrib, h5pt_writestepattrib, h5pt_close, h5pt_openr, h5pt_openw, & h5pt_setnpoints, h5pt_setstep, h5pt_writedata, h5pt_readdata implicit none private public :: pif_obj !> 粒子二进制文件对象 !> @note 粒子集数组分区：| 实弹性粒子 | 实刚体粒子 | 虚粒子 | type pif_h5part_type integer(8), private :: fid !! 输出文件句柄 integer, private :: fid_stat !! stat.csv 文件句柄 integer, private :: fid_float !! float.csv 文件句柄 integer, pointer :: numbers(:) !! 弹性实、刚性浮体实、边界虚、总粒子数，刚性浮体数量 integer, pointer :: float_objects(:, :) !! 浮体对象 integer :: nreal !! 实粒子数 integer(8) :: istep !! 步数 type(particle_type), pointer :: ptc !! 粒子集 real(rk), allocatable :: z(:) contains procedure :: read, write, open, close end type pif_h5part_type type(pif_h5part_type) :: pif_obj !! 粒子二进制文件对象 contains !> 读取粒子二进制文件 subroutine read (self, region) class(pif_h5part_type), intent(inout) :: self type(particle_type), intent(inout), target :: region real(rk), allocatable :: z(:, :) integer(8) :: istat, fid integer :: i self%numbers => pnl_obj%numbers self%float_objects => pnl_obj%float_objects self%ptc => region ! 重定向粒子集 self%ptc%hsml = pnl_obj%hsml self%nreal = self%numbers(1) + self%numbers(2) fid = h5pt_openr(cli_obj%working_directory.join.ioc_obj%input_file) if (fid == 0_8) then call terminal_obj%error('open pif h5part file failed') end if allocate (self%ptc%loc(2, self%numbers(4)), & self%ptc%vel(2, self%numbers(4)), & self%ptc%acc(2, self%numbers(4)), & self%ptc%rho(self%numbers(4)), & self%ptc%p(self%numbers(4)), & self%ptc%mass(self%numbers(4)), & self%ptc%u(self%numbers(4)), & self%ptc%drho(self%numbers(4)), & self%ptc%du(self%numbers(4)), & self%ptc%c(self%numbers(4)), & source=0.0_rk) allocate (self%ptc%itype(self%numbers(4)), source=0) allocate (z(self%numbers(4), 2), source=0.0_rk) ! if (cfg_obj%is_go_on) then ! istat = h5pt_setstep(fid, int(cfg_obj%go_on_step, 8)) ! else istat = h5pt_setstep(fid, 1_8) ! end if istat = h5pt_readdata(fid, "X", z(:, 1)) istat = h5pt_readdata(fid, "Y", z(:, 2)) do i = 1, self%numbers(4) self%ptc%loc(:, i) = z(i, :) end do istat = h5pt_readdata(fid, "Vx", z(:, 1)) istat = h5pt_readdata(fid, "Vy", z(:, 2)) do i = 1, self%numbers(4) self%ptc%vel(:, i) = z(i, :) end do self%ptc%acc = 0.0_rk istat = h5pt_readdata(fid, "Mass", self%ptc%mass) istat = h5pt_readdata(fid, "Rho", self%ptc%rho) istat = h5pt_readdata(fid, "P", self%ptc%p) istat = h5pt_readdata(fid, "U", self%ptc%u) istat = h5pt_readdata(fid, "Itype", self%ptc%itype) istat = h5pt_close(fid) end subroutine read !> 打开 subroutine open (self) class(pif_h5part_type), intent(inout) :: self integer(8) :: istat allocate (self%z(self%numbers(4)), source=0.0_rk) self%istep = 0_8 self%fid = h5pt_openw(cli_obj%working_directory.join.ioc_obj%output_file) if (self%fid == 0_8) then call terminal_obj%error('open pif--out h5part file failed') end if istat = h5pt_writefileattrib(self%fid, "DateTime", nowtime()) istat = h5pt_writefileattrib(self%fid, "Creator", "SPX") istat = h5pt_writefileattrib(self%fid, "Project", "default") istat = h5pt_writefileattrib(self%fid, "Nreal, Nvirt, Ntotal, Dimension", & [self%numbers(1), self%numbers(3), self%numbers(4), 2], 4_8) istat = h5pt_writefileattrib(self%fid, "HSML", [self%ptc%hsml], 1_8) istat = h5pt_setnpoints(self%fid, int(self%numbers(4), 8)) open (newunit=self%fid_stat, file=cli_obj%working_directory.join."STAT.csv", action='write') write (self%fid_stat, "(*(A10,:,','))") "Time/s", "X_MIN/m", "Y_MIN/m", "X_MAX/m", "Y_MAX/m" if (pif_obj%numbers(5) /= 0) then open (newunit=self%fid_float, file=cli_obj%working_directory.join."FLOAT.csv", action='write') write (self%fid_float, "(*(A10,:,','))") "Time/s", "X/m", "Y/m", "Fx/N", "Fy/N", "M/N*m" end if end subroutine open !> 关闭 subroutine close (self) class(pif_h5part_type), intent(inout) :: self integer(8) :: istat istat = h5pt_close(self%fid) close (self%fid_stat) if (pif_obj%numbers(5) /= 0) then close (self%fid_float) end if end subroutine close !> 写入 subroutine write (self, t) class(pif_h5part_type), intent(inout) :: self real(rk), intent(in) :: t integer(8) :: istat pif_obj%ptc%acc(:, pif_obj%nreal + 1:) = 0.0_rk ! 边界粒子加速度清零输出, TOCHECK: 是否需要为浮体粒子加速度清零 self%istep = self%istep + 1_8 istat = h5pt_setstep(self%fid, self%istep) istat = h5pt_writestepattrib(self%fid, "RealTime", [t], 1_8) istat = h5pt_writedata(self%fid, "X", self%ptc%loc(1, :)) istat = h5pt_writedata(self%fid, "Y", self%ptc%loc(2, :)) istat = h5pt_writedata(self%fid, "Z", self%z) istat = h5pt_writedata(self%fid, "Vx", self%ptc%vel(1, :)) istat = h5pt_writedata(self%fid, "Vy", self%ptc%vel(2, :)) istat = h5pt_writedata(self%fid, "Fx", self%ptc%acc(1, :)) istat = h5pt_writedata(self%fid, "Fy", self%ptc%acc(2, :)) istat = h5pt_writedata(self%fid, "Rho", self%ptc%rho) istat = h5pt_writedata(self%fid, "P", self%ptc%p) istat = h5pt_writedata(self%fid, "Mass", self%ptc%mass) istat = h5pt_writedata(self%fid, "Itype", self%ptc%itype) write (self%fid_stat, "(*(ES10.3,:,','))") t, stat%loc_min, stat%loc_max if (pif_obj%numbers(5) /= 0) then write (self%fid_float, "(*(ES10.3,:,','))") t, float_obj(1)%loc, & float_obj(1)%acc, float_obj(1)%ang_acc end if end subroutine write end module spx_pif_h5part