From e43ae9b1d460f234113bf908f6596b1627400a75 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 23 Aug 2021 09:39:58 -0400 Subject: [PATCH] Consolidated interaction step into a single parallelized subroutine instead of having a different one for SyMBA. --- src/io/io.f90 | 6 +- src/kick/kick.f90 | 148 ++++++++++++++++++------------- src/modules/swiftest_classes.f90 | 41 +++++---- src/symba/symba_kick.f90 | 46 ++-------- 4 files changed, 124 insertions(+), 117 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index f8a0d112a..a4b619847 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1054,7 +1054,7 @@ module function io_read_frame_body(self, iu, param) result(ierr) read(iu, err = 667, iomsg = errmsg) pl%Gmass(:) pl%mass(:) = pl%Gmass(:) / param%GU if (param%lrhill_present) read(iu, err = 667, iomsg = errmsg) pl%rhill(:) - if (param%lclose) read(iu, err = 667, iomsg = errmsg) pl%radius(:) + read(iu, err = 667, iomsg = errmsg) pl%radius(:) if (param%lrotation) then read(iu, err = 667, iomsg = errmsg) pl%Ip(1, :) read(iu, err = 667, iomsg = errmsg) pl%Ip(2, :) @@ -1081,7 +1081,7 @@ module function io_read_frame_body(self, iu, param) result(ierr) end if self%Gmass(i) = real(val, kind=DP) self%mass(i) = real(val / param%GU, kind=DP) - if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i) + read(iu, *, err = 667, iomsg = errmsg) self%radius(i) class is (swiftest_tp) read(iu, *, err = 667, iomsg = errmsg) self%id(i) end select @@ -1508,7 +1508,7 @@ module subroutine io_write_frame_body(self, iu, param) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body write(iu, err = 667, iomsg = errmsg) pl%Gmass(1:n) if (param%lrhill_present) write(iu, err = 667, iomsg = errmsg) pl%rhill(1:n) - if (param%lclose) write(iu, err = 667, iomsg = errmsg) pl%radius(1:n) + write(iu, err = 667, iomsg = errmsg) pl%radius(1:n) if (param%lrotation) then write(iu, err = 667, iomsg = errmsg) pl%Ip(1, 1:n) write(iu, err = 667, iomsg = errmsg) pl%Ip(2, 1:n) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index c351d0656..19d717b0d 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -12,37 +12,8 @@ module subroutine kick_getacch_int_pl(self) implicit none ! Arguments class(swiftest_pl), intent(inout) :: self - ! Internals - integer(I8B) :: k, nplpl - real(DP) :: rji2 - real(DP) :: dx, dy, dz - integer(I4B) :: i, j - real(DP), dimension(:,:), pointer :: ah, xh - real(DP), dimension(NDIM,self%nbody) :: ahi, ahj - integer(I4B), dimension(:,:), pointer :: k_plpl - logical, dimension(:), pointer :: lmask - real(DP), dimension(:), pointer :: Gmass - - associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass) - nplpl = self%nplpl - ahi(:,:) = 0.0_DP - ahj(:,:) = 0.0_DP - !$omp parallel do default(shared)& - !$omp private(k, i, j, dx, dy, dz, rji2) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do k = 1_I8B, nplpl - i = k_plpl(1,k) - j = k_plpl(2,k) - dx = xh(1, j) - xh(1, i) - 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)) - end do - !$omp end parallel do - ah(:,1:self%nbody) = ah(:,1:self%nbody) + ahi(:,1:self%nbody) + ahj(:,1:self%nbody) - end associate + + call kick_getacch_int_all_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) return end subroutine kick_getacch_int_pl @@ -63,7 +34,6 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) integer(I4B), intent(in) :: npl !! Number of active massive bodies ! Internals integer(I4B) :: i, j - !real(DP), dimension(:,:), allocatable :: aht if ((self%nbody == 0) .or. (npl == 0)) return @@ -73,60 +43,116 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) do j = 1, npl block real(DP) :: rji2 - real(DP) :: dx, dy, dz - dx = tp%xh(1, i) - xhp(1, j) - dy = tp%xh(2, i) - xhp(1, j) - dz = tp%xh(3, i) - xhp(1, j) - rji2 = dx**2 + dy**2 + dz**2 - call kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i)) + real(DP) :: xr, yr, zr + xr = tp%xh(1, i) - xhp(1, j) + yr = tp%xh(2, i) - xhp(1, j) + zr = tp%xh(3, i) - xhp(1, j) + rji2 = xr**2 + yr**2 + zr**2 + call kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i)) end block end do end if end do - !call move_alloc(aht, tp%ah) end associate return end subroutine kick_getacch_int_tp - module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) + + module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization + !! + !! 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 - real(DP), intent(in) :: rji2 - real(DP), intent(in) :: dx, dy, dz - real(DP), intent(in) :: Gmi - real(DP), intent(in) :: Gmj - real(DP), intent(inout) :: axi, ayi, azi - real(DP), intent(inout) :: axj, ayj, azj + 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) :: 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 + ! Internals + integer(I8B) :: k + real(DP) :: rji2, rlim2 + real(DP) :: xr, yr, zr + integer(I4B) :: i, j + real(DP), dimension(NDIM,npl) :: ahi, ahj + + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + !$omp parallel do default(private)& + !$omp shared(nplpl, k_plpl, x, Gmass, radius) & + !$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 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 + acc(:,1:npl) = acc(:,1:npl) + ahi(:,1:npl) + ahj(:,1:npl) + return + end subroutine kick_getacch_int_all_pl + + + module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations for a single pair of massive bodies + !! + !! 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 + real(DP), intent(in) :: rji2 !! Square of distance between the two bodies + real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions + real(DP), intent(in) :: Gmi !! G*mass of body i + real(DP), intent(in) :: Gmj !! G*mass of body j + real(DP), intent(inout) :: axi, ayi, azi !! Acceleration vector components of body i + real(DP), intent(inout) :: axj, ayj, azj !! Acceleration vector components of body j ! Internals real(DP) :: faci, facj, irij3 irij3 = 1.0_DP / (rji2 * sqrt(rji2)) faci = Gmi * irij3 facj = Gmj * irij3 - axi = axi + facj * dx - ayi = ayi + facj * dy - azi = azi + facj * dz - axj = axj - faci * dx - ayj = ayj - faci * dy - azj = azj - faci * dz + axi = axi + facj * xr + ayi = ayi + facj * yr + azi = azi + facj * zr + axj = axj - faci * xr + ayj = ayj - faci * yr + azj = azj - faci * zr + return end subroutine kick_getacch_int_one_pl - !module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az) - module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az) + module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, az) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations of a single test particle massive body pair. + !! + !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int_tp.f90 implicit none - real(DP), intent(in) :: rji2 - real(DP), intent(in) :: dx, dy, dz - real(DP), intent(in) :: GMpl - real(DP), intent(inout) :: ax, ay, az + real(DP), intent(in) :: rji2 !! Square of distance between the test particle and massive body + real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions + real(DP), intent(in) :: Gmpl !! G*mass of massive body + real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle ! Internals real(DP) :: fac fac = GMpl / (rji2 * sqrt(rji2)) - ax = ax - fac * dx - ay = ay - fac * dy - az = az - fac * dz + ax = ax - fac * xr + ay = ay - fac * yr + az = az - fac * zr return end subroutine kick_getacch_int_one_tp diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 383091b8c..19c7967bd 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -709,29 +709,40 @@ module subroutine io_write_frame_system(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_frame_system - module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) - implicit none - real(DP), intent(in) :: rji2 - real(DP), intent(in) :: dx, dy, dz - real(DP), intent(in) :: Gmi - real(DP), intent(in) :: Gmj - real(DP), intent(inout) :: axi, ayi, azi - real(DP), intent(inout) :: axj, ayj, azj + module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, 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) :: x !! 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 kick_getacch_int_all_pl + + module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) + !$omp declare simd(kick_getacch_int_one_pl) + implicit none + real(DP), intent(in) :: rji2 !! Square of distance between the two bodies + real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions + real(DP), intent(in) :: Gmi !! G*mass of body i + real(DP), intent(in) :: Gmj !! G*mass of body j + real(DP), intent(inout) :: axi, ayi, azi !! Acceleration vector components of body i + real(DP), intent(inout) :: axj, ayj, azj !! Acceleration vector components of body j end subroutine kick_getacch_int_one_pl - !module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, Gmpl, ax, ay, az) - module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, Gmpl, ax, ay, az) + module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, ax, ay, az) !$omp declare simd(kick_getacch_int_one_tp) implicit none - real(DP), intent(in) :: rji2 - real(DP), intent(in) :: dx, dy, dz - real(DP), intent(in) :: Gmpl - real(DP), intent(inout) :: ax, ay, az + real(DP), intent(in) :: rji2 !! Square of distance between the test particle and massive body + real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions + real(DP), intent(in) :: Gmpl !! G*mass of massive body + real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle end subroutine kick_getacch_int_one_tp module subroutine kick_getacch_int_pl(self) implicit none - class(swiftest_pl), intent(inout) :: self + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object end subroutine kick_getacch_int_pl module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index d282f9c53..af440eada 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -12,38 +12,8 @@ module subroutine symba_kick_getacch_int_pl(self) implicit none ! Arguments class(symba_pl), intent(inout) :: self - ! Internals - integer(I8B) :: k, nplplm - real(DP) :: rji2, rlim2 - real(DP) :: dx, dy, dz - integer(I4B) :: i, j - real(DP), dimension(:,:), pointer :: ah, xh - real(DP), dimension(NDIM,self%nbody) :: ahi, ahj - integer(I4B), dimension(:,:), pointer :: k_plpl - logical, dimension(:), pointer :: lmask - real(DP), dimension(:), pointer :: Gmass, radius - 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, rlim2) & - !$omp reduction(+:ahi) & - !$omp reduction(-:ahj) - do k = 1_I8B, nplplm - i = k_plpl(1,k) - j = k_plpl(2,k) - dx = xh(1, j) - xh(1, i) - dy = xh(2, j) - xh(2, i) - dz = xh(3, j) - xh(3, i) - rji2 = dx**2 + dy**2 + dz**2 - 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(:,1:self%nbody) = ah(:,1:self%nbody) + ahi(:,1:self%nbody) + ahj(:,1:self%nbody) - end associate + call kick_getacch_int_all_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) return end subroutine symba_kick_getacch_int_pl @@ -66,7 +36,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) ! Internals integer(I4B) :: i, j integer(I8B) :: k, nplplenc - real(DP) :: rji2, rlim2, dx, dy, dz + 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 @@ -82,18 +52,18 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) ahi(:,:) = 0.0_DP ahj(:,:) = 0.0_DP !$omp parallel do default(shared)& - !$omp private(k, i, j, dx, dy, dz, rji2, rlim2) & + !$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) - dx = xh(1, j) - xh(1, i) - dy = xh(2, j) - xh(2, i) - dz = xh(3, j) - xh(3, i) - rji2 = dx**2 + dy**2 + dz**2 + 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, 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)) + 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)