Merge branch 'OOPSymba' into OOPTides

examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb

@@ -81,8 +81,8 @@
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x2b768e2c9d10>,\n",
-       " <matplotlib.lines.Line2D at 0x2b768e2b9c90>]"
+       "[<matplotlib.lines.Line2D at 0x2af2269e12d0>,\n",
+       " <matplotlib.lines.Line2D at 0x2af2269eeb10>]"
       ]
      },
      "execution_count": 6,

src/gr/gr.f90

@@ -1,7 +1,7 @@
 submodule(swiftest_classes) s_gr
    use swiftest
 contains
-   module pure subroutine gr_getaccb_ns_body(self, system, param) 
+   module pure subroutine gr_kick_getaccb_ns_body(self, system, param) 
       !! author: David A. Minton
       !!
       !! Add relativistic correction acceleration for non-symplectic integrators.
@@ -11,7 +11,7 @@ module pure subroutine gr_getaccb_ns_body(self, system, param)
       !!    Quinn, T.R., Tremaine, S., Duncan, M., 1991. A three million year integration of the earth’s orbit. 
       !!       AJ 101, 2287–2305. https://doi.org/10.1086/115850
       !!
-      !! Adapted from David A. Minton's Swifter routine routine gr_getaccb_ns.f90
+      !! Adapted from David A. Minton's Swifter routine routine gr_kick_getaccb_ns.f90
       implicit none
       ! Arguments
       class(swiftest_body),         intent(inout) :: self   !! Swiftest generic body object
@@ -41,7 +41,7 @@ module pure subroutine gr_getaccb_ns_body(self, system, param)
 
       return
 
-   end subroutine gr_getaccb_ns_body
+   end subroutine gr_kick_getaccb_ns_body
 
    module pure subroutine gr_p4_pos_kick(param, x, v, dt)
       !! author: David A. Minton

src/helio/helio_drift.f90

@@ -71,7 +71,7 @@ module subroutine helio_drift_tp(self, system, param, dt, mask)
       return
    end subroutine helio_drift_tp
 
-   module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
+   module subroutine helio_drift_linear_pl(self, cb, dt, mask, lbeg)
       !! author: David A. Minton
       !!
       !! Perform linear drift of massive bodies due to barycentric momentum of Sun
@@ -80,21 +80,24 @@ module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
       !! Adapted from Hal Levison's Swift routine helio_lindrift.f
       implicit none
       ! Arguments
-      class(helio_pl), intent(inout) :: self   !! Helio massive body object
-      class(helio_cb), intent(inout) :: cb   !! Helio central bod
-      real(DP),        intent(in)    :: dt     !! Stepsize
-      logical,         intent(in)    :: lbeg   !! Argument that determines whether or not this is the beginning or end of the step
+      class(helio_pl),               intent(inout) :: self !! Helio massive body object
+      class(helio_cb),               intent(inout) :: cb   !! Helio central body
+      real(DP),                      intent(in)    :: dt   !! Stepsize
+      logical,         dimension(:), intent(in)    :: mask !! Mask that determines which bodies to kick
+      logical,                       intent(in)    :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
       ! Internals
-      real(DP), dimension(NDIM)      :: pt     !! negative barycentric velocity of the central body
+      real(DP), dimension(NDIM) :: pt     !! negative barycentric velocity of the central body
+      integer(I4B)              :: i    
 
       associate(pl => self, npl => self%nbody)
-         pt(1) = sum(pl%Gmass(1:npl) * pl%vb(1,1:npl))
-         pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl))
-         pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl))
+         if (npl == 0) return
+         pt(1) = sum(pl%Gmass(1:npl) * pl%vb(1,1:npl), mask)
+         pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl), mask)
+         pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl), mask)
          pt(:) = pt(:) / cb%Gmass
-         pl%xh(1,1:npl) = pl%xh(1,1:npl) + pt(1) * dt
-         pl%xh(2,1:npl) = pl%xh(2,1:npl) + pt(2) * dt
-         pl%xh(3,1:npl) = pl%xh(3,1:npl) + pt(3) * dt
+         do concurrent(i = 1:npl, mask(i))
+            pl%xh(:,i) = pl%xh(:,i) + pt(:) * dt
+         end do
 
          if (lbeg) then
             cb%ptbeg = pt(:)
@@ -106,7 +109,7 @@ module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
       return
    end subroutine helio_drift_linear_pl
 
-   module subroutine helio_drift_linear_tp(self, cb, dt, lbeg)
+   module subroutine helio_drift_linear_tp(self, cb, dt, mask, lbeg)
       !! author: David A. Minton
       !!
       !! Perform linear drift of test particles due to barycentric momentum of Sun
@@ -116,20 +119,22 @@ module subroutine helio_drift_linear_tp(self, cb, dt, lbeg)
       !! Adapted from Hal Levison's Swift routine helio_lindrift_tp.f
       implicit none
       ! Arguments
-      class(helio_tp), intent(inout) :: self !! Helio test particleb object
-      class(helio_cb), intent(in)    :: cb   !! Helio central body
-      real(DP),        intent(in)    :: dt   !! Stepsize
-      logical,         intent(in)    :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
+      class(helio_tp),               intent(inout) :: self !! Helio test particleb object
+      class(helio_cb),               intent(in)    :: cb   !! Helio central body
+      real(DP),                      intent(in)    :: dt   !! Stepsize
+      logical,         dimension(:), intent(in)    :: mask !! Mask that determines which bodies to kick
+      logical,                       intent(in)    :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
       ! Internals
       real(DP), dimension(NDIM)      :: pt     !! negative barycentric velocity of the central body
 
       associate(tp => self, ntp => self%nbody)
+         if (ntp == 0) return
          if (lbeg) then
             pt(:) = cb%ptbeg
          else
             pt(:) = cb%ptend
          end if
-         where (tp%status(1:ntp) == ACTIVE)
+         where (mask(1:ntp))
             tp%xh(1, 1:ntp) = tp%xh(1, 1:ntp) + pt(1) * dt
             tp%xh(2, 1:ntp) = tp%xh(2, 1:ntp) + pt(2) * dt
             tp%xh(3, 1:ntp) = tp%xh(3, 1:ntp) + pt(3) * dt

src/helio/helio_kick.f90

@@ -1,53 +1,141 @@
 submodule(helio_classes) s_helio_kick
    use swiftest
 contains
-   module subroutine helio_kickvb_pl(self, dt)
+module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg)
+   !! author: David A. Minton
+   !!
+   !! Compute heliocentric accelerations of massive bodies
+   !!
+   !! Adapted from David E. Kaufmann's Swifter routine helio_kick_getacch.f90
+   !! Adapted from Hal Levison's Swift routine helio_kick_getacch.f
+   implicit none
+   ! Arguments
+   class(helio_pl),              intent(inout) :: self   !! Helio massive body particle data structure
+   class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+   class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+   real(DP),                     intent(in)    :: t      !! Current simulation time
+   logical,                      intent(in)    :: lbeg   !! Logical flag that determines whether or not this is the beginning or end of the step
+
+   associate(cb => system%cb, pl => self, npl => self%nbody)
+      call pl%accel_int()
+      if (param%loblatecb) then 
+         call pl%accel_obl(system)
+         if (lbeg) then
+            cb%aoblbeg = cb%aobl
+         else
+            cb%aoblend = cb%aobl
+         end if
+         if (param%ltides) then
+            call pl%accel_tides(system)
+            if (lbeg) then
+               cb%atidebeg = cb%atide
+            else
+               cb%atideend = cb%atide
+            end if
+         end if
+      end if
+      if (param%lextra_force) call pl%accel_user(system, param, t, lbeg)
+      !if (param%lgr) call pl%gr_accel(param)
+   end associate
+
+   return
+   end subroutine helio_kick_getacch_pl
+
+   module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg)
+      !! author: David A. Minton
+      !!
+      !! Compute heliocentric accelerations of test particles
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine helio_kick_getacch_tp.f90
+      !! Adapted from Hal Levison's Swift routine helio_kick_getacch_tp.f
+      implicit none
+      ! Arguments
+      class(helio_tp),              intent(inout) :: self   !! Helio test particle data structure
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+      real(DP),                     intent(in)    :: t      !! Current time
+      logical,                      intent(in)    :: lbeg   !! Logical flag that determines whether or not this is the beginning or end of the step
+
+      associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody)
+         system%lbeg = lbeg
+         if (system%lbeg) then
+            call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl)
+         else
+            call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl)
+         end if
+         if (param%loblatecb) call tp%accel_obl(system)
+         if (param%lextra_force) call tp%accel_user(system, param, t, lbeg)
+         !if (param%lgr) call tp%gr_accel(param)
+      end associate
+      return
+   end subroutine helio_kick_getacch_tp
+
+   module subroutine helio_kick_vb_pl(self, system, param, t, dt, mask, lbeg)
       !! author: David A. Minton
       !!
       !! Kick barycentric velocities of bodies
       !!
       !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f 
-      !! Adapted from David E. Kaufmann's Swifter routine helio_kickvb.f90
+      !! Adapted from David E. Kaufmann's Swifter routine helio_kick_vb.f90
       implicit none
       ! Arguments
-      class(helio_pl), intent(inout) :: self !! Swiftest generic body object
-      real(DP),        intent(in)    :: dt   !! Stepsize
+      class(helio_pl),              intent(inout) :: self   !! Swiftest generic body object
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+      real(DP),                     intent(in)    :: t      !! Current time
+      real(DP),                     intent(in)    :: dt     !! Stepsize
+      logical, dimension(:),        intent(in)    :: mask   !! Mask that determines which bodies to kick
+      logical,                      intent(in)    :: lbeg   !! Logical flag indicating whether this is the beginning of the half step or not. 
       ! Internals
       integer(I4B) :: i
 
       associate(pl => self, npl => self%nbody)
          if (npl ==0) return
-         do concurrent(i = 1:npl, pl%status(i) == ACTIVE) 
+         pl%ah(:,:) = 0.0_DP
+         call pl%accel(system, param, t, lbeg)
+         if (lbeg) then
+            call pl%set_beg_end(xbeg = pl%xh)
+         else
+            call pl%set_beg_end(xend = pl%xh)
+         end if
+         do concurrent(i = 1:npl, mask(i)) 
             pl%vb(:, i) = pl%vb(:, i) + pl%ah(:, i) * dt
          end do
       end associate
 
       return
 
-   end subroutine helio_kickvb_pl
+   end subroutine helio_kick_vb_pl
 
-   module subroutine helio_kickvb_tp(self, dt)
+   module subroutine helio_kick_vb_tp(self, system, param, t, dt, mask, lbeg)
       !! author: David A. Minton
       !!
       !! Kick barycentric velocities of bodies
       !!
       !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh_tp.f
-      !! Adapted from David E. Kaufmann's Swifter routine helio_kickvb_tp.f90
+      !! Adapted from David E. Kaufmann's Swifter routine helio_kick_vb_tp.f90
       implicit none
       ! Arguments
-      class(helio_tp), intent(inout) :: self !! Swiftest generic body object
-      real(DP),        intent(in)    :: dt   !! Stepsize
+      class(helio_tp),              intent(inout) :: self !! Swiftest generic body object
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param !! Current run configuration parameters 
+      real(DP),                     intent(in)    :: t     !! Current time
+      real(DP),                     intent(in)    :: dt    !! Stepsize
+      logical, dimension(:),        intent(in)    :: mask  !! Mask that determines which bodies to kick
+      logical,                      intent(in)    :: lbeg  !! Logical flag indicating whether this is the beginning of the half step or not. 
       ! Internals
       integer(I4B) :: i
 
       associate(tp => self, ntp => self%nbody)
          if (ntp ==0) return
-         do concurrent(i = 1:ntp, tp%status(i) == ACTIVE) 
+         tp%ah(:,:) = 0.0_DP
+         call tp%accel(system, param, t, lbeg)
+         do concurrent(i = 1:ntp, mask(i)) 
             tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt
          end do
       end associate
 
       return
 
-   end subroutine helio_kickvb_tp
+   end subroutine helio_kick_vb_tp
 end submodule s_helio_kick