Restructured GR routines to prepare for testing in Helio/SyMBA. Built…

… additional templates for SyMBA.

src/gr/gr.f90

@@ -1,50 +1,180 @@
 submodule(swiftest_classes) s_gr
    use swiftest
 contains
-   subroutine gr_getaccb_ns_body(self, cb, param, agr, agr0) 
-
+   module pure subroutine gr_getaccb_ns_body(self, system, param) 
       !! author: David A. Minton
       !!
-      !! Add relativistic correction acceleration for non-symplectic integrators
-      !!    Based on Quinn, Tremaine, & Duncan 1998
+      !! Add relativistic correction acceleration for non-symplectic integrators.
+      !! Based on Quinn et al. (1991) eq. 5
+      !!
+      !! Reference:  
+      !!    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
       implicit none
       ! Arguments
-      class(swiftest_body),       intent(inout) :: self
-      class(swiftest_cb),         intent(inout) :: cb
-      class(swiftest_parameters), intent(in)    :: param
-      real(DP), dimension(:, :),  intent(inout) :: agr
-      real(DP), dimension(NDIM),  intent(out)   :: agr0
+      class(swiftest_body),         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 
       ! Internals
-      real(DP), dimension(NDIM) :: xh, vh
       real(DP)                  :: rmag, rdotv, vmag2
       integer(I4B)              :: i
 
-      associate(n => self%nbody, msun => cb%Gmass, vbsun => cb%vb, xbsun => cb%xb, mu => self%mu, c2 => param%inv_c2, &
-                xb => self%xb, vb => self%vb)
+      associate(n => self%nbody, cb => system%cb,  inv_c2 => param%inv_c2)
          if (n == 0) return
          do i = 1, n
-            xh(:) = xb(:, i) - xbsun(:)
-            vh(:) = vb(:, i) - vbsun(:)
-            rmag = norm2(xh(:))
-            vmag2 = dot_product(vh(:), vh(:))
-            rdotv = dot_product(xh(:), vh(:))
-            agr(:, i) =  mu * c2 / rmag**3 * ((4 * mu(i) / rmag - vmag2) * xh(:) + 4 * rdotv * vh(:))
+            rmag = norm2(self%xh(:,i))
+            vmag2 = dot_product(self%vh(:,i), self%vh(:,i))
+            rdotv = dot_product(self%xh(:,i), self%vh(:,i))
+            self%agr(:, i) = self%mu * inv_c2 / rmag**3 * ((4 * self%mu(i) / rmag - vmag2) * self%xh(:,i) + 4 * rdotv * self%vh(:,i))
          end do
 
-         agr0 =  0.0_DP
          select type(self)
          class is (swiftest_pl)
             do i = 1, NDIM
-               agr0(i) = -sum(self%Gmass(1:n) * agr(1:n, i) / msun)
+               cb%agr(i) = -sum(self%Gmass(1:n) * self%agr(1:n, i) / cb%Gmass)
             end do
          end select
 
       end associate 
 
       return
 
-      end subroutine gr_getaccb_ns_body
+   end subroutine gr_getaccb_ns_body
+
+   module pure subroutine gr_p4_pos_kick(param, x, v, dt)
+      !! author: David A. Minton
+      !!
+      !! Position kick due to p**4 term in the post-Newtonian correction
+      !!    Based on Saha & Tremaine (1994) Eq. 28
+      !!
+      !! Reference:
+      !!    Saha, P., Tremaine, S., 1994. Long-term planetary integration with individual time steps. 
+      !!       AJ 108, 1962–1969. https://doi.org/10.1086/117210
+      !!
+      !! Adapted from David A. Minton's Swifter routine gr_whm_p4.f90  
+      implicit none
+      ! Arguments
+      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
+      real(DP), dimension(:),     intent(inout) :: x     !! Position vector
+      real(DP), dimension(:),     intent(in)    :: v     !! Velocity vector
+      real(DP),                   intent(in)    :: dt    !! Step size
+      ! Internals
+      real(DP), dimension(NDIM)             :: dr
+      real(DP)                              :: vmag2
+
+      vmag2 = dot_product(v(:), v(:)) 
+      dr(:) = -2 * param%inv_c2 * vmag2 * v(:)
+      x(:) = x(:) + dr(:) * dt
+
+      return
+   end subroutine gr_p4_pos_kick
+
+   module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh) 
+      !! author: David A. Minton
+      !!
+      !! Converts the relativistic pseudovelocity back into a veliocentric velocity
+      !!    Based on Saha & Tremaine (1994) Eq. 32
+      !! 
+      !! Reference:
+      !!    Saha, P., Tremaine, S., 1994. Long-term planetary integration with individual time steps. 
+      !!       AJ 108, 1962–1969. https://doi.org/10.1086/117210
+      !!
+      !! Adapted from David A. Minton's Swifter routine gr_pseudovel2vel.f90 
+      implicit none
+      ! Arguments
+      class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
+      real(DP),                   intent(in)  :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
+      real(DP), dimension(:),     intent(in)  :: xh    !! Heliocentric position vector 
+      real(DP), dimension(:),     intent(in)  :: pv    !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
+      real(DP), dimension(:),     intent(out) :: vh    !! Heliocentric velocity vector 
+      ! Internals
+      real(DP) :: vmag2, rmag, grterm
+
+      associate(c2 => param%inv_c2)
+         vmag2 = dot_product(pv(:), pv(:))
+         rmag  = norm2(xh(:))
+         grterm = 1.0_DP - c2 * (0.5_DP * vmag2 + 3 * mu / rmag)
+         vh(:) = pv(:) * grterm
+      end associate
+      return
+   end subroutine gr_pseudovel2vel
+
+   module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
+      !! author: David A. Minton
+      !!
+      !! Converts the heliocentric velocity into a pseudovelocity with relativistic corrections. 
+      !! Uses Newton-Raphson method with direct inversion of the Jacobian (yeah, it's slow, but 
+      !! this is only done once per run).
+      !!
+      !! Reference:
+      !!    Saha, P., Tremaine, S., 1994. Long-term planetary integration with individual time steps. 
+      !!       AJ 108, 1962–1969. https://doi.org/10.1086/117210
+      !!
+      !! Adapted from David A. Minton's Swifter routine gr_vel2pseudovel.f90
+      implicit none
+      ! Arguments
+      class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
+      real(DP),                   intent(in)  :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
+      real(DP), dimension(:),     intent(in)  :: xh    !! Heliocentric position vector 
+      real(DP), dimension(:),     intent(in)  :: vh    !! Heliocentric velocity vector 
+      real(DP), dimension(:),     intent(out) :: pv    !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
+      ! Internals
+      real(DP)                       :: v2, G, pv2, rterm, det
+      real(DP), dimension(NDIM,NDIM) :: J,Jinv
+      real(DP), dimension(NDIM)      :: F
+      integer(I4B)                   :: n,i,k
+      integer(I4B), parameter        :: MAXITER = 50
+      real(DP),parameter             :: TOL = 1.0e-12_DP
+
+      associate (c2 => param%inv_c2)
+         pv(1:NDIM) = vh(1:NDIM) ! Initial guess
+         rterm = 3 * mu / norm2(xh(:))
+         v2 = dot_product(vh(:), vh(:))
+         do n = 1, MAXITER
+            pv2 = dot_product(pv(:), pv(:))
+            G = 1.0_DP - c2 * (0.5_DP * pv2 + rterm)
+            F(:) = pv(:) * G - vh(:)
+            if (abs(sum(F) / v2 ) < TOL) exit ! Root found
+
+            ! Calculate the Jacobian
+            do k = 1, NDIM
+               do i = 1, NDIM
+                  if (i == k) then
+                     J(i,k) = G - c2 * pv(k)
+                  else
+                     J(i,k) = - c2 * pv(k)
+                  end if
+               end do
+            end do
+
+            ! Inverse of the Jacobian
+            det = J(1,1) * (J(3,3) * J(2,2) - J(3,2) * J(2,3))
+            det = det - J(2,1) * (J(3,3) * J(1,2)-J(3,2) * J(1,3))
+            det = det + J(3,1) * (J(2,3) * J(1,2)-J(2,2) * J(1,3))
+
+            Jinv(1,1) =   J(3,3) * J(2,2) - J(3,2) * J(2,3)
+            Jinv(1,2) = -(J(3,3) * J(1,2) - J(3,2) * J(1,3))
+            Jinv(1,3) =   J(2,3) * J(1,2) - J(2,2) * J(1,3)
+
+            Jinv(2,1) = -(J(3,3) * J(2,1) - J(3,1) * J(2,3))
+            Jinv(2,2) =   J(3,3) * J(1,1) - J(3,1) * J(1,3)
+            Jinv(2,3) = -(J(2,3) * J(1,1) - J(2,1) * J(1,3))
+
+            Jinv(3,1) =   J(3,2) * J(2,1) - J(3,1) * J(2,2)
+            Jinv(3,2) = -(J(3,2) * J(1,1) - J(3,1) * J(1,2))
+            Jinv(3,3) =   J(2,2) * J(1,1) - J(2,1) * J(1,2)
+
+            Jinv = Jinv * det
+
+            do i = 1, NDIM
+               pv(i) = pv(i) - dot_product(Jinv(i,:) ,F(:))
+            end do
+         end do 
+
+      end associate
+      return
+   end subroutine gr_vel2pseudovel
 
 end submodule s_gr

src/modules/swiftest.f90

@@ -14,9 +14,6 @@ module swiftest
    use module_nrutil
    !use advisor_annotate
    !$ use omp_lib
-
-
-
    implicit none
    public
 

src/modules/swiftest_classes.f90

@@ -9,6 +9,7 @@ module swiftest_classes
    public :: discard_pl, discard_system, discard_tp 
    public :: drift_one
    public :: eucl_dist_index_plpl, eucl_dist_index_pltp, eucl_irij3_plpl
+   public :: gr_getaccb_ns_body, gr_p4_pos_kick, gr_pseudovel2vel, gr_vel2pseudovel
    public :: io_dump_param, io_dump_swiftest, io_dump_system, io_get_args, io_get_token, io_param_reader, io_param_writer, io_read_body_in, &
              io_read_cb_in, io_read_param_in, io_read_frame_body, io_read_frame_cb, io_read_frame_system, io_read_initialize_system, &
              io_write_discard, io_write_encounter, io_write_frame_body, io_write_frame_cb, io_write_frame_system
@@ -111,6 +112,7 @@ module swiftest_classes
       real(DP), dimension(NDIM) :: aoblend = 0.0_DP !! Barycentric acceleration due to central body oblatenes at end of step
       real(DP), dimension(NDIM) :: xb      = 0.0_DP !! Barycentric position (units DU)
       real(DP), dimension(NDIM) :: vb      = 0.0_DP !! Barycentric velocity (units DU / TU)
+      real(DP), dimension(NDIM) :: agr     = 0.0_DP !! Acceleration due to post-Newtonian correction
    contains
       private
       procedure, public         :: initialize  => io_read_cb_in        !! I/O routine for reading in central body data
@@ -137,6 +139,7 @@ module swiftest_classes
       real(DP),              dimension(:,:), allocatable :: vb              !! Barycentric velocity
       real(DP),              dimension(:,:), allocatable :: ah              !! Total heliocentric acceleration
       real(DP),              dimension(:,:), allocatable :: aobl            !! Barycentric accelerations of bodies due to central body oblatenes
+      real(DP),              dimension(:,:), allocatable :: agr             !! Acceleration due to post-Newtonian correction
       real(DP),              dimension(:),   allocatable :: ir3h            !! Inverse heliocentric radius term (1/rh**3)
       real(DP),              dimension(:),   allocatable :: a               !! Semimajor axis (pericentric distance for a parabolic orbit)
       real(DP),              dimension(:),   allocatable :: e               !! Eccentricity
@@ -380,6 +383,39 @@ module subroutine eucl_irij3_plpl(self)
          class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
       end subroutine eucl_irij3_plpl
 
+      module pure subroutine gr_getaccb_ns_body(self, system, param)
+         implicit none
+         class(swiftest_body),         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 
+      end subroutine gr_getaccb_ns_body
+
+      module pure subroutine gr_p4_pos_kick(param, x, v, dt)
+         implicit none
+         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
+         real(DP), dimension(:),     intent(inout) :: x     !! Position vector
+         real(DP), dimension(:),     intent(in)    :: v     !! Velocity vector
+         real(DP),                   intent(in)    :: dt    !! Step size
+      end subroutine gr_p4_pos_kick
+
+      module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh) 
+         implicit none
+         class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
+         real(DP),                   intent(in)  :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
+         real(DP), dimension(:),     intent(in)  :: xh    !! Heliocentric position vector 
+         real(DP), dimension(:),     intent(in)  :: pv    !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
+         real(DP), dimension(:),     intent(out) :: vh    !! Heliocentric velocity vector 
+      end subroutine gr_pseudovel2vel
+
+      module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
+         implicit none
+         class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
+         real(DP),                   intent(in)  :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
+         real(DP), dimension(:),     intent(in)  :: xh    !! Heliocentric position vector 
+         real(DP), dimension(:),     intent(in)  :: vh    !! Heliocentric velocity vector 
+         real(DP), dimension(:),     intent(out) :: pv    !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
+      end subroutine gr_vel2pseudovel
+
       module subroutine io_dump_param(self, param_file_name)
          implicit none
          class(swiftest_parameters),intent(in)    :: self            !! Output collection of parameters

src/modules/symba_classes.f90

@@ -271,7 +271,6 @@ module subroutine symba_io_read_frame_pl(self, iu, param, form, ierr)
          integer(I4B),               intent(out)   :: ierr    !! Error code
       end subroutine symba_io_read_frame_pl
 
-
       module subroutine symba_io_read_pl_in(self, param) 
          use swiftest_classes, only : swiftest_parameters
          implicit none

src/rmvs/rmvs_setup.f90

@@ -24,7 +24,6 @@ module subroutine rmvs_setup_pl(self,n)
          if (.not.pl%lplanetocentric) then
             allocate(pl%nenc(n))
             pl%nenc(:) = 0
-
             ! Set up inner and outer planet interpolation vector storage containers
             do i = 0, NTENC
                allocate(pl%outer(i)%x(NDIM, n))

src/setup/setup.f90

@@ -46,7 +46,19 @@ module subroutine setup_construct_system(system, param)
             allocate(rmvs_tp :: system%tp_discards)
          end select
       case (SYMBA)
-         write(*,*) 'SyMBA integrator not yet enabled'
+         allocate(symba_nbody_system :: system)
+         select type(system)
+         class is (symba_nbody_system)
+            allocate(symba_cb :: system%cb)
+            allocate(symba_pl :: system%pl)
+            allocate(symba_tp :: system%tp)
+            allocate(symba_pl :: system%pl_discards)
+            allocate(symba_tp :: system%tp_discards)
+            allocate(symba_pl :: system%mergeadd_list)
+            allocate(symba_pl :: system%mergesub_list)
+            allocate(symba_plplenc :: system%plplenc_list)
+            allocate(symba_pltpenc :: system%pltpenc_list)
+         end select
       case (RINGMOONS)
          write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled'
       case default
@@ -73,6 +85,7 @@ module subroutine setup_body(self,n)
 
       !write(*,*) 'Allocating the basic Swiftest particle'
       allocate(self%id(n))
+      allocate(self%name(n))
       allocate(self%status(n))
       allocate(self%ldiscard(n))
       allocate(self%xh(NDIM, n))
@@ -81,6 +94,7 @@ module subroutine setup_body(self,n)
       allocate(self%vb(NDIM, n))
       allocate(self%ah(NDIM, n))
       allocate(self%aobl(NDIM, n))
+      allocate(self%agr(NDIM, n))
       allocate(self%ir3h(n))
       allocate(self%a(n))
       allocate(self%e(n))
@@ -91,6 +105,7 @@ module subroutine setup_body(self,n)
       allocate(self%mu(n))
 
       self%id(:)   = 0
+      self%name(:) = "UNNAMED"
       self%status(:) = INACTIVE
       self%ldiscard(:) = .false.
       self%xh(:,:)   = 0.0_DP

src/symba/symba_discard.f90

@@ -0,0 +1,25 @@
+submodule (symba_classes) s_symba_discard
+   use swiftest
+contains
+
+   module subroutine symba_discard_pl(self, system, param)
+      implicit none
+      ! Arguments
+      class(symba_pl),              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 
+
+      return
+   end subroutine symba_discard_pl
+
+   module subroutine symba_discard_tp(self, system, param)
+      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 
+
+      return
+   end subroutine symba_discard_tp
+
+end submodule s_symba_discard

src/symba/symba_encounter_check.f90

@@ -0,0 +1,30 @@
+submodule (symba_classes) s_symba_encounter_check
+   use swiftest
+contains
+   module function symba_encounter_check_pl(self, system, dt) result(lencounter)
+      implicit none
+      ! Arguments
+      class(symba_pl),           intent(inout) :: self       !! SyMBA test particle object  
+      class(symba_nbody_system), intent(inout) :: system     !! SyMBA nbody system object
+      real(DP),                  intent(in)    :: dt         !! step size
+      ! Result
+      logical                                  :: lencounter !! Returns true if there is at least one close encounter      
+
+      lencounter = .false.
+      return
+   end function symba_encounter_check_pl
+
+   module function symba_encounter_check_tp(self, system, dt) result(lencounter)
+      implicit none
+      ! Arguments
+      class(symba_tp),           intent(inout) :: self       !! SyMBA test particle object  
+      class(symba_nbody_system), intent(inout) :: system     !! SyMBA nbody system object
+      real(DP),                  intent(in)    :: dt         !! step size
+      ! Result
+      logical                                  :: lencounter !! Returns true if there is at least one close encounter      
+
+      lencounter = .false.
+      return
+   end function symba_encounter_check_tp
+
+end submodule s_symba_encounter_check