Streamlined kick step in SyMBA by saving the k_plpl index values of e…

…ach encoutering pair. Then the subtraction step can re-use the main parallelized interaction subroutine without needing to duplicate effort.

src/modules/swiftest_classes.f90

@@ -318,6 +318,7 @@ module swiftest_classes
       integer(I4B)                              :: nenc   !! Total number of encounters
       logical,      dimension(:),   allocatable :: lvdotr !! relative vdotr flag
       integer(I4B), dimension(:),   allocatable :: status !! status of the interaction
+      integer(I8B), dimension(:),   allocatable :: kidx   !! index value of the encounter from the master k_plpl encounter list
       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
@@ -1412,6 +1413,14 @@ module subroutine util_spill_arr_I4B(keeps, 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_arr_I4B
 
+      module subroutine util_spill_arr_I8B(keeps, discards, lspill_list, ldestructive)
+         implicit none
+         integer(I8B), dimension(:), allocatable, intent(inout) :: keeps        !! Array of values to keep 
+         integer(I8B), dimension(:), allocatable, intent(inout) :: discards     !! Array of discards
+         logical,      dimension(:),              intent(in)    :: lspill_list  !! Logical array of bodies to spill into the discardss
+         logical,                                 intent(in)    :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not
+      end subroutine util_spill_arr_I8B
+
       module subroutine util_spill_arr_logical(keeps, discards, lspill_list, ldestructive)
          implicit none
          logical, dimension(:), allocatable, intent(inout) :: keeps        !! Array of values to keep 

src/setup/setup.f90

@@ -86,6 +86,7 @@ module subroutine setup_encounter(self, n)
 
       if (allocated(self%lvdotr)) deallocate(self%lvdotr)
       if (allocated(self%status)) deallocate(self%status)
+      if (allocated(self%kidx)) deallocate(self%kidx)
       if (allocated(self%index1)) deallocate(self%index1)
       if (allocated(self%index2)) deallocate(self%index2)
       if (allocated(self%id1)) deallocate(self%id1)
@@ -98,6 +99,7 @@ module subroutine setup_encounter(self, n)
 
       allocate(self%lvdotr(n))
       allocate(self%status(n))
+      allocate(self%kidx(n))
       allocate(self%index1(n))
       allocate(self%index2(n))
       allocate(self%id1(n))
@@ -110,6 +112,7 @@ module subroutine setup_encounter(self, n)
 
       self%lvdotr(:) = .false.
       self%status(:) = INACTIVE
+      self%kidx(:) = 0_I8B
       self%index1(:) = 0
       self%index2(:) = 0
       self%id1(:) = 0

src/symba/symba_collision.f90

@@ -300,11 +300,9 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad
             end do
          end do
 
-
          status = MERGED
          call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status)
 
-
       end select
 
       return 

src/symba/symba_encounter_check.f90

@@ -54,8 +54,9 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
             associate(plplenc_list => system%plplenc_list)
                call plplenc_list%resize(nenc)
                plplenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(1:nplplm), lencounter(1:nplplm))
-               plplenc_list%index1(1:nenc) = pack(pl%k_plpl(1,1:nplplm), lencounter(1:nplplm))
-               plplenc_list%index2(1:nenc) = pack(pl%k_plpl(2,1:nplplm), lencounter(1:nplplm))
+               plplenc_list%kidx(1:nenc) = pack([(k, k = 1_I8B, nplplm)], lencounter(1:nplplm))
+               plplenc_list%index1(1:nenc) = pl%k_plpl(1,plplenc_list%kidx(1:nenc))
+               plplenc_list%index2(1:nenc) = pl%k_plpl(2,plplenc_list%kidx(1:nenc))
                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

src/symba/symba_kick.f90

@@ -37,36 +37,24 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
       integer(I4B)              :: i, j
       integer(I8B)              :: k, nplplenc
       real(DP)                  :: rji2, rlim2, xr, yr, zr
-      real(DP), dimension(NDIM,self%nbody) :: ahi, ahj
-      class(symba_plplenc), pointer :: plplenc_list
-      real(DP), dimension(:), pointer :: Gmass, radius
-      real(DP), dimension(:,:), pointer :: ah, xh
+      real(DP), dimension(NDIM,self%nbody) :: ah_enc
+      integer(I4B), dimension(:,:), allocatable :: k_plpl_enc
 
       if (self%nbody == 0) return
       select type(system)
       class is (symba_nbody_system)
-         associate(pl => self, xh => self%xh, ah => self%ah, Gmass => self%Gmass, plplenc_list => system%plplenc_list, radius => self%radius)
+         associate(pl => self, npl => self%nbody, plplenc_list => system%plplenc_list, radius => self%radius)
+            ! Apply kicks to all bodies (including those in the encounter list)
             call helio_kick_getacch_pl(pl, system, param, t, lbeg)
-            ! Remove accelerations from encountering pairs
+
+            ! Remove kicks from bodies involved currently in the encounter list, as these are dealt with separately.
             nplplenc = int(plplenc_list%nenc, kind=I8B)
-            ahi(:,:) = 0.0_DP
-            ahj(:,:) = 0.0_DP
-            !$omp parallel do default(shared)&
-            !$omp private(k, i, j, xr, yr, zr, rji2, rlim2)  &
-            !$omp reduction(+:ahi) &
-            !$omp reduction(-:ahj) 
-            do k = 1_I8B, nplplenc
-               i = plplenc_list%index1(k)
-               j = plplenc_list%index2(k)
-               xr = xh(1, j) - xh(1, i)
-               yr = xh(2, j) - xh(2, i)
-               zr = xh(3, j) - xh(3, i)
-               rji2 = xr**2 + yr**2 + zr**2
-               rlim2 = (radius(i)+radius(j))**2
-               if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
-            end do
-            !$omp end parallel do
-            ah(:,1:self%nbody) = ah(:,1:self%nbody) - ahi(:,1:self%nbody) - ahj(:,1:self%nbody)
+            allocate(k_plpl_enc(2,nplplenc))
+            k_plpl_enc(:,1:nplplenc) = pl%k_plpl(:,plplenc_list%kidx(1:nplplenc))
+            ah_enc(:,:) = 0.0_DP
+            call kick_getacch_int_all_pl(npl, nplplenc, k_plpl_enc, pl%xh, pl%Gmass, pl%radius, ah_enc)
+            pl%ah(:,1:npl) = pl%ah(:,1:npl) - ah_enc(:,1:npl)
+
          end associate
       end select
 

src/util/util_copy.f90

@@ -15,6 +15,7 @@ module subroutine util_copy_encounter(self, source)
       self%nenc = n
       self%lvdotr(1:n) = source%lvdotr(1:n) 
       self%status(1:n) = source%status(1:n) 
+      self%kidx(1:n) = source%kidx(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)

src/util/util_spill.f90

@@ -160,6 +160,45 @@ module subroutine util_spill_arr_I4B(keeps, discards, lspill_list, ldestructive)
 
       return
    end subroutine util_spill_arr_I4B
+
+
+   module subroutine util_spill_arr_I8B(keeps, discards, lspill_list, ldestructive)
+      !! author: David A. Minton
+      !!
+      !! Performs a spill operation on a single array of type I4B
+      !! This is the inverse of a spill operation
+      implicit none
+      ! Arguments
+      integer(I8B), dimension(:), allocatable, intent(inout) :: keeps        !! Array of values to keep 
+      integer(I8B), dimension(:), allocatable, intent(inout) :: discards     !! Array of discards
+      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 the keeps array or not
+      ! Internals
+      integer(I4B) :: nspill, nkeep, nlist
+
+      nkeep = count(.not.lspill_list(:))
+      nspill = count(lspill_list(:))
+      nlist = size(lspill_list(:))
+
+      if (.not.allocated(keeps) .or. nspill == 0) return
+      if (.not.allocated(discards)) then
+         allocate(discards(nspill))
+      else if (size(discards) /= nspill) then
+         deallocate(discards)
+         allocate(discards(nspill))
+      end if
+
+      discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist))
+      if (ldestructive) then
+         if (nkeep > 0) then
+            keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist))
+         else
+            deallocate(keeps)
+         end if
+      end if
+
+      return
+   end subroutine util_spill_arr_I8B
 
 
    module subroutine util_spill_arr_logical(keeps, discards, lspill_list, ldestructive)
@@ -262,6 +301,7 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive
       associate(keeps => self)
          call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive)
          call util_spill(keeps%status, discards%status, lspill_list, ldestructive)
+         call util_spill(keeps%kidx, discards%kidx, 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)