Fixed bug causing encounter history frames to not be stored properly

examples/Fragmentation/Fragmentation_Movie.py

@@ -127,7 +127,7 @@ def encounter_combiner(sim):
     ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time")
 
     # Interpolate in time to make a smooth, constant time step dataset 
-    smooth_time =  np.linspace(start=tgood.isel(time=0), stop=tgood.isel(time=-1), num=2400)
+    smooth_time =  np.linspace(start=tgood.isel(time=0), stop=ds.time[-1], num=2400)
     ds = ds.interp(time=smooth_time)
 
     return ds

python/swiftest/swiftest/simulation_class.py

@@ -2775,6 +2775,9 @@ def _preprocess(ds, param):
         # Remove any overlapping time values
         tgood,tid = np.unique(self.encounters.time,return_index=True)
         self.encounters = self.encounters.isel(time=tid)
+        # Remove any NaN values
+        tgood=self.encounters.time.where(~np.isnan(self.encounters.time),drop=True)
+        self.encounters = self.encounters.sel(time=tgood)
 
         return
 

src/collision/collision_check.f90

@@ -79,9 +79,7 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l
       real(DP), dimension(NDIM)                 :: xr, vr
       integer(I4B)                              :: i, j, k, nenc
       real(DP)                                  :: rlim, Gmtot
-      logical                                   :: isplpl, lany_closest
-      character(len=STRMAX)                     :: timestr, idstri, idstrj, message
-      class(collision_list_plpl), allocatable   :: tmp       
+      logical                                   :: lany_closest
 
       lany_collision = .false.
       if (self%nenc == 0) return

src/encounter/encounter_io.f90

@@ -17,8 +17,8 @@ module subroutine encounter_io_netcdf_dump(self, param)
       !! Dumps the time history of an encounter to file.
       implicit none
       ! Arguments
-      class(encounter_storage(*)),  intent(inout)        :: self   !! Encounter storage object
-      class(base_parameters),   intent(inout)        :: param  !! Current run configuration parameters 
+      class(encounter_storage(*)), intent(inout) :: self   !! Encounter storage object
+      class(base_parameters),      intent(inout) :: param  !! Current run configuration parameters 
       ! Internals
       integer(I4B) :: i
 
@@ -33,7 +33,7 @@ module subroutine encounter_io_netcdf_dump(self, param)
             write(nc%file_name, '("encounter_",I0.6,".nc")') nc%file_number
             call nc%initialize(param)
 
-            do i = 1, self%nframes
+            do i = 1, self%iframe
                if (allocated(self%frame(i)%item)) then
                   select type(snapshot => self%frame(i)%item)
                   class is (encounter_snapshot)

src/encounter/encounter_util.f90

@@ -332,8 +332,6 @@ module subroutine encounter_util_setup_list(self, n)
       ! Arguments
       class(encounter_list), intent(inout) :: self !! Swiftest encounter structure
       integer(I8B),          intent(in)    :: n    !! Number of encounters to allocate space for
-      ! Internals
-      integer(I8B) :: i
 
       if (n < 0) return
       call self%dealloc()
@@ -450,7 +448,6 @@ subroutine encounter_util_save_snapshot(encounter_history, snapshot)
          end do
          deallocate(encounter_history)
          call move_alloc(tmp,encounter_history)
-         nnew = nbig
       end if
 
       ! Find out which time slot this belongs in by searching for an existing slot

src/fraggle/fraggle_generate.f90

@@ -131,7 +131,6 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
 
          lfailure_local = .false.
 
-         call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message)))
 
          ! Use the disruption collision model to generate initial conditions
          ! Compute the "before" energy/momentum and the budgets
@@ -140,6 +139,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
          do try = 1, MAXTRY
             self%fail_scale = (fail_scale_initial)**try
             write(message,*) try
+            call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message)))
             call fraggle_generate_pos_vec(self)
             call fraggle_generate_rot_vec(self)
             call fraggle_generate_vel_vec(self,lfailure_local)
@@ -187,7 +187,6 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t)
       logical                                 :: lpure 
       character(len=STRMAX) :: message
 
-
       select type(nbody_system)
       class is (swiftest_nbody_system)
       select type(pl => nbody_system%pl)
@@ -250,29 +249,26 @@ module subroutine fraggle_generate_pos_vec(collider)
       ! Internals
       real(DP)  :: dis
       real(DP), dimension(NDIM,2) :: fragment_cloud_center
-      real(DP), dimension(2) :: fragment_cloud_radius
+      real(DP), dimension(2) :: fragment_cloud_radius, rdistance
       logical, dimension(collider%fragments%nbody) :: loverlap
       integer(I4B) :: i, j, loop
       logical :: lcat, lhitandrun
       integer(I4B), parameter :: MAXLOOP = 10000
-      real(DP) :: rdistance
+      real(DP), parameter :: rdistance_scale_factor = 0.01_DP ! Scale factor to apply to distance scaling in the event of a fail
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
          lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) 
          lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) 
 
          ! We will treat the first two fragments of the list as special cases. 
          ! Place the first two bodies at the centers of the two fragment clouds, but be sure they are sufficiently far apart to avoid overlap
-         rdistance = .mag. (impactors%rc(:,2) - impactors%rc(:,1)) - sum(fragments%radius(1:2))
-         rdistance = 0.5_DP*rdistance
-
          fragment_cloud_radius(:) = impactors%radius(:)
-
          loverlap(:) = .true.
+         rdistance(:) = rdistance_scale_factor * .mag.impactors%vc(:,:)
          do loop = 1, MAXLOOP
             if (.not.any(loverlap(:))) exit
-            fragment_cloud_center(:,1) = impactors%rc(:,1) !+ rdistance * impactors%bounce_unit(:)
-            fragment_cloud_center(:,2) = impactors%rc(:,2) !- rdistance * impactors%bounce_unit(:)
+            fragment_cloud_center(:,1) = impactors%rc(:,1) - rdistance(1) * impactors%bounce_unit(:)
+            fragment_cloud_center(:,2) = impactors%rc(:,2) + rdistance(2) * impactors%bounce_unit(:)
             do concurrent(i = 1:nfrag, loverlap(i))
                if (i < 3) then
                   fragments%rc(:,i) = fragment_cloud_center(:,i)
@@ -301,7 +297,7 @@ module subroutine fraggle_generate_pos_vec(collider)
                   loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) 
                end do
             end do
-            rdistance = rdistance * collider%fail_scale
+            rdistance(:) = rdistance(:) * collider%fail_scale
             fragment_cloud_radius(:) = fragment_cloud_radius(:) * collider%fail_scale
          end do
 

src/symba/symba_step.f90

@@ -34,19 +34,17 @@ module subroutine symba_step_system(self, param, t, dt)
       class is (symba_tp)
       select type(param)
       class is (swiftest_parameters)
-         associate(encounter_history => self%encounter_history)
-            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%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t, "trajectory") 
-               call self%interp(param, t, dt)
-               if (param%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t+dt, "trajectory") 
-            else
-               self%irec = -1
-               call helio_step_system(self, param, t, dt)
-            end if
-            param%lfirstkick = pl%lfirst
-         end associate
+         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%lenc_save_trajectory) call self%encounter_history%take_snapshot(param, self, t, "trajectory") 
+            call self%interp(param, t, dt)
+            if (param%lenc_save_trajectory) call self%encounter_history%take_snapshot(param, self, t+dt, "trajectory") 
+         else
+            self%irec = -1
+            call helio_step_system(self, param, t, dt)
+         end if
+         param%lfirstkick = pl%lfirst
       end select
       end select
       end select
@@ -192,7 +190,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci)
       class is (symba_list_plpl)
       select type(pltp_encounter => self%pltp_encounter)
       class is (symba_list_pltp)
-         associate(nbody_system => self, lplpl_collision => plpl_encounter%lcollision, lpltp_collision => pltp_encounter%lcollision,  encounter_history => self%encounter_history)
+         associate(nbody_system => self, lplpl_collision => plpl_encounter%lcollision, lpltp_collision => pltp_encounter%lcollision)
             nbody_system%irec = ireci
             dtl = param%dt / (NTENC**ireci)
             dth = 0.5_DP * dtl
@@ -249,7 +247,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci)
                   if (lplpl_collision) call plpl_encounter%resolve_collision(nbody_system, param, t+j*dtl, dtl, ireci)
                   if (lpltp_collision) call pltp_encounter%resolve_collision(nbody_system, param, t+j*dtl, dtl, ireci)
                end if
-               if (param%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t+j*dtl, "trajectory") 
+               if (param%lenc_save_trajectory) call self%encounter_history%take_snapshot(param, self, t+j*dtl, "trajectory") 
 
                call self%set_recur_levels(ireci)