From 567e5cff5d4cffb75b37532b307085b39624f5fc Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 10 Sep 2021 18:09:54 -0400 Subject: [PATCH] Restructured expensive loops to try to profile and get better performance --- src/drift/drift.f90 | 17 +++--- src/gr/gr.f90 | 2 +- src/helio/helio_gr.f90 | 4 +- src/helio/helio_kick.f90 | 4 +- src/kick/kick.f90 | 82 +++++++++++++++++++---------- src/modules/helio_classes.f90 | 4 +- src/modules/rmvs_classes.f90 | 12 +++-- src/modules/swiftest_classes.f90 | 43 +++++++++------ src/modules/symba_classes.f90 | 3 +- src/modules/whm_classes.f90 | 4 +- src/rmvs/rmvs_encounter_check.f90 | 66 +++++++++++++---------- src/rmvs/rmvs_kick.f90 | 1 - src/symba/symba_encounter_check.f90 | 48 ++++++++++------- src/symba/symba_kick.f90 | 4 +- src/whm/whm_gr.f90 | 5 +- 15 files changed, 183 insertions(+), 116 deletions(-) diff --git a/src/drift/drift.f90 b/src/drift/drift.f90 index 79744c0f3..e929a95fa 100644 --- a/src/drift/drift.f90 +++ b/src/drift/drift.f90 @@ -43,7 +43,7 @@ module subroutine drift_body(self, system, param, dt) end subroutine drift_body - module pure subroutine drift_all(mu, x, v, n, param, dt, mask, iflag) + module subroutine drift_all(mu, x, v, n, param, dt, lmask, iflag) !! author: David A. Minton !! !! Loop bodies and call Danby drift routine on all bodies for the given position and velocity vector. @@ -55,9 +55,9 @@ module pure subroutine drift_all(mu, x, v, n, param, dt, mask, iflag) real(DP), dimension(:), intent(in) :: mu !! Vector of gravitational constants real(DP), dimension(:,:), intent(inout) :: x, v !! Position and velocity vectors integer(I4B), intent(in) :: n !! number of bodies - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: dt !! Stepsize - logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift. + logical, dimension(:), intent(in) :: lmask !! Logical mask of size self%nbody that determines which bodies to drift. integer(I4B), dimension(:), intent(out) :: iflag !! Vector of error flags. 0 means no problem ! Internals integer(I4B) :: i @@ -68,18 +68,21 @@ module pure subroutine drift_all(mu, x, v, n, param, dt, mask, iflag) allocate(dtp(n)) if (param%lgr) then - do concurrent(i = 1:n, mask(i)) + do concurrent(i = 1:n, lmask(i)) rmag = norm2(x(:, i)) vmag2 = dot_product(v(:, i), v(:, i)) energy = 0.5_DP * vmag2 - mu(i) / rmag dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy) end do else - where(mask(1:n)) dtp(1:n) = dt + where(lmask(1:n)) dtp(1:n) = dt end if - do concurrent(i = 1:n, mask(i)) - call drift_one(mu(i), x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), dtp(i), iflag(i)) + !$omp parallel do default(private) & + !$omp shared(n, lmask, mu, x, v, dtp, iflag) + do i = 1,n + if (lmask(i)) call drift_one(mu(i), x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), dtp(i), iflag(i)) end do + !$omp end parallel do return end subroutine drift_all diff --git a/src/gr/gr.f90 b/src/gr/gr.f90 index 0c0333907..71d50f448 100644 --- a/src/gr/gr.f90 +++ b/src/gr/gr.f90 @@ -43,7 +43,7 @@ module pure subroutine gr_kick_getaccb_ns_body(self, system, param) end subroutine gr_kick_getaccb_ns_body - module subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) + module pure subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) !! author: David A. Minton !! !! Compute relativisitic accelerations of massive bodies diff --git a/src/helio/helio_gr.f90 b/src/helio/helio_gr.f90 index 4902c45b8..5d2936324 100644 --- a/src/helio/helio_gr.f90 +++ b/src/helio/helio_gr.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine helio_gr_kick_getacch_pl(self, param) + module pure subroutine helio_gr_kick_getacch_pl(self, param) !! author: David A. Minton !! !! Compute relativisitic accelerations of massive bodies @@ -30,7 +30,7 @@ module subroutine helio_gr_kick_getacch_pl(self, param) end subroutine helio_gr_kick_getacch_pl - module subroutine helio_gr_kick_getacch_tp(self, param) + module pure subroutine helio_gr_kick_getacch_tp(self, param) !! author: David A. Minton !! !! Compute relativisitic accelerations of test particles diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index 8e7571a0b..b06abd581 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -107,7 +107,9 @@ module subroutine helio_kick_vb_pl(self, system, param, t, dt, lbeg) call pl%set_beg_end(xend = pl%xh) end if do concurrent(i = 1:npl, pl%lmask(i)) - pl%vb(:, i) = pl%vb(:, i) + pl%ah(:, i) * dt + pl%vb(1, i) = pl%vb(1, i) + pl%ah(1, i) * dt + pl%vb(2, i) = pl%vb(2, i) + pl%ah(2, i) * dt + pl%vb(3, i) = pl%vb(3, i) + pl%ah(3, i) * dt end do end associate diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 5c57918e2..9f678dcae 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -11,7 +11,7 @@ module subroutine kick_getacch_int_pl(self) !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f90 implicit none ! Arguments - class(swiftest_pl), intent(inout) :: self + 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) @@ -19,7 +19,7 @@ module subroutine kick_getacch_int_pl(self) end subroutine kick_getacch_int_pl - module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) + module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies @@ -28,32 +28,14 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int_tp.f90 implicit none ! Arguments - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors integer(I4B), intent(in) :: npl !! Number of active massive bodies - ! Internals - integer(I4B) :: i, j if ((self%nbody == 0) .or. (npl == 0)) return - associate(tp => self, ntp => self%nbody) - do i = 1, ntp - if (tp%lmask(i)) then - do j = 1, npl - block - real(DP) :: rji2 - real(DP) :: xr, yr, zr - xr = tp%xh(1, i) - xhp(1, j) - yr = tp%xh(2, i) - xhp(1, j) - zr = tp%xh(3, i) - xhp(1, j) - rji2 = xr**2 + yr**2 + zr**2 - call kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i)) - end block - end do - end if - end do - end associate + call kick_getacch_int_all_tp(self%nbody, npl, self%xh, xhp, GMpl, self%lmask, self%ah) return end subroutine kick_getacch_int_tp @@ -76,10 +58,11 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array ! Internals integer(I8B) :: k - real(DP) :: rji2, rlim2 - real(DP) :: xr, yr, zr - integer(I4B) :: i, j real(DP), dimension(NDIM,npl) :: ahi, ahj + integer(I4B) :: i, j + real(DP) :: rji2, rlim2 + real(DP) :: xr, yr, zr + ahi(:,:) = 0.0_DP ahj(:,:) = 0.0_DP @@ -98,12 +81,54 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) end do !$omp end parallel do - - acc(:,1:npl) = acc(:,1:npl) + ahi(:,1:npl) + ahj(:,1:npl) + + do concurrent(i = 1:npl) + acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) + end do + return end subroutine kick_getacch_int_all_pl + module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies with parallelisim + !! + !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int_tp.f99 + implicit none + integer(I4B), intent(in) :: ntp !! Number of test particles + integer(I4B), intent(in) :: npl !! Number of massive bodies + real(DP), dimension(:,:), intent(in) :: xtp !! Test particle position vector array + real(DP), dimension(:,:), intent(in) :: xpl !! Massive body particle position vector array + real(DP), dimension(:), intent(in) :: GMpl !! Array of massive body G*mass + logical, dimension(:), intent(in) :: lmask !! Logical mask indicating which test particles should be computed + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + ! Internals + real(DP) :: rji2 + real(DP) :: xr, yr, zr + integer(I4B) :: i, j + + !$omp parallel do default(private)& + !$omp shared(npl, ntp, lmask, xtp, xpl, acc) + do i = 1, ntp + if (lmask(i)) then + do j = 1, npl + xr = xtp(1, i) - xpl(1, j) + yr = xtp(2, i) - xpl(2, j) + zr = xtp(3, i) - xpl(3, j) + rji2 = xr**2 + yr**2 + zr**2 + call kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl(i), acc(1,i), acc(2,i), acc(3,i)) + end do + end if + end do + !$omp end parallel do + + return + end subroutine kick_getacch_int_all_tp + + module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) !! author: David A. Minton !! @@ -150,10 +175,11 @@ module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, a ! Internals real(DP) :: fac - fac = GMpl / (rji2 * sqrt(rji2)) + fac = GMpl * sqrt(rji2**(-3)) ax = ax - fac * xr ay = ay - fac * yr az = az - fac * zr + return end subroutine kick_getacch_int_one_tp diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 1c93bdf6b..3a7de3814 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -105,14 +105,14 @@ module subroutine helio_drift_linear_tp(self, cb, dt, lbeg) logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step end subroutine helio_drift_linear_tp - module subroutine helio_gr_kick_getacch_pl(self, param) + module pure subroutine helio_gr_kick_getacch_pl(self, param) use swiftest_classes, only : swiftest_parameters implicit none class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine helio_gr_kick_getacch_pl - module subroutine helio_gr_kick_getacch_tp(self, param) + module pure subroutine helio_gr_kick_getacch_tp(self, param) use swiftest_classes, only : swiftest_parameters implicit none class(helio_tp), intent(inout) :: self !! Helio massive body particle data structure diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 0b837ec0a..9ec0c8b86 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -103,11 +103,15 @@ module rmvs_classes end type rmvs_pl interface - module elemental function rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) result(lflag) + module pure subroutine rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit, lencounter, lvdotr) implicit none - real(DP), intent(in) :: r2, v2, vdotr, dt, r2crit - logical :: lflag - end function rmvs_chk_ind + real(DP), intent(in) :: xr, yr, zr !! Relative distance vector components + real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components + real(DP), intent(in) :: dt !! Step size + real(DP), intent(in) :: r2crit !! Square of the critical encounter distance + logical, intent(out) :: lencounter !! Flag indicating that an encounter has occurred + logical, intent(out) :: lvdotr !! Logical flag indicating the direction of the v .dot. r vector + end subroutine rmvs_chk_ind module subroutine rmvs_discard_tp(self, system, param) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 82d684f1a..8c7376c88 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -534,14 +534,14 @@ module subroutine discard_tp(self, system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine discard_tp - module pure subroutine drift_all(mu, x, v, n, param, dt, mask, iflag) + module subroutine drift_all(mu, x, v, n, param, dt, lmask, iflag) implicit none real(DP), dimension(:), intent(in) :: mu !! Vector of gravitational constants real(DP), dimension(:,:), intent(inout) :: x, v !! Position and velocity vectors integer(I4B), intent(in) :: n !! number of bodies class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: dt !! Stepsize - logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift. + logical, dimension(:), intent(in) :: lmask !! Logical mask of size self%nbody that determines which bodies to drift. integer(I4B), dimension(:), intent(out) :: iflag !! Vector of error flags. 0 means no problem end subroutine drift_all @@ -581,7 +581,7 @@ module pure subroutine gr_kick_getaccb_ns_body(self, system, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine gr_kick_getaccb_ns_body - module subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) + module pure subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) implicit none real(DP), dimension(:), intent(in) :: mu !! Gravitational constant real(DP), dimension(:,:), intent(in) :: x !! Position vectors @@ -823,6 +823,19 @@ module subroutine io_write_hdr_system(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_hdr_system + module subroutine kick_getacch_int_pl(self) + implicit none + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + end subroutine kick_getacch_int_pl + + module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) + implicit none + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle + real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses + real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors + 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) implicit none integer(I4B), intent(in) :: npl !! Number of massive bodies @@ -834,6 +847,17 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array end subroutine kick_getacch_int_all_pl + module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc) + implicit none + integer(I4B), intent(in) :: ntp !! Number of test particles + integer(I4B), intent(in) :: npl !! Number of massive bodies + real(DP), dimension(:,:), intent(in) :: xtp !! Test particle position vector array + real(DP), dimension(:,:), intent(in) :: xpl !! Massive body particle position vector array + real(DP), dimension(:), intent(in) :: GMpl !! Array of massive body G*mass + 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 kick_getacch_int_all_tp + module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) implicit none real(DP), intent(in) :: rji2 !! Square of distance between the two bodies @@ -852,19 +876,6 @@ module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, ax, ay, a real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle end subroutine kick_getacch_int_one_tp - module subroutine kick_getacch_int_pl(self) - implicit none - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - end subroutine kick_getacch_int_pl - - module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) - implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle - real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses - real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors - integer(I4B), intent(in) :: npl !! Number of active massive bodies - end subroutine kick_getacch_int_tp - module subroutine netcdf_close(self, param) implicit none class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 21986b920..bd4b14486 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -6,7 +6,6 @@ module symba_classes use swiftest_globals use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_encounter, swiftest_particle_info, netcdf_parameters use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system - use rmvs_classes, only : rmvs_chk_ind use fraggle_classes, only : fraggle_colliders, fraggle_fragments implicit none public @@ -257,7 +256,7 @@ module subroutine symba_drift_tp(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine symba_drift_tp - module pure elemental subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) + module pure subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) implicit none real(DP), intent(in) :: xr, yr, zr, vxr, vyr, vzr real(DP), intent(in) :: rhill1, rhill2, dt diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index 9675f289d..bd17f98e0 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -159,14 +159,14 @@ module subroutine whm_kick_vh_tp(self, system, param, t, dt, lbeg) logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. end subroutine whm_kick_vh_tp - module subroutine whm_gr_kick_getacch_pl(self, param) + module pure subroutine whm_gr_kick_getacch_pl(self, param) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine whm_gr_kick_getacch_pl - module subroutine whm_gr_kick_getacch_tp(self, param) + module pure subroutine whm_gr_kick_getacch_tp(self, param) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none class(whm_tp), intent(inout) :: self !! WHM test particle data structure diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 4cb119e48..71465870d 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -18,10 +18,12 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) logical :: lencounter !! Returns true if there is at least one close encounter ! Internals integer(I4B) :: i, j - real(DP) :: r2, v2, vdotr - real(DP), dimension(NDIM) :: xr, vr + real(DP) :: xr, yr, zr, vxr, vyr, vzr real(DP), dimension(system%pl%nbody) :: r2crit - logical :: lflag + logical :: lflag, lvdotr + + ! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we + ! can simply fail the attempt and try again. So we need to turn off any floating point exception halting modes temporarily lencounter = .false. if (self%nbody == 0) return @@ -31,27 +33,32 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) associate(tp => self, ntp => self%nbody, npl => pl%nbody, rts => system%rts) r2crit(1:npl) = (rts * pl%rhill(1:npl))**2 tp%plencP(1:ntp) = 0 + !$omp parallel do default(private)& + !$omp shared(npl, ntp, tp, pl, dt, r2crit) do j = 1, npl do i = 1, ntp if ((.not.tp%lmask(i)).or.(tp%plencP(i) /= 0)) cycle - xr(:) = tp%xh(:, i) - pl%xbeg(:, j) - vr(:) = tp%vh(:, i) - pl%vbeg(:, j) - r2 = dot_product(xr(:), xr(:)) - v2 = dot_product(vr(:), vr(:)) - vdotr = dot_product(vr(:), xr(:)) - lflag = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit(j)) + xr = tp%xh(1, i) - pl%xbeg(1, j) + yr = tp%xh(2, i) - pl%xbeg(2, j) + zr = tp%xh(3, i) - pl%xbeg(3, j) + vxr = tp%vh(1, i) - pl%vbeg(1, j) + vyr = tp%vh(2, i) - pl%vbeg(2, j) + vzr = tp%vh(3, i) - pl%vbeg(3, j) + call rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit(j), lflag, lvdotr) if (lflag) tp%plencP(i) = j end do pl%nenc(j) = count(tp%plencP(1:ntp) == j) end do + !$omp end parallel do lencounter = any(pl%nenc(1:npl) > 0) end associate end select + return end function rmvs_encounter_check_tp - module elemental function rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) result(lflag) + module pure subroutine rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit, lencounter, lvdotr) !! author: David A. Minton !! !! Determine whether a test particle and planet are having or will have an encounter within the next time step @@ -60,28 +67,33 @@ module elemental function rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) result(lflag) !! Adapted from Hal Levison's Swift routine rmvs_chk_ind.f implicit none ! Arguments - real(DP), intent(in) :: r2, v2, vdotr, dt, r2crit - logical :: lflag + real(DP), intent(in) :: xr, yr, zr !! Relative distance vector components + real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components + real(DP), intent(in) :: dt !! Step size + real(DP), intent(in) :: r2crit !! Square of the critical encounter distance + logical, intent(out) :: lencounter !! Flag indicating that an encounter has occurred + logical, intent(out) :: lvdotr !! Logical flag indicating the direction of the v .dot. r vector ! Internals - real(DP) :: tmin, r2min + real(DP) :: r2min, r2, v2, vdotr + + r2 = xr**2 + yr**2 + zr**2 + lencounter = (r2 < r2crit) + if (lencounter) return + + vdotr = vxr * xr + vyr * yr + vzr * zr + lvdotr = (vdotr < 0.0_DP) + if (.not.lvdotr) return + + v2 = vxr**2 + vyr**2 + vzr**2 - lflag = .false. - if (r2 < r2crit) then - lflag = .true. + if (-vdotr < v2 * dt) then + r2min = r2 - vdotr**2 / v2 else - if (vdotr < 0.0_DP) then - tmin = -vdotr / v2 - if (tmin < dt) then - r2min = r2 - vdotr**2 / v2 - else - r2min = r2 + 2 * vdotr * dt + v2 * dt**2 - end if - r2min = min(r2min, r2) - lflag = (r2min <= r2crit) - end if + r2min = r2 + 2 * vdotr * dt + v2 * dt**2 end if + lencounter = (r2min <= r2crit) return - end function rmvs_chk_ind + end subroutine rmvs_chk_ind end submodule s_rmvs_chk diff --git a/src/rmvs/rmvs_kick.f90 b/src/rmvs/rmvs_kick.f90 index 5d9514763..6a8c0786e 100644 --- a/src/rmvs/rmvs_kick.f90 +++ b/src/rmvs/rmvs_kick.f90 @@ -3,7 +3,6 @@ contains module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) - !! author: David A. Minton !! !! Compute the oblateness acceleration in the inner encounter region with planets diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 6f6010047..b99184465 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -32,8 +32,6 @@ subroutine symba_encounter_check_all(nplplm, k_plpl, x, v, rhill, dt, irec, len vzr = v(3, j) - v(3, i) rhill1 = rhill(i) rhill2 = rhill(j) - lencounter(k) = .false. - loc_lvdotr(k) = .false. call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter(k), loc_lvdotr(k)) end do !$omp end parallel do @@ -58,7 +56,8 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc ! Internals integer(I8B) :: k, nplplm integer(I4B) :: i, j, nenc - logical, dimension(:), allocatable :: lencounter, loc_lvdotr + logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr + integer(I8B), dimension(:), allocatable :: kidx if (self%nbody == 0) return @@ -76,9 +75,13 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc lany_encounter = nenc > 0 if (lany_encounter) then associate(plplenc_list => system%plplenc_list) + allocate(lvdotr(nenc)) + allocate(kidx(nenc)) + lvdotr(:) = pack(loc_lvdotr(:), lencounter(:)) + kidx(:) = pack([(k, k = 1_I8B, nplplm)], lencounter(:)) + call move_alloc(lvdotr, plplenc_list%lvdotr) + call move_alloc(kidx, plplenc_list%kidx) call plplenc_list%resize(nenc) - plplenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(1:nplplm), lencounter(1:nplplm)) - plplenc_list%kidx(1:nenc) = pack([(k, k = 1_I8B, nplplm)], lencounter(1:nplplm)) deallocate(lencounter, loc_lvdotr) plplenc_list%index1(1:nenc) = pl%k_plpl(1,plplenc_list%kidx(1:nenc)) plplenc_list%index2(1:nenc) = pl%k_plpl(2,plplenc_list%kidx(1:nenc)) @@ -120,11 +123,12 @@ module function symba_encounter_check(self, system, dt, irec) result(lany_encoun integer(I4B), intent(in) :: irec !! Current recursion level logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals - integer(I4B) :: i, j,k + integer(I4B) :: i, j, k, lidx, nenc_enc real(DP), dimension(NDIM) :: xr, vr logical :: lencounter, isplpl real(DP) :: rlim2, rji2 logical, dimension(:), allocatable :: lencmask + integer(I4B), dimension(:), allocatable :: encidx lany_encounter = .false. if (self%nenc == 0) return @@ -142,8 +146,14 @@ module function symba_encounter_check(self, system, dt, irec) result(lany_encoun class is (symba_tp) allocate(lencmask(self%nenc)) lencmask(:) = (self%status(1:self%nenc) == ACTIVE) .and. (self%level(1:self%nenc) == irec - 1) - if (.not.any(lencmask(:))) return - do concurrent(k = 1:self%nenc, lencmask(k)) + nenc_enc = count(lencmask(:)) + if (nenc_enc == 0) return + + allocate(encidx(nenc_enc)) + encidx(:) = pack([(k, k = 1, self%nenc)], lencmask(:)) + + do lidx = 1, nenc_enc + k = encidx(lidx) i = self%index1(k) j = self%index2(k) if (isplpl) then @@ -211,8 +221,12 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc do j = 1, npl do i = 1, ntp - xr(:) = tp%xh(:, i) - pl%xh(:, j) - vr(:) = tp%vh(:, i) - pl%vh(:, j) + xr(1) = tp%xh(1, i) - pl%xh(1, j) + xr(2) = tp%xh(2, i) - pl%xh(2, j) + xr(3) = tp%xh(3, i) - pl%xh(3, j) + vr(1) = tp%vh(1, i) - pl%vh(1, j) + vr(2) = tp%vh(2, i) - pl%vh(2, j) + vr(3) = tp%vh(3, i) - pl%vh(3, j) call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(j), 0.0_DP, dt, irec, lencounter(i,j), loc_lvdotr(i,j)) end do end do @@ -252,7 +266,7 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc end function symba_encounter_check_tp - module pure elemental subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) + module pure subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) !! author: David A. Minton !! !! Check for an encounter. @@ -266,15 +280,11 @@ module pure elemental subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, integer(I4B), intent(in) :: irec logical, intent(out) :: lencounter, lvdotr ! Internals - real(DP) :: r2, v2, rcrit, r2crit, vdotr + real(DP) :: r2crit - rcrit = (rhill1 + rhill2)*RHSCALE*(RSHELL**(irec)) - r2crit = rcrit**2 - r2 = xr**2 + yr**2 + zr**2 - v2 = vxr**2 + vyr**2 + vzr**2 - vdotr = xr * vxr + yr * vyr + zr * vzr - lencounter = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) - lvdotr = (vdotr < 0.0_DP) + r2crit = (rhill1 + rhill2)*RHSCALE*(RSHELL**(irec)) + r2crit = r2crit**2 + call rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit, lencounter, lvdotr) return end subroutine symba_encounter_check_one diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 2c074cb0d..8a70e3c21 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -90,7 +90,7 @@ module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) do k = 1, npltpenc i = pltpenc_list%index1(k) j = pltpenc_list%index2(k) - if (tp%lmask(j)) THEN + if (tp%lmask(j)) then if (lbeg) then dx(:) = tp%xh(:,j) - pl%xbeg(:,i) else @@ -99,7 +99,7 @@ module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) rjj = dot_product(dx(:), dx(:)) fac = pl%Gmass(i) / (rjj * sqrt(rjj)) tp%ah(:,j) = tp%ah(:,j) + fac * dx(:) - end IF + end if end do end associate end select diff --git a/src/whm/whm_gr.f90 b/src/whm/whm_gr.f90 index bfba5c6a2..bbe6ca30b 100644 --- a/src/whm/whm_gr.f90 +++ b/src/whm/whm_gr.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine whm_gr_kick_getacch_pl(self, param) + module pure subroutine whm_gr_kick_getacch_pl(self, param) !! author: David A. Minton !! !! Compute relativisitic accelerations of massive bodies @@ -35,7 +35,7 @@ module subroutine whm_gr_kick_getacch_pl(self, param) end subroutine whm_gr_kick_getacch_pl - module subroutine whm_gr_kick_getacch_tp(self, param) + module pure subroutine whm_gr_kick_getacch_tp(self, param) !! author: David A. Minton !! !! Compute relativisitic accelerations of test particles @@ -86,6 +86,7 @@ module pure subroutine whm_gr_p4_pl(self, param, dt) return end subroutine whm_gr_p4_pl + module pure subroutine whm_gr_p4_tp(self, param, dt) !! author: David A. Minton !!