From 321da65b476a96668b973ecbc43f76e6e00a4c2f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 22 Aug 2021 20:56:19 -0400 Subject: [PATCH] Put distance limit check back in --- src/kick/kick.f90 | 2 +- src/symba/symba_kick.f90 | 20 ++++++++++++-------- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index c47c9f174..49da7b14c 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -14,7 +14,7 @@ module subroutine kick_getacch_int_pl(self) class(swiftest_pl), intent(inout) :: self ! Internals integer(I8B) :: k, nplpl - real(DP) :: rji2, rlim2 + real(DP) :: rji2 real(DP) :: dx, dy, dz integer(I4B) :: i, j real(DP), dimension(:,:), pointer :: ah, xh diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index e03c12816..68ce0086a 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -21,14 +21,14 @@ module subroutine symba_kick_getacch_int_pl(self) real(DP), dimension(NDIM,self%nbody) :: ahi, ahj integer(I4B), dimension(:,:), pointer :: k_plpl logical, dimension(:), pointer :: lmask - real(DP), dimension(:), pointer :: Gmass + real(DP), dimension(:), pointer :: Gmass, radius - associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass) + associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass, radius => self%radius) nplplm = self%nplplm ahi(:,:) = 0.0_DP ahj(:,:) = 0.0_DP !$omp parallel do default(shared)& - !$omp private(k, i, j, dx, dy, dz, rji2) & + !$omp private(k, i, j, dx, dy, dz, rji2, rlim2) & !$omp reduction(+:ahi) & !$omp reduction(-:ahj) do k = 1_I8B, nplplm @@ -38,7 +38,8 @@ module subroutine symba_kick_getacch_int_pl(self) dy = xh(2, j) - xh(2, i) dz = xh(3, j) - xh(3, i) rji2 = dx**2 + dy**2 + dz**2 - if (lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) + rlim2 = (radius(i)+radius(j))**2 + if ((rji2 > rlim2) .and. lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, 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(:,:) = ah(:,:) + ahi(:,:) + ahj(:,:) @@ -65,21 +66,23 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) ! Internals integer(I4B) :: i, j integer(I8B) :: k, nplplenc - real(DP) :: rji2, dx, dy, dz + real(DP) :: rji2, rlim2, dx, dy, dz 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 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) + associate(pl => self, xh => self%xh, ah => self%ah, Gmass => self%Gmass, plplenc_list => system%plplenc_list, radius => self%radius) call helio_kick_getacch_pl(pl, system, param, t, lbeg) ! Remove accelerations from encountering pairs 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, dx, dy, dz, rji2) & + !$omp private(k, i, j, dx, dy, dz, rji2, rlim2) & !$omp reduction(+:ahi) & !$omp reduction(-:ahj) do k = 1_I8B, nplplenc @@ -89,7 +92,8 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) dy = xh(2, j) - xh(2, i) dz = xh(3, j) - xh(3, i) rji2 = dx**2 + dy**2 + dz**2 - call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) + rlim2 = (radius(i)+radius(j))**2 + if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, dx, dy, dz, 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(:,:) = ah(:,:) - ahi(:,:) - ahj(:,:)