From d28560c60ebe00e7e654c1da8b289e160b053b0e Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 2 Apr 2021 20:19:11 -0400 Subject: [PATCH] Simplified the encounter check code --- .../mars_ejecta/profmaker.sh | 2 +- src/modules/rmvs_classes.f90 | 5 +- src/rmvs/rmvs_encounter_check.f90 | 56 ++++++++----------- 3 files changed, 25 insertions(+), 38 deletions(-) diff --git a/examples/rmvs_swifter_comparison/mars_ejecta/profmaker.sh b/examples/rmvs_swifter_comparison/mars_ejecta/profmaker.sh index cb31cb5b4..9b1adcd8c 100755 --- a/examples/rmvs_swifter_comparison/mars_ejecta/profmaker.sh +++ b/examples/rmvs_swifter_comparison/mars_ejecta/profmaker.sh @@ -1,2 +1,2 @@ #!/bin/bash -gprof ./swifter_symba_ringmoons | /home/daminton/git/gprof2dot/gprof2dot.py | dot -Tpng -o output.png +gprof ./swiftest_driver | /home/daminton/git/gprof2dot/gprof2dot.py | dot -Tpng -o swiftest_profile.png diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 70330f1b7..c0b51660d 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -221,10 +221,9 @@ module subroutine rmvs_step_out2(cb, pl, tp, t, dt, index, config) class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters end subroutine rmvs_step_out2 - module function rmvs_chk_ind(xr, vr, dt, r2crit) result(lflag) + module elemental function rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) result(lflag) implicit none - real(DP), intent(in) :: dt, r2crit - real(DP), dimension(:), intent(in) :: xr, vr + real(DP), intent(in) :: r2, v2, vdotr, dt, r2crit logical :: lflag end function rmvs_chk_ind diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index ebfc2095c..63ead1338 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -17,44 +17,36 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter real(DP), intent(in) :: rts !! fraction of Hill's sphere radius to use as radius of encounter regio logical :: lencounter !! Returns true if there is at least one close encounter ! Internals - integer(I4B) :: i, j, k, nenc - real(DP) :: r2crit + integer(I4B) :: i, j + real(DP) :: r2, v2, vdotr real(DP), dimension(NDIM) :: xr, vr - integer(I4B) :: tpencPindex - logical :: lflag - logical, save :: lfirst = .true. + real(DP), dimension(pl%nbody) :: r2crit associate(tp => self, ntp => self%nbody, npl => pl%nbody, rhill => pl%rhill, & xht => self%xh, vht => self%vh, xbeg => self%xbeg, vbeg => self%vbeg) if (.not.allocated(pl%encmask)) allocate(pl%encmask(ntp, npl)) pl%encmask(:,:) = .false. - lencounter = .false. - pl%nenc(:) = 0 + r2crit(:) = (rts * rhill(:))**2 do i = 1, ntp - if (tp%status(i) == ACTIVE) then - tp%plencP(i) = 0 - lflag = .false. - - do j = 1, npl - r2crit = (rts * rhill(j))**2 - xr(:) = xht(:, i) - xbeg(:, j) - vr(:) = vht(:, i) - vbeg(:, j) - lflag = rmvs_chk_ind(xr(:), vr(:), dt, r2crit) - if (lflag) then - lencounter = .true. - pl%encmask(i,j) = .true. - pl%nenc(j) = pl%nenc(j) + 1 - tp%plencP(i) = j - exit - end if - end do - end if + if (tp%status(i) /= ACTIVE) cycle + tp%plencP(i) = 0 + do j = 1, npl + xr(:) = xht(:, i) - xbeg(:, j) + vr(:) = vht(:, i) - vbeg(:, j) + r2 = dot_product(xr(:), xr(:)) + v2 = dot_product(vr(:), vr(:)) + vdotr = dot_product(vr(:), xr(:)) + pl%encmask(i,j) = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit(j)) + if (pl%encmask(i,j)) tp%plencP(i) = j + end do end do + pl%nenc(:) = count(pl%encmask(:,:), dim=1) + lencounter = any(pl%nenc(:) > 0) end associate return end function rmvs_encounter_check_tp - module function rmvs_chk_ind(xr, vr, dt, r2crit) result(lflag) + module elemental function rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) result(lflag) !! author: David A. Minton !! !! Determine whether a test particle and planet are having or will have an encounter within the next time step @@ -63,20 +55,16 @@ module function rmvs_chk_ind(xr, vr, dt, r2crit) result(lflag) !! Adapted from Hal Levison's Swift routine rmvs_chk_ind.f implicit none ! Arguments - real(DP), intent(in) :: dt, r2crit - real(DP), dimension(:), intent(in) :: xr, vr - logical :: lflag + real(DP), intent(in) :: r2, v2, vdotr, dt, r2crit + logical :: lflag ! Internals - real(DP) :: r2, v2, vdotr, tmin, r2min + real(DP) :: tmin, r2min lflag = .false. - r2 = dot_product(xr(:), xr(:)) if (r2 < r2crit) then lflag = .true. else - vdotr = dot_product(vr(:), xr(:)) if (vdotr < 0.0_DP) then - v2 = dot_product(vr(:), vr(:)) tmin = -vdotr / v2 if (tmin < dt) then r2min = r2 - vdotr**2 / v2 @@ -84,7 +72,7 @@ module function rmvs_chk_ind(xr, vr, dt, r2crit) result(lflag) r2min = r2 + 2 * vdotr * dt + v2 * dt**2 end if r2min = min(r2min, r2) - if (r2min <= r2crit) lflag = .true. + lflag = (r2min <= r2crit) end if end if