By carefully reconstructing the plplenc_list object in rearray, and k…

…eeping track of recursion levels, I've enabled the ability to resolve fragments at the moment they are identified, and not wait until the recursion steps are done.

examples/symba_mars_disk/mars.in

@@ -1,4 +1,16 @@
 1500        ! Mars System in SI units
+727 1.71022032e+06 2.13948145e+04       ! particle number mass Rhill
+1.24108926e+04                    !particle radius in m
+-8.12608230e+06 -4.37306608e+06 -9.62736144e+03      ! x y z 
+9.87984575e+02 -1.88769371e+03 1.06882012e+01       ! vx vy vz 
+0.4 0.4 0.4 ! Ip
+0.0 0.0 0.0 ! rot
+231 6.25152932e+05 1.58916481e+04       ! particle number mass Rhill
+8.87389776e+03                    !particle radius in m
+-8.21586374e+06 -4.28792953e+06 2.41010139e+04      ! x y z 
+1.01581225e+03 -1.90933511e+03 -2.60449634e+00       ! vx vy vz 
+0.4 0.4 0.4 ! Ip
+0.0 0.0 0.0 ! rot
 2 9.90685589e+04 8.35558297e+03       ! particle number mass Rhill
 7.07643092e+03                    !particle radius in m
 -2.35807426e+06 8.60445552e+06 1.25224401e+04      ! x y z 
@@ -1373,12 +1385,6 @@
 -1.24261883e+03 1.71209694e+03 -6.95777672e+00       ! vx vy vz 
 0.4 0.4 0.4 ! Ip
 0.0 0.0 0.0 ! rot
-231 6.25152932e+05 1.58916481e+04       ! particle number mass Rhill
-8.87389776e+03                    !particle radius in m
--8.21586374e+06 -4.28792953e+06 2.41010139e+04      ! x y z 
-1.01581225e+03 -1.90933511e+03 -2.60449634e+00       ! vx vy vz 
-0.4 0.4 0.4 ! Ip
-0.0 0.0 0.0 ! rot
 232 1.02687634e+05 1.54462392e+04       ! particle number mass Rhill
 4.85987440e+03                    !particle radius in m
 1.62736579e+07 2.82256969e+06 3.66384128e+04      ! x y z 
@@ -4349,12 +4355,6 @@
 -1.87648989e+03 1.17295601e+03 -5.20044045e-01       ! vx vy vz 
 0.4 0.4 0.4 ! Ip
 0.0 0.0 0.0 ! rot
-727 1.71022032e+06 2.13948145e+04       ! particle number mass Rhill
-1.24108926e+04                    !particle radius in m
--8.12608230e+06 -4.37306608e+06 -9.62736144e+03      ! x y z 
-9.87984575e+02 -1.88769371e+03 1.06882012e+01       ! vx vy vz 
-0.4 0.4 0.4 ! Ip
-0.0 0.0 0.0 ! rot
 728 7.56089690e+05 4.18623276e+04       ! particle number mass Rhill
 9.45460600e+03                    !particle radius in m
 -2.34951111e+05 2.32053308e+07 -5.22036528e+04      ! x y z 

examples/symba_mars_disk/param.in

@@ -22,9 +22,9 @@ CHK_QMIN_RANGE 3389500.0 338950000000.0
 EXTRA_FORCE    no
 BIG_DISCARD    no
 RHILL_PRESENT  yes
-GMTINY          1000.0
+GMTINY         1000.0
 ENERGY         yes
-FRAGMENTATION  yes
+FRAGMENTATION  no  
 ROTATION       yes
 MU2KG          1.0
 DU2M           1.0

examples/symba_swifter_comparison/1pl_1pl_encounter/param.swifter.in

@@ -24,3 +24,4 @@ ENC_OUT        enc.swifter.dat
 EXTRA_FORCE    no
 BIG_DISCARD    no
 RHILL_PRESENT  yes
+ENERGY         yes

src/modules/swiftest_classes.f90

@@ -317,6 +317,8 @@ module swiftest_classes
       integer(I4B), dimension(:),   allocatable :: status !! status of the interaction
       integer(I4B), dimension(:),   allocatable :: index1 !! position of the first body in the encounter
       integer(I4B), dimension(:),   allocatable :: index2 !! position of the second body in the encounter
+      integer(I4B), dimension(:),   allocatable :: id1    !! id of the first body in the encounter
+      integer(I4B), dimension(:),   allocatable :: id2    !! id of the second body in the encounter
       real(DP),     dimension(:,:), allocatable :: x1     !! the position of body 1 in the encounter
       real(DP),     dimension(:,:), allocatable :: x2     !! the position of body 2 in the encounter
       real(DP),     dimension(:,:), allocatable :: v1     !! the velocity of body 1 in the encounter

src/modules/symba_classes.f90

@@ -711,7 +711,7 @@ end subroutine symba_util_spill_pl
       module subroutine symba_util_spill_encounter(self, discards, lspill_list, ldestructive)
          use swiftest_classes, only : swiftest_encounter
          implicit none
-         class(symba_encounter),      intent(inout) :: self         !! SyMBA pl-tp encounter list
+         class(symba_encounter),    intent(inout) :: self         !! SyMBA pl-tp encounter list
          class(swiftest_encounter), intent(inout) :: discards     !! Discarded object 
          logical, dimension(:),     intent(in)    :: lspill_list  !! Logical array of bodies to spill into the discards
          logical,                   intent(in)    :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list

src/setup/setup.f90

@@ -88,6 +88,8 @@ module subroutine setup_encounter(self, n)
       if (allocated(self%status)) deallocate(self%status)
       if (allocated(self%index1)) deallocate(self%index1)
       if (allocated(self%index2)) deallocate(self%index2)
+      if (allocated(self%id1)) deallocate(self%id1)
+      if (allocated(self%id2)) deallocate(self%id2)
       if (allocated(self%x1)) deallocate(self%x1)
       if (allocated(self%x2)) deallocate(self%x2)
       if (allocated(self%v1)) deallocate(self%v1)
@@ -98,6 +100,8 @@ module subroutine setup_encounter(self, n)
       allocate(self%status(n))
       allocate(self%index1(n))
       allocate(self%index2(n))
+      allocate(self%id1(n))
+      allocate(self%id2(n))
       allocate(self%x1(NDIM,n))
       allocate(self%x2(NDIM,n))
       allocate(self%v1(NDIM,n))
@@ -108,6 +112,8 @@ module subroutine setup_encounter(self, n)
       self%status(:) = INACTIVE
       self%index1(:) = 0
       self%index2(:) = 0
+      self%id1(:) = 0
+      self%id2(:) = 0
       self%x1(:,:) = 0.0_DP
       self%x2(:,:) = 0.0_DP
       self%v1(:,:) = 0.0_DP

src/symba/symba_encounter_check.f90

@@ -43,6 +43,8 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
                plplenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(:), lencounter(:))
                plplenc_list%index1(1:nenc) = pack(pl%k_plpl(1,:), lencounter(:))
                plplenc_list%index2(1:nenc) = pack(pl%k_plpl(2,:), lencounter(:))
+               plplenc_list%id1(1:nenc) = pl%id(plplenc_list%index1(1:nenc))
+               plplenc_list%id2(1:nenc) = pl%id(plplenc_list%index2(1:nenc))
                do k = 1, nenc
                   plplenc_list%status(k) = ACTIVE
                   plplenc_list%level(k) = irec
@@ -178,6 +180,8 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc
                pltpenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(:,:), lencounter(:,:))
                pltpenc_list%index1(1:nenc) = pack(spread([(i, i = 1, npl)], dim=1, ncopies=ntp), lencounter(:,:)) 
                pltpenc_list%index2(1:nenc) = pack(spread([(i, i = 1, ntp)], dim=2, ncopies=npl), lencounter(:,:))
+               pltpenc_list%id1(1:nenc) = pl%id(pltpenc_list%index1(1:nenc))
+               pltpenc_list%id2(1:nenc) = tp%id(pltpenc_list%index2(1:nenc))
                select type(pl)
                class is (symba_pl)
                   pl%lencounter(:) = .false.

src/symba/symba_util.f90

@@ -395,12 +395,13 @@ module subroutine symba_util_rearray_pl(self, system, param)
       class(symba_parameters),   intent(in)    :: param  !! Current run configuration parameters
       ! Internals
       class(symba_pl), allocatable :: tmp !! The discarded body list.
-      integer(I4B) :: i
+      integer(I4B) :: i, j, k
       logical, dimension(:), allocatable :: lmask
       class(symba_plplenc), allocatable :: plplenc_old
       logical :: lencounter
 
       associate(pl => self, pl_adds => system%pl_adds)
+
          allocate(tmp, mold=pl)
          ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere
          allocate(lmask, source=pl%ldiscard(:))
@@ -424,15 +425,47 @@ module subroutine symba_util_rearray_pl(self, system, param)
 
             allocate(lmask(pl%nbody)) 
             lmask(:) = pl%status(1:pl%nbody) == NEW_PARTICLE
+
             call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, pl%nbody)], lmask))
             where(pl%status(:) /= INACTIVE) pl%status(:) = ACTIVE
 
             call pl%sort("mass", ascending=.false.)
             pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
             pl%nplm = count(pl%lmtiny(:))
 
-            ! Reindex
+            ! Reindex the bodies and calculate the level 0 encounter list
             call pl%eucl_index()
+            lencounter = pl%encounter_check(system, param%dt, 0) 
+            select type(tp => system%tp)
+            class is (symba_tp)
+               lencounter = tp%encounter_check(system, param%dt, 0)
+            end select
+
+            associate(idnew1 => system%plplenc_list%id1, idnew2 => system%plplenc_list%id2, idold1 => plplenc_old%id1, idold2 => plplenc_old%id2)
+               do k = 1, system%plplenc_list%nenc
+                  if ((idnew1(k) == idold1(k)) .and. (idnew2(k) == idold2(k))) then
+                     ! This is an encounter we already know about, so save the old information
+                     system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) 
+                     system%plplenc_list%status(k) = plplenc_old%status(k) 
+                     system%plplenc_list%x1(:,k) = plplenc_old%x1(:,k)
+                     system%plplenc_list%x2(:,k) = plplenc_old%x2(:,k)
+                     system%plplenc_list%v1(:,k) = plplenc_old%v1(:,k)
+                     system%plplenc_list%v2(:,k) = plplenc_old%v2(:,k)
+                     system%plplenc_list%t(k) = plplenc_old%t(k)
+                     system%plplenc_list%level(k) = plplenc_old%level(k)
+                  else if (((idnew1(k) == idold2(k)) .and. (idnew2(k) == idold1(k)))) then
+                     ! This is an encounter we already know about, but with the order reversed, so save the old information
+                     system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) 
+                     system%plplenc_list%status(k) = plplenc_old%status(k) 
+                     system%plplenc_list%x1(:,k) = plplenc_old%x2(:,k)
+                     system%plplenc_list%x2(:,k) = plplenc_old%x1(:,k)
+                     system%plplenc_list%v1(:,k) = plplenc_old%v2(:,k)
+                     system%plplenc_list%v2(:,k) = plplenc_old%v1(:,k)
+                     system%plplenc_list%t(k) = plplenc_old%t(k)
+                     system%plplenc_list%level(k) = plplenc_old%level(k)
+                  end if
+               end do
+            end associate
 
          end if
 
@@ -956,4 +989,4 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive)
       return
    end subroutine symba_util_spill_tp
 
-end submodule s_symba_util
+end submodule s_symba_util

src/util/util_append.f90

@@ -200,6 +200,8 @@ module subroutine util_append_encounter(self, source, lsource_mask)
          call util_append(self%status, source%status, nold, nsrc, lsource_mask)
          call util_append(self%index1, source%index1, nold, nsrc, lsource_mask)
          call util_append(self%index2, source%index2, nold, nsrc, lsource_mask)
+         call util_append(self%id1, source%id1, nold, nsrc, lsource_mask)
+         call util_append(self%id2, source%id2, nold, nsrc, lsource_mask)
          call util_append(self%x1, source%x1, nold, nsrc, lsource_mask)
          call util_append(self%x2, source%x2, nold, nsrc, lsource_mask)
          call util_append(self%v1, source%v1, nold, nsrc, lsource_mask)

src/util/util_copy.f90

@@ -17,6 +17,8 @@ module subroutine util_copy_encounter(self, source)
       self%status(1:n) = source%status(1:n) 
       self%index1(1:n) = source%index1(1:n)
       self%index2(1:n) = source%index2(1:n)
+      self%id1(1:n) = source%id1(1:n)
+      self%id2(1:n) = source%id2(1:n)
       self%x1(:,1:n) = source%x1(:,1:n)
       self%x2(:,1:n) = source%x2(:,1:n)
       self%v1(:,1:n) = source%v1(:,1:n)

src/util/util_spill.f90

@@ -260,6 +260,8 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive
          call util_spill(keeps%status, discards%status, lspill_list, ldestructive)
          call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive)
          call util_spill(keeps%index2, discards%index2, lspill_list, ldestructive)
+         call util_spill(keeps%id1, discards%id1, lspill_list, ldestructive)
+         call util_spill(keeps%id2, discards%id2, lspill_list, ldestructive)
          call util_spill(keeps%x1, discards%x1, lspill_list, ldestructive)
          call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive)
          call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive)