From 0fe3f1f90cfa394f70dd49a17968be31b754e1b7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 11 Dec 2022 09:31:08 -0500 Subject: [PATCH] Another major restructuring aimed toward improving how the storage objects handle the indexing of bodies and time steps --- src/CMakeLists.txt | 1 + src/encounter/encounter_io.f90 | 3 +- src/encounter/encounter_util.f90 | 23 ++++++++++---- src/io/io.f90 | 8 +++-- src/main/swiftest_driver.f90 | 8 ++--- src/modules/encounter_classes.f90 | 12 +++---- src/modules/swiftest_classes.f90 | 50 ++++++++++++++++++++--------- src/netcdf/netcdf.f90 | 4 +-- src/symba/symba_io.f90 | 3 +- src/symba/symba_util.f90 | 17 +--------- src/util/util_index.f90 | 42 ++++++++++++++++++++++++ src/util/util_reset.f90 | 7 ++-- src/util/util_snapshot.f90 | 33 +++++++++++++++++++ src/util/util_unique.f90 | 53 +++++++++++++++++++++++++++---- 14 files changed, 198 insertions(+), 66 deletions(-) create mode 100644 src/util/util_snapshot.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6f382cd8e..594850a50 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -85,6 +85,7 @@ SET(FAST_MATH_FILES ${SRC}/util/util_reset.f90 ${SRC}/util/util_resize.f90 ${SRC}/util/util_set.f90 + ${SRC}/util/util_snapshot.f90 ${SRC}/util/util_solve.f90 ${SRC}/util/util_sort.f90 ${SRC}/util/util_spill.f90 diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index c5e458e57..0f821f49b 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -50,12 +50,11 @@ module subroutine encounter_io_dump_storage(self, param) ! Internals integer(I4B) :: i - call self%mapid() do i = 1, self%nframes if (allocated(self%frame(i)%item)) then select type(snapshot => self%frame(i)%item) class is (encounter_snapshot) - param%ioutput = self%tslot(i) + param%ioutput = self%tmap(i) call snapshot%write_frame(self%nc,param) end select else diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 9385ad2c0..ef89b3cf3 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -189,11 +189,14 @@ module subroutine encounter_util_index_map_storage(self) ! Arguments class(encounter_storage(*)), intent(inout) :: self !! Swiftest storage object ! Internals - integer(I4B) :: i, n, nold - integer(I4B), dimension(:), allocatable :: idlist + ! Internals + integer(I4B) :: i, n, nold, nt + integer(I4B), dimension(:), allocatable :: idvals + real(DP), dimension(:), allocatable :: tvals if (self%nid == 0) return - allocate(idlist(self%nid)) + allocate(idvals(self%nid)) + allocate(tvals(self%nframes)) n = 0 nold = 1 @@ -201,23 +204,29 @@ module subroutine encounter_util_index_map_storage(self) 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 - idlist(nold:n) = snapshot%pl%id(:) + idvals(nold:n) = snapshot%pl%id(:) nold = n+1 - end if + end if if (allocated(snapshot%tp)) then n = n + snapshot%tp%nbody - idlist(nold:n) = snapshot%tp%id(:) + idvals(nold:n) = snapshot%tp%id(:) nold = n+1 end if end select else + nt = i-1 exit end if end do - call util_unique(idlist,self%idmap) + call util_unique(idvals,self%idvals,self%idmap) + self%nid = size(self%idvals) + + call util_unique(tvals(1:nt),self%tvals,self%tmap) + self%nt = size(self%tvals) return end subroutine encounter_util_index_map_storage diff --git a/src/io/io.f90 b/src/io/io.f90 index 74116d5b0..3d694e027 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -294,9 +294,11 @@ module subroutine io_dump_storage(self, param) integer(I4B) :: i integer(I8B) :: iloop_start - iloop_start = max(param%iloop - int(param%istep_out * param%dump_cadence, kind=I8B),1) + if (self%iframe == 0) return + iloop_start = param%iloop - int(param%istep_out * param%dump_cadence, kind=I8B) + 1 + call self%make_index_map() do i = 1, param%dump_cadence - param%ioutput = max(int(iloop_start / param%istep_out, kind=I4B),1) + i + param%ioutput = iloop_start + self%tmap(i) if (allocated(self%frame(i)%item)) then select type(system => self%frame(i)%item) class is (swiftest_nbody_system) @@ -305,7 +307,7 @@ module subroutine io_dump_storage(self, param) deallocate(self%frame(i)%item) end if end do - + call self%reset() return end subroutine io_dump_storage diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index bfc0b38c6..803ca6849 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -88,15 +88,15 @@ program swiftest_driver call system%initialize(param) - ! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file. if (param%lrestart) then if (param%lenergy) call system%conservation_report(param, lterminal=.true.) else if (param%lenergy) call system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum - call system%write_frame(param) - call system%dump(param) + call system_history%take_snapshot(system) + call system_history%dump(param) end if + call system%dump(param) write(display_unit, *) " *************** Main Loop *************** " @@ -130,7 +130,7 @@ program swiftest_driver if (iout == istep_out) then iout = 0 idump = idump + 1 - system_history%frame(idump) = system ! Store a snapshot of the system for posterity + call system_history%take_snapshot(system) if (idump == dump_cadence) then idump = 0 diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 49065fc2a..b1b7bb079 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -68,18 +68,18 @@ module encounter_classes type, extends(swiftest_storage) :: encounter_storage class(encounter_io_parameters), allocatable :: nc !! NetCDF parameter object containing the details about the file attached to this storage object contains - procedure :: dump => encounter_io_dump_storage !! Dumps contents of encounter history to file - procedure :: mapid => encounter_util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - final :: encounter_util_final_storage + procedure :: dump => encounter_io_dump_storage !! Dumps contents of encounter history to file + procedure :: make_index_map => encounter_util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id + final :: encounter_util_final_storage end type encounter_storage !> A class that that is used to store simulation history data between file output type, extends(swiftest_storage) :: collision_storage class(encounter_io_parameters), allocatable :: nc !! NetCDF parameter object containing the details about the file attached to this storage object contains - procedure :: dump => encounter_io_dump_collision_storage !! Dumps contents of encounter history to file - procedure :: mapid => encounter_util_index_map_collision_storage !! 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 + procedure :: dump => encounter_io_dump_collision_storage !! Dumps contents of encounter history to file + procedure :: make_index_map => encounter_util_index_map_collision_storage !! 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 type encounter_bounding_box_1D diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 6fb35383e..4894e305d 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -148,7 +148,7 @@ module swiftest_classes contains procedure :: store => util_copy_store !! Stores a snapshot of the nbody system so that later it can be retrieved for saving to file. generic :: assignment(=) => store - final :: util_final_storage_frame + final :: util_final_storage_frame end type type :: swiftest_storage(nframes) @@ -156,15 +156,18 @@ module swiftest_classes integer(I4B), len :: nframes = 4096 !! Total number of frames that can be stored type(swiftest_storage_frame), dimension(nframes) :: frame !! Array of stored frames integer(I4B) :: iframe = 0 !! Index of the last frame stored in the system - integer(I4B), dimension(nframes) :: tslot !! The value of the time dimension index associated with each frame - real(DP), dimension(nframes) :: tvals !! Stored time values for snapshots - integer(I4B), dimension(:), allocatable :: idmap !! The id value -> index map integer(I4B) :: nid !! Number of unique id values in all saved snapshots + integer(I4B), dimension(:), allocatable :: idvals !! The set of unique id values contained in the snapshots + integer(I4B), dimension(:), allocatable :: idmap !! The id value -> index map + integer(I4B) :: nt !! Number of unique time values in all saved snapshots + real(DP), dimension(:), allocatable :: tvals !! The set of unique time values contained in the snapshots + integer(I4B), dimension(:), allocatable :: tmap !! The t value -> index map contains - procedure :: dump => io_dump_storage !! Dumps storage object contents to file - procedure :: mapid => util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id - procedure :: reset => util_reset_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 - final :: util_final_storage + procedure :: dump => io_dump_storage !! Dumps storage object contents to file + procedure :: make_index_map => util_index_map_storage !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id + procedure :: reset => util_reset_storage !! Resets a storage object by deallocating all items and resetting the frame counter to 0 + procedure :: take_snapshot => util_snapshot_system !! Takes a snapshot of the system for later file storage + final :: util_final_storage end type swiftest_storage !******************************************************************************************************************************** @@ -550,7 +553,7 @@ module swiftest_classes procedure :: finalize => setup_finalize_system !! Runs any finalization subroutines when ending the simulation. procedure :: initialize => setup_initialize_system !! Initialize the system from input files procedure :: init_particle_info => setup_initialize_particle_info_system !! Initialize the system from input files - ! procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. + ! 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 :: rescale => util_rescale_system !! Rescales the system into a new set of units @@ -611,11 +614,9 @@ subroutine abstract_step_system(self, param, t, dt) real(DP), intent(in) :: t !! Simulation time real(DP), intent(in) :: dt !! Current stepsize end subroutine abstract_step_system - end interface interface - module subroutine check(status, call_identifier) implicit none integer, intent (in) :: status !! The status code returned by a NetCDF function @@ -1690,6 +1691,12 @@ module subroutine util_set_rhill_approximate(self,cb) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine util_set_rhill_approximate + + module subroutine util_snapshot_system(self, system) + implicit none + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object to store + end subroutine util_snapshot_system end interface interface util_solve_linear_system @@ -1957,12 +1964,25 @@ module subroutine util_spill_tp(self, discards, lspill_list, ldestructive) logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not end subroutine util_spill_tp - module subroutine util_unique(input_array, output_array) + end interface + + interface util_unique + module subroutine util_unique_DP(input_array, output_array, index_map) implicit none - integer(I4B), dimension(:), intent(in) :: input_array - integer(I4B), dimension(:), allocatable, intent(out) :: output_array - end subroutine util_unique + real(DP), dimension(:), intent(in) :: input_array !! Unsorted input array + real(DP), dimension(:), allocatable, intent(out) :: output_array !! Sorted array of unique values + integer(I4B), dimension(:), allocatable, intent(out) :: index_map !! An array of the same size as input_array that such that any for any index i, output_array(index_map(i)) = input_array(i) + end subroutine util_unique_DP + module subroutine util_unique_I4B(input_array, output_array, index_map) + implicit none + integer(I4B), dimension(:), intent(in) :: input_array !! Unsorted input array + integer(I4B), dimension(:), allocatable, intent(out) :: output_array !! Sorted array of unique values + integer(I4B), dimension(:), allocatable, intent(out) :: index_map !! An array of the same size as input_array that such that any for any index i, output_array(index_map(i)) = input_array(i) + end subroutine util_unique_I4B + end interface util_unique + + interface module subroutine util_valid_id_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 8cfc6432f..9bc945227 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -825,8 +825,8 @@ module subroutine netcdf_read_particle_info_system(self, nc, param, plmask, tpma pl%id(:) = pack(itemp, plmask) tp%id(:) = pack(itemp, tpmask) cb%id = 0 - pl%id(:) = pack([(i,i=1,idmax)],plmask) - tp%id(:) = pack([(i,i=1,idmax)],tpmask) + pl%id(:) = pack([(i,i=0,idmax-1)],plmask) + tp%id(:) = pack([(i,i=0,idmax-1)],tpmask) call check( nf90_get_var(nc%id, nc%name_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_read_particle_info_system nf90_getvar name_varid" ) call cb%info%set_value(name=ctemp(1)) diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index cd8693598..d1f36abca 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -30,7 +30,7 @@ module subroutine symba_io_dump_encounter(self, param) if (num_enc_frames > 0) then ! Create and save the output files for this encounter and fragmentation nce%file_number = nce%file_number + 1 - nce%time_dimsize = maxval(encounter_history%tslot(:)) + call encounter_history%make_index_map() write(nce%file_name, '("encounter_",I0.6,".nc")') nce%file_number call nce%initialize(param) call encounter_history%dump(param) @@ -44,6 +44,7 @@ module subroutine symba_io_dump_encounter(self, param) if (num_coll_frames > 0) then ncc%file_number = ncc%file_number + 1 ncc%event_dimsize = num_coll_frames + call collision_history%make_index_map() write(ncc%file_name, '("collision_",I0.6,".nc")') ncc%file_number call ncc%initialize(param) call collision_history%dump(param) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index f7711e6bf..4381ce93b 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -897,10 +897,6 @@ subroutine symba_util_save_collision(system, snapshot) nbig = nbig * 2 end do allocate(collision_storage(nbig) :: tmp) - tmp%tvals(1:nold) = system%collision_history%tvals(1:nold) - tmp%tvals(nold+1:nbig) = huge(1.0_DP) - tmp%tslot(1:nold) = system%collision_history%tslot(1:nold) - tmp%tslot(nold+1:nbig) = 0 tmp%iframe = system%collision_history%iframe call move_alloc(system%collision_history%nc, tmp%nc) @@ -947,10 +943,6 @@ subroutine symba_util_save_encounter(system, snapshot, t) nbig = nbig * 2 end do allocate(encounter_storage(nbig) :: tmp) - tmp%tvals(1:nold) = system%encounter_history%tvals(1:nold) - tmp%tvals(nold+1:nbig) = huge(1.0_DP) - tmp%tslot(1:nold) = system%encounter_history%tslot(1:nold) - tmp%tslot(nold+1:nbig) = 0 tmp%iframe = system%encounter_history%iframe call move_alloc(system%encounter_history%nc, tmp%nc) @@ -964,14 +956,7 @@ subroutine symba_util_save_encounter(system, snapshot, t) ! Find out which time slot this belongs in by searching for an existing slot ! with the same value of time or the first available one - do i = 1, nnew - if (t <= system%encounter_history%tvals(i)) then - system%encounter_history%tvals(i) = t - system%encounter_history%tslot(nnew) = i - system%encounter_history%frame(nnew) = snapshot - exit - end if - end do + system%encounter_history%frame(nnew) = snapshot return end subroutine symba_util_save_encounter diff --git a/src/util/util_index.f90 b/src/util/util_index.f90 index 20772a2a6..ae4a80ce8 100644 --- a/src/util/util_index.f90 +++ b/src/util/util_index.f90 @@ -53,6 +53,48 @@ module subroutine util_index_map_storage(self) ! Arguments 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_unique(idvals,self%idvals,self%idmap) + self%nid = size(self%idvals) + + call util_unique(tvals(1:nt),self%tvals,self%tmap) + self%nt = size(self%tvals) return end subroutine util_index_map_storage diff --git a/src/util/util_reset.f90 b/src/util/util_reset.f90 index a588c65fe..9b37f7d15 100644 --- a/src/util/util_reset.f90 +++ b/src/util/util_reset.f90 @@ -24,11 +24,12 @@ module subroutine util_reset_storage(self) do i = 1, self%nframes if (allocated(self%frame(i)%item)) deallocate(self%frame(i)%item) end do - self%tslot(:) = 0 - self%tvals(:) = huge(1.0_DP) - self%iframe = 0 + if (allocated(self%idmap)) deallocate(self%idmap) + if (allocated(self%tmap)) deallocate(self%tmap) self%nid = 0 + self%nt = 0 + self%iframe = 0 return end subroutine util_reset_storage diff --git a/src/util/util_snapshot.f90 b/src/util/util_snapshot.f90 new file mode 100644 index 000000000..1c67bc7f8 --- /dev/null +++ b/src/util/util_snapshot.f90 @@ -0,0 +1,33 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule(swiftest_classes) s_util_snapshot + use swiftest +contains + + module subroutine util_snapshot_system(self, system) + !! author: David A. Minton + !! + !! Takes a snapshot of the system for later file storage + implicit none + ! Arguments + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest storage object + class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object to store + + self%iframe = self%iframe + 1 + self%nt = self%iframe + self%frame(self%iframe) = system ! Store a snapshot of the system for posterity + self%nid = self%nid + 1 ! Central body + if (allocated(system%pl)) self%nid = self%nid + system%pl%nbody + if (allocated(system%tp)) self%nid = self%nid + system%tp%nbody + + return + end subroutine util_snapshot_system + +end submodule s_util_snapshot \ No newline at end of file diff --git a/src/util/util_unique.f90 b/src/util/util_unique.f90 index 9cf77536c..19eb4ba78 100644 --- a/src/util/util_unique.f90 +++ b/src/util/util_unique.f90 @@ -11,31 +11,70 @@ use swiftest contains - module subroutine util_unique(input_array, output_array) + module subroutine util_unique_DP(input_array, output_array, index_map) !! author: David A. Minton !! - !! Takes an input unsorted integer array and returns a new array of sorted, unique values + !! Takes an input unsorted integer array and returns a new array of sorted, unique values (DP version) implicit none ! Arguments - integer(I4B), dimension(:), intent(in) :: input_array - integer(I4B), dimension(:), allocatable, intent(out) :: output_array + real(DP), dimension(:), intent(in) :: input_array !! Unsorted input array + real(DP), dimension(:), allocatable, intent(out) :: output_array !! Sorted array of unique values + integer(I4B), dimension(:), allocatable, intent(out) :: index_map !! An array of the same size as input_array that such that any for any index i, output_array(index_map(i)) = input_array(i) + ! Internals + real(DP), dimension(:), allocatable :: unique_array + integer(I4B) :: n + real(DP) :: lo, hi + + allocate(unique_array, mold=input_array) + allocate(index_map(size(input_array))) + lo = minval(input_array) - 1 + hi = maxval(input_array) + + n = 0 + do + n = n + 1 + lo = minval(input_array(:), mask=input_array(:) > lo) + unique_array(n) = lo + where(input_array(:) == lo) index_map(:) = n + if (lo >= hi) exit + enddo + allocate(output_array(n), source=unique_array(1:n)) + + return + end subroutine util_unique_DP + + + module subroutine util_unique_I4B(input_array, output_array, index_map) + !! author: David A. Minton + !! + !! Takes an input unsorted integer array and returns a new array of sorted, unique values (I4B version) + implicit none + ! Arguments + integer(I4B), dimension(:), intent(in) :: input_array !! Unsorted input array + integer(I4B), dimension(:), allocatable, intent(out) :: output_array !! Sorted array of unique values + integer(I4B), dimension(:), allocatable, intent(out) :: index_map !! An array of the same size as input_array that such that any for any index i, output_array(index_map(i)) = input_array(i) ! Internals integer(I4B), dimension(:), allocatable :: unique_array integer(I4B) :: n, lo, hi allocate(unique_array, mold=input_array) + allocate(index_map, mold=input_array) lo = minval(input_array) - 1 hi = maxval(input_array) n = 0 - do while (lo < hi) + do n = n + 1 - lo = minval(input_array, mask=input_array > lo) + lo = minval(input_array(:), mask=input_array(:) > lo) unique_array(n) = lo + where(input_array(:) == lo) index_map(:) = n + if (lo >= hi) exit enddo allocate(output_array(n), source=unique_array(1:n)) return - end subroutine util_unique + end subroutine util_unique_I4B + + end submodule s_util_unique \ No newline at end of file