Improved the stability of the encounter snapshot

examples/Fragmentation/Fragmentation_Movie.py

@@ -41,8 +41,8 @@
 movie_titles = dict(zip(available_movie_styles, movie_title_list))
 
 # These initial conditions were generated by trial and error
-pos_vectors = {"disruption_headon"         : [np.array([1.0, -5.0e-05, 0.0]),
-                                              np.array([1.0,  5.0e-05 ,0.0])],
+pos_vectors = {"disruption_headon"         : [np.array([1.0, -5.0e-04, 0.0]),
+                                              np.array([1.0,  5.0e-04 ,0.0])],
                "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05,      0.0]),
                                               np.array([1.0,  4.2e-05,      0.0])],
                "hitandrun"                 : [np.array([1.0, -2.0e-05,      0.0]),
@@ -203,7 +203,7 @@ def data_stream(self, frame=0):
         minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body
         gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades
         sim.set_parameter(fragmentation=True, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
-        sim.run(dt=1e-5, tstop=4.0e-3, istep_out=10, dump_cadence=0)
+        sim.run(dt=1e-4, tstop=1.0e-4, istep_out=1, dump_cadence=0)
 
         print("Generating animation")
         anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1)

src/modules/symba_classes.f90

@@ -229,11 +229,15 @@ module symba_classes
    end type symba_nbody_system
 
 
-   type, extends(symba_nbody_system) :: symba_encounter_snapshot
-      integer(I4B)                   :: tslot !! The index for the time array in the final NetCDF file
+   type :: symba_encounter_snapshot
+      type(symba_pl) :: pl         !! Massive body data structure
+      type(symba_tp) :: tp         !! Test particle data structure
+      real(DP)       :: t = 0.0_DP !! Time at the snapshot 
+      integer(I4B)   :: tslot = 0  !! The index for the time array in the final NetCDF file
    contains
-      procedure :: write_encounter_frame => symba_io_encounter_write_frame    !! Writes a frame of encounter data to file 
-      generic   :: write_frame           => write_encounter_frame
+      procedure :: write_encounter_frame => symba_io_encounter_write_frame !! Writes a frame of encounter data to file 
+      procedure :: dealloc               => symba_util_dealloc_snapshot    !! Deallocates all allocatable arrays
+      generic   :: write_frame           => write_encounter_frame          !! Writes a snaphot frame to file
       final     :: symba_util_final_encounter_snapshot
    end type symba_encounter_snapshot
 
@@ -657,6 +661,11 @@ module subroutine symba_util_dealloc_merger(self)
          class(symba_merger),  intent(inout) :: self !! SyMBA body merger object
       end subroutine symba_util_dealloc_merger
 
+      module subroutine symba_util_dealloc_snapshot(self)
+         implicit none
+         class(symba_encounter_snapshot),  intent(inout) :: self !! SyMBA nbody system object
+      end subroutine symba_util_dealloc_snapshot
+
       module subroutine symba_util_dealloc_system(self)
          implicit none
          class(symba_nbody_system),  intent(inout) :: self !! SyMBA nbody system object

src/symba/symba_io.f90

@@ -95,7 +95,8 @@ module subroutine symba_io_encounter_initialize_output(self, param)
          call check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%id_dimid], nc%ptype_varid), "symba_io_encounter_initialize_output 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), "symba_io_encounter_initialize_output 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), "symba_io_encounter_initialize_output 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), "symba_io_encounter_initialize_output nf90_def_var Gmass_varid"  )
+         call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "symba_io_encounter_initialize_output 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), "symba_io_encounter_initialize_output nf90_def_var level_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), "symba_io_encounter_initialize_output nf90_def_var radius_varid"  )
          end if
@@ -153,15 +154,15 @@ module subroutine symba_io_encounter_write_frame(self, nc, param)
       tslot = self%tslot
       call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "symba_io_encounter_write_frame nf90_put_var time_varid"  )
 
-      select type(pl => self%pl)
-      class is (symba_pl)
+      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]), "symba_io_encounter_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]), "symba_io_encounter_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]), "symba_io_encounter_write_frame nf90_put_var pl vh_varid"  )
             call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl Gmass_varid"  )
+            call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl level_varid"  )
 
             if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "symba_io_encounter_write_frame nf90_put_var pl radius_varid"  )
 
@@ -175,10 +176,7 @@ module subroutine symba_io_encounter_write_frame(self, nc, param)
             charstring = trim(adjustl(pl%info(i)%particle_type))
             call check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "symba_io_encounter_write_frame nf90_put_var pl particle_type_varid"  )
          end do
-
-      end select
 
-      associate(tp => self%tp) 
          ntp = tp%nbody
          do i = 1, ntp
             idslot = tp%id(i)

src/symba/symba_util.f90

@@ -238,10 +238,28 @@ module subroutine symba_util_dealloc_merger(self)
    end subroutine symba_util_dealloc_merger
 
 
+   module subroutine symba_util_dealloc_snapshot(self)
+      !! author: David A. Minton
+      !!
+      !! Deallocates allocatable arrays in an encounter snapshot
+      implicit none
+      ! Arguments
+      class(symba_encounter_snapshot),  intent(inout) :: self !! SyMBA nbody system object
+
+
+      call self%pl%dealloc()
+      call self%tp%dealloc()
+      self%t = 0.0_DP
+      self%tslot = 0
+
+      return
+   end subroutine symba_util_dealloc_snapshot
+
+
    module subroutine symba_util_dealloc_system(self)
       !! author: David A. Minton
       !!
-      !! Deallocates all allocatabale arrays
+      !! Deallocates all allocatabale arrays in a SyMBA system object
       implicit none
       ! Arguments
       class(symba_nbody_system),  intent(inout) :: self !! SyMBA nbody system object
@@ -937,6 +955,8 @@ module subroutine symba_util_resize_storage(self, nnew)
             nbig = nbig * 2
          end do
          allocate(symba_encounter_storage(nbig) :: tmp) 
+         tmp%tvals(1:nold) = self%encounter_history%tvals(1:nold)
+         tmp%tvals(nold+1:nbig) = huge(1.0_DP)
          if (lmalloc) then
             do i = 1, nold
                if (allocated(self%encounter_history%frame(i)%item)) call move_alloc(self%encounter_history%frame(i)%item, tmp%frame(i)%item)
@@ -1328,8 +1348,6 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t)
          allocate(symba_encounter_snapshot :: snapshot)
          snapshot%t = t
 
-         allocate(symba_pl :: snapshot%pl)
-         allocate(symba_tp :: snapshot%tp)
          if (npl + ntp == 0) return
          npl_snap = npl
          ntp_snap = ntp
@@ -1353,6 +1371,8 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t)
                   allocate(snapshot%pl%id(npl_snap))
                   allocate(snapshot%pl%info(npl_snap))
                   allocate(snapshot%pl%Gmass(npl_snap))
+                  allocate(snapshot%pl%levelg(npl_snap))
+                  snapshot%pl%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl))
                   snapshot%pl%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl))
                   snapshot%pl%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl))
                   snapshot%pl%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl))
@@ -1394,15 +1414,15 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t)
                end if
 
                ! Save the snapshot
-               self%encounter_history%iframe = self%encounter_history%iframe + 1
-               call self%resize_storage(self%encounter_history%iframe)
 
                ! 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, self%encounter_history%nframes
                   if (t <= self%encounter_history%tvals(i)) then
                      snapshot%tslot = i
                      self%encounter_history%tvals(i) = t
+                     self%encounter_history%iframe = i
+                     call self%resize_storage(i)
                      exit
                   end if
                end do