From 4e4712ef049e33c56e377e4dd25ee2d2e62681ee Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 7 Mar 2023 11:41:23 -0500 Subject: [PATCH] Refactored kick and also fixed issue that was causing semi-interacting bodies to not be computed properly --- src/swiftest/swiftest_kick.f90 | 372 ++++++++++++++++++------------- src/swiftest/swiftest_module.f90 | 37 ++- src/symba/symba_kick.f90 | 36 +-- 3 files changed, 260 insertions(+), 185 deletions(-) diff --git a/src/swiftest/swiftest_kick.f90 b/src/swiftest/swiftest_kick.f90 index 87d58ce64..b3a5f7d39 100644 --- a/src/swiftest/swiftest_kick.f90 +++ b/src/swiftest/swiftest_kick.f90 @@ -21,41 +21,22 @@ module subroutine swiftest_kick_getacch_int_pl(self, param) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters ! Internals - ! type(interaction_timer), save :: itimer logical, save :: lfirst = .true. - ! if (param%ladaptive_interactions) then - ! if (self%nplpl > 0) then - ! if (lfirst) then - ! write(itimer%loopname, *) "kick_getacch_int_pl" - ! write(itimer%looptype, *) "INTERACTION" - ! call itimer%time_this_loop(param, self%nplpl, self) - ! lfirst = .false. - ! else - ! if (itimer%netcdf_io_check(param, self%nplpl)) call itimer%time_this_loop(param, self%nplpl, self) - ! end if - ! else - ! param%lflatten_interactions = .false. - ! end if - ! end if - - ! if (param%lflatten_interactions) then - ! if (param%lclose) then - ! call swiftest_kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%rh, self%Gmass, self%radius, self%ah) - ! else - ! call swiftest_kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%rh, self%Gmass, acc=self%ah) - ! end if - ! else + + if (param%lflatten_interactions) then if (param%lclose) then - call swiftest_kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%rh, self%Gmass, self%radius, self%ah) + call swiftest_kick_getacch_int_all(self%nbody, self%nplpl, self%k_plpl, self%rh, self%Gmass, self%radius, self%ah) else - call swiftest_kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%rh, self%Gmass, acc=self%ah) + call swiftest_kick_getacch_int_all(self%nbody, self%nplpl, self%k_plpl, self%rh, self%Gmass, self%ah) end if - ! end if - - ! if (param%ladaptive_interactions .and. self%nplpl > 0) then - ! if (itimer%is_on) call itimer%adapt(param, self%nplpl, self) - ! end if + else + if (param%lclose) then + call swiftest_kick_getacch_int_all(self%nbody, self%nbody, self%rh, self%Gmass, self%radius, self%ah) + else + call swiftest_kick_getacch_int_all(self%nbody, self%nbody, self%rh, self%Gmass, self%ah) + end if + end if return end subroutine swiftest_kick_getacch_int_pl @@ -84,7 +65,7 @@ module subroutine swiftest_kick_getacch_int_tp(self, param, GMpl, rhp, npl) end subroutine swiftest_kick_getacch_int_tp - module subroutine swiftest_kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + module subroutine swiftest_kick_getacch_int_all_flat_rad_pl(npl, nplpl, k_plpl, r, Gmass, radius, acc) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization. @@ -96,64 +77,93 @@ module subroutine swiftest_kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, G integer(I4B), intent(in) :: npl !! Number of massive bodies integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix) - real(DP), dimension(:,:), intent(in) :: x !! Position vector array + real(DP), dimension(:,:), intent(in) :: r !! Position vector array real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass - real(DP), dimension(:), intent(in), optional :: radius !! Array of massive body radii + real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array ! Internals integer(I8B) :: k real(DP), dimension(NDIM,npl) :: ahi, ahj integer(I4B) :: i, j real(DP) :: rji2, rlim2 - real(DP) :: xr, yr, zr + real(DP) :: rx, ry, rz ahi(:,:) = 0.0_DP ahj(:,:) = 0.0_DP - if (present(radius)) then - !$omp parallel do default(private) schedule(static)& - !$omp shared(nplpl, k_plpl, x, Gmass, radius) & - !$omp lastprivate(i, j, rji2, rlim2, xr, yr, zr) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do k = 1_I8B, nplpl - i = k_plpl(1, k) - j = k_plpl(2, k) - xr = x(1, j) - x(1, i) - yr = x(2, j) - x(2, i) - zr = x(3, j) - x(3, i) - rji2 = xr**2 + yr**2 + zr**2 - rlim2 = (radius(i) + radius(j))**2 - if (rji2 > rlim2) call swiftest_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 - else - !$omp parallel do default(private) schedule(static)& - !$omp shared(nplpl, k_plpl, x, Gmass, radius) & - !$omp lastprivate(i, j, rji2, xr, yr, zr) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do k = 1_I8B, nplpl - i = k_plpl(1, k) - j = k_plpl(2, k) - xr = x(1, j) - x(1, i) - yr = x(2, j) - x(2, i) - zr = x(3, j) - x(3, i) - rji2 = xr**2 + yr**2 + zr**2 - call swiftest_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 - end if + !$omp parallel do default(private) schedule(static)& + !$omp shared(nplpl, k_plpl, r, Gmass, radius) & + !$omp lastprivate(i, j, rji2, rlim2, rx, ry, rz) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do k = 1_I8B, nplpl + i = k_plpl(1, k) + j = k_plpl(2, k) + rx = r(1, j) - r(1, i) + ry = r(2, j) - r(2, i) + rz = r(3, j) - r(3, i) + rji2 = rx**2 + ry**2 + rz**2 + rlim2 = (radius(i) + radius(j))**2 + if (rji2 > rlim2) call swiftest_kick_getacch_int_one_pl(rji2, rx, ry, rz, 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 + + acc(:,:) = acc(:,:) + ahi(:,:) + ahj(:,:) + + return + end subroutine swiftest_kick_getacch_int_all_flat_rad_pl + + + module subroutine swiftest_kick_getacch_int_all_flat_norad_pl(npl, nplpl, k_plpl, r, Gmass, acc) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization. + !! This is the flattened (single loop) version that uses the k_plpl interaction pair index array + !! + !! Adapted from Hal Levison's Swift routine getacch_ah3.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9 + implicit none + integer(I4B), intent(in) :: npl !! Number of massive bodies + integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute + integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix) + real(DP), dimension(:,:), intent(in) :: r !! Position vector array + real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + ! Internals + integer(I8B) :: k + real(DP), dimension(NDIM,npl) :: ahi, ahj + integer(I4B) :: i, j + real(DP) :: rji2, rlim2 + real(DP) :: rx, ry, rz + + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + + !$omp parallel do default(private) schedule(static)& + !$omp shared(nplpl, k_plpl, r, Gmass) & + !$omp lastprivate(i, j, rji2, rx, ry, rz) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do k = 1_I8B, nplpl + i = k_plpl(1, k) + j = k_plpl(2, k) + rx = r(1, j) - r(1, i) + ry = r(2, j) - r(2, i) + rz = r(3, j) - r(3, i) + rji2 = rx**2 + ry**2 + rz**2 + call swiftest_kick_getacch_int_one_pl(rji2, rx, ry, rz, 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 acc(:,:) = acc(:,:) + ahi(:,:) + ahj(:,:) return - end subroutine swiftest_kick_getacch_int_all_flat_pl + end subroutine swiftest_kick_getacch_int_all_flat_norad_pl - module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, r, Gmass, radius, acc) + module subroutine swiftest_kick_getacch_int_all_tri_rad_pl(npl, nplm, r, Gmass, radius, acc) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization. @@ -166,7 +176,7 @@ module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, r, Gmas integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies real(DP), dimension(:,:), intent(in) :: r !! Position vector array real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass - real(DP), dimension(:), intent(in), optional :: radius !! Array of massive body radii + real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array ! Internals integer(I4B) :: i, j, nplt @@ -177,92 +187,156 @@ module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, r, Gmas nplt = npl - nplm lmtiny = (nplt > nplm) - if (present(radius)) then - if (lmtiny) then - ahi(:,:) = 0.0_DP - ahj(:,:) = 0.0_DP - !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, nplm, r, Gmass, radius) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do i = 1, nplm - do concurrent(j = i+1:npl) - rx = r(1, j) - r(1, i) - ry = r(2, j) - r(2, i) - rz = r(3, j) - r(3, i) - rji2 = rx**2 + ry**2 + rz**2 - rlim2 = (radius(i) + radius(j))**2 - if (rji2 > rlim2) call swiftest_kick_getacch_int_one_pl(rji2, rx, ry, rz, Gmass(i), Gmass(j), & - ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) - end do - end do - !$omp end parallel do - do concurrent(i = 1:npl) - acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) - end do - else - !$omp parallel do default(private) schedule(static)& - !$omp shared(npl,nplm, r, Gmass, radius, acc) - do i = 1, nplm - do concurrent(j = 1:npl, j/=i) - rx = r(1,j) - r(1,i) - ry = r(2,j) - r(2,i) - rz = r(3,j) - r(3,i) - rji2 = rx**2 + ry**2 + rz**2 - rlim2 = (radius(i) + radius(j))**2 - if (rji2 > rlim2) then - fac = Gmass(j) / (rji2 * sqrt(rji2)) - acc(1,i) = acc(1,i) + fac * rx - acc(2,i) = acc(2,i) + fac * ry - acc(3,i) = acc(3,i) + fac * rz - end if - end do - end do - !$omp end parallel do - end if - else - if (lmtiny) then - ahi(:,:) = 0.0_DP - ahj(:,:) = 0.0_DP - !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, nplm, r, Gmass) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do i = 1, nplm - do concurrent(j = i+1:npl) - rx = r(1, j) - r(1, i) - ry = r(2, j) - r(2, i) - rz = r(3, j) - r(3, i) - rji2 = rx**2 + ry**2 + rz**2 - call swiftest_kick_getacch_int_one_pl(rji2, rx, ry, rz, Gmass(i), Gmass(j), & - ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) - end do + if (lmtiny) then + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, nplm, r, Gmass, radius) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do i = 1, nplm + do concurrent(j = i+1:npl) + rx = r(1, j) - r(1, i) + ry = r(2, j) - r(2, i) + rz = r(3, j) - r(3, i) + rji2 = rx**2 + ry**2 + rz**2 + rlim2 = (radius(i) + radius(j))**2 + if (rji2 > rlim2) call swiftest_kick_getacch_int_one_pl(rji2, rx, ry, rz, 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 - do concurrent(i = 1:npl) - acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) + end do + !$omp end parallel do + do concurrent(i = 1:npl) + acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) + end do + else + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl,nplm, r, Gmass, radius, acc) + do i = 1, nplm + do concurrent(j = 1:npl, j/=i) + rx = r(1,j) - r(1,i) + ry = r(2,j) - r(2,i) + rz = r(3,j) - r(3,i) + rji2 = rx**2 + ry**2 + rz**2 + rlim2 = (radius(i) + radius(j))**2 + if (rji2 > rlim2) then + fac = Gmass(j) / (rji2 * sqrt(rji2)) + acc(1,i) = acc(1,i) + fac * rx + acc(2,i) = acc(2,i) + fac * ry + acc(3,i) = acc(3,i) + fac * rz + end if end do - else - !$omp parallel do default(private) schedule(static)& - !$omp shared(npl,nplm, r, Gmass, acc) - do i = 1, nplm - do concurrent(j = 1:npl, j/=i) - rx = r(1,j) - r(1,i) - ry = r(2,j) - r(2,i) - rz = r(3,j) - r(3,i) - rji2 = rx**2 + ry**2 + rz**2 + end do + !$omp end parallel do + + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl,nplm, r, Gmass, radius, acc) + do i = nplm+1,npl + do concurrent(j = 1:nplm, j/=i) + rx = r(1,j) - r(1,i) + ry = r(2,j) - r(2,i) + rz = r(3,j) - r(3,i) + rji2 = rx**2 + ry**2 + rz**2 + rlim2 = (radius(i) + radius(j))**2 + if (rji2 > rlim2) then fac = Gmass(j) / (rji2 * sqrt(rji2)) acc(1,i) = acc(1,i) + fac * rx acc(2,i) = acc(2,i) + fac * ry acc(3,i) = acc(3,i) + fac * rz - end do + end if end do - !$omp end parallel do - end if + end do + !$omp end parallel do + + end if + + + return + end subroutine swiftest_kick_getacch_int_all_tri_rad_pl + + + module subroutine swiftest_kick_getacch_int_all_tri_norad_pl(npl, nplm, r, Gmass, acc) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization. + !! This is the upper triangular matrix (double loop) version. + !! + !! Adapted from Hal Levison's Swift routine getacch_ah3.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9 + implicit none + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vector array + real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + ! Internals + integer(I4B) :: i, j, nplt + real(DP) :: rji2, rlim2, fac, rx, ry, rz + real(DP), dimension(NDIM,npl) :: ahi, ahj + logical :: lmtiny + + nplt = npl - nplm + lmtiny = (nplt > nplm) + + if (lmtiny) then + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, nplm, r, Gmass) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do i = 1, nplm + do concurrent(j = i+1:npl) + rx = r(1, j) - r(1, i) + ry = r(2, j) - r(2, i) + rz = r(3, j) - r(3, i) + rji2 = rx**2 + ry**2 + rz**2 + call swiftest_kick_getacch_int_one_pl(rji2, rx, ry, rz, Gmass(i), Gmass(j), & + ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) + end do + end do + !$omp end parallel do + do concurrent(i = 1:npl) + acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) + end do + else + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl,nplm, r, Gmass, acc) + do i = 1, nplm + do concurrent(j = 1:npl, j/=i) + rx = r(1,j) - r(1,i) + ry = r(2,j) - r(2,i) + rz = r(3,j) - r(3,i) + rji2 = rx**2 + ry**2 + rz**2 + fac = Gmass(j) / (rji2 * sqrt(rji2)) + acc(1,i) = acc(1,i) + fac * rx + acc(2,i) = acc(2,i) + fac * ry + acc(3,i) = acc(3,i) + fac * rz + end do + end do + !$omp end parallel do + + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl,nplm, r, Gmass, acc) + do i = nplm+1,npl + do concurrent(j = 1:nplm, j/=i) + rx = r(1,j) - r(1,i) + ry = r(2,j) - r(2,i) + rz = r(3,j) - r(3,i) + rji2 = rx**2 + ry**2 + rz**2 + fac = Gmass(j) / (rji2 * sqrt(rji2)) + acc(1,i) = acc(1,i) + fac * rx + acc(2,i) = acc(2,i) + fac * ry + acc(3,i) = acc(3,i) + fac * rz + end do + end do + !$omp end parallel do + end if return - end subroutine swiftest_kick_getacch_int_all_triangular_pl + end subroutine swiftest_kick_getacch_int_all_tri_norad_pl + module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, xtp, rpl, GMpl, lmask, acc) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 78b307977..d360db56b 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -890,27 +890,48 @@ module subroutine swiftest_kick_getacch_int_tp(self, param, GMpl, rhp, npl) real(DP), dimension(:,:), intent(in) :: rhp !! Massive body position vectors integer(I4B), intent(in) :: npl !! Number of active massive bodies end subroutine swiftest_kick_getacch_int_tp + end interface + + interface swiftest_kick_getacch_int_all + module subroutine swiftest_kick_getacch_int_all_flat_rad_pl(npl, nplpl, k_plpl, r, Gmass, radius, acc) + implicit none + integer(I4B), intent(in) :: npl !! Number of massive bodies + integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute + integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix) + real(DP), dimension(:,:), intent(in) :: r !! Position vector array + real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass + real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + end subroutine swiftest_kick_getacch_int_all_flat_rad_pl - module subroutine swiftest_kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + module subroutine swiftest_kick_getacch_int_all_flat_norad_pl(npl, nplpl, k_plpl, r, Gmass, acc) implicit none integer(I4B), intent(in) :: npl !! Number of massive bodies integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix) - real(DP), dimension(:,:), intent(in) :: x !! Position vector array + real(DP), dimension(:,:), intent(in) :: r !! Position vector array real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass - real(DP), dimension(:), intent(in), optional :: radius !! Array of massive body radii real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array - end subroutine swiftest_kick_getacch_int_all_flat_pl + end subroutine swiftest_kick_getacch_int_all_flat_norad_pl - module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, r, Gmass, radius, acc) + module subroutine swiftest_kick_getacch_int_all_tri_rad_pl(npl, nplm, r, Gmass, radius, acc) implicit none integer(I4B), intent(in) :: npl !! Total number of massive bodies integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies real(DP), dimension(:,:), intent(in) :: r !! Position vector array real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass - real(DP), dimension(:), intent(in), optional :: radius !! Array of massive body radii + real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array - end subroutine swiftest_kick_getacch_int_all_triangular_pl + end subroutine swiftest_kick_getacch_int_all_tri_rad_pl + + module subroutine swiftest_kick_getacch_int_all_tri_norad_pl(npl, nplm, r, Gmass, acc) + implicit none + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vector array + real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + end subroutine swiftest_kick_getacch_int_all_tri_norad_pl module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, xtp, rpl, GMpl, lmask, acc) implicit none @@ -922,7 +943,9 @@ module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, xtp, rpl, GMpl, lma logical, dimension(:), intent(in) :: lmask !! Logical mask indicating which test particles should be computed real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array end subroutine swiftest_kick_getacch_int_all_tp + end interface + interface pure module subroutine swiftest_kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) !$omp declare simd(swiftest_kick_getacch_int_one_pl) implicit none diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 05256309a..0c75e1054 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -22,34 +22,12 @@ module subroutine symba_kick_getacch_int_pl(self, param) ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter - ! Internals - ! type(interaction_timer), save :: itimer - ! logical, save :: lfirst = .true. - - ! if (param%ladaptive_interactions) then - ! if (self%nplplm > 0) then - ! if (lfirst) then - ! write(itimer%loopname, *) "symba_kick_getacch_int_pl" - ! write(itimer%looptype, *) "INTERACTION" - ! call itimer%time_this_loop(param, self%nplplm, self) - ! lfirst = .false. - ! else - ! if (itimer%netcdf_io_check(param, self%nplplm)) call itimer%time_this_loop(param, self%nplplm, self) - ! end if - ! else - ! param%lflatten_interactions = .false. - ! end if - ! end if - - ! if (param%lflatten_interactions) then - ! call swiftest_kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%rh, self%Gmass, self%radius, self%ah) - ! else - call swiftest_kick_getacch_int_all_triangular_pl(self%nbody, self%nplm, self%rh, self%Gmass, self%radius, self%ah) - ! end if - - ! if (param%ladaptive_interactions .and. self%nplplm > 0) then - ! if (itimer%is_on) call itimer%adapt(param, self%nplplm, self) - ! end if + + if (param%lflatten_interactions) then + call swiftest_kick_getacch_int_all(self%nbody, self%nplplm, self%k_plpl, self%rh, self%Gmass, self%radius, self%ah) + else + call swiftest_kick_getacch_int_all(self%nbody, self%nplm, self%rh, self%Gmass, self%radius, self%ah) + end if return end subroutine symba_kick_getacch_int_pl @@ -87,7 +65,7 @@ module subroutine symba_kick_getacch_pl(self, nbody_system, param, t, lbeg) allocate(k_plpl_enc(2,nplplenc)) k_plpl_enc(1,1:nplplenc) = plpl_encounter%index1(1:nplplenc) k_plpl_enc(2,1:nplplenc) = plpl_encounter%index2(1:nplplenc) - call swiftest_kick_getacch_int_all_flat_pl(npl, nplplenc, k_plpl_enc, pl%rh, pl%Gmass, pl%radius, ah_enc) + call swiftest_kick_getacch_int_all(npl, nplplenc, k_plpl_enc, pl%rh, pl%Gmass, pl%radius, ah_enc) pl%ah(:,1:npl) = pl%ah(:,1:npl) - ah_enc(:,1:npl) end if