Simplified the encounter check code

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

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
 

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,28 +55,24 @@ 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
             else
                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