From d09d8c5e90b68a85e041fe72334a0292736ca076 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Dec 2022 11:46:09 -0500 Subject: [PATCH] Made a number of improvements to the handling of encounters. time and id dimensions are now limited, which significantly reduces the file sizes of the individual encounter data files --- src/modules/symba_classes.f90 | 12 +++--- src/symba/symba_io.f90 | 79 ++++++++++++++++++++++++----------- src/symba/symba_util.f90 | 11 +++-- 3 files changed, 66 insertions(+), 36 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 68b94c373..6c3aaae12 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -182,12 +182,12 @@ module symba_classes !! NetCDF dimension and variable names for the enounter save object type, extends(netcdf_parameters) :: symba_io_encounter_parameters - integer(I4B) :: COLLIDER_DIM_SIZE = 2 !! Size of collider dimension - integer(I4B) :: ienc_frame = 1 !! Current frame number for the encounter history - character(STRMAX) :: enc_file = "encounter.nc" !! Encounter output file name - - character(NAMELEN) :: level_varname = "level" !! Recursion depth + integer(I4B) :: ienc_frame = 1 !! Current frame number for the encounter history + character(STRMAX) :: enc_file = "encounter.nc" !! Encounter output file name + character(NAMELEN) :: level_varname = "level" !! Recursion depth integer(I4B) :: level_varid !! ID for the recursion level variable + integer(I4B) :: time_dimsize = 0 !! Number of time values in snapshot + integer(I4B) :: id_dimsize = 0 !! Number of potential id values in snapshot contains procedure :: initialize => symba_io_encounter_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object end type symba_io_encounter_parameters @@ -230,7 +230,7 @@ module symba_classes type, extends(symba_nbody_system) :: symba_encounter_snapshot - integer(I4B) :: tslot !! The index for the time array in the final NetCDF file + integer(I4B) :: tslot !! The index for the time array in the final NetCDF file contains procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file generic :: write_frame => write_encounter_frame diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 0fceaa724..bc0f1b985 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -31,6 +31,8 @@ module subroutine symba_io_encounter_dump(self, param) self%nc%ienc_frame = i call snapshot%write_frame(self%nc,param) end select + else + exit end if end do @@ -57,7 +59,6 @@ module subroutine symba_io_encounter_initialize_output(self, param) character(len=STRMAX) :: errmsg integer(I4B) :: ndims - associate(nc => self) dfill = ieee_value(dfill, IEEE_QUIET_NAN) sfill = ieee_value(sfill, IEEE_QUIET_NAN) @@ -69,7 +70,6 @@ module subroutine symba_io_encounter_initialize_output(self, param) self%out_type = NF90_DOUBLE end select - ! Check if the file exists, and if it does, delete it inquire(file=nc%enc_file, exist=fileExists) if (fileExists) then @@ -80,9 +80,9 @@ module subroutine symba_io_encounter_initialize_output(self, param) call check( nf90_create(nc%enc_file, NF90_NETCDF4, nc%id), "symba_io_encounter_initialize_output nf90_create" ) ! Dimensions - call check( nf90_def_dim(nc%id, nc%time_dimname, NF90_UNLIMITED, nc%time_dimid), "symba_io_encounter_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension - call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "symba_io_encounter_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call check( nf90_def_dim(nc%id, nc%id_dimname, NF90_UNLIMITED, nc%id_dimid), "symba_io_encounter_initialize_output nf90_def_dim id_dimid" ) ! dimension to store particle id numbers + call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "symba_io_encounter_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension + call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM , nc%space_dimid), "symba_io_encounter_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call check( nf90_def_dim(nc%id, nc%id_dimname, param%maxid, nc%id_dimid), "symba_io_encounter_initialize_output nf90_def_dim id_dimid" ) ! dimension to store particle id numbers call check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "symba_io_encounter_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) ! Dimension coordinates @@ -141,37 +141,58 @@ module subroutine symba_io_encounter_write_frame(self, nc, param) use netcdf implicit none ! Arguments - class(symba_encounter_snapshot), intent(in) :: self !! Swiftest encounter structure - class(symba_io_encounter_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(symba_encounter_snapshot), intent(in) :: self !! Swiftest encounter structure + class(symba_io_encounter_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i, tslot, idslot, old_mode, n - character(len=NAMELEN) :: charstring + integer(I4B) :: i, tslot, idslot, old_mode, npl, ntp + character(len=NAMELEN) :: charstring - tslot = self%tslot call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "symba_io_encounter_write_frame nf90_set_fill" ) + + tslot = self%tslot + call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "symba_io_encounter_write_frame nf90_put_var time_varid" ) + select type(pl => self%pl) class is (symba_pl) - n = size(pl%id(:)) - do i = 1, n + npl = pl%nbody + do i = 1, npl idslot = pl%id(i) - call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "symba_io_encounter_write_frame nf90_put_var time_varid" ) - call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "symba_io_encounter_write_frame nf90_put_var id_varid" ) - call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var rh_varid" ) - call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var vh_varid" ) - call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var body Gmass_varid" ) - if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var body radius_varid" ) + call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "symba_io_encounter_write_frame nf90_put_var pl id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl vh_varid" ) + call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl Gmass_varid" ) + + if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl radius_varid" ) + if (param%lrotation) then - call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var body Ip_varid" ) - call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var body rotx_varid" ) + call check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl Ip_varid" ) + call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1,idslot, tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var pl rotx_varid" ) end if + charstring = trim(adjustl(pl%info(i)%name)) - call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var name_varid" ) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var pl name_varid" ) charstring = trim(adjustl(pl%info(i)%particle_type)) - call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var particle_type_varid" ) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var pl particle_type_varid" ) end do + end select + associate(tp => self%tp) + ntp = tp%nbody + do i = 1, ntp + idslot = tp%id(i) + call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "symba_io_encounter_write_frame nf90_put_var tp id_varid" ) + call check( nf90_put_var(nc%id, nc%rh_varid, tp%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var tp rh_varid" ) + call check( nf90_put_var(nc%id, nc%vh_varid, tp%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "symba_io_encounter_write_frame nf90_put_var tp vh_varid" ) + + charstring = trim(adjustl(tp%info(i)%name)) + call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var tp name_varid" ) + charstring = trim(adjustl(tp%info(i)%particle_type)) + call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var tp particle_type_varid" ) + end do + end associate + call check( nf90_set_fill(nc%id, old_mode, old_mode) ) return @@ -388,10 +409,20 @@ module subroutine symba_io_stop_encounter(self, param, t) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time ! Internals - !character(STRMAX) + integer(I4B) :: i ! Create and save the output file for this encounter + + ! Figure out how many time slots we need + do i = 1, self%encounter_history%nframes + if (self%t + param%dt <= self%encounter_history%tvals(i)) then + self%encounter_history%nc%time_dimsize = i + exit + end if + end do + write(self%encounter_history%nc%enc_file, '("encounter_",I0.6,".nc")') param%iloop + call self%encounter_history%nc%initialize(param) call self%encounter_history%dump(param) call self%encounter_history%nc%close() diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 20d295b9c..b96270fc9 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -1326,9 +1326,10 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) associate(npl => self%pl%nbody, ntp => self%tp%nbody) allocate(symba_encounter_snapshot :: snapshot) + snapshot%t = t - if (npl > 0) allocate(symba_pl :: snapshot%pl) - if (ntp > 0) allocate(symba_tp :: snapshot%tp) + allocate(symba_pl :: snapshot%pl) + allocate(symba_tp :: snapshot%tp) if (npl + ntp == 0) return npl_snap = npl ntp_snap = ntp @@ -1345,6 +1346,7 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == self%irec ntp_snap = count(tp%lmask(1:ntp)) end if + snapshot%pl%nbody = npl_snap ! Take snapshot of the currently encountering massive bodies if (npl_snap > 0) then @@ -1377,6 +1379,7 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) end if ! Take snapshot of the currently encountering test particles + snapshot%tp%nbody = ntp_snap if (ntp_snap > 0) then allocate(snapshot%tp%id(ntp_snap)) allocate(snapshot%tp%info(ntp_snap)) @@ -1390,10 +1393,6 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) end do end if - if (npl_snap + ntp_snap == 0) return - - snapshot%t = t - ! Save the snapshot self%encounter_history%iframe = self%encounter_history%iframe + 1 call self%resize_storage(self%encounter_history%iframe)