Restructured expensive loops to try to profile and get better perform…

…ance

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

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

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

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
 

src/kick/kick.f90

@@ -11,15 +11,15 @@ 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)
 
       return
    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
 

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

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

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

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

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