From 1b82c228d5f5829c7b83e0e9fd399743fd2cfb55 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 12 Dec 2022 12:39:48 -0500 Subject: [PATCH] Lots of fixes to indexing of encounters and collisions --- src/encounter/encounter_io.f90 | 160 ++++++++++++------------ src/encounter/encounter_util.f90 | 144 ++++++++++++++++----- src/fraggle/fraggle_io.f90 | 200 +++++++++++++++--------------- src/fraggle/fraggle_util.f90 | 36 +++++- src/modules/encounter_classes.f90 | 17 ++- src/modules/fraggle_classes.f90 | 15 ++- src/modules/swiftest_classes.f90 | 7 ++ src/util/util_index.f90 | 130 +++++++++++++------ 8 files changed, 447 insertions(+), 262 deletions(-) diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index bc3bfaac2..09c2dd676 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -27,15 +27,17 @@ module subroutine encounter_io_dump_collision(self, param) class is (fraggle_io_parameters) if (self%iframe > 0) then nc%file_number = nc%file_number + 1 - nc%event_dimsize = self%iframe call self%make_index_map() + nc%event_dimsize = self%nt + nc%id_dimsize = self%nid + write(nc%file_name, '("collision_",I0.6,".nc")') nc%file_number call nc%initialize(param) do i = 1, self%nframes if (allocated(self%frame(i)%item)) then select type(snapshot => self%frame(i)%item) - class is (fraggle_collision_snapshot) + class is (fraggle_snapshot) param%ioutput = i call snapshot%write_frame(nc,param) end select @@ -70,6 +72,8 @@ module subroutine encounter_io_dump_encounter(self, param) ! Create and save the output files for this encounter and fragmentation nc%file_number = nc%file_number + 1 call self%make_index_map() + nc%time_dimsize = self%nt + nc%id_dimsize = self%nid write(nc%file_name, '("encounter_",I0.6,".nc")') nc%file_number call nc%initialize(param) @@ -113,77 +117,72 @@ module subroutine encounter_io_initialize(self, param) integer(I4B) :: ndims, i associate(nc => self) - select type(param) - class is (symba_parameters) - dfill = ieee_value(dfill, IEEE_QUIET_NAN) - sfill = ieee_value(sfill, IEEE_QUIET_NAN) - - select case (param%out_type) - case("NETCDF_FLOAT") - self%out_type = NF90_FLOAT - case("NETCDF_DOUBLE") - self%out_type = NF90_DOUBLE - end select + dfill = ieee_value(dfill, IEEE_QUIET_NAN) + sfill = ieee_value(sfill, IEEE_QUIET_NAN) + + select case (param%out_type) + case("NETCDF_FLOAT") + self%out_type = NF90_FLOAT + case("NETCDF_DOUBLE") + self%out_type = NF90_DOUBLE + end select - ! Check if the file exists, and if it does, delete it - inquire(file=nc%file_name, exist=fileExists) - if (fileExists) then - open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg) - close(unit=LUN, status="delete") - end if - - call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "encounter_io_initialize nf90_create" ) - - ! Dimensions - call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize nf90_def_dim time_dimid" ) ! Simulation time dimension - call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM , nc%space_dimid), "encounter_io_initialize nf90_def_dim space_dimid" ) ! 3D space dimension - call check( nf90_def_dim(nc%id, nc%id_dimname, param%encounter_history%nid, nc%id_dimid), "encounter_io_initialize 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), "encounter_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - - ! Dimension coordinates - call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize nf90_def_var time_varid" ) - call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_io_initialize nf90_def_var space_varid" ) - call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "encounter_io_initialize nf90_def_var id_varid" ) - - ! Variables - call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%name_varid), "encounter_io_initialize nf90_def_var name_varid" ) - call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "encounter_io_initialize nf90_def_var ptype_varid" ) - call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_initialize nf90_def_var rh_varid" ) - call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "encounter_io_initialize nf90_def_var vh_varid" ) - call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize nf90_def_var Gmass_varid" ) - call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, [nc%time_dimid], nc%loop_varid), "encounter_io_initialize nf90_def_var loop_varid" ) - if (param%lclose) then - call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize nf90_def_var radius_varid" ) - end if - if (param%lrotation) then - call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%Ip_varid), "encounter_io_initialize nf90_def_var Ip_varid" ) - call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rot_varid), "encounter_io_initialize nf90_def_var rot_varid" ) - end if - - call check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize nf90_inquire nVariables" ) - do varid = 1, nvar - call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize nf90_inquire_variable" ) - select case(vartype) - case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" ) - case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" ) - case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" ) - case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" ) - end select - end do + ! Check if the file exists, and if it does, delete it + inquire(file=nc%file_name, exist=fileExists) + if (fileExists) then + open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg) + close(unit=LUN, status="delete") + end if - ! Take the file out of define mode - call check( nf90_enddef(nc%id), "encounter_io_initialize nf90_enddef" ) + call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "encounter_io_initialize nf90_create" ) + + ! Dimensions + call check( nf90_def_dim(nc%id, nc%time_dimname, nc%time_dimsize, nc%time_dimid), "encounter_io_initialize nf90_def_dim time_dimid" ) ! Simulation time dimension + call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM , nc%space_dimid), "encounter_io_initialize nf90_def_dim space_dimid" ) ! 3D space dimension + call check( nf90_def_dim(nc%id, nc%id_dimname, nc%id_dimsize, nc%id_dimid), "encounter_io_initialize 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), "encounter_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + + ! Dimension coordinates + call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), "encounter_io_initialize nf90_def_var time_varid" ) + call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "encounter_io_initialize nf90_def_var space_varid" ) + call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "encounter_io_initialize nf90_def_var id_varid" ) + + ! Variables + call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%name_varid), "encounter_io_initialize nf90_def_var name_varid" ) + call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "encounter_io_initialize nf90_def_var ptype_varid" ) + call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_initialize nf90_def_var rh_varid" ) + call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "encounter_io_initialize nf90_def_var vh_varid" ) + call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize nf90_def_var Gmass_varid" ) + call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, [nc%time_dimid], nc%loop_varid), "encounter_io_initialize nf90_def_var loop_varid" ) + if (param%lclose) then + call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize nf90_def_var radius_varid" ) + end if + if (param%lrotation) then + call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%Ip_varid), "encounter_io_initialize nf90_def_var Ip_varid" ) + call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rot_varid), "encounter_io_initialize nf90_def_var rot_varid" ) + end if - ! Add in the space dimension coordinates - call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize nf90_put_var space" ) + call check( nf90_inquire(nc%id, nVariables=nvar), "encounter_io_initialize nf90_inquire nVariables" ) + do varid = 1, nvar + call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "encounter_io_initialize nf90_inquire_variable" ) + select case(vartype) + case(NF90_INT) + call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "encounter_io_initialize nf90_def_var_fill NF90_INT" ) + case(NF90_FLOAT) + call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "encounter_io_initialize nf90_def_var_fill NF90_FLOAT" ) + case(NF90_DOUBLE) + call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "encounter_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + case(NF90_CHAR) + call check( nf90_def_var_fill(nc%id, varid, 0, 0), "encounter_io_initialize nf90_def_var_fill NF90_CHAR" ) + end select + end do + + ! Take the file out of define mode + call check( nf90_enddef(nc%id), "encounter_io_initialize nf90_enddef" ) + + ! Add in the space dimension coordinates + call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize nf90_put_var space" ) - ! Pre-fill name slots with ids - !call check( nf90_put_var(nc%id, nc%id_varid, [(-1,i=1,param%maxid+1)], start=[1], count=[param%maxid+1]), "encounter_io_initialize nf90_put_var pl id_varid" ) - end select end associate return @@ -205,16 +204,14 @@ module subroutine encounter_io_write_frame(self, nc, param) class(netcdf_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, npl, ntp + integer(I4B) :: i, idslot, old_mode, npl, ntp character(len=:), allocatable :: charstring - tslot = param%ioutput - associate(pl => self%pl, tp => self%tp) - select type (nc) - class is (encounter_io_parameters) - select type (param) - class is (symba_parameters) - + select type (nc) + class is (encounter_io_parameters) + select type (param) + class is (symba_parameters) + associate(pl => self%pl, tp => self%tp, encounter_history => param%encounter_history, tslot => param%ioutput) call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame nf90_set_fill" ) call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" ) @@ -222,8 +219,7 @@ module subroutine encounter_io_write_frame(self, nc, param) npl = pl%nbody do i = 1, npl - idslot = pl%id(i) + 1 - idslot = findloc(param%encounter_history%idvals,pl%id(i),dim=1) + idslot = findloc(encounter_history%idvals,pl%id(i),dim=1) call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_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]), "encounter_io_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]), "encounter_io_write_frame nf90_put_var pl vh_varid" ) @@ -244,7 +240,7 @@ module subroutine encounter_io_write_frame(self, nc, param) ntp = tp%nbody do i = 1, ntp - idslot = tp%id(i) + 1 + idslot = findloc(param%encounter_history%idvals,tp%id(i),dim=1) call check( nf90_put_var(nc%id, nc%id_varid, tp%id(i), start=[idslot]), "encounter_io_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]), "encounter_io_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]), "encounter_io_write_frame nf90_put_var tp vh_varid" ) @@ -256,9 +252,9 @@ module subroutine encounter_io_write_frame(self, nc, param) end do call check( nf90_set_fill(nc%id, old_mode, old_mode) ) - end select + end associate end select - end associate + end select return end subroutine encounter_io_write_frame diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 5c4212253..3a9f4b062 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -180,6 +180,97 @@ module subroutine encounter_util_final_storage(self) end subroutine encounter_util_final_storage + module subroutine encounter_util_get_idvalues_snapshot(self, idvals) + !! author: David A. Minton + !! + !! Returns an array of all id values saved in this snapshot + implicit none + ! Arguments + class(encounter_snapshot), intent(in) :: self !! Encounter snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + ! Internals + integer(I4B) :: npl, ntp + + if (allocated(self%pl)) then + npl = self%pl%nbody + else + npl = 0 + end if + if (allocated(self%tp)) then + ntp = self%tp%nbody + else + ntp = 0 + end if + + if (npl + ntp == 0) return + allocate(idvals(npl+ntp)) + + if (npl > 0) idvals(1:npl) = self%pl%id(:) + if (ntp >0) idvals(npl+1:npl+ntp) = self%tp%id(:) + + return + + end subroutine encounter_util_get_idvalues_snapshot + + + subroutine encounter_util_get_vals_storage(storage, idvals, tvals) + !! author: David A. Minton + !! + !! Gets the id values in a storage object, regardless of whether it is encounter of collision + ! Argument + class(swiftest_storage(*)), intent(in) :: storage !! Swiftest storage object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values in all snapshots + real(DP), dimension(:), allocatable, intent(out) :: tvals !! Array of all time values in all snapshots + ! Internals + integer(I4B) :: i, n, nlo, nhi, ntotal + integer(I4B), dimension(:), allocatable :: itmp + + associate(nsnaps => storage%iframe) + + allocate(tvals(nsnaps)) + + tvals(:) = 0.0_DP + + ! First pass to get total number of ids + ntotal = 0 + do i = 1, nsnaps + if (allocated(storage%frame(i)%item)) then + select type(snapshot => storage%frame(i)%item) + class is (encounter_snapshot) + tvals(i) = snapshot%t + call snapshot%get_idvals(itmp) + if (allocated(itmp)) then + n = size(itmp) + ntotal = ntotal + n + end if + end select + end if + end do + + allocate(idvals(ntotal)) + nlo = 1 + ! Second pass to store all ids get all of the ids stored + do i = 1, nsnaps + if (allocated(storage%frame(i)%item)) then + select type(snapshot => storage%frame(i)%item) + class is (encounter_snapshot) + tvals(i) = snapshot%t + call snapshot%get_idvals(itmp) + if (allocated(itmp)) then + n = size(itmp) + nhi = nlo + n - 1 + idvals(nlo:nhi) = itmp(1:n) + nlo = nhi + 1 + end if + end select + end if + end do + + end associate + return + end subroutine encounter_util_get_vals_storage + + module subroutine encounter_util_index_map_encounter(self) !! author: David A. Minton !! @@ -189,43 +280,17 @@ module subroutine encounter_util_index_map_encounter(self) ! Arguments class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object ! Internals - ! Internals - integer(I4B) :: i, n, nold, nt - integer(I4B), dimension(:), allocatable :: idvals, tmp + integer(I4B), dimension(:), allocatable :: idvals real(DP), dimension(:), allocatable :: tvals - if (self%nid == 0) return - allocate(idvals(self%nid)) - allocate(tvals(self%nframes)) - - n = 0 - nold = 1 - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) then - select type(snapshot => self%frame(i)%item) - class is (encounter_snapshot) - tvals(i) = snapshot%t - if (allocated(snapshot%pl)) then - n = n + snapshot%pl%nbody - idvals(nold:n) = snapshot%pl%id(:) - nold = n+1 - end if - if (allocated(snapshot%tp)) then - n = n + snapshot%tp%nbody - idvals(nold:n) = snapshot%tp%id(:) - nold = n+1 - end if - end select - else - nt = i-1 - exit - end if - end do + call encounter_util_get_vals_storage(self, idvals, tvals) + ! Consolidate ids to only unique values call util_unique(idvals,self%idvals,self%idmap) self%nid = size(self%idvals) - call util_unique(tvals(1:nt),self%tvals,self%tmap) + ! Consolidate time values to only unique values + call util_unique(tvals,self%tvals,self%tmap) self%nt = size(self%tvals) return @@ -239,8 +304,19 @@ module subroutine encounter_util_index_map_collision(self) !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id implicit none ! Arguments - class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object + class(collision_storage(*)), intent(inout) :: self !! Swiftest storage object ! Internals + integer(I4B), dimension(:), allocatable :: idvals + real(DP), dimension(:), allocatable :: tvals + + call encounter_util_get_vals_storage(self, idvals, tvals) + + ! Consolidate ids to only unique values + call util_unique(idvals,self%idvals,self%idmap) + self%nid = size(self%idvals) + + ! Don't consolidate time values (multiple collisions can happen in a single time step) + self%nt = size(self%tvals) return end subroutine encounter_util_index_map_collision @@ -428,7 +504,7 @@ module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) real(DP), intent(in), optional :: t !! Time of snapshot if different from system time character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) ! Arguments - class(fraggle_collision_snapshot), allocatable :: snapshot + class(fraggle_snapshot), allocatable :: snapshot type(symba_pl) :: pl character(len=:), allocatable :: stage integer(I4B) :: i,j @@ -460,7 +536,7 @@ module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) allocate(system%colliders%pl, source=pl) end associate case("after") - allocate(fraggle_collision_snapshot :: snapshot) + allocate(fraggle_snapshot :: snapshot) allocate(snapshot%colliders, source=system%colliders) allocate(snapshot%fragments, source=system%fragments) select type (param) diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index d4a8d3d9e..ff20394e7 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -30,115 +30,116 @@ module subroutine fraggle_io_initialize_output(self, param) character(len=STRMAX) :: errmsg integer(I4B) :: i, ndims - associate(nc => self) - dfill = ieee_value(dfill, IEEE_QUIET_NAN) - sfill = ieee_value(sfill, IEEE_QUIET_NAN) - - select case (param%out_type) - case("NETCDF_FLOAT") - self%out_type = NF90_FLOAT - case("NETCDF_DOUBLE") - self%out_type = NF90_DOUBLE - end select + select type(param) + class is (symba_parameters) + associate(nc => self, collision_history => param%collision_history) + dfill = ieee_value(dfill, IEEE_QUIET_NAN) + sfill = ieee_value(sfill, IEEE_QUIET_NAN) + + select case (param%out_type) + case("NETCDF_FLOAT") + self%out_type = NF90_FLOAT + case("NETCDF_DOUBLE") + self%out_type = NF90_DOUBLE + end select - ! Check if the file exists, and if it does, delete it - inquire(file=nc%file_name, exist=fileExists) - if (fileExists) then - open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg) - close(unit=LUN, status="delete") - end if - - call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "fraggle_io_initialize nf90_create" ) - - ! Dimensions - call check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "fraggle_io_initialize nf90_def_dim event_dimid" ) ! Dimension to store individual collision events - call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "fraggle_io_initialize nf90_def_dim space_dimid" ) ! 3D space dimension - call check( nf90_def_dim(nc%id, nc%id_dimname, param%maxid+1, nc%id_dimid), "fraggle_io_initialize 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), "fraggle_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - call check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "fraggle_io_initialize nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" - - ! Dimension coordinates - call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "fraggle_io_initialize nf90_def_var space_varid" ) - call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "fraggle_io_initialize nf90_def_var id_varid" ) - call check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "fraggle_io_initialize nf90_def_var stage_varid" ) - - ! Variables - call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & - nc%event_dimid, nc%time_varid), "fraggle_io_initialize nf90_def_var time_varid" ) - call check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & - [nc%str_dimid, nc%event_dimid], nc%regime_varid), "fraggle_io_initialize nf90_def_var regime_varid") - call check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & - [ nc%event_dimid], nc%Qloss_varid), "fraggle_io_initialize nf90_def_var Qloss_varid") - call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, & - [nc%str_dimid, nc%id_dimid ], nc%name_varid), "fraggle_io_initialize nf90_def_var name_varid") + ! Check if the file exists, and if it does, delete it + inquire(file=nc%file_name, exist=fileExists) + if (fileExists) then + open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg) + close(unit=LUN, status="delete") + end if + + call check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "fraggle_io_initialize nf90_create" ) + + ! Dimensions + call check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "fraggle_io_initialize nf90_def_dim event_dimid" ) ! Dimension to store individual collision events + call check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "fraggle_io_initialize nf90_def_dim space_dimid" ) ! 3D space dimension + call check( nf90_def_dim(nc%id, nc%id_dimname, nc%id_dimsize, nc%id_dimid), "fraggle_io_initialize 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), "fraggle_io_initialize nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + call check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "fraggle_io_initialize nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" + + ! Dimension coordinates + call check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "fraggle_io_initialize nf90_def_var space_varid" ) + call check( nf90_def_var(nc%id, nc%id_dimname, NF90_INT, nc%id_dimid, nc%id_varid), "fraggle_io_initialize nf90_def_var id_varid" ) + call check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "fraggle_io_initialize nf90_def_var stage_varid" ) - call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & - [nc%str_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%ptype_varid), "fraggle_io_initialize nf90_def_var ptype_varid") + ! Variables + call check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & + nc%event_dimid, nc%time_varid), "fraggle_io_initialize nf90_def_var time_varid" ) + call check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & + [nc%str_dimid, nc%event_dimid], nc%regime_varid), "fraggle_io_initialize nf90_def_var regime_varid") + call check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & + [ nc%event_dimid], nc%Qloss_varid), "fraggle_io_initialize nf90_def_var Qloss_varid") + call check( nf90_def_var(nc%id, nc%name_varname, NF90_CHAR, & + [nc%str_dimid, nc%id_dimid ], nc%name_varid), "fraggle_io_initialize nf90_def_var name_varid") + + call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & + [nc%str_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%ptype_varid), "fraggle_io_initialize nf90_def_var ptype_varid") - call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, & - [ nc%event_dimid], nc%loop_varid), "fraggle_io_initialize nf90_def_var loop_varid") + call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, & + [ nc%event_dimid], nc%loop_varid), "fraggle_io_initialize nf90_def_var loop_varid") - call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& - [ nc%space_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%rh_varid), "fraggle_io_initialize nf90_def_var rh_varid") + call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& + [ nc%space_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%rh_varid), "fraggle_io_initialize nf90_def_var rh_varid") - call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& - [ nc%space_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%vh_varid), "fraggle_io_initialize nf90_def_var vh_varid") + call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& + [ nc%space_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%vh_varid), "fraggle_io_initialize nf90_def_var vh_varid") - call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& - [ nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%Gmass_varid), "fraggle_io_initialize nf90_def_var Gmass_varid") + call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& + [ nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%Gmass_varid), "fraggle_io_initialize nf90_def_var Gmass_varid") - call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& - [ nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%radius_varid), "fraggle_io_initialize nf90_def_var radius_varid") + call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& + [ nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%radius_varid), "fraggle_io_initialize nf90_def_var radius_varid") - call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& - [ nc%space_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%Ip_varid), "fraggle_io_initialize nf90_def_var Ip_varid") + call check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& + [ nc%space_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%Ip_varid), "fraggle_io_initialize nf90_def_var Ip_varid") - call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& - [ nc%space_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%rot_varid), "fraggle_io_initialize nf90_def_var rot_varid") + call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& + [ nc%space_dimid, nc%id_dimid, nc%stage_dimid, nc%event_dimid], nc%rot_varid), "fraggle_io_initialize nf90_def_var rot_varid") - call check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "fraggle_io_initialize_output nf90_def_var KE_orb_varid") + call check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "fraggle_io_initialize_output nf90_def_var KE_orb_varid") - call check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "fraggle_io_initialize_output nf90_def_var KE_spin_varid" ) + call check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "fraggle_io_initialize_output nf90_def_var KE_spin_varid" ) - call check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "fraggle_io_initialize_output nf90_def_var PE_varid" ) + call check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "fraggle_io_initialize_output nf90_def_var PE_varid" ) - call check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "fraggle_io_initialize_output nf90_def_var L_orb_varid" ) + call check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & + [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "fraggle_io_initialize_output nf90_def_var L_orb_varid" ) - call check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_spin_varid), "fraggle_io_initialize_output nf90_def_var L_spin_varid" ) + call check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& + [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_spin_varid), "fraggle_io_initialize_output nf90_def_var L_spin_varid" ) - call check( nf90_inquire(nc%id, nVariables=nvar), "fraggle_io_initialize nf90_inquire nVariables" ) - do varid = 1, nvar - call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "fraggle_io_initialize nf90_inquire_variable" ) - select case(vartype) - case(NF90_INT) - call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "fraggle_io_initialize nf90_def_var_fill NF90_INT" ) - case(NF90_FLOAT) - call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "fraggle_io_initialize nf90_def_var_fill NF90_FLOAT" ) - case(NF90_DOUBLE) - call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "fraggle_io_initialize nf90_def_var_fill NF90_DOUBLE" ) - case(NF90_CHAR) - call check( nf90_def_var_fill(nc%id, varid, 0, 0), "fraggle_io_initialize nf90_def_var_fill NF90_CHAR" ) - end select - end do - ! Take the file out of define mode - call check( nf90_enddef(nc%id), "fraggle_io_initialize nf90_enddef" ) + call check( nf90_inquire(nc%id, nVariables=nvar), "fraggle_io_initialize nf90_inquire nVariables" ) + do varid = 1, nvar + call check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), "fraggle_io_initialize nf90_inquire_variable" ) + select case(vartype) + case(NF90_INT) + call check( nf90_def_var_fill(nc%id, varid, 0, NF90_FILL_INT), "fraggle_io_initialize nf90_def_var_fill NF90_INT" ) + case(NF90_FLOAT) + call check( nf90_def_var_fill(nc%id, varid, 0, sfill), "fraggle_io_initialize nf90_def_var_fill NF90_FLOAT" ) + case(NF90_DOUBLE) + call check( nf90_def_var_fill(nc%id, varid, 0, dfill), "fraggle_io_initialize nf90_def_var_fill NF90_DOUBLE" ) + case(NF90_CHAR) + call check( nf90_def_var_fill(nc%id, varid, 0, 0), "fraggle_io_initialize nf90_def_var_fill NF90_CHAR" ) + end select + end do + ! Take the file out of define mode + call check( nf90_enddef(nc%id), "fraggle_io_initialize nf90_enddef" ) - ! Add in the space and stage dimension coordinates - call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "fraggle_io_initialize nf90_put_var space" ) - call check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(1), start=[1,1], count=[len(nc%stage_coords(1)),1]), "fraggle_io_initialize nf90_put_var stage 1" ) - call check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(2), start=[1,2], count=[len(nc%stage_coords(2)),1]), "fraggle_io_initialize nf90_put_var stage 2" ) + ! Add in the space and stage dimension coordinates + call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "fraggle_io_initialize nf90_put_var space" ) + call check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(1), start=[1,1], count=[len(nc%stage_coords(1)),1]), "fraggle_io_initialize nf90_put_var stage 1" ) + call check( nf90_put_var(nc%id, nc%stage_varid, nc%stage_coords(2), start=[1,2], count=[len(nc%stage_coords(2)),1]), "fraggle_io_initialize nf90_put_var stage 2" ) - ! Pre-fill id slots with ids - call check( nf90_put_var(nc%id, nc%id_varid, [(-1,i=1,param%maxid+1)], start=[1], count=[param%maxid+1]), "fraggle_io_initialize nf90_put_varid_varid" ) - end associate + end associate + end select return @@ -155,20 +156,19 @@ module subroutine fraggle_io_write_frame(self, nc, param) use netcdf implicit none ! Arguments - class(fraggle_collision_snapshot), intent(in) :: self !! Swiftest encounter structure + class(fraggle_snapshot), intent(in) :: self !! Swiftest encounter structure class(netcdf_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, eslot, idslot, old_mode, npl, stage + integer(I4B) :: i, idslot, old_mode, npl, stage character(len=:), allocatable :: charstring class(swiftest_pl), allocatable :: pl - associate(colliders => self%colliders, fragments => self%fragments) - select type(nc) - class is (encounter_io_parameters) - eslot = param%ioutput - select type(nc) - class is (fraggle_io_parameters) + select type(nc) + class is (fraggle_io_parameters) + select type (param) + class is (symba_parameters) + associate(colliders => self%colliders, fragments => self%fragments, collision_history => param%collision_history, eslot => param%ioutput) call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "fraggle_io_write_frame nf90_set_fill" ) call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[eslot]), "fraggle_io_write_frame nf90_put_var time_varid" ) @@ -188,7 +188,7 @@ module subroutine fraggle_io_write_frame(self, nc, param) end select npl = pl%nbody do i = 1, npl - idslot = pl%id(i) + 1 + idslot = findloc(collision_history%idvals,pl%id(i),dim=1) call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[ idslot ]), "fraggle_io_write_frame nf90_put_var id_varid" ) charstring = trim(adjustl(pl%info(i)%name)) call check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot ], count=[len(charstring), 1]), "fraggle_io_write_frame nf90_put_var name_varid" ) @@ -212,9 +212,9 @@ module subroutine fraggle_io_write_frame(self, nc, param) call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid after" ) call check( nf90_set_fill(nc%id, old_mode, old_mode) ) - end select + end associate end select - end associate + end select return end subroutine fraggle_io_write_frame diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 2f708859c..22ca5cd55 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -174,7 +174,7 @@ module subroutine fraggle_util_final_snapshot(self) !! Finalizer will deallocate all allocatables implicit none ! Arguments - type(fraggle_collision_snapshot), intent(inout) :: self !! Fraggle encountar storage object + type(fraggle_snapshot), intent(inout) :: self !! Fraggle encountar storage object call encounter_util_final_snapshot(self%encounter_snapshot) @@ -257,6 +257,40 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para end subroutine fraggle_util_get_energy_momentum + module subroutine fraggle_util_get_idvalues_snapshot(self, idvals) + !! author: David A. Minton + !! + !! Returns an array of all id values saved in this snapshot + implicit none + ! Arguments + class(fraggle_snapshot), intent(in) :: self !! Fraggle snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + ! Internals + integer(I4B) :: ncoll, nfrag + + if (allocated(self%colliders)) then + ncoll = self%colliders%pl%nbody + else + ncoll = 0 + end if + + if (allocated(self%fragments)) then + nfrag = self%fragments%pl%nbody + else + nfrag = 0 + end if + + if (ncoll + nfrag == 0) return + allocate(idvals(ncoll+nfrag)) + + if (ncoll > 0) idvals(1:ncoll) = self%colliders%pl%id(:) + if (nfrag > 0) idvals(ncoll+1:ncoll+nfrag) = self%fragments%pl%id(:) + + return + + end subroutine fraggle_util_get_idvalues_snapshot + + module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_start) !! Author: David A. Minton !! diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index f9ed5b896..f5429fd12 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -49,7 +49,8 @@ module encounter_classes real(DP) :: t !! Simulation time when snapshot was taken integer(I8B) :: iloop !! Loop number at time of snapshot contains - procedure :: write_frame => encounter_io_write_frame !! Writes a frame of encounter data to file + procedure :: write_frame => encounter_io_write_frame !! Writes a frame of encounter data to file + procedure :: get_idvals => encounter_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot final :: encounter_util_final_snapshot end type encounter_snapshot @@ -68,8 +69,8 @@ module encounter_classes type, extends(swiftest_storage) :: collision_storage contains procedure :: dump => encounter_io_dump_collision !! Dumps contents of encounter history to file - procedure :: make_index_map => encounter_util_index_map_collision !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - procedure :: take_snapshot => encounter_util_snapshot_collision !! Take a minimal snapshot of the system through an encounter + procedure :: take_snapshot => encounter_util_snapshot_collision !! Take a minimal snapshot of the system through an encounter + procedure :: make_index_map => encounter_util_index_map_collision !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id final :: encounter_util_final_collision_storage end type collision_storage @@ -300,14 +301,20 @@ module subroutine encounter_util_final_storage(self) type(encounter_storage(*)), intent(inout) :: self !! SyMBA nbody system object end subroutine encounter_util_final_storage + module subroutine encounter_util_get_idvalues_snapshot(self, idvals) + implicit none + class(encounter_snapshot), intent(in) :: self !! Encounter snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + end subroutine encounter_util_get_idvalues_snapshot + module subroutine encounter_util_index_map_collision(self) implicit none - class(collision_storage(*)), intent(inout) :: self !! E + class(collision_storage(*)), intent(inout) :: self !! Collision storage object end subroutine encounter_util_index_map_collision module subroutine encounter_util_index_map_encounter(self) implicit none - class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object + class(encounter_storage(*)), intent(inout) :: self !! Encounter storage object end subroutine encounter_util_index_map_encounter module subroutine encounter_util_resize_list(self, nnew) diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 2bd120b83..989399f94 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -132,14 +132,15 @@ module fraggle_classes procedure :: initialize => fraggle_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object end type fraggle_io_parameters - type, extends(encounter_snapshot) :: fraggle_collision_snapshot + type, extends(encounter_snapshot) :: fraggle_snapshot logical :: lcollision !! Indicates that this snapshot contains at least one collision class(fraggle_colliders), allocatable :: colliders !! Colliders object at this snapshot class(fraggle_fragments), allocatable :: fragments !! Fragments object at this snapshot contains procedure :: write_frame => fraggle_io_write_frame !! Writes a frame of encounter data to file + procedure :: get_idvals => fraggle_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot final :: fraggle_util_final_snapshot - end type fraggle_collision_snapshot + end type fraggle_snapshot interface module subroutine fraggle_generate_fragments(self, colliders, system, param, lfailure) @@ -160,7 +161,7 @@ end subroutine fraggle_io_initialize_output module subroutine fraggle_io_write_frame(self, nc, param) implicit none - class(fraggle_collision_snapshot), intent(in) :: self !! Swiftest encounter structure + class(fraggle_snapshot), intent(in) :: self !! Swiftest encounter structure class(netcdf_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine fraggle_io_write_frame @@ -295,7 +296,7 @@ end subroutine fraggle_util_final_fragments module subroutine fraggle_util_final_snapshot(self) implicit none - type(fraggle_collision_snapshot), intent(inout) :: self !! Fraggle storage snapshot object + type(fraggle_snapshot), intent(inout) :: self !! Fraggle storage snapshot object end subroutine fraggle_util_final_snapshot module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore) @@ -308,6 +309,12 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa end subroutine fraggle_util_get_energy_momentum + module subroutine fraggle_util_get_idvalues_snapshot(self, idvals) + implicit none + class(fraggle_snapshot), intent(in) :: self !! Fraggle snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + end subroutine fraggle_util_get_idvalues_snapshot + module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_start) implicit none class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index b89e9eee8..6dc6a9877 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -556,6 +556,7 @@ module swiftest_classes ! procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies. procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum + procedure :: get_idvals => util_get_idvalues_system !! Returns an array of all id values in use in the system procedure :: rescale => util_rescale_system !! Rescales the system into a new set of units procedure :: validate_ids => util_valid_id_system !! Validate the numerical ids passed to the system and save the maximum value generic :: write_frame => write_frame_system, write_frame_netcdf !! Generic method call for reading a frame of output data @@ -1620,6 +1621,12 @@ module subroutine util_get_energy_momentum_system(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine util_get_energy_momentum_system + module subroutine util_get_idvalues_system(self, idvals) + implicit none + class(swiftest_nbody_system), intent(in) :: self !! Encounter snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + end subroutine util_get_idvalues_system + module subroutine util_set_beg_end_pl(self, xbeg, xend, vbeg) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object diff --git a/src/util/util_index.f90 b/src/util/util_index.f90 index ae4a80ce8..0fbd40319 100644 --- a/src/util/util_index.f90 +++ b/src/util/util_index.f90 @@ -45,55 +45,113 @@ module subroutine util_index_array(ind_arr, n) end subroutine util_index_array + module subroutine util_get_idvalues_system(self, idvals) + !! author: David A. Minton + !! + !! Returns an array of all id values saved in this snapshot + implicit none + ! Arguments + class(swiftest_nbody_system), intent(in) :: self !! Encounter snapshot object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values saved in this snapshot + ! Internals + integer(I4B) :: npl, ntp + + if (allocated(self%pl)) then + npl = self%pl%nbody + else + npl = 0 + end if + if (allocated(self%tp)) then + ntp = self%tp%nbody + else + ntp = 0 + end if + + allocate(idvals(1 + npl+ntp)) + + idvals(1) = self%cb%id + if (npl > 0) idvals(2:npl+1) = self%pl%id(:) + if (ntp > 0) idvals(npl+2:npl+ntp+1) = self%tp%id(:) + + return + + end subroutine util_get_idvalues_system + + + subroutine util_get_vals_storage(storage, idvals, tvals) + !! author: David A. Minton + !! + !! Gets the id values in a storage object, regardless of whether it is encounter of collision + ! Argument + class(swiftest_storage(*)), intent(in) :: storage !! Swiftest storage object + integer(I4B), dimension(:), allocatable, intent(out) :: idvals !! Array of all id values in all snapshots + real(DP), dimension(:), allocatable, intent(out) :: tvals !! Array of all time values in all snapshots + ! Internals + integer(I4B) :: i, n, nlo, nhi, ntotal + integer(I4B), dimension(:), allocatable :: itmp + + associate(nsnaps => storage%iframe) + + allocate(tvals(nsnaps)) + tvals(:) = 0.0_DP + + ! First pass to get total number of ids + ntotal = 0 + do i = 1, nsnaps + if (allocated(storage%frame(i)%item)) then + select type(snapshot => storage%frame(i)%item) + class is (swiftest_nbody_system) + tvals(i) = snapshot%t + call snapshot%get_idvals(itmp) + if (allocated(itmp)) then + n = size(itmp) + ntotal = ntotal + n + end if + end select + end if + end do + + allocate(idvals(ntotal)) + nlo = 1 + ! Second pass to store all ids get all of the ids stored + do i = 1, nsnaps + if (allocated(storage%frame(i)%item)) then + select type(snapshot => storage%frame(i)%item) + class is (swiftest_nbody_system) + tvals(i) = snapshot%t + call snapshot%get_idvals(itmp) + if (allocated(itmp)) then + n = size(itmp) + nhi = nlo + n - 1 + idvals(nlo:nhi) = itmp(1:n) + nlo = nhi + 1 + end if + end select + end if + end do + + end associate + return + end subroutine util_get_vals_storage + + module subroutine util_index_map_storage(self) !! author: David A. Minton !! !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id implicit none ! Arguments - class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object ! Internals - integer(I4B) :: i, n, nold, nt integer(I4B), dimension(:), allocatable :: idvals real(DP), dimension(:), allocatable :: tvals - - if (self%nid == 0) return - allocate(idvals(self%nid)) - allocate(tvals(self%nframes)) - - n = 0 - nold = 1 - nt = 0 - do i = 1, self%nframes - if (allocated(self%frame(i)%item)) then - nt = i - select type(snapshot => self%frame(i)%item) - class is (swiftest_nbody_system) - tvals(i) = snapshot%t - ! Central body - n = n + 1 - idvals(n) = snapshot%cb%id - nold = n + 1 - if (allocated(snapshot%pl)) then - n = n + snapshot%pl%nbody - idvals(nold:n) = snapshot%pl%id(:) - nold = n+1 - end if - if (allocated(snapshot%tp)) then - n = n + snapshot%tp%nbody - idvals(nold:n) = snapshot%tp%id(:) - nold = n+1 - end if - end select - else - exit - end if - end do + + call util_get_vals_storage(self, idvals, tvals) call util_unique(idvals,self%idvals,self%idmap) self%nid = size(self%idvals) - call util_unique(tvals(1:nt),self%tvals,self%tmap) + call util_unique(tvals,self%tvals,self%tmap) self%nt = size(self%tvals) return