From f4d1d7bfe9416d96933073fcfbef8a01529043d9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 4 Mar 2023 12:19:41 -0500 Subject: [PATCH 1/9] Update axis labels --- .../parallel-performance-plots.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/swiftest_performance/symba-performance-comparison_JOSS/parallel-performance-plots.py b/examples/swiftest_performance/symba-performance-comparison_JOSS/parallel-performance-plots.py index dd18e1cdc..a2ec2a85a 100755 --- a/examples/swiftest_performance/symba-performance-comparison_JOSS/parallel-performance-plots.py +++ b/examples/swiftest_performance/symba-performance-comparison_JOSS/parallel-performance-plots.py @@ -36,8 +36,8 @@ ax.set_yticks([1, 2, 4, 6, 8, 12, 16, 20, 24]) ax.set_xticks([1, 2, 4, 6, 8, 12, 16, 20, 24]) ax.grid(True) - ax.set_xlabel("Number of CPUs", fontsize=axes_fontsize) - ax.set_ylabel("Speedup (relative to Swift-SyMBA)", fontsize=axes_fontsize) + ax.set_xlabel("Number of cores", fontsize=axes_fontsize) + ax.set_ylabel("Speedup (relative to Swift)", fontsize=axes_fontsize) ax.set_facecolor("white") plt.plot(swest[n]['N cores'], swift[n]['wall time(s)'][0] / swest[n]['wall time(s)'], alpha=0.5, linewidth=6, label="Swiftest") plt.plot(swomp[n]['N cores'], swift[n]['wall time(s)'][0] / swomp[n]['wall time(s)'], c="darkgreen", alpha=0.5, linewidth=6, label="Swifter-OMP") From 4e4712ef049e33c56e377e4dd25ee2d2e62681ee Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 7 Mar 2023 11:41:23 -0500 Subject: [PATCH 2/9] 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 From 5688e850cb13d0a04b082f0c9bad69afb867ebcb Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 7 Mar 2023 12:36:43 -0500 Subject: [PATCH 3/9] removed redundant mask --- src/swiftest/swiftest_kick.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/swiftest/swiftest_kick.f90 b/src/swiftest/swiftest_kick.f90 index b3a5f7d39..e7d454c12 100644 --- a/src/swiftest/swiftest_kick.f90 +++ b/src/swiftest/swiftest_kick.f90 @@ -232,7 +232,7 @@ module subroutine swiftest_kick_getacch_int_all_tri_rad_pl(npl, nplm, r, Gmass, !$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) + do concurrent(j = 1:nplm) rx = r(1,j) - r(1,i) ry = r(2,j) - r(2,i) rz = r(3,j) - r(3,i) @@ -319,7 +319,7 @@ module subroutine swiftest_kick_getacch_int_all_tri_norad_pl(npl, nplm, r, Gmass !$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) + do concurrent(j = 1:nplm) rx = r(1,j) - r(1,i) ry = r(2,j) - r(2,i) rz = r(3,j) - r(3,i) From 9351997477042163a936bf58bed6de8463cbcd0f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 7 Mar 2023 17:30:23 -0500 Subject: [PATCH 4/9] Some minor refinements to kick and stripped away adaptive code for encounter checks --- src/encounter/encounter_check.f90 | 102 ++---------------------------- src/swiftest/swiftest_kick.f90 | 72 +++++++++++---------- 2 files changed, 43 insertions(+), 131 deletions(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 46b712910..917fcae45 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -34,38 +34,8 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, ind logical, save :: skipit = .false. ! This will be used to ensure that the sort & sweep subroutine gets called at least once before timing it so that the extent array is nearly sorted when it is timed integer(I8B) :: nplpl = 0_I8B - ! if (param%ladaptive_encounters_plpl .and. (.not. skipit)) then - ! nplpl = (npl * (npl - 1) / 2) - ! if (nplpl > 0) then - ! if (lfirst) then - ! write(itimer%loopname, *) "encounter_check_all_plpl" - ! write(itimer%looptype, *) "ENCOUNTER_PLPL" - ! lfirst = .false. - ! itimer%step_counter = INTERACTION_TIMER_CADENCE - ! else - ! if (itimer%netcdf_io_check(param, nplpl)) call itimer%time_this_loop(param, nplpl) - ! end if - ! else - ! param%lencounter_sas_plpl = .false. - ! end if - ! end if - - ! if (param%lencounter_sas_plpl) then - ! call encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) - ! else - call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) - ! end if - - ! if (skipit) then - ! skipit = .false. - ! else - ! if (param%ladaptive_encounters_plpl .and. nplpl > 0) then - ! if (itimer%is_on) then - ! call itimer%adapt(param, nplpl) - ! skipit = .true. - ! end if - ! end if - ! end if + + call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_plpl @@ -107,23 +77,6 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, integer(I4B), dimension(:), allocatable :: itmp logical, dimension(:), allocatable :: ltmp - ! if (param%ladaptive_encounters_plpl .and. (.not. skipit)) then - ! npl = nplm + nplt - ! nplplm = nplm * npl - nplm * (nplm + 1) / 2 - ! if (nplplm > 0) then - ! if (lfirst) then - ! write(itimer%loopname, *) "encounter_check_all_plpl" - ! write(itimer%looptype, *) "ENCOUNTER_PLPL" - ! lfirst = .false. - ! itimer%step_counter = INTERACTION_TIMER_CADENCE - ! else - ! if (itimer%netcdf_io_check(param, nplplm)) call itimer%time_this_loop(param, nplplm) - ! end if - ! else - ! param%lencounter_sas_plpl = .false. - ! end if - ! end if - allocate(tmp_param, source=param) ! Turn off adaptive encounter checks for the pl-pl group @@ -132,23 +85,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, ! Start with the pl-pl group call encounter_check_all_plpl(tmp_param, nplm, rplm, vplm, rencm, dt, nenc, index1, index2, lvdotr) - ! if (param%lencounter_sas_plpl) then - ! call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, & - ! plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) - ! else - call encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) - ! end if - - ! if (skipit) then - ! skipit = .false. - ! else - ! if (param%ladaptive_encounters_plpl .and. nplplm > 0) then - ! if (itimer%is_on) then - ! call itimer%adapt(param, nplplm) - ! skipit = .true. - ! end if - ! end if - ! end if + call encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) if (plmplt_nenc > 0) then ! Consolidate the two lists allocate(itmp(nenc+plmplt_nenc)) @@ -202,37 +139,7 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, logical, save :: lsecond = .false. integer(I8B) :: npltp = 0_I8B - ! if (param%ladaptive_encounters_pltp) then - ! npltp = npl * ntp - ! if (npltp > 0) then - ! if (lfirst) then - ! write(itimer%loopname, *) "encounter_check_all_pltp" - ! write(itimer%looptype, *) "ENCOUNTER_PLTP" - ! lfirst = .false. - ! lsecond = .true. - ! else - ! if (lsecond) then ! This ensures that the encounter check methods are run at least once prior to timing. Sort and sweep improves on the second pass due to the bounding box extents needing to be nearly sorted - ! call itimer%time_this_loop(param, npltp) - ! lsecond = .false. - ! else if (itimer%netcdf_io_check(param, npltp)) then - ! lsecond = .true. - ! itimer%is_on = .false. - ! end if - ! end if - ! else - ! param%lencounter_sas_pltp = .false. - ! end if - ! end if - - ! if (param%lencounter_sas_pltp) then - ! call encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) - ! else - call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) - ! end if - - ! if (.not.lfirst .and. param%ladaptive_encounters_pltp .and. npltp > 0) then - ! if (itimer%is_on) call itimer%adapt(param, npltp) - ! end if + call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_pltp @@ -576,6 +483,7 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr) + nenc = 0 return end subroutine encounter_check_all_triangular_plpl diff --git a/src/swiftest/swiftest_kick.f90 b/src/swiftest/swiftest_kick.f90 index e7d454c12..16bb2ac31 100644 --- a/src/swiftest/swiftest_kick.f90 +++ b/src/swiftest/swiftest_kick.f90 @@ -211,9 +211,9 @@ module subroutine swiftest_kick_getacch_int_all_tri_rad_pl(npl, nplm, r, Gmass, end do else !$omp parallel do default(private) schedule(static)& - !$omp shared(npl,nplm, r, Gmass, radius, acc) + !$omp shared(npl, nplm, r, Gmass, radius, acc) do i = 1, nplm - do concurrent(j = 1:npl, j/=i) + do concurrent(j = 1:npl, i/=j) rx = r(1,j) - r(1,i) ry = r(2,j) - r(2,i) rz = r(3,j) - r(3,i) @@ -229,24 +229,26 @@ module subroutine swiftest_kick_getacch_int_all_tri_rad_pl(npl, nplm, r, Gmass, 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) - 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 + if (nplt > 0) then + !$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) + 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 - end do - !$omp end parallel do + !$omp end parallel do + end if end if @@ -301,7 +303,7 @@ module subroutine swiftest_kick_getacch_int_all_tri_norad_pl(npl, nplm, r, Gmass end do else !$omp parallel do default(private) schedule(static)& - !$omp shared(npl,nplm, r, Gmass, acc) + !$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) @@ -316,21 +318,23 @@ module subroutine swiftest_kick_getacch_int_all_tri_norad_pl(npl, nplm, r, Gmass 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) - 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 + if (nplt > 0) then + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, nplm, r, Gmass, acc) + do i = nplm+1,npl + do concurrent(j = 1:nplm) + 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 - end do - !$omp end parallel do + !$omp end parallel do + end if end if From 5e6eeb0b2bb256ccceefc1213c325ec0b2515f18 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 8 Mar 2023 08:54:45 -0500 Subject: [PATCH 5/9] Simplified sort and sweep to use radial distance and set up for testing --- src/encounter/encounter_check.f90 | 270 ++++++++++++---------------- src/encounter/encounter_module.f90 | 16 +- src/encounter/encounter_util.f90 | 32 ++-- src/rmvs/rmvs_encounter_check.f90 | 2 +- src/swiftest/swiftest_module.f90 | 4 +- src/swiftest/swiftest_orbel.f90 | 4 +- src/symba/symba_encounter_check.f90 | 2 +- src/symba/symba_step.f90 | 2 +- 8 files changed, 141 insertions(+), 191 deletions(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 917fcae45..aa8da45a5 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -11,18 +11,18 @@ use swiftest contains - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plpl(param, npl, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. Choose between the standard triangular or the Sort & Sweep method based on user inputs !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies - real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), intent(in) :: dt !! Step size integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter @@ -35,7 +35,8 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, ind integer(I8B) :: nplpl = 0_I8B - call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + !call encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) + call encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_plpl @@ -48,7 +49,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: nplm !! Total number of fully interacting massive bodies integer(I4B), intent(in) :: nplt !! Total number of partially interacting masive bodies (GM < GMTINY) real(DP), dimension(:,:), intent(in) :: rplm !! Position vectors of fully interacting massive bodies @@ -113,19 +114,19 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, end subroutine encounter_check_all_plplm - module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. Choose between the standard triangular or the Sort & Sweep method based on user inputs !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies - real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of test particlse + real(DP), dimension(:,:), intent(in) :: rtp !! Position vectors of test particlse real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of test particles real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -139,13 +140,13 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, logical, save :: lsecond = .false. integer(I8B) :: npltp = 0_I8B - call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_pltp - subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. @@ -155,7 +156,7 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in implicit none ! Arguments integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -164,10 +165,11 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals - integer(I4B) :: dim, n + integer(I4B) :: i, n integer(I4B), save :: npl_last = 0 type(encounter_bounding_box), save :: boundingbox - integer(I2B), dimension(npl) :: vshift_min, vshift_max + real(DP), dimension(npl) :: rmin,rmax + real(DP) :: rmag if (npl == 0) return @@ -178,23 +180,15 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in npl_last = npl end if - !$omp parallel do default(private) schedule(static) & - !$omp shared(x, v, renc, boundingbox) & - !$omp firstprivate(dt, npl, n) - do dim = 1, SWEEPDIM - where(v(dim,1:npl) < 0.0_DP) - vshift_min(1:npl) = 1 - vshift_max(1:npl) = 0 - elsewhere - vshift_min(1:npl) = 0 - vshift_max(1:npl) = 1 - end where - call boundingbox%aabb(dim)%sort(npl, [x(dim,1:npl) - renc(1:npl) + vshift_min(1:npl) * v(dim,1:npl) * dt, & - x(dim,1:npl) + renc(1:npl) + vshift_max(1:npl) * v(dim,1:npl) * dt]) + do concurrent (i = 1:npl) + rmag = .mag.r(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * renc(i) + rmin(i) = rmag - RSWEEP_FACTOR * renc(i) end do - !$omp end parallel do - call boundingbox%sweep(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + call boundingbox%aabb%sort(npl, [rmin,rmax]) + + call boundingbox%sweep(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_sort_and_sweep_plpl @@ -224,10 +218,10 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(encounter_bounding_box), save :: boundingbox - integer(I4B) :: dim, n, ntot + integer(I4B) :: i, n, ntot integer(I4B), save :: ntot_last = 0 - integer(I2B), dimension(nplm) :: vplmshift_min, vplmshift_max - integer(I2B), dimension(nplt) :: vpltshift_min, vpltshift_max + real(DP), dimension(nplm+nplt) :: rmin,rmax + real(DP) :: rmag ! If this is the first time through, build the index lists if ((nplm == 0) .or. (nplt == 0)) return @@ -239,32 +233,18 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt ntot_last = ntot end if - !$omp parallel do default(private) schedule(static) & - !$omp shared(rplm, rplt, vplm, vplt, rencm, renct, boundingbox) & - !$omp firstprivate(dt, nplm, nplt, ntot) - do dim = 1, SWEEPDIM - where(vplm(dim,1:nplm) < 0.0_DP) - vplmshift_min(1:nplm) = 1 - vplmshift_max(1:nplm) = 0 - elsewhere - vplmshift_min(1:nplm) = 0 - vplmshift_max(1:nplm) = 1 - end where - - where(vplt(dim,1:nplt) < 0.0_DP) - vpltshift_min(1:nplt) = 1 - vpltshift_max(1:nplt) = 0 - elsewhere - vpltshift_min(1:nplt) = 0 - vpltshift_max(1:nplt) = 1 - end where - - call boundingbox%aabb(dim)%sort(ntot, [rplm(dim,1:nplm) - rencm(1:nplm) + vplmshift_min(1:nplm) * vplm(dim,1:nplm) * dt, & - rplt(dim,1:nplt) - renct(1:nplt) + vpltshift_min(1:nplt) * vplt(dim,1:nplt) * dt, & - rplm(dim,1:nplm) + rencm(1:nplm) + vplmshift_max(1:nplm) * vplm(dim,1:nplm) * dt, & - rplt(dim,1:nplt) + renct(1:nplt) + vpltshift_max(1:nplt) * vplt(dim,1:nplt) * dt]) + do concurrent (i = 1:nplm) + rmag = .mag.rplm(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * rencm(i) + rmin(i) = rmag - RSWEEP_FACTOR * rencm(i) end do - !$omp end parallel do + do concurrent (i = 1:nplt) + rmag = .mag.rplt(:,i) + rmax(nplm+i) = rmag + RSWEEP_FACTOR * renct(i) + rmin(nplm+i) = rmag - RSWEEP_FACTOR * renct(i) + end do + + call boundingbox%aabb%sort(ntot, [rmin, rmax]) call boundingbox%sweep(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) @@ -272,7 +252,7 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt end subroutine encounter_check_all_sort_and_sweep_plplm - subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, rencpl, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, rtp, vtp, rencpl, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. @@ -285,7 +265,7 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies - real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: rtp !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: rencpl !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -295,11 +275,11 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(encounter_bounding_box), save :: boundingbox - integer(I4B) :: dim, n, ntot + integer(I4B) :: i, n, ntot integer(I4B), save :: ntot_last = 0 - integer(I2B), dimension(npl) :: vplshift_min, vplshift_max - integer(I2B), dimension(ntp) :: vtpshift_min, vtpshift_max - real(DP), dimension(ntp) :: renctp + real(DP), dimension(npl+ntp) :: rmin,rmax + real(DP), dimension(ntp) :: renctp + real(DP) :: rmag ! If this is the first time through, build the index lists if ((ntp == 0) .or. (npl == 0)) return @@ -312,35 +292,21 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, end if renctp(:) = 0.0_DP - - !$omp parallel do default(private) schedule(static) & - !$omp shared(rpl, xtp, vpl, vtp, rencpl, renctp, boundingbox) & - !$omp firstprivate(dt, npl, ntp, ntot) - do dim = 1, SWEEPDIM - where(vpl(dim,1:npl) < 0.0_DP) - vplshift_min(1:npl) = 1 - vplshift_max(1:npl) = 0 - elsewhere - vplshift_min(1:npl) = 0 - vplshift_max(1:npl) = 1 - end where - - where(vtp(dim,1:ntp) < 0.0_DP) - vtpshift_min(1:ntp) = 1 - vtpshift_max(1:ntp) = 0 - elsewhere - vtpshift_min(1:ntp) = 0 - vtpshift_max(1:ntp) = 1 - end where - - call boundingbox%aabb(dim)%sort(ntot, [rpl(dim,1:npl) - rencpl(1:npl) + vplshift_min(1:npl) * vpl(dim,1:npl) * dt, & - xtp(dim,1:ntp) - renctp(1:ntp) + vtpshift_min(1:ntp) * vtp(dim,1:ntp) * dt, & - rpl(dim,1:npl) + rencpl(1:npl) + vplshift_max(1:npl) * vpl(dim,1:npl) * dt, & - xtp(dim,1:ntp) + renctp(1:ntp) + vtpshift_max(1:ntp) * vtp(dim,1:ntp) * dt]) + + do concurrent (i = 1:npl) + rmag = .mag.rpl(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * rencpl(i) + rmin(i) = rmag - RSWEEP_FACTOR * rencpl(i) + end do + do concurrent (i = 1:ntp) + rmag = .mag.rtp(:,i) + rmax(npl+i) = rmag + RSWEEP_FACTOR * renctp(i) + rmin(npl+i) = rmag - RSWEEP_FACTOR * renctp(i) end do - !$omp end parallel do - call boundingbox%sweep(npl, ntp, rpl, vpl, xtp, vtp, rencpl, renctp, dt, nenc, index1, index2, lvdotr) + call boundingbox%aabb%sort(ntot, [rmin, rmax]) + + call boundingbox%sweep(npl, ntp, rpl, vpl, rtp, vtp, rencpl, renctp, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_sort_and_sweep_pltp @@ -445,7 +411,7 @@ subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x end subroutine encounter_check_all_triangular_one - subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. Split off from the main subroutine for performance @@ -453,9 +419,9 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 implicit none ! Arguments integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies - real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter @@ -469,12 +435,12 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 call swiftest_util_index_array(ind_arr, npl) !$omp parallel do default(private) schedule(static)& - !$omp shared(x, v, renc, lenc, ind_arr) & + !$omp shared(r, v, renc, lenc, ind_arr) & !$omp firstprivate(npl, dt) do i = 1,npl - call encounter_check_all_triangular_one(i, npl, x(1,i), x(2,i), x(3,i), & + call encounter_check_all_triangular_one(i, npl, r(1,i), r(2,i), r(3,i), & v(1,i), v(2,i), v(3,i), & - x(1,:), x(2,:), x(3,:), & + r(1,:), r(2,:), r(3,:), & v(1,:), v(2,:), v(3,:), & renc(i), renc(:), dt, ind_arr(:), lenc(i)) if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i @@ -483,7 +449,6 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr) - nenc = 0 return end subroutine encounter_check_all_triangular_plpl @@ -820,9 +785,8 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r ! Internals integer(I4B) :: ii, i, ntot, nbox, dim logical, dimension(n1+n2) :: loverlap - logical, dimension(SWEEPDIM,n1+n2) :: loverlap_by_dimension - logical, dimension(SWEEPDIM,2*(n1+n2)) :: llist1 - integer(I4B), dimension(SWEEPDIM,2*(n1+n2)) :: ext_ind + logical, dimension(2*(n1+n2)) :: llist1 + integer(I4B), dimension(2*(n1+n2)) :: ext_ind type(collision_list_pltp), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I8B) :: ibeg, iend @@ -831,39 +795,33 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r ntot = n1 + n2 call swiftest_util_index_array(ind_arr, ntot) - do concurrent(dim = 1:SWEEPDIM) - loverlap_by_dimension(dim,:) = (self%aabb(dim)%ibeg(:) + 1_I8B) < (self%aabb(dim)%iend(:) - 1_I8B) - where(self%aabb(dim)%ind(:) > ntot) - ext_ind(dim,:) = self%aabb(dim)%ind(:) - ntot - elsewhere - ext_ind(dim,:) = self%aabb(dim)%ind(:) - endwhere - end do - llist1(:,:) = ext_ind(:,:) <= n1 - where(.not.llist1(:,:)) ext_ind(:,:) = ext_ind(:,:) - n1 + loverlap(:) = (self%aabb%ibeg(:) + 1_I8B) < (self%aabb%iend(:) - 1_I8B) + where(self%aabb%ind(:) > ntot) + ext_ind(:) = self%aabb%ind(:) - ntot + elsewhere + ext_ind(:) = self%aabb%ind(:) + endwhere + llist1(:) = ext_ind(:) <= n1 + where(.not.llist1(:)) ext_ind(:) = ext_ind(:) - n1 - loverlap(:) = loverlap_by_dimension(1,:) - do dim = 2, SWEEPDIM - loverlap(:) = loverlap(:) .and. loverlap_by_dimension(dim,:) - end do dim = 1 - where(llist1(dim,:)) - xind(:) = r1(1,ext_ind(dim,:)) - yind(:) = r1(2,ext_ind(dim,:)) - zind(:) = r1(3,ext_ind(dim,:)) - vxind(:) = v1(1,ext_ind(dim,:)) - vyind(:) = v1(2,ext_ind(dim,:)) - vzind(:) = v1(3,ext_ind(dim,:)) - rencind(:) = renc1(ext_ind(dim,:)) + where(llist1(:)) + xind(:) = r1(1,ext_ind(:)) + yind(:) = r1(2,ext_ind(:)) + zind(:) = r1(3,ext_ind(:)) + vxind(:) = v1(1,ext_ind(:)) + vyind(:) = v1(2,ext_ind(:)) + vzind(:) = v1(3,ext_ind(:)) + rencind(:) = renc1(ext_ind(:)) elsewhere - xind(:) = r2(1,ext_ind(dim,:)) - yind(:) = r2(2,ext_ind(dim,:)) - zind(:) = r2(3,ext_ind(dim,:)) - vxind(:) = v2(1,ext_ind(dim,:)) - vyind(:) = v2(2,ext_ind(dim,:)) - vzind(:) = v2(3,ext_ind(dim,:)) - rencind(:) = renc2(ext_ind(dim,:)) + xind(:) = r2(1,ext_ind(:)) + yind(:) = r2(2,ext_ind(:)) + zind(:) = r2(3,ext_ind(:)) + vxind(:) = v2(1,ext_ind(:)) + vyind(:) = v2(2,ext_ind(:)) + vzind(:) = v2(3,ext_ind(:)) + rencind(:) = renc2(ext_ind(:)) endwhere where(.not.loverlap(:)) lenc(:)%nenc = 0 @@ -875,14 +833,14 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r !$omp do schedule(static) do i = 1, n1 if (loverlap(i)) then - ibeg = self%aabb(dim)%ibeg(i) + 1_I8B - iend = self%aabb(dim)%iend(i) - 1_I8B + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = iend - ibeg + 1 call encounter_check_all_sweep_one(i, nbox, r1(1,i), r1(2,i), r1(3,i), v1(1,i), v1(2,i), v1(3,i), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & - renc1(i), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & - .not.llist1(dim,ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) + renc1(i), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & + .not.llist1(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) end if end do !$omp end do nowait @@ -891,15 +849,15 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r !$omp do schedule(static) do i = n1+1, ntot if (loverlap(i)) then - ibeg = self%aabb(dim)%ibeg(i) + 1_I8B - iend = self%aabb(dim)%iend(i) - 1_I8B + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = iend - ibeg + 1 ii = i - n1 call encounter_check_all_sweep_one(ii, nbox, r1(1,ii), r1(2,ii), r1(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & - renc2(ii), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & - llist1(dim,ibeg:iend), lenc(i)%nenc, lenc(i)%index2, lenc(i)%index1, lenc(i)%lvdotr) + renc2(ii), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & + llist1(ibeg:iend), lenc(i)%nenc, lenc(i)%index2, lenc(i)%index1, lenc(i)%lvdotr) end if end do !$omp end do nowait @@ -914,7 +872,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r end subroutine encounter_check_sweep_aabb_double_list - module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_sweep_aabb_single_list(self, n, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Sweeps the sorted bounding box extents and returns the true encounters (combines broad and narrow phases) @@ -923,7 +881,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt ! Arguments class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n !! Number of bodies - real(DP), dimension(:,:), intent(in) :: x, v !! Array of position and velocity vectors + real(DP), dimension(:,:), intent(in) :: r, v !! Array of position and velocity vectors real(DP), dimension(:), intent(in) :: renc !! Radius of encounter regions of bodies 1 real(DP), intent(in) :: dt !! Step size integer(I8B), intent(out) :: nenc !! Total number of encounter candidates @@ -936,7 +894,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt logical, dimension(n) :: loverlap logical, dimension(2*n) :: lencounteri real(DP), dimension(2*n) :: xind, yind, zind, vxind, vyind, vzind, rencind - integer(I4B), dimension(SWEEPDIM,2*n) :: ext_ind + integer(I4B), dimension(2*n) :: ext_ind type(collision_list_plpl), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I8B) :: ibeg, iend @@ -946,36 +904,36 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt ! Sweep the intervals for each of the massive bodies along one dimension ! This will build a ragged pair of index lists inside of the lenc data structure - where(self%aabb(dim)%ind(:) > n) - ext_ind(dim,:) = self%aabb(1)%ind(:) - n + where(self%aabb%ind(:) > n) + ext_ind(:) = self%aabb%ind(:) - n elsewhere - ext_ind(dim,:) = self%aabb(1)%ind(:) + ext_ind(:) = self%aabb%ind(:) endwhere - xind(:) = x(1,ext_ind(dim,:)) - yind(:) = x(2,ext_ind(dim,:)) - zind(:) = x(3,ext_ind(dim,:)) - vxind(:) = v(1,ext_ind(dim,:)) - vyind(:) = v(2,ext_ind(dim,:)) - vzind(:) = v(3,ext_ind(dim,:)) - rencind(:) = renc(ext_ind(dim,:)) + xind(:) = r(1,ext_ind(:)) + yind(:) = r(2,ext_ind(:)) + zind(:) = r(3,ext_ind(:)) + vxind(:) = v(1,ext_ind(:)) + vyind(:) = v(2,ext_ind(:)) + vzind(:) = v(3,ext_ind(:)) + rencind(:) = renc(ext_ind(:)) - loverlap(:) = (self%aabb(dim)%ibeg(:) + 1_I8B) < (self%aabb(dim)%iend(:) - 1_I8B) + loverlap(:) = (self%aabb%ibeg(:) + 1_I8B) < (self%aabb%iend(:) - 1_I8B) where(.not.loverlap(:)) lenc(:)%nenc = 0 !$omp parallel do default(private) schedule(static)& - !$omp shared(self, ext_ind, lenc, loverlap, x, v, renc, xind, yind, zind, vxind, vyind, vzind, rencind) & + !$omp shared(self, ext_ind, lenc, loverlap, r, v, renc, xind, yind, zind, vxind, vyind, vzind, rencind) & !$omp firstprivate(n, dt, dim) do i = 1, n if (loverlap(i)) then - ibeg = self%aabb(dim)%ibeg(i) + 1_I8B - iend = self%aabb(dim)%iend(i) - 1_I8B + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = int(iend - ibeg + 1, kind=I4B) lencounteri(ibeg:iend) = .true. - call encounter_check_all_sweep_one(i, nbox, x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), & + call encounter_check_all_sweep_one(i, nbox, r(1,i), r(2,i), r(3,i), v(1,i), v(2,i), v(3,i), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & - renc(i), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & + renc(i), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & lencounteri(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) end if end do diff --git a/src/encounter/encounter_module.f90 b/src/encounter/encounter_module.f90 index 720627485..f3c24c763 100644 --- a/src/encounter/encounter_module.f90 +++ b/src/encounter/encounter_module.f90 @@ -17,8 +17,8 @@ module encounter implicit none public - integer(I4B), parameter :: SWEEPDIM = 3 character(len=*), parameter :: ENCOUNTER_OUTFILE = 'encounters.nc' !! Name of NetCDF output file for encounter information + real(DP), parameter :: RSWEEP_FACTOR = 1.1_DP type, abstract :: encounter_list integer(I8B) :: nenc = 0 !! Total number of encounters @@ -95,7 +95,7 @@ module encounter type encounter_bounding_box - type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb + type(encounter_bounding_box_1D) :: aabb contains procedure :: dealloc => encounter_util_dealloc_bounding_box !! Deallocates all allocatables procedure :: setup => encounter_util_setup_aabb !! Setup a new axis-aligned bounding box structure @@ -107,12 +107,12 @@ module encounter interface - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plpl(param, npl, r, v, renc, dt, nenc, index1, index2, lvdotr) use base, only: base_parameters implicit none class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -142,7 +142,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x end subroutine encounter_check_all_plplm - module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) use base, only: base_parameters implicit none class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s @@ -150,7 +150,7 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies - real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: rtp !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -205,11 +205,11 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_sweep_aabb_double_list - module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_sweep_aabb_single_list(self, n, r, v, renc, dt, nenc, index1, index2, lvdotr) implicit none class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n !! Number of bodies - real(DP), dimension(:,:), intent(in) :: x, v !! Array of position and velocity vectors + real(DP), dimension(:,:), intent(in) :: r, v !! Array of position and velocity vectors real(DP), dimension(:), intent(in) :: renc !! Radius of encounter regions of bodies 1 real(DP), intent(in) :: dt !! Step size integer(I8B), intent(out) :: nenc !! Total number of encounter candidates diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 941ed6e16..01848d571 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -104,9 +104,7 @@ module subroutine encounter_util_dealloc_bounding_box(self) ! Internals integer(I4B) :: i - do i = 1,NDIM - call self%aabb(i)%dealloc() - end do + call self%aabb%dealloc() return end subroutine encounter_util_dealloc_bounding_box @@ -351,26 +349,20 @@ module subroutine encounter_util_setup_aabb(self, n, n_last) next_last = 2 * n_last if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies - do dim = 1, SWEEPDIM - allocate(itmp(next)) - if (n_last > 0) itmp(1:next_last) = self%aabb(dim)%ind(1:next_last) - call move_alloc(itmp, self%aabb(dim)%ind) - self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)] - end do + allocate(itmp(next)) + if (n_last > 0) itmp(1:next_last) = self%aabb%ind(1:next_last) + call move_alloc(itmp, self%aabb%ind) + self%aabb%ind(next_last+1:next) = [(k, k = next_last+1, next)] else ! The number of bodies has gone down. Resize and chop of the old indices - do dim = 1, SWEEPDIM - allocate(itmp(next)) - itmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) - call move_alloc(itmp, self%aabb(dim)%ind) - end do + allocate(itmp(next)) + itmp(1:next) = pack(self%aabb%ind(1:next_last), self%aabb%ind(1:next_last) <= next) + call move_alloc(itmp, self%aabb%ind) end if - do dim = 1, SWEEPDIM - if (allocated(self%aabb(dim)%ibeg)) deallocate(self%aabb(dim)%ibeg) - allocate(self%aabb(dim)%ibeg(n)) - if (allocated(self%aabb(dim)%iend)) deallocate(self%aabb(dim)%iend) - allocate(self%aabb(dim)%iend(n)) - end do + if (allocated(self%aabb%ibeg)) deallocate(self%aabb%ibeg) + allocate(self%aabb%ibeg(n)) + if (allocated(self%aabb%iend)) deallocate(self%aabb%iend) + allocate(self%aabb%iend(n)) return end subroutine encounter_util_setup_aabb diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 00aafd1fb..0880cb136 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -40,7 +40,7 @@ module function rmvs_encounter_check_tp(self, param, nbody_system, dt) result(le select type(pl => nbody_system%pl) class is (rmvs_pl) - associate(tp => self, ntp => self%nbody, npl => pl%nbody) + associate(tp => self, ntp => self%nbody, npl => pl%nbody, cb => nbody_system%cb) tp%plencP(1:ntp) = 0 call encounter_check_all_pltp(param, npl, ntp, pl%rbeg, pl%vbeg, tp%rh, tp%vh, pl%renc, dt, & nenc, index1, index2, lvdotr) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index d360db56b..7913fd1ae 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -140,8 +140,8 @@ module swiftest procedure :: write_info => swiftest_io_netcdf_write_info_body !! Dump contents of particle information metadata to file procedure :: accel_obl => swiftest_obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: el2xv => swiftest_orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors - procedure :: xv2el => swiftest_orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements - procedure :: setup => swiftest_util_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays + procedure :: xv2el => swiftest_orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements + procedure :: setup => swiftest_util_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays procedure :: accel_user => swiftest_user_kick_getacch_body !! Add user-supplied heliocentric accelerations to planets procedure :: append => swiftest_util_append_body !! Appends elements from one structure to another procedure :: dealloc => swiftest_util_dealloc_body !! Deallocates all allocatable arrays diff --git a/src/swiftest/swiftest_orbel.f90 b/src/swiftest/swiftest_orbel.f90 index bf72c91b2..8e351e911 100644 --- a/src/swiftest/swiftest_orbel.f90 +++ b/src/swiftest/swiftest_orbel.f90 @@ -680,7 +680,7 @@ real(DP) pure function swiftest_orbel_fhybrid(e,n) abn = n if(n < 0._DP) abn = -abn - if(abn < 0.636_DP * e -0.6_DP) then + if(abn < 0.636_DP * e - 0.6_DP) then swiftest_orbel_fhybrid = swiftest_orbel_flon(e,n) else swiftest_orbel_fhybrid = swiftest_orbel_fget(e,n) @@ -715,7 +715,7 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v q = 0.0_DP x = [px, py, pz] v = [vx, vy, vz] - r = sqrt(dot_product(x(:), x(:))) + r = .mag.x(:) v2 = dot_product(v(:), v(:)) hvec(:) = x(:) .cross. v(:) h2 = dot_product(hvec(:), hvec(:)) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index ca7761ff0..895a9f34d 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -249,7 +249,7 @@ module function symba_encounter_check_tp(self, param, nbody_system, dt, irec) re lany_encounter = .false. if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody) + associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb) call pl%set_renc(irec) call encounter_check_all_pltp(param, npl, ntp, pl%rh, pl%vb, tp%rh, tp%vb, pl%renc, dt, nenc, index1, index2, lvdotr) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index c811b1968..6e72c3e95 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -27,6 +27,7 @@ module subroutine symba_step_system(self, param, t, dt) real(DP), intent(in) :: dt !! Current stepsize ! Internals logical :: lencounter + type(walltimer) :: timer1,timer2,timer3 !! Object used for computing elapsed wall time select type(pl => self%pl) class is (symba_pl) @@ -48,7 +49,6 @@ module subroutine symba_step_system(self, param, t, dt) end select end select end select - return end subroutine symba_step_system From 048904481d123fb31a2dd1c7f4480b5cff123e17 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 8 Mar 2023 09:38:52 -0500 Subject: [PATCH 6/9] Removed adaptive option, as I plan to determine the best option algorithmically. --- python/swiftest/swiftest/simulation_class.py | 36 +++++--------------- 1 file changed, 8 insertions(+), 28 deletions(-) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index c908d322d..c6342fa6a 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -283,29 +283,19 @@ def __init__(self,read_param: bool = False, execute. If false, will start a new run. If the file given by `output_file_name` exists, it will be replaced when the run is executed. Parameter input file equivalent: `OUT_STAT` - interaction_loops : {"TRIANGULAR","FLAT","ADAPTIVE"}, default "TRIANGULAR" + interaction_loops : {"TRIANGULAR","FLAT"}, default "TRIANGULAR" > *Swiftest Experimental feature* Specifies which algorithm to use for the computation of body-body gravitational forces. * "TRIANGULAR" : Upper-triangular double-loops . * "FLAT" : Body-body interation pairs are flattened into a 1-D array. - * "ADAPTIVE" : Periodically times the TRIANGULAR and FLAT methods and determines which one to use based on - the wall time to complete the loop. *Note:* Using ADAPTIVE means that bit-identical repeatability cannot - be assured, as the choice of algorithm depends on possible external factors that can influence the wall - time calculation. The exact floating-point results of the interaction will be different between the two - algorithm types. Parameter input file equivalent: `INTERACTION_LOOPS` - encounter_check_loops : {"TRIANGULAR","SORTSWEEP","ADAPTIVE"}, default "TRIANGULAR" + encounter_check_loops : {"TRIANGULAR","SORTSWEEP"}, default "TRIANGULAR" > *Swiftest Experimental feature* Specifies which algorithm to use for checking whether bodies are in a close encounter state or not. * "TRIANGULAR" : Upper-triangular double-loops. * "SORTSWEEP" : A Sort-Sweep algorithm is used to reduce the population of potential close encounter bodies. This algorithm is still in development, and does not necessarily speed up the encounter checking. Use with caution. - * "ADAPTIVE" : Periodically times the TRIANGULAR and SORTSWEEP methods and determines which one to use based - on the wall time to complete the loop. *Note:* Using ADAPTIVE means that bit-identical repeatability cannot - be assured, as the choice of algorithm depends on possible external factors that can influence the wall - time calculation. The exact floating-point results of the interaction will be different between the two - algorithm types. Parameter input file equivalent: `ENCOUNTER_CHECK` dask : bool, default False Use Dask to lazily load data (useful for very large datasets) @@ -1072,8 +1062,8 @@ def set_feature(self, rhill_present: bool | None = None, restart: bool | None = None, tides: bool | None = None, - interaction_loops: Literal["TRIANGULAR", "FLAT", "ADAPTIVE"] | None = None, - encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] | None = None, + interaction_loops: Literal["TRIANGULAR", "FLAT"] | None = None, + encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP"] | None = None, encounter_save: Literal["NONE", "TRAJECTORY", "CLOSEST", "BOTH"] | None = None, verbose: bool | None = None, simdir: str | os.PathLike = None, @@ -1128,28 +1118,18 @@ def set_feature(self, Includes big bodies when performing a discard (Swifter only) rhill_present: bool, optional Include the Hill's radius with the input files. - interaction_loops : {"TRIANGULAR","FLAT","ADAPTIVE"}, default "TRIANGULAR" + interaction_loops : {"TRIANGULAR","FLAT"}, default "TRIANGULAR" *Swiftest Experimental feature* Specifies which algorithm to use for the computation of body-body gravitational forces. * "TRIANGULAR" : Upper-triangular double-loops . * "FLAT" : Body-body interation pairs are flattened into a 1-D array. - * "ADAPTIVE" : Periodically times the TRIANGULAR and FLAT methods and determines which one to use based on - the wall time to complete the loop. *Note:* Using ADAPTIVE means that bit-identical repeatability cannot - be assured, as the choice of algorithm depends on possible external factors that can influence the wall - time calculation. The exact floating-point results of the interaction will be different between the two - algorithm types. - encounter_check_loops : {"TRIANGULAR","SORTSWEEP","ADAPTIVE"}, default "TRIANGULAR" + encounter_check_loops : {"TRIANGULAR","SORTSWEEP"}, default "TRIANGULAR" *Swiftest Experimental feature* Specifies which algorithm to use for checking whether bodies are in a close encounter state or not. * "TRIANGULAR" : Upper-triangular double-loops. * "SORTSWEEP" : A Sort-Sweep algorithm is used to reduce the population of potential close encounter bodies. This algorithm is still in development, and does not necessarily speed up the encounter checking. Use with caution. - * "ADAPTIVE" : Periodically times the TRIANGULAR and SORTSWEEP methods and determines which one to use based - on the wall time to complete the loop. *Note:* Using ADAPTIVE means that bit-identical repeatability cannot - be assured, as the choice of algorithm depends on possible external factors that can influence the wall - time calculation. The exact floating-point results of the interaction will be different between the two - algorithm types. tides: bool, optional Turns on tidal model (IN DEVELOPMENT - IGNORED) Yarkovsky: bool, optional @@ -1250,7 +1230,7 @@ def set_feature(self, update_list.append("restart") if interaction_loops is not None: - valid_vals = ["TRIANGULAR", "FLAT", "ADAPTIVE"] + valid_vals = ["TRIANGULAR", "FLAT"] if interaction_loops not in valid_vals: msg = f"{interaction_loops} is not a valid option for interaction loops." msg += f"\nMust be one of {valid_vals}" @@ -1262,7 +1242,7 @@ def set_feature(self, update_list.append("interaction_loops") if encounter_check_loops is not None: - valid_vals = ["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] + valid_vals = ["TRIANGULAR", "SORTSWEEP"] if encounter_check_loops not in valid_vals: msg = f"{encounter_check_loops} is not a valid option for interaction loops." msg += f"\nMust be one of {valid_vals}" From f1eab7c9bf20b08e07af52be05d4454fefd777cc Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 8 Mar 2023 09:46:30 -0500 Subject: [PATCH 7/9] Upper cased the interaction_loops and encounter_check_loops argument values so that the user doesn't have to --- python/swiftest/swiftest/simulation_class.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index c6342fa6a..991758441 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -1231,6 +1231,7 @@ def set_feature(self, if interaction_loops is not None: valid_vals = ["TRIANGULAR", "FLAT"] + interaction_loops = interaction_loops.upper() if interaction_loops not in valid_vals: msg = f"{interaction_loops} is not a valid option for interaction loops." msg += f"\nMust be one of {valid_vals}" @@ -1243,6 +1244,7 @@ def set_feature(self, if encounter_check_loops is not None: valid_vals = ["TRIANGULAR", "SORTSWEEP"] + encounter_check_loops = encounter_check_loops.upper() if encounter_check_loops not in valid_vals: msg = f"{encounter_check_loops} is not a valid option for interaction loops." msg += f"\nMust be one of {valid_vals}" From 2ed8bc06e01e9d3e6ac06364e9c525dd7c0c1e3f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 8 Mar 2023 11:07:12 -0500 Subject: [PATCH 8/9] Added back user control over enounter check algorithm --- src/base/base_module.f90 | 3 -- src/encounter/encounter_check.f90 | 35 ++++++++++---------- src/swiftest/swiftest_io.f90 | 54 +++++++------------------------ 3 files changed, 29 insertions(+), 63 deletions(-) diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index 480587b04..0a9429dbc 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -70,11 +70,8 @@ module base ! The following are used internally, and are not set by the user, but instead are determined by the input value of INTERACTION_LOOPS logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interaction loops - logical :: ladaptive_interactions = .false. !! Adaptive interaction loop is turned on (choose between TRIANGULAR and FLAT based on periodic timing tests) logical :: lencounter_sas_plpl = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters logical :: lencounter_sas_pltp = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters - logical :: ladaptive_encounters_plpl = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests) - logical :: ladaptive_encounters_pltp = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests) ! Logical flags to turn on or off various features of the code logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index aa8da45a5..9be5da7cc 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -28,15 +28,12 @@ module subroutine encounter_check_all_plpl(param, npl, r, v, renc, dt, nenc, ind integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x - ! Internals - !type(interaction_timer), save :: itimer - logical, save :: lfirst = .true. - logical, save :: skipit = .false. ! This will be used to ensure that the sort & sweep subroutine gets called at least once before timing it so that the extent array is nearly sorted when it is timed - integer(I8B) :: nplpl = 0_I8B - - !call encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) - call encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) + if (param%lencounter_sas_plpl) then + call encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) + else + call encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) + end if return end subroutine encounter_check_all_plpl @@ -80,13 +77,16 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, allocate(tmp_param, source=param) - ! Turn off adaptive encounter checks for the pl-pl group - tmp_param%ladaptive_encounters_plpl = .false. - ! Start with the pl-pl group call encounter_check_all_plpl(tmp_param, nplm, rplm, vplm, rencm, dt, nenc, index1, index2, lvdotr) - call encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) + if (param%lencounter_sas_plpl) then + call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, & + plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) + else + call encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, & + plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) + end if if (plmplt_nenc > 0) then ! Consolidate the two lists allocate(itmp(nenc+plmplt_nenc)) @@ -134,13 +134,12 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, vtp, integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x - ! Internals - ! type(interaction_timer), save :: itimer - logical, save :: lfirst = .true. - logical, save :: lsecond = .false. - integer(I8B) :: npltp = 0_I8B - call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + if (param%lencounter_sas_pltp) then + call encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + else + call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + end if return end subroutine encounter_check_all_pltp diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 3e5efa606..9406e44c2 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -2270,72 +2270,42 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i end if select case(trim(adjustl(param%interaction_loops))) - case("ADAPTIVE") - param%ladaptive_interactions = .true. - param%lflatten_interactions = .true. - call swiftest_io_log_start(param, INTERACTION_TIMER_LOG_OUT, "Interaction loop timer logfile") - call swiftest_io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") case("TRIANGULAR") - param%ladaptive_interactions = .false. param%lflatten_interactions = .false. case("FLAT") - param%ladaptive_interactions = .false. param%lflatten_interactions = .true. case default write(*,*) "Unknown value for parameter INTERACTION_LOOPS: -> ",trim(adjustl(param%interaction_loops)) - write(*,*) "Must be one of the following: TRIANGULAR, FLAT, or ADAPTIVE" - write(*,*) "Using default value of ADAPTIVE" - param%interaction_loops = "ADAPTIVE" - param%ladaptive_interactions = .true. - param%lflatten_interactions = .true. - call swiftest_io_log_start(param, INTERACTION_TIMER_LOG_OUT, "Interaction loop timer logfile") - call swiftest_io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") + write(*,*) "Must be one of the following: TRIANGULAR or FLAT" + write(*,*) "Using default value of TRIANGULAR" + param%interaction_loops = "TRIANGULAR" + param%lflatten_interactions = .false. end select select case(trim(adjustl(param%encounter_check_plpl))) - case("ADAPTIVE") - param%ladaptive_encounters_plpl = .true. - param%lencounter_sas_plpl = .true. - call swiftest_io_log_start(param, ENCOUNTER_PLPL_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call swiftest_io_log_one_message(ENCOUNTER_PLPL_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") case("TRIANGULAR") - param%ladaptive_encounters_plpl = .false. param%lencounter_sas_plpl = .false. case("SORTSWEEP") - param%ladaptive_encounters_plpl = .false. param%lencounter_sas_plpl = .true. case default write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLPL: -> ",trim(adjustl(param%encounter_check_plpl)) - write(*,*) "Must be one of the following: TRIANGULAR, SORTSWEEP, or ADAPTIVE" - write(*,*) "Using default value of ADAPTIVE" - param%encounter_check_plpl = "ADAPTIVE" - param%ladaptive_encounters_plpl = .true. - param%lencounter_sas_plpl = .true. - call swiftest_io_log_start(param, ENCOUNTER_PLPL_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call swiftest_io_log_one_message(ENCOUNTER_PLPL_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") + write(*,*) "Must be one of the following: TRIANGULAR or SORTSWEEP" + write(*,*) "Using default value of TRIANGULAR" + param%encounter_check_plpl = "TRIANGULAR" + param%lencounter_sas_plpl = .false. end select select case(trim(adjustl(param%encounter_check_pltp))) - case("ADAPTIVE") - param%ladaptive_encounters_pltp = .true. - param%lencounter_sas_pltp = .true. - call swiftest_io_log_start(param, ENCOUNTER_PLTP_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call swiftest_io_log_one_message(ENCOUNTER_PLTP_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, npltp, metric") case("TRIANGULAR") - param%ladaptive_encounters_pltp = .false. param%lencounter_sas_pltp = .false. case("SORTSWEEP") - param%ladaptive_encounters_pltp = .false. param%lencounter_sas_pltp = .true. case default write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLTP: -> ",trim(adjustl(param%encounter_check_pltp)) - write(*,*) "Must be one of the following: TRIANGULAR, SORTSWEEP, or ADAPTIVE" - write(*,*) "Using default value of ADAPTIVE" - param%encounter_check_pltp = "ADAPTIVE" - param%ladaptive_encounters_pltp = .true. - param%lencounter_sas_pltp = .true. - call swiftest_io_log_start(param, ENCOUNTER_PLTP_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call swiftest_io_log_one_message(ENCOUNTER_PLTP_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, npltp, metric") + write(*,*) "Must be one of the following: TRIANGULAR or SORTSWEEP" + write(*,*) "Using default value of TRIANGULAR" + param%encounter_check_pltp = "TRIANGULAR" + param%lencounter_sas_pltp = .false. end select iostat = 0 From d6e5d088ca134c239b6638022233847d784c8127 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 8 Mar 2023 17:05:38 -0500 Subject: [PATCH 9/9] Fixed bad radius variable name --- src/encounter/encounter_check.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 9be5da7cc..105aaf7fe 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -852,7 +852,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r iend = self%aabb%iend(i) - 1_I8B nbox = iend - ibeg + 1 ii = i - n1 - call encounter_check_all_sweep_one(ii, nbox, r1(1,ii), r1(2,ii), r1(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), & + call encounter_check_all_sweep_one(ii, nbox, r2(1,ii), r2(2,ii), r2(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & renc2(ii), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), &