diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 19c7967bd..f93f8d714 100644 --- a/src/modules/swiftest_classes.f90 +++ b/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 diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 4a253a263..070308838 100644 --- a/src/setup/setup.f90 +++ b/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 diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 5d95ba44f..3672f959b 100644 --- a/src/symba/symba_collision.f90 +++ b/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 diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 42b67d079..c1b83f9dd 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/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 diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index af440eada..9a5b11229 100644 --- a/src/symba/symba_kick.f90 +++ b/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 diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index c8407416d..73195df91 100644 --- a/src/util/util_copy.f90 +++ b/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) diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 77e00cf66..ccc494474 100644 --- a/src/util/util_spill.f90 +++ b/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)