Added in override to gr_pos_kick for SyMBA, which allows the position…

… kick operation to be masked during recursive steps so only bodies at the current recursion depth are kicked

src/helio/helio_gr.f90

@@ -56,7 +56,7 @@ module pure subroutine helio_gr_kick_getacch_tp(self, param)
    end subroutine helio_gr_kick_getacch_tp
 
 
-   module pure subroutine helio_gr_p4_pl(self, param, dt)
+   module pure subroutine helio_gr_p4_pl(self, system, param, dt)
       !! author: David A. Minton
       !!
       !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction
@@ -65,11 +65,12 @@ module pure subroutine helio_gr_p4_pl(self, param, dt)
       !! Adapted from David A. Minton's Swifter routine routine gr_helio_p4.f90
       implicit none
       ! Arguments
-      class(helio_pl),            intent(inout) :: self   !! Swiftest particle object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
-      real(DP),                   intent(in)    :: dt     !! Step size
+      class(helio_pl),              intent(inout) :: self   !! Swiftest particle 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)    :: dt     !! Step size
       ! Internals
-      integer(I4B)                                 :: i
+      integer(I4B) :: i
 
       if (self%nbody == 0) return
 
@@ -82,7 +83,8 @@ module pure subroutine helio_gr_p4_pl(self, param, dt)
      return
    end subroutine helio_gr_p4_pl
 
-   module pure subroutine helio_gr_p4_tp(self, param, dt)
+
+   module pure subroutine helio_gr_p4_tp(self, system, param, dt)
       !! author: David A. Minton
       !!
       !! Position kick to test particles due to p**4 term in the post-Newtonian correction
@@ -91,11 +93,12 @@ module pure subroutine helio_gr_p4_tp(self, param, dt)
       !! Adapted from David A. Minton's Swifter routine routine gr_helio_p4.f90
       implicit none
       ! Arguments
-      class(helio_tp),              intent(inout) :: self  !! Swiftest particle object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
-      real(DP),                   intent(in)    :: dt    !! Step size
+      class(helio_tp),              intent(inout) :: self   !! Swiftest particle 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)    :: dt     !! Step size
       ! Internals
-      integer(I4B)                              :: i
+      integer(I4B) :: i
 
       if (self%nbody == 0) return
 

src/helio/helio_step.f90

@@ -35,7 +35,7 @@ module subroutine helio_step_pl(self, system, param, t, dt)
       implicit none
       ! Arguments
       class(helio_pl),              intent(inout) :: self   !! Helio massive body particle data structure
-      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system
       class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
       real(DP),                     intent(in)    :: t      !! Current simulation time
       real(DP),                     intent(in)    :: dt     !! Stepsize
@@ -54,10 +54,10 @@ module subroutine helio_step_pl(self, system, param, t, dt)
             end if
             call pl%lindrift(cb, dth, lbeg=.true.)
             call pl%kick(system, param, t, dth, lbeg=.true.)
-            if (param%lgr) call pl%gr_pos_kick(param, dth)
+            if (param%lgr) call pl%gr_pos_kick(system, param, dth)
             call pl%drift(system, param, dt)
             call pl%kick(system, param, t + dt, dth, lbeg=.false.)
-            if (param%lgr) call pl%gr_pos_kick(param, dth)
+            if (param%lgr) call pl%gr_pos_kick(system, param, dth)
             call pl%lindrift(cb, dth, lbeg=.false.)
             call pl%vb2vh(cb)
          end select
@@ -78,7 +78,7 @@ module subroutine helio_step_tp(self, system, param, t, dt)
       implicit none
       ! Arguments
       class(helio_tp),              intent(inout) :: self    !! Helio test particle data structure
-      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system
       class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
       real(DP),                     intent(in)    :: t      !! Current simulation time
       real(DP),                     intent(in)    :: dt     !! Stepsize
@@ -97,10 +97,10 @@ module subroutine helio_step_tp(self, system, param, t, dt)
             end if
             call tp%lindrift(cb, dth, lbeg=.true.)
             call tp%kick(system, param, t, dth, lbeg=.true.)
-            if (param%lgr) call tp%gr_pos_kick(param, dth)
+            if (param%lgr) call tp%gr_pos_kick(system, param, dth)
             call tp%drift(system, param, dt)
             call tp%kick(system, param, t + dt, dth, lbeg=.false.)
-            if (param%lgr) call tp%gr_pos_kick(param, dth)
+            if (param%lgr) call tp%gr_pos_kick(system, param, dth)
             call tp%lindrift(cb, dth, lbeg=.false.)
             call tp%vb2vh(vbcb = -cb%ptend)
          end select

src/modules/fraggle_classes.f90

@@ -150,7 +150,7 @@ module subroutine fraggle_placeholder_step(self, system, param, t, dt)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
          implicit none
          class(fraggle_fragments),     intent(inout) :: self   !! Helio massive body particle object
-         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system
          class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
          real(DP),                     intent(in)    :: t      !! Current simulation time
          real(DP),                     intent(in)    :: dt     !! Stepsiz

src/modules/helio_classes.f90

@@ -122,20 +122,22 @@ module pure subroutine helio_gr_kick_getacch_tp(self, param)
          class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters 
       end subroutine helio_gr_kick_getacch_tp
 
-      module pure subroutine helio_gr_p4_pl(self, param, dt)
-         use swiftest_classes, only : swiftest_parameters
+      module pure subroutine helio_gr_p4_pl(self, system, param, dt)
+         use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
          implicit none
-         class(helio_pl),            intent(inout) :: self   !! Swiftest particle object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
-         real(DP),                   intent(in)    :: dt     !! Step size
+         class(helio_pl),              intent(inout) :: self   !! Swiftest particle 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)    :: dt     !! Step size
       end subroutine helio_gr_p4_pl
 
-      module pure subroutine helio_gr_p4_tp(self, param, dt)
-         use swiftest_classes, only : swiftest_parameters
+      module pure subroutine helio_gr_p4_tp(self, system, param, dt)
+         use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
          implicit none
-         class(helio_tp),            intent(inout) :: self  !! Swiftest particle object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
-         real(DP),                   intent(in)    :: dt    !! Step size
+         class(helio_tp),              intent(inout) :: self   !! Swiftest particle 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)    :: dt     !! Step size
       end subroutine helio_gr_p4_tp
 
       module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg)
@@ -191,7 +193,7 @@ module subroutine helio_step_pl(self, system, param, t, dt)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
          implicit none
          class(helio_pl),              intent(inout) :: self   !! Helio massive body particle object
-         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system
          class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
          real(DP),                     intent(in)    :: t      !! Current simulation time
          real(DP),                     intent(in)    :: dt     !! Stepsize

src/modules/symba_classes.f90

@@ -75,6 +75,7 @@ module symba_classes
       procedure :: discard         => symba_discard_pl                  !! Process massive body discards
       procedure :: drift           => symba_drift_pl                    !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level
       procedure :: encounter_check => symba_encounter_check_pl          !! Checks if massive bodies are going through close encounters with each other
+      procedure :: gr_pos_kick     => symba_gr_p4_pl                    !! Position kick due to p**4 term in the post-Newtonian correction
       procedure :: accel_int       => symba_kick_getacch_int_pl         !! Compute direct cross (third) term heliocentric accelerations of massive bodiess, with no mutual interactions between bodies below GMTINY
       procedure :: accel           => symba_kick_getacch_pl             !! Compute heliocentric accelerations of massive bodies
       procedure :: setup           => symba_setup_pl                    !! Constructor method - Allocates space for the input number of bodies
@@ -113,6 +114,7 @@ module symba_classes
    contains
       procedure :: drift           => symba_drift_tp               !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level
       procedure :: encounter_check => symba_encounter_check_tp     !! Checks if any test particles are undergoing a close encounter with a massive body
+      procedure :: gr_pos_kick     => symba_gr_p4_tp               !! Position kick due to p**4 term in the post-Newtonian correction
       procedure :: accel           => symba_kick_getacch_tp        !! Compute heliocentric accelerations of test particles
       procedure :: setup           => symba_setup_tp               !! Constructor method - Allocates space for the input number of bodies
       procedure :: append          => symba_util_append_tp         !! Appends elements from one structure to another
@@ -272,6 +274,7 @@ module subroutine symba_drift_tp(self, system, param, dt)
       end subroutine symba_drift_tp
 
       module function symba_encounter_check_pl(self, param, system, dt, irec) result(lany_encounter)
+         use swiftest_classes, only : swiftest_nbody_system
          implicit none
          class(symba_pl),            intent(inout) :: self           !! SyMBA test particle object  
          class(swiftest_parameters), intent(inout) :: param          !! Current swiftest run configuration parameters
@@ -282,6 +285,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
       end function symba_encounter_check_pl
 
       module function symba_encounter_check(self, param, system, dt, irec) result(lany_encounter)
+         use swiftest_classes, only : swiftest_parameters
          implicit none
          class(symba_encounter),     intent(inout) :: self           !! SyMBA pl-pl encounter list object
          class(swiftest_parameters), intent(inout) :: param          !! Current swiftest run configuration parameters
@@ -292,6 +296,7 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany
       end function symba_encounter_check
 
       module function symba_encounter_check_tp(self, param, system, dt, irec) result(lany_encounter)
+         use swiftest_classes, only : swiftest_parameters
          implicit none
          class(symba_tp),            intent(inout) :: self           !! SyMBA test particle object  
          class(swiftest_parameters), intent(inout) :: param          !! Current swiftest run configuration parameters
@@ -301,6 +306,24 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l
          logical                                   :: lany_encounter !! Returns true if there is at least one close encounter      
       end function symba_encounter_check_tp
 
+      module pure subroutine symba_gr_p4_pl(self, system, param, dt)
+         use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
+         implicit none
+         class(symba_pl),              intent(inout) :: self   !! SyMBA massive 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)    :: dt     !! Step size
+      end subroutine symba_gr_p4_pl
+
+      module pure subroutine symba_gr_p4_tp(self, system, param, dt)
+         use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
+         implicit none
+         class(symba_tp),              intent(inout) :: self   !! SyMBA test particle 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)    :: dt     !! Step size
+      end subroutine symba_gr_p4_tp
+
       module function symba_collision_casedisruption(system, param, colliders, frag) result(status)
          use fraggle_classes, only : fraggle_colliders, fraggle_fragments
          implicit none

src/modules/whm_classes.f90

@@ -170,20 +170,22 @@ module pure subroutine whm_gr_kick_getacch_tp(self, param)
          class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
       end subroutine whm_gr_kick_getacch_tp
 
-      module pure subroutine whm_gr_p4_pl(self, param, dt)
+      module pure subroutine whm_gr_p4_pl(self, system, param, dt)
          use swiftest_classes, only : swiftest_parameters
          implicit none
-         class(whm_pl),              intent(inout) :: self  !! WHM massive body object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
-         real(DP),                   intent(in)    :: dt    !! Step size
+         class(whm_pl),                intent(inout) :: self  !! WHM massive body object
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system
+         class(swiftest_parameters),   intent(in)    :: param !! Current run configuration parameters 
+         real(DP),                     intent(in)    :: dt    !! Step size
       end subroutine whm_gr_p4_pl
 
-      module pure subroutine whm_gr_p4_tp(self, param, dt)
+      module pure subroutine whm_gr_p4_tp(self, system, param, dt)
          use swiftest_classes, only : swiftest_parameters
          implicit none
-         class(whm_tp),              intent(inout) :: self  !! WHM test particle object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
-         real(DP),                    intent(in)   :: dt    !! Step size
+         class(whm_tp),                intent(inout) :: self   !! WHM test particle object
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system
+         class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+         real(DP),                     intent(in)    :: dt     !! Step size
       end subroutine whm_gr_p4_tp
 
       !> Reads WHM massive body object in from file

src/symba/symba_gr.f90

@@ -0,0 +1,66 @@
+submodule(symba_classes) s_symba_gr
+   use swiftest
+contains
+
+   module pure subroutine symba_gr_p4_pl(self, system, param, dt)
+      !! author: David A. Minton
+      !!
+      !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction
+      !!    Based on Saha & Tremaine (1994) Eq. 28
+      !!
+      !! Adapted from David A. Minton's Swifter routine routine gr_symba_p4.f90
+      implicit none
+      ! Arguments
+      class(symba_pl),              intent(inout) :: self   !! SyMBA massive 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)    :: dt     !! Step size
+      ! Internals
+      integer(I4B)                                 :: i
+
+      if (self%nbody == 0) return
+
+      associate(pl => self, npl => self%nbody)
+         select type(system)
+         class is (symba_nbody_system)
+            do concurrent(i = 1:npl, pl%lmask(i) .and. pl%levelg(i) == system%irec )
+               call gr_p4_pos_kick(param, pl%xh(:, i), pl%vb(:, i), dt)
+            end do
+         end select
+      end associate
+
+   return
+   end subroutine symba_gr_p4_pl
+
+
+   module pure subroutine symba_gr_p4_tp(self, system, param, dt)
+      !! author: David A. Minton
+      !!
+      !! Position kick to test particles due to p**4 term in the post-Newtonian correction
+      !!    Based on Saha & Tremaine (1994) Eq. 28
+      !!
+      !! Adapted from David A. Minton's Swifter routine routine gr_symba_p4.f90
+      implicit none
+      ! Arguments
+      class(symba_tp),              intent(inout) :: self   !! SyMBA test particle 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)    :: dt     !! Step size
+      ! Internals
+      integer(I4B)                              :: i
+
+      if (self%nbody == 0) return
+
+      associate(tp => self, ntp => self%nbody)
+         select type(system)
+         class is (symba_nbody_system)
+            do concurrent(i = 1:ntp, tp%lmask(i) .and. tp%levelg(i) == system%irec)
+               call gr_p4_pos_kick(param, tp%xh(:, i), tp%vb(:, i), dt)
+            end do
+         end select
+      end associate
+
+   return
+   end subroutine symba_gr_p4_tp
+
+end submodule s_symba_gr