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)