Fixed the code so that it only outputs encounter data on the dump steps

python/swiftest/swiftest/io.py

@@ -64,7 +64,7 @@
 # handles strings differently than Python's Xarray.
 string_varnames = ["name", "particle_type", "status", "origin_type"]
 char_varnames = ["space"]
-int_varnames = ["id", "ntp", "npl", "nplm", "discard_body_id", "collision_id"]
+int_varnames = ["id", "ntp", "npl", "nplm", "discard_body_id", "collision_id", "loopnum"]
 
 def bool2yesno(boolval):
     """
@@ -816,8 +816,8 @@ def process_netcdf_input(ds, param):
     -------
     ds : xarray dataset
     """
-
-    ds = ds.where(~np.isnan(ds.id) ,drop=True)
+    #
+    ds = ds.where(ds.id >=0,drop=True)
     if param['OUT_TYPE'] == "NETCDF_DOUBLE":
         ds = fix_types(ds,ftype=np.float64)
     elif param['OUT_TYPE'] == "NETCDF_FLOAT":

python/swiftest/swiftest/simulation_class.py

@@ -2780,6 +2780,11 @@ def _preprocess(ds, param):
         tgood,tid = np.unique(self.enc.time,return_index=True)
         self.enc = self.enc.isel(time=tid)
 
+        # Reduce the dimensionality of variables that got expanded in the combine process
+        self.enc['loopnum'] = self.enc['loopnum'].max(dim="name")
+        self.enc['id'] = self.enc['id'].max(dim="time")
+        self.enc['particle_type'] = self.enc['particle_type'].max(dim="time")
+
         return
 
 

src/encounter/encounter_io.f90

@@ -53,7 +53,7 @@ module subroutine encounter_io_initialize(self, param)
       class(encounter_io_parameters), intent(inout) :: self    !! Parameters used to identify a particular NetCDF dataset
       class(swiftest_parameters),     intent(in)    :: param   !! Current run configuration parameters
       ! Internals
-      integer(I4B) :: nvar, varid, vartype
+      integer(I4B) :: i, nvar, varid, vartype
       real(DP) :: dfill
       real(SP) :: sfill
       logical :: fileExists
@@ -97,7 +97,7 @@ module subroutine encounter_io_initialize(self, param)
          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%level_varname, NF90_INT, [nc%id_dimid, nc%time_dimid], nc%level_varid), "encounter_io_initialize nf90_def_var level_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
@@ -126,6 +126,10 @@ module subroutine encounter_io_initialize(self, param)
 
          ! 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 id slots with ids
+
+         call check( nf90_put_var(nc%id, nc%id_varid, [(-1,i=1,param%maxid)], start=[1], count=[param%maxid]), "encounter_io_initialize nf90_put_var pl id_varid"  )
       end associate
 
       return
@@ -145,7 +149,7 @@ module subroutine encounter_io_write_frame(self, nc, param)
       ! Arguments
       class(encounter_snapshot),      intent(in)    :: self  !! Swiftest encounter structure
       class(encounter_io_parameters), intent(inout) :: nc    !! Parameters used to identify a particular encounter io NetCDF dataset
-      class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
+      class(swiftest_parameters),     intent(inout) :: param !! Current run configuration parameters
       ! Internals
       integer(I4B)           :: i, tslot, idslot, old_mode, npl, ntp
       character(len=NAMELEN) :: charstring
@@ -155,19 +159,16 @@ module subroutine encounter_io_write_frame(self, nc, param)
       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"  )
+      call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_write_frame nf90_put_var pl loop_varid"  )
 
       associate(pl => self%pl, tp => self%tp)
          npl = pl%nbody
          do i = 1, npl
             idslot = pl%id(i)
-            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%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"  )
             call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl Gmass_varid"  )
-            select type(pl)
-            class is (symba_pl)
-               call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl level_varid"  )
-            end select
 
             if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl radius_varid"  )
 

src/io/io.f90

@@ -264,6 +264,17 @@ module subroutine io_dump_system(self, param)
       idx = idx + 1
       if (idx > NDUMPFILES) idx = 1
 
+      ! Dump the encounter history if necessary
+      select type(param)
+      class is (symba_parameters)
+         if (param%lencounter_save) then
+            select type(self)
+            class is (symba_nbody_system)
+               call self%dump_encounter(param)
+            end select
+         end if
+      end select
+
       return
    end subroutine io_dump_system
 

src/main/swiftest_driver.f90

@@ -51,6 +51,8 @@ program swiftest_driver
    end select
    param%integrator = trim(adjustl(integrator))
    call param%set_display(display_style)
+   call param%read_in(param_file_name)
+   call setup_construct_system(nbody_system, param)
 
    !> Define the maximum number of threads
    nthreads = 1            ! In the *serial* case
@@ -60,8 +62,6 @@ program swiftest_driver
    !$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads 
    !$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads
 
-   call setup_construct_system(nbody_system, param)
-   call param%read_in(param_file_name)
 
    associate(t               => nbody_system%t, &
              t0              => param%t0, &

src/modules/encounter_classes.f90

@@ -43,29 +43,30 @@ module encounter_classes
 
    type :: encounter_snapshot
       !! A simplified version of a SyMBA nbody system object for storing minimal snapshots of the system state during encounters
-      class(swiftest_pl), allocatable :: pl  !! Massive body data structure
-      class(swiftest_tp), allocatable :: tp  !! Test particle data structure
-      real(DP)                        :: t   !! Simulation time when snapshot was taken
+      class(swiftest_pl), allocatable :: pl    !! Massive body data structure
+      class(swiftest_tp), allocatable :: tp    !! Test particle data structure
+      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 
       final     ::                encounter_util_final_snapshot
    end type encounter_snapshot
 
    !> NetCDF dimension and variable names for the enounter save object
    type, extends(netcdf_parameters) :: encounter_io_parameters
-      integer(I4B)       :: ienc_frame    = 1       !! Current frame number for the encounter history
-      character(STRMAX)  :: enc_file                !! Encounter output file name
-      character(NAMELEN) :: level_varname = "level" !! Recursion depth
-      integer(I4B)       :: level_varid             !! ID for the recursion level variable
-      integer(I4B)       :: time_dimsize = 0        !! Number of time values in snapshot
-      integer(I4B)       :: id_dimsize   = 0        !! Number of potential id values in snapshot
+      integer(I4B)       :: ienc_frame    = 1        !! Current frame number for the encounter history
+      character(STRMAX)  :: enc_file                 !! Encounter output file name
+      character(NAMELEN) :: loop_varname = "loopnum" !! Loop number for encounter
+      integer(I4B)       :: loop_varid               !! ID for the recursion level variable
+      integer(I4B)       :: time_dimsize = 0         !! Number of time values in snapshot
+      integer(I4B)       :: id_dimsize   = 0         !! Number of potential id values in snapshot
    contains
       procedure :: initialize => encounter_io_initialize !! Initialize a set of parameters used to identify a NetCDF output object
    end type encounter_io_parameters
 
    !> A class that that is used to store simulation history data between file output
    type, extends(swiftest_storage) :: encounter_storage
-      type(encounter_io_parameters) :: nc  !! NetCDF parameter object containing the details about the file attached to this storage object 
+      type(encounter_io_parameters) :: nc  !! NetCDF parameter object containing the details about the file attached to this storage object
    contains
       procedure :: dump   => encounter_io_dump !! Dumps contents of encounter history to file
       final     ::           encounter_util_final_storage

src/modules/symba_classes.f90

@@ -199,8 +199,7 @@ module symba_classes
       procedure :: recursive_step   => symba_step_recur_system            !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current recursion level, if applicable, and descend to the next deeper level if necessary
       procedure :: reset            => symba_step_reset_system            !! Resets pl, tp,and encounter structures at the start of a new step 
       procedure :: snapshot         => symba_util_take_encounter_snapshot !! Take a minimal snapshot of the system through an encounter
-      procedure :: start_encounter  => symba_io_start_encounter           !! Initializes the new encounter history
-      procedure :: stop_encounter   => symba_io_stop_encounter            !! Saves the encounter and/or fragmentation data to file(s)   
+      procedure :: dump_encounter   => symba_io_dump_encounter            !! Saves the encounter and/or fragmentation data to file(s)   
       final     ::                     symba_util_final_system            !! Finalizes the SyMBA nbody system object - deallocates all allocatables
    end type symba_nbody_system
 
@@ -409,19 +408,11 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms
          character(len=*),       intent(inout) :: iomsg     !! Message to pass if iostat /= 0
       end subroutine symba_io_param_writer
 
-      module subroutine symba_io_start_encounter(self, param, t)
+      module subroutine symba_io_dump_encounter(self, param)
          implicit none
          class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
          class(symba_parameters),    intent(inout) :: param !! Current run configuration parameters 
-         real(DP),                   intent(in)    :: t     !! Current simulation time
-      end subroutine symba_io_start_encounter
-
-      module subroutine symba_io_stop_encounter(self, param, t)
-         implicit none
-         class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
-         class(symba_parameters),    intent(inout) :: param !! Current run configuration parameters 
-         real(DP),                   intent(in)    :: t     !! Current simulation time
-      end subroutine symba_io_stop_encounter
+      end subroutine symba_io_dump_encounter
 
       module subroutine symba_io_write_discard(self, param)
          use swiftest_classes, only : swiftest_parameters

src/setup/setup.f90

@@ -68,6 +68,14 @@ module subroutine setup_construct_system(system, param)
             allocate(symba_pltpenc :: system%pltpenc_list)
             allocate(symba_plplenc :: system%plplenc_list)
             allocate(symba_plplenc :: system%plplcollision_list)
+
+            select type(param)
+            class is (symba_parameters)
+               if (param%lencounter_save) then
+                  if (.not. allocated(system%encounter_history)) allocate(encounter_storage :: system%encounter_history)
+                  call system%encounter_history%reset()
+               end if
+            end select
          end select
       case (RINGMOONS)
          write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled'

src/symba/symba_io.f90

@@ -11,6 +11,30 @@
    use swiftest
 contains
 
+
+   module subroutine symba_io_dump_encounter(self, param)
+      !! author: David A. Minton
+      !!
+      !! Saves the encounter and/or fragmentation data to file with the name of the current output interval number attached
+      implicit none
+      ! Arguments
+      class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
+      class(symba_parameters),    intent(inout) :: param !! Current run configuration parameters 
+
+      if (self%encounter_history%iframe == 0) return ! No enounters in this interval
+
+      ! Create and save the output file for this encounter
+      self%encounter_history%nc%time_dimsize = maxval(self%encounter_history%tslot(:))
+      write(self%encounter_history%nc%enc_file, '("encounter_",I0.6,".nc")') param%iloop / param%dump_cadence
+      call self%encounter_history%nc%initialize(param)
+      call self%encounter_history%dump(param)
+      call self%encounter_history%nc%close()
+      call self%encounter_history%reset()
+
+      return
+   end subroutine symba_io_dump_encounter
+
+
    module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) 
       !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
       !!
@@ -188,55 +212,6 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms
    end subroutine symba_io_param_writer
 
 
-   module subroutine symba_io_start_encounter(self, param, t)
-      !! author: David A. Minton
-      !!
-      !! Initializes the new encounter and/or fragmentation history
-      implicit none
-      ! Arguments
-      class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
-      class(symba_parameters),    intent(inout) :: param !! Current run configuration parameters 
-      real(DP),                   intent(in)    :: t     !! Current simulation time
-
-      if (.not. allocated(self%encounter_history)) then
-         allocate(encounter_storage :: self%encounter_history)
-      end if
-      call self%encounter_history%reset()
-
-      ! Take the snapshot at the start of the encounter
-      call self%snapshot(param, t) 
-
-      return
-   end subroutine symba_io_start_encounter
-
-
-   module subroutine symba_io_stop_encounter(self, param, t)
-      !! author: David A. Minton
-      !!
-      !! Saves the encounter and/or fragmentation data to file(s)  
-      implicit none
-      ! Arguments
-      class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
-      class(symba_parameters),    intent(inout) :: param !! Current run configuration parameters 
-      real(DP),                   intent(in)    :: t     !! Current simulation time
-      ! Internals
-      integer(I4B) :: i
-
-      ! Create and save the output file for this encounter
-
-      self%encounter_history%nc%time_dimsize = maxval(self%encounter_history%tslot(:))
-
-      write(self%encounter_history%nc%enc_file, '("encounter_",I0.6,".nc")') param%iloop
-
-      call self%encounter_history%nc%initialize(param)
-      call self%encounter_history%dump(param)
-      call self%encounter_history%nc%close()
-      call self%encounter_history%reset()
-
-      return
-   end subroutine symba_io_stop_encounter
-
-
    module subroutine symba_io_write_discard(self, param)
       !! author: David A. Minton
       !!

src/symba/symba_step.f90

@@ -37,9 +37,9 @@ module subroutine symba_step_system(self, param, t, dt)
                call self%reset(param)
                lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0)
                if (lencounter) then
-                  if (param%lencounter_save) call self%start_encounter(param, t)
+                  if (param%lencounter_save) call self%snapshot(param, t) 
                   call self%interp(param, t, dt)
-                  if (param%lencounter_save) call self%stop_encounter(param, t+dt)
+                  if (param%lencounter_save) call self%snapshot(param, t+dt) 
                else
                   self%irec = -1
                   call helio_step_system(self, param, t, dt)

src/symba/symba_util.f90

@@ -1302,6 +1302,7 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t)
       associate(npl => self%pl%nbody,  ntp => self%tp%nbody)
 
          snapshot%t = t
+         snapshot%iloop = param%iloop
 
          if (npl + ntp == 0) return
          npl_snap = npl