diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 45da5cb31..05603143b 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -13,7 +13,7 @@ module subroutine kick_getacch_int_pl(self) ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - call kick_getacch_int_all_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) return end subroutine kick_getacch_int_pl @@ -41,10 +41,11 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) end subroutine kick_getacch_int_tp - module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + module subroutine kick_getacch_int_all_flat_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 + !! 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 @@ -88,7 +89,55 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, end do return - end subroutine kick_getacch_int_all_pl + end subroutine kick_getacch_int_all_flat_pl + + + module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, 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 !! Number of massive bodies + 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 + 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 + + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, x, Gmass, radius) & + !$omp lastprivate(rji2, rlim2, xr, yr, zr) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do i = 1, npl + 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 + 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 + end do + !$omp end parallel do + + do concurrent(i = 1:npl) + acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) + end do + + return + end subroutine kick_getacch_int_all_triangular_pl module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 598cb4811..61f0de8ef 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -853,7 +853,7 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) integer(I4B), intent(in) :: npl !! Number of active massive bodies end subroutine kick_getacch_int_tp - module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + module subroutine kick_getacch_int_all_flat_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 @@ -862,7 +862,16 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, 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 + end subroutine kick_getacch_int_all_flat_pl + + module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) + implicit none + integer(I4B), intent(in) :: npl !! Number of massive bodies + 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_triangular_pl module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc) implicit none diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index e20e51ea9..fa5913507 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -2,10 +2,11 @@ use swiftest contains - subroutine symba_encounter_check_all(nplplm, k_plpl, x, v, rhill, dt, irec, lencounter, loc_lvdotr) + subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec, lencounter, loc_lvdotr) !! author: David A. Minton !! - !! Check for encounters between massive bodies. Split off from the main subroutine for performance + !! Check for encounters between massive bodies. Split off from the main subroutine for performance. + !! This is the flat (single loop) version. implicit none integer(I8B), intent(in) :: nplplm integer(I4B), dimension(:,:), intent(in) :: k_plpl @@ -38,7 +39,80 @@ subroutine symba_encounter_check_all(nplplm, k_plpl, x, v, rhill, dt, irec, len !$omp end parallel do simd return - end subroutine symba_encounter_check_all + end subroutine symba_encounter_check_all_flat + + + subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, irec, loc_lvdotr, k_plplenc, nenc) + !! author: David A. Minton + !! + !! Check for encounters between massive bodies. Split off from the main subroutine for performance + !! This is the upper triangular (double loop) version. + implicit none + integer(I4B), intent(in) :: npl, nplm + real(DP), dimension(:,:), intent(in) :: x, v + real(DP), dimension(:), intent(in) :: rhill + real(DP), intent(in) :: dt + integer(I4B), intent(in) :: irec + logical, dimension(:), allocatable, intent(out) :: loc_lvdotr + integer(I4B), dimension(:,:), allocatable, intent(out) :: k_plplenc + integer(I4B), intent(out) :: nenc + ! Internals + integer(I4B) :: i, j, nenci, j0, j1 + real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2 + logical, dimension(npl) :: lencounteri, loc_lvdotri + integer(I4B), dimension(npl) :: ind_arr + type lenctype + logical, dimension(:), allocatable :: loc_lvdotr + integer(I4B), dimension(:), allocatable :: index2 + integer(I4B) :: nenc + end type + type(lenctype), dimension(nplm) :: lenc + + ind_arr(:) = [(i, i = 1, npl)] + !$omp parallel do simd default(private) schedule(static)& + !$omp shared(npl, nplm, x, v, rhill, dt, irec, lenc) & + !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, loc_lvdotri) + 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) + vxr = v(1, j) - v(1, i) + vyr = v(2, j) - v(2, i) + vzr = v(3, j) - v(3, i) + rhill1 = rhill(i) + rhill2 = rhill(j) + call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounteri(j), loc_lvdotri(j)) + end do + lenc(i)%nenc = count(lencounteri(i+1:npl)) + if (lenc(i)%nenc > 0) then + allocate(lenc(i)%loc_lvdotr(nenci), lenc(i)%index2(nenci)) + lenc(i)%loc_lvdotr(:) = pack(loc_lvdotri(i+1:npl), lencounteri(i+1:npl)) + lenc(i)%index2(:) = pack(ind_arr(i+1:npl), lencounteri(i+1:npl)) + end if + end do + !$omp end parallel do simd + + associate(nenc_arr => lenc(:)%nenc) + nenc = sum(nenc_arr(1:nplm)) + end associate + if (nenc > 0) then + allocate(loc_lvdotr(nenc)) + allocate(k_plplenc(2,nenc)) + j0 = 1 + do i = 1, nplm + if (lenc(i)%nenc > 0) then + j1 = j0 + lenc(i)%nenc - 1 + loc_lvdotr(j0:j1) = lenc(i)%loc_lvdotr(:) + k_plplenc(1,j0:j1) = i + k_plplenc(2,j0:j1) = lenc(i)%index2(:) + j0 = j1 + 1 + end if + end do + end if + + return + end subroutine symba_encounter_check_all_triangular module function symba_encounter_check_pl(self, system, dt, irec) result(lany_encounter) @@ -68,7 +142,7 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc allocate(lencounter(nplplm)) allocate(loc_lvdotr(nplplm)) - call symba_encounter_check_all(nplplm, pl%k_plpl, pl%xh, pl%vh, pl%rhill, dt, irec, lencounter, loc_lvdotr) + call symba_encounter_check_all_flat(nplplm, pl%k_plpl, pl%xh, pl%vh, pl%rhill, dt, irec, lencounter, loc_lvdotr) nenc = count(lencounter(:)) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 325ef5a45..97b52e96e 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -13,7 +13,7 @@ module subroutine symba_kick_getacch_int_pl(self) ! Arguments class(symba_pl), intent(inout) :: self - call kick_getacch_int_all_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) return end subroutine symba_kick_getacch_int_pl @@ -52,7 +52,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) 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) + call kick_getacch_int_all_flat_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 if