From 0a4fae0b0f9216b01901730672e0d305bb339aa6 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 26 Feb 2023 16:20:26 -0500 Subject: [PATCH] Improved the parallel speedup of the pl-pl kick subroutine --- src/swiftest/swiftest_kick.f90 | 61 ++++++++++++++------------------ src/swiftest/swiftest_module.f90 | 4 +-- 2 files changed, 29 insertions(+), 36 deletions(-) diff --git a/src/swiftest/swiftest_kick.f90 b/src/swiftest/swiftest_kick.f90 index 23740432a..76819ac59 100644 --- a/src/swiftest/swiftest_kick.f90 +++ b/src/swiftest/swiftest_kick.f90 @@ -153,7 +153,7 @@ module subroutine swiftest_kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, G end subroutine swiftest_kick_getacch_int_all_flat_pl - module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc) + module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, r, Gmass, radius, acc) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization. @@ -164,60 +164,53 @@ module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmas 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) :: 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 ! Internals - real(DP), dimension(NDIM,npl) :: ahi, ahj integer(I4B) :: i, j - real(DP) :: rji2, rlim2 - real(DP) :: xr, yr, zr - - ahi(:,:) = 0.0_DP - ahj(:,:) = 0.0_DP + real(DP) :: rji2, rlim2, fac, rx, ry, rz if (present(radius)) then !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, nplm, x, Gmass, radius) & - !$omp lastprivate(rji2, rlim2, xr, yr, zr) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) + !$omp shared(npl, nplm, r, Gmass, radius) & + !$omp reduction(-:acc) do i = 1, nplm - do concurrent(j = i+1:npl) - 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 + 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) 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)) + if (rji2 > rlim2) then + fac = Gmass(i) / (rji2 * sqrt(rji2)) + acc(1,j) = acc(1,j) - fac * rx + acc(2,j) = acc(2,j) - fac * ry + acc(3,j) = acc(3,j) - fac * rz + end if end do end do !$omp end parallel do else !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, nplm, x, Gmass) & - !$omp lastprivate(rji2, xr, yr, zr) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) + !$omp shared(npl, nplm, r, Gmass) & + !$omp reduction(-:acc) do i = 1, nplm - do concurrent(j = i+1:npl) - 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)) + 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(i) / (rji2 * sqrt(rji2)) + acc(1,j) = acc(1,j) - fac * rx + acc(2,j) = acc(2,j) - fac * ry + acc(3,j) = acc(3,j) - fac * rz end do end do !$omp end parallel do end if - do concurrent(i = 1:npl) - acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) - end do - return end subroutine swiftest_kick_getacch_int_all_triangular_pl diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 7633e6e44..78b307977 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -902,11 +902,11 @@ module subroutine swiftest_kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, G real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array end subroutine swiftest_kick_getacch_int_all_flat_pl - module subroutine swiftest_kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc) + module subroutine swiftest_kick_getacch_int_all_triangular_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) :: 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