Numerous small changes and refactoring for clarity and hunting for source of divergence with original Swifter RMVS

src/driftkick/drift.f90

@@ -9,8 +9,7 @@
    integer(I2B), parameter :: NLAG2    = 40
 
 contains
-   module pure subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) 
-      !$omp declare simd
+   module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) 
       !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
       !!
       !! Perform Danby drift for one body, redoing drift with smaller substeps if original accuracy is insufficient
@@ -36,7 +35,7 @@ module pure subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag)
          dttmp = 0.1_DP * dt
          do i = 1, 10
             call drift_dan(mu, x(:), v(:), dttmp, iflag)
-            if (iflag /= 0) return
+            if (iflag /= 0) exit
          end do
       end if
       px = x(1); py = x(2); pz = x(3)

src/helio/helio_drift.f90

@@ -16,9 +16,10 @@ module subroutine helio_drift_pl(self, cb, config, dt)
       class(swiftest_configuration), intent(in)    :: config !! Input collection of 
       real(DP),                      intent(in)    :: dt     !! Stepsize
       ! Internals
-      integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag
       integer(I4B) :: i !! Loop counter
-      real(DP) :: rmag, vmag2, energy, dtp
+      real(DP) :: rmag, vmag2, energy
+      integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag
+      real(DP), dimension(:), allocatable          :: dtp
 
       associate(npl    => self%nbody, &
          xh     => self%xh, &
@@ -30,23 +31,28 @@ module subroutine helio_drift_pl(self, cb, config, dt)
 
          allocate(iflag(npl))
          iflag(:) = 0
-
-         do i = 1, npl
-            if (config%lgr) then
+         allocate(dtp(npl))
+
+         if (config%lgr) then
+            do i = 1,npl
                rmag = norm2(xh(:, i))
                vmag2 = dot_product(vb(:, i),  vb(:, i))
                energy = 0.5_DP * vmag2 - mu(i) / rmag
-               dtp = dt * (1.0_DP + 3 * config%inv_c2 * energy)
-            else
-               dtp = dt
-            end if
-            call drift_one(mu(i), xh(1, i), xh(2, i), xh(3, i), vb(1, i), vb(2, i), vb(3, i), dtp, iflag(i))
-         end do 
+               dtp(i) = dt * (1.0_DP + 3 * config%inv_c2 * energy)
+            end do
+         else
+            dtp(:) = dt
+         end if 
+
+         !!$omp simd
+         call drift_one(mu(1:npl), xh(1,1:npl), xh(2,1:npl), xh(3,1:npl), &
+                                   vb(1,1:npl), vb(2,1:npl), vb(3,1:npl), &
+                                   dtp(1:npl), iflag(1:npl))
          if (any(iflag(1:npl) /= 0)) then
             do i = 1, npl
                write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!"
-               write(*, *) xh
-               write(*, *) vb
+               write(*, *) xh(:,i)
+               write(*, *) vb(:,i)
                write(*, *) " stopping "
                call util_exit(FAILURE)
             end do
@@ -103,33 +109,39 @@ module subroutine helio_drift_tp(self, cb, config, dt)
       class(swiftest_configuration), intent(in)    :: config !! Input collection of 
       real(DP),                      intent(in)    :: dt     !! Stepsize
       ! Internals
-      integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag
       integer(I4B) :: i !! Loop counter
-      real(DP) :: rmag, vmag2, energy, dtp
+      real(DP) :: rmag, vmag2, energy
+      real(DP), dimension(:), allocatable          :: dtp
+      integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag
 
       associate(ntp    => self%nbody, &
          xh     => self%xh, &
-         vb     => self%vb, &
+         vh     => self%vh, &
          status => self%status,&
          mu     => self%mu)
          if (ntp == 0) return
          allocate(iflag(ntp))
+         allocate(dtp(ntp))
          iflag(:) = 0
-         do i = 1,ntp
-            if (status(i) == ACTIVE) then
-               if (config%lgr) then
-                  rmag = norm2(xh(:, i))
-                  vmag2 = dot_product(vb(:, i), vb(:, i))
-                  energy = 0.5_DP * vmag2 - cb%Gmass / rmag
-                  dtp = dt * (1.0_DP + 3 * config%inv_c2 * energy)
-               else
-                  dtp = dt
-               end if
-               call drift_one(mu(i), xh(1, i), xh(2, i), xh(3, i), vb(1, i), vb(2, i), vb(3, i), dtp, iflag(i))
-               if (iflag(i) /= 0) status = DISCARDED_DRIFTERR
-            end if
-         end do
+
+         iflag(:) = 0
+         allocate(dtp(ntp))
+
+         if (config%lgr) then
+            do i = 1,ntp
+               rmag = norm2(xh(:, i))
+               vmag2 = dot_product(vh(:, i),  vh(:, i))
+               energy = 0.5_DP * vmag2 - mu(i) / rmag
+               dtp(i) = dt * (1.0_DP + 3 * config%inv_c2 * energy)
+            end do
+         else
+            dtp(:) = dt
+         end if 
+         call drift_one(mu(1:ntp), xh(1,1:ntp), xh(2,1:ntp), xh(3,1:ntp), &
+                                   vh(1,1:ntp), vh(2,1:ntp), vh(3,1:ntp), &
+                                   dtp(1:ntp), iflag(1:ntp))
          if (any(iflag(1:ntp) /= 0)) then
+            status = DISCARDED_DRIFTERR
             do i = 1, ntp
                if (iflag(i) /= 0) write(*, *) "Particle ", self%name(i), " lost due to error in Danby drift"
             end do

src/modules/rmvs_classes.f90

@@ -29,14 +29,20 @@ module rmvs_classes
       procedure, public :: step          => rmvs_step_system
    end type rmvs_nbody_system
 
+   type, private :: rmvs_interp
+      real(DP), dimension(:, :), allocatable :: x    !! interpolated heliocentric planet position for outer encounter
+      real(DP), dimension(:, :), allocatable :: v    !! interpolated heliocentric planet velocity for outer encounter
+      real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value
+   end type rmvs_interp
+
    !********************************************************************************************************************************
    ! rmvs_cb class definitions and method interfaces
    !*******************************************************************************************************************************
    !> RMVS central body particle class 
    type, public, extends(whm_cb) :: rmvs_cb
       logical                                   :: lplanetocentric = .false.  !! Flag that indicates that the object is a planetocentric set of test particles used for close encounter calculations
-      real(DP), dimension(NDIM) :: xin = [0.0_DP, 0.0_DP, 0.0_DP]
-      real(DP), dimension(NDIM) :: vin = [0.0_DP, 0.0_DP, 0.0_DP]
+      type(rmvs_interp),           dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters
+      type(rmvs_interp),           dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters
    end type rmvs_cb
 
    !********************************************************************************************************************************
@@ -73,12 +79,7 @@ module rmvs_classes
    !********************************************************************************************************************************
    !                                    rmvs_pl class definitions and method interfaces
    !*******************************************************************************************************************************
-
-   type, private :: rmvs_interp
-      real(DP), dimension(:, :), allocatable :: x    !! interpolated heliocentric planet position for outer encounter
-      real(DP), dimension(:, :), allocatable :: v    !! interpolated heliocentric planet velocity for outer encounter
-      real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value
-   end type rmvs_interp 
+
 
    !> RMVS massive body particle class
    type, private, extends(whm_pl) :: rmvs_base_pl

src/modules/swiftest_classes.f90

@@ -356,8 +356,7 @@ module subroutine discard_system(self, config)
          class(swiftest_configuration), intent(in)    :: config  !! Input collection of  configuration parameters 
       end subroutine discard_system
 
-      module pure subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag)
-         !$omp declare simd(drift_one)
+      module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag)
          implicit none
          real(DP), intent(in)       :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift
          real(DP),  intent(inout)   :: px, py, pz, vx, vy, vz  !! Position and velocity of body to drift

src/rmvs/rmvs_getacch.f90

@@ -28,34 +28,37 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh)
          if (tp%lplanetocentric) then  ! This is a close encounter step, so any accelerations requiring heliocentric position values
                                        ! must be handeled outside the normal WHM method call
             select type(pl)
-               class is (rmvs_pl)
+            class is (rmvs_pl)
+            select type (cb)
+            class is (rmvs_cb)
+               associate(xpc => pl%xh, xpct => self%xh)
                   allocate(xh_original, source=tp%xh)
                   config_planetocen = config
                   ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter
                   config_planetocen%loblatecb = .false.
                   config_planetocen%lextra_force = .false.
                   config_planetocen%lgr = .false.
                   ! Now compute the planetocentric values of acceleration
-                  call whm_getacch_tp(tp, cb, pl%planetocentric(ipleP)%pl, config_planetocen, t, pl%xh)
+                  call whm_getacch_tp(tp, cb, pl, config_planetocen, t, pl%xh)
 
                   ! Now compute any heliocentric values of acceleration 
                   if (tp%lfirst) then
                      do i = 1, ntp
-                        tp%xheliocentric(:,i) = tp%xh(:,i) + pl%inner(enc_index - 1)%x(:,ipleP)
+                        tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(enc_index - 1)%x(:,1)
                      end do
                   else
                      do i = 1, ntp
-                        tp%xheliocentric(:,i) = tp%xh(:,i) + pl%inner(enc_index    )%x(:,ipleP)
+                        tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(enc_index    )%x(:,1)
                      end do
                   end if
                   ! Swap the planetocentric and heliocentric position vectors
                   tp%xh(:,:) = tp%xheliocentric(:,:)
                   if (config%loblatecb) then
                      ! Put in the current encountering planet's oblateness acceleration as the central body's
                      if (tp%lfirst) then
-                        cb_heliocentric%aobl(:) = pl%inner(enc_index - 1)%aobl(:, ipleP)
+                        cb_heliocentric%aobl(:) = cb%inner(enc_index - 1)%aobl(:,1)
                      else
-                        cb_heliocentric%aobl(:) = pl%inner(enc_index    )%aobl(:, ipleP)
+                        cb_heliocentric%aobl(:) = cb%inner(enc_index    )%aobl(:,1)
                      end if
                      call tp%obl_acc(cb_heliocentric)
                   end if
@@ -64,6 +67,8 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh)
                   if (config%lgr) call tp%gr_getacch(cb_heliocentric, config)
 
                   call move_alloc(xh_original, tp%xh)
+               end associate
+            end select
             end select
          else ! Not a close encounter, so just proceded with the standard WHM method
             call whm_getacch_tp(tp, cb, pl, config, t, pl%xh)

src/rmvs/rmvs_setup.f90

@@ -20,17 +20,17 @@ module subroutine rmvs_setup_pl(self,n)
          call whm_setup_pl(pl, n) 
          if (n <= 0) return
 
+         allocate(pl%outer(0:NTENC))
+         allocate(pl%inner(0:NTPHENC))
          if (.not.pl%lplanetocentric) then
             allocate(pl%nenc(n))
             pl%nenc(:)         = 0
 
             ! Set up inner and outer planet interpolation vector storage containers
-            allocate(pl%outer(0:NTENC))
             do i = 0, NTENC
                allocate(pl%outer(i)%x(NDIM, n))
                allocate(pl%outer(i)%v(NDIM, n))
             end do
-            allocate(pl%inner(0:NTPHENC))
             do i = 0, NTPHENC
                allocate(pl%inner(i)%x(NDIM, n))
                allocate(pl%inner(i)%v(NDIM, n))