Another major restructuring aimed toward improving how the storage ob…

…jects handle the indexing of bodies and time steps

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

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

src/encounter/encounter_util.f90

@@ -189,35 +189,44 @@ 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
       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
-                  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

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
 

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

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

src/modules/swiftest_classes.f90

@@ -148,23 +148,26 @@ 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)
       !! An class that establishes the pattern for various storage objects
       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

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))

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)

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

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

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