From 2b2027d1f8d1a212c6d1d3d61fbb4cd43faad1d9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 28 Aug 2021 15:34:18 -0400 Subject: [PATCH] Fixed issues related to gaps in the id record that were causing the processing code to fail because of empty string values if idslots were not filled between writes. Now particle info is dumped to the bin file in NetCDF mode after a rearray operation to catch all newly created bodies and fill all id slots. --- src/io/io.f90 | 70 ++++++++++--------- src/modules/swiftest_classes.f90 | 18 +++-- src/modules/symba_classes.f90 | 2 +- src/netcdf/netcdf.f90 | 111 ++++++++++++++++++++++--------- src/setup/setup.f90 | 11 +-- src/symba/symba_collision.f90 | 1 - src/symba/symba_util.f90 | 40 ++++++----- 7 files changed, 154 insertions(+), 99 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index b3cd08b59..c3417f094 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -145,7 +145,7 @@ module subroutine io_dump_particle_info_base(self, param, idx) implicit none ! Arguments class(swiftest_base), intent(inout) :: self !! Swiftest base object (can be cb, pl, or tp) - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters integer(I4B), dimension(:), optional, intent(in) :: idx !! Array of test particle indices to append to the particle file ! Internals @@ -154,41 +154,47 @@ module subroutine io_dump_particle_info_base(self, param, idx) integer(I4B) :: i character(STRMAX) :: errmsg - if (lfirst) then - select case(param%out_stat) - case('APPEND') - open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - case('NEW', 'UNKNOWN', 'REPLACE') - open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', err = 667, iomsg = errmsg) - case default - write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) - call util_exit(FAILURE) - end select - - lfirst = .false. - else - open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - end if + if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then + if (lfirst) then + select case(param%out_stat) + case('APPEND') + open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + case('NEW', 'UNKNOWN', 'REPLACE') + open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', err = 667, iomsg = errmsg) + case default + write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) + call util_exit(FAILURE) + end select - select type(self) - class is (swiftest_cb) - write(LUN, err = 667, iomsg = errmsg) self%id - call self%info%dump(LUN) - class is (swiftest_body) - if (present(idx)) then - do i = 1, size(idx) - write(LUN, err = 667, iomsg = errmsg) self%id(idx(i)) - call self%info(idx(i))%dump(LUN) - end do + lfirst = .false. else - do i = 1, self%nbody - write(LUN, err = 667, iomsg = errmsg) self%id(i) - call self%info(i)%dump(LUN) - end do + open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) end if - end select - close(unit = LUN, err = 667, iomsg = errmsg) + select type(self) + class is (swiftest_cb) + write(LUN, err = 667, iomsg = errmsg) self%id + call self%info%dump(LUN) + class is (swiftest_body) + if (present(idx)) then + do i = 1, size(idx) + write(LUN, err = 667, iomsg = errmsg) self%id(idx(i)) + call self%info(idx(i))%dump(LUN) + end do + else + do i = 1, self%nbody + write(LUN, err = 667, iomsg = errmsg) self%id(i) + call self%info(i)%dump(LUN) + end do + end if + end select + + close(unit = LUN, err = 667, iomsg = errmsg) + else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + call param%nciu%open(param) + call self%write_particle_info(param%nciu) + call param%nciu%close(param) + end if return diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 45ea97f15..75b792b57 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -179,10 +179,12 @@ module swiftest_classes !! An abstract superclass for a generic Swiftest object contains !! The minimal methods that all systems must have - procedure :: dump => io_dump_base !! Dump contents to file - procedure :: dump_particle_info => io_dump_particle_info_base !! Dump contents of particle information metadata to file - procedure :: write_frame_netcdf => netcdf_write_frame_base !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format - generic :: write_frame => write_frame_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments + procedure :: dump => io_dump_base !! Dump contents to file + procedure :: dump_particle_info => io_dump_particle_info_base !! Dump contents of particle information metadata to file + procedure :: write_frame_netcdf => netcdf_write_frame_base !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format + procedure :: write_particle_info_netcdf => netcdf_write_particle_info_base !! Writes out the particle information metadata to NetCDF file + generic :: write_frame => write_frame_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments + generic :: write_particle_info => write_particle_info_netcdf end type swiftest_base !******************************************************************************************************************************** @@ -652,7 +654,7 @@ end subroutine io_dump_param module subroutine io_dump_particle_info_base(self, param, idx) implicit none class(swiftest_base), intent(inout) :: self !! Swiftest base object (can be cb, pl, or tp) - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters integer(I4B), dimension(:), optional, intent(in) :: idx !! Array of test particle indices to append to the particle file end subroutine io_dump_particle_info_base @@ -913,6 +915,12 @@ module subroutine netcdf_write_hdr_system(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine netcdf_write_hdr_system + module subroutine netcdf_write_particle_info_base(self, iu) + implicit none + class(swiftest_base), intent(in) :: self !! Swiftest particle object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to identify a particular NetCDF dataset + end subroutine netcdf_write_particle_info_base + module subroutine obl_acc_body(self, system) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 3a9e37fce..60dbd4feb 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -571,7 +571,7 @@ module subroutine symba_util_rearray_pl(self, system, param) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions end subroutine symba_util_rearray_pl end interface diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index d2f7a0f80..e37703782 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -238,26 +238,22 @@ module subroutine netcdf_write_frame_base(self, iu, param) integer(I4B), dimension(:), allocatable :: ind character(len=:), allocatable :: charstring + call self%write_particle_info(iu) + tslot = int(param%ioutput, kind=I4B) + 1 select type(self) class is (swiftest_body) associate(n => self%nbody) if (n == 0) return + + allocate(ind(n)) call util_sort(self%id(1:n), ind) do i = 1, n j = ind(i) idslot = self%id(j) + 1 - call check( nf90_put_var(iu%ncid, iu%id_varid, self%id(j), start=[idslot]) ) - charstring = trim(adjustl(self%info(j)%name)) - strlen = len(charstring) - call check( nf90_put_var(iu%ncid, iu%name_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) - - charstring = trim(adjustl(self%info(j)%particle_type)) - strlen = len(charstring) - call check( nf90_put_var(iu%ncid, iu%ptype_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) if ((param%out_form == XV) .or. (param%out_form == XVEL)) then call check( nf90_put_var(iu%ncid, iu%xhx_varid, self%xh(1, j), start=[idslot, tslot]) ) @@ -277,19 +273,6 @@ module subroutine netcdf_write_frame_base(self, iu, param) call check( nf90_put_var(iu%ncid, iu%capm_varid, self%capm(j), start=[idslot, tslot]) ) end if - if (iu%ltrack_origin) then - charstring = trim(adjustl(self%info(j)%origin_type)) - strlen = len(charstring) - call check( nf90_put_var(iu%ncid, iu%origin_type_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) - call check( nf90_put_var(iu%ncid, iu%origin_time_varid, self%info(j)%origin_time, start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_xhx_varid, self%info(j)%origin_xh(1), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_xhy_varid, self%info(j)%origin_xh(2), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_xhz_varid, self%info(j)%origin_xh(3), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_vhx_varid, self%info(j)%origin_vh(1), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_vhy_varid, self%info(j)%origin_vh(2), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_vhz_varid, self%info(j)%origin_vh(3), start=[idslot]) ) - end if - select type(self) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]) ) @@ -342,23 +325,85 @@ module subroutine netcdf_write_frame_base(self, iu, param) call check( nf90_put_var(iu%ncid, iu%Q_varid, self%Q, start=[idslot, tslot]) ) end if - if (iu%ltrack_origin) then - charstring = trim(adjustl(self%info%origin_type)) - strlen = len(charstring) - call check( nf90_put_var(iu%ncid, iu%origin_type_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) - call check( nf90_put_var(iu%ncid, iu%origin_time_varid, self%info%origin_time, start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_xhx_varid, self%info%origin_xh(1), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_xhy_varid, self%info%origin_xh(2), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_xhz_varid, self%info%origin_xh(3), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_vhx_varid, self%info%origin_vh(1), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_vhy_varid, self%info%origin_vh(2), start=[idslot]) ) - call check( nf90_put_var(iu%ncid, iu%origin_vhz_varid, self%info%origin_vh(3), start=[idslot]) ) - end if end select return end subroutine netcdf_write_frame_base + + module subroutine netcdf_write_particle_info_base(self, iu) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Write all current particle to file + implicit none + ! Arguments + class(swiftest_base), intent(in) :: self !! Swiftest particle object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to identify a particular NetCDF dataset + ! Internals + integer(I4B) :: i, j, tslot, strlen, idslot + integer(I4B), dimension(:), allocatable :: ind + character(len=:), allocatable :: charstring + + + select type(self) + class is (swiftest_body) + associate(n => self%nbody) + if (n == 0) return + allocate(ind(n)) + call util_sort(self%id(1:n), ind) + + do i = 1, n + j = ind(i) + idslot = self%id(j) + 1 + call check( nf90_put_var(iu%ncid, iu%id_varid, self%id(j), start=[idslot]) ) + charstring = trim(adjustl(self%info(j)%name)) + strlen = len(charstring) + call check( nf90_put_var(iu%ncid, iu%name_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + + charstring = trim(adjustl(self%info(j)%particle_type)) + strlen = len(charstring) + call check( nf90_put_var(iu%ncid, iu%ptype_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + + charstring = trim(adjustl(self%info(j)%origin_type)) + strlen = len(charstring) + call check( nf90_put_var(iu%ncid, iu%origin_type_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + call check( nf90_put_var(iu%ncid, iu%origin_time_varid, self%info(j)%origin_time, start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_xhx_varid, self%info(j)%origin_xh(1), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_xhy_varid, self%info(j)%origin_xh(2), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_xhz_varid, self%info(j)%origin_xh(3), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_vhx_varid, self%info(j)%origin_vh(1), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_vhy_varid, self%info(j)%origin_vh(2), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_vhz_varid, self%info(j)%origin_vh(3), start=[idslot]) ) + end do + end associate + + class is (swiftest_cb) + idslot = self%id + 1 + call check( nf90_put_var(iu%ncid, iu%id_varid, self%id, start=[idslot]) ) + + charstring = trim(adjustl(self%info%name)) + strlen = len(charstring) + call check( nf90_put_var(iu%ncid, iu%name_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + + charstring = trim(adjustl(self%info%particle_type)) + strlen = len(charstring) + call check( nf90_put_var(iu%ncid, iu%ptype_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + + charstring = trim(adjustl(self%info%origin_type)) + strlen = len(charstring) + call check( nf90_put_var(iu%ncid, iu%origin_type_varid, charstring, start=[1, idslot], count=[strlen, 1]) ) + call check( nf90_put_var(iu%ncid, iu%origin_time_varid, self%info%origin_time, start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_xhx_varid, self%info%origin_xh(1), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_xhy_varid, self%info%origin_xh(2), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_xhz_varid, self%info%origin_xh(3), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_vhx_varid, self%info%origin_vh(1), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_vhy_varid, self%info%origin_vh(2), start=[idslot]) ) + call check( nf90_put_var(iu%ncid, iu%origin_vhz_varid, self%info%origin_vh(3), start=[idslot]) ) + end select + + return + end subroutine netcdf_write_particle_info_base + module subroutine netcdf_write_hdr_system(self, iu, param) !! author: David A. Minton diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 342818672..11bd69fd1 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -175,15 +175,8 @@ module subroutine setup_initialize_particle_info_system(self, param) end if cb%info%particle_type = CB_TYPE_NAME - call cb%dump_particle_info(param) - if (npl > 0) then - pl%info(1:npl)%particle_type = PL_TYPE_NAME - call pl%dump_particle_info(param) - end if - if (ntp > 0) then - tp%info(1:ntp)%particle_type = TP_TYPE_NAME - call tp%dump_particle_info(param) - end if + if (npl > 0) pl%info(1:npl)%particle_type = PL_TYPE_NAME + if (ntp > 0) tp%info(1:ntp)%particle_type = TP_TYPE_NAME end associate diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index cf4fdb131..7a68a56cd 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -907,7 +907,6 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, ! Copy over identification, information, and physical properties of the new bodies from the fragment list plnew%id(1:nfrag) = id_frag(1:nfrag) - param%maxid = param%maxid + nfrag plnew%xb(:, 1:nfrag) = xb_frag(:, 1:nfrag) plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) do i = 1, nfrag diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 51af8159d..150bf7825 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -291,15 +291,6 @@ module subroutine symba_util_index_eucl_plpl(self, param) npl = int(self%nbody, kind=I8B) call pl%sort("mass", ascending=.false.) - select type(param) - class is (symba_parameters) - pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY - where(pl%lmtiny(1:npl)) - pl%info(1:npl)%particle_type = PL_TINY_TYPE_NAME - elsewhere - pl%info(1:npl)%particle_type = PL_TYPE_NAME - end where - end select nplm = count(.not. pl%lmtiny(1:npl)) pl%nplm = int(nplm, kind=I4B) @@ -415,11 +406,11 @@ module subroutine symba_util_rearray_pl(self, system, param) ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(symba_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals class(symba_pl), allocatable :: tmp !! The discarded body list. integer(I4B) :: i, j, k, npl, nencmin - logical, dimension(:), allocatable :: lmask + logical, dimension(:), allocatable :: lmask, ldump_mask class(symba_plplenc), allocatable :: plplenc_old logical :: lencounter integer(I4B), dimension(:), allocatable :: levelg_orig_pl, levelm_orig_pl, levelg_orig_tp, levelm_orig_tp, nplenc_orig_pl, nplenc_orig_tp, ntpenc_orig_pl @@ -450,16 +441,16 @@ module subroutine symba_util_rearray_pl(self, system, param) ! Append the adds to the main pl object call pl%append(pl_adds, lsource_mask=[(.true., i=1, nadd)]) - allocate(lmask(npl+nadd)) - lmask(1:npl) = .false. - lmask(npl+1:npl+nadd) = pl%status(npl+1:npl+nadd) == NEW_PARTICLE - + allocate(ldump_mask(npl+nadd)) ! This mask is used only to append the original Fortran binary particle.dat file with new bodies. This is ignored for NetCDF output + ldump_mask(1:npl) = .false. + ldump_mask(npl+1:npl+nadd) = pl%status(npl+1:npl+nadd) == NEW_PARTICLE npl = pl%nbody - call pl%dump_particle_info(param, idx=pack([(i, i=1, npl)], lmask)) - - deallocate(lmask) + else + allocate(ldump_mask(npl)) + ldump_mask(:) = .false. end if + ! Reset all of the status flags for this body where(pl%status(1:npl) /= INACTIVE) pl%status(1:npl) = ACTIVE @@ -471,6 +462,19 @@ module subroutine symba_util_rearray_pl(self, system, param) pl%lmask(1:npl) = .false. end where + select type(param) + class is (symba_parameters) + pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY + where(pl%lmtiny(1:npl)) + pl%info(1:npl)%particle_type = PL_TINY_TYPE_NAME + elsewhere + pl%info(1:npl)%particle_type = PL_TYPE_NAME + end where + end select + + call pl%dump_particle_info(param, idx=pack([(i, i=1, npl)], ldump_mask)) + deallocate(ldump_mask) + ! Reindex the new list of bodies call pl%index(param)