Restructured gr pseudovelocity call to be part of the base swiftest_b…

…ody class

examples/whm_gr_test/init_cond.py

@@ -8,7 +8,7 @@
 sim.param['DU2M'] = swiftest.AU2M
 sim.param['T0'] = 0.0
 sim.param['DT'] = 0.25 * swiftest.JD2S / swiftest.YR2S
-sim.param['TSTOP'] = 1000.0
+sim.param['TSTOP'] = 100.0
 sim.param['ISTEP_OUT']  = 1461
 sim.param['ISTEP_DUMP'] = 1461
 sim.param['CHK_QMIN_COORD'] = "HELIO"

examples/whm_gr_test/param.swifter.in

@@ -1,6 +1,6 @@
 ! VERSION        Swifter parameter file converted from Swiftest
 T0               0.0
-TSTOP            1000.0
+TSTOP            100.0
 DT               0.0006844626967830253
 ISTEP_OUT        1461
 ISTEP_DUMP       1461

examples/whm_gr_test/param.swiftest.in

@@ -1,6 +1,6 @@
 ! VERSION        Swiftest parameter input
 T0               0.0
-TSTOP            1000.0
+TSTOP            100.0
 DT               0.0006844626967830253
 ISTEP_OUT        1461
 ISTEP_DUMP       1461

src/gr/gr.f90

@@ -101,6 +101,29 @@ module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh)
       return
    end subroutine gr_pseudovel2vel
 
+   module pure subroutine gr_pv2vh_body(self, param)
+      !! author: David A. Minton
+      !!
+      !! Wrapper function that converts from pseudovelocity to heliocentric velocity for swiftest bodies
+      implicit none
+      ! Arguments
+      class(swiftest_body),       intent(inout) :: self  !! Swiftest particle object
+      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
+      ! Internals
+      integer(I4B)                              :: i
+      real(DP), dimension(:,:), allocatable     :: vh    !! Temporary holder of pseudovelocity for in-place conversion
+
+      associate(n => self%nbody)
+         if (n == 0) return
+         allocate(vh, mold = self%vh)
+         do i = 1, n
+            call gr_pseudovel2vel(param, self%mu(i), self%xh(:, i), self%vh(:, i), vh(:, i))
+         end do
+         call move_alloc(vh, self%vh)
+      end associate
+      return
+   end subroutine gr_pv2vh_body
+
    module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
       !! author: David A. Minton
       !!
@@ -177,4 +200,28 @@ module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
       return
    end subroutine gr_vel2pseudovel
 
+   module pure subroutine gr_vh2pv_body(self, param)
+      !! author: David A. Minton
+      !!
+      !! Wrapper function that converts from heliocentric velocity to pseudovelocity for Swiftest bodies
+      implicit none
+      ! Arguments
+      class(swiftest_body),       intent(inout) :: self  !! Swiftest particle object
+      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
+      ! Internals
+      integer(I4B)                                 :: i
+      real(DP), dimension(:,:), allocatable        :: pv !! Temporary holder of pseudovelocity for in-place conversion
+
+      associate(n => self%nbody)
+         if (n == 0) return
+         allocate(pv, mold = self%vh)
+         do i = 1, n
+            call gr_vel2pseudovel(param, self%mu(i), self%xh(:, i), self%vh(:, i), pv(:, i))
+         end do
+         call move_alloc(pv, self%vh)
+      end associate
+      return
+   end subroutine gr_vh2pv_body
+
+
 end submodule s_gr

src/modules/swiftest_classes.f90

@@ -159,6 +159,8 @@ module swiftest_classes
       procedure(abstract_step_body),    public, deferred :: step
       procedure(abstract_accel),        public, deferred :: accel
       ! These are concrete because the implementation is the same for all types of particles
+      procedure, public :: vh2pv          => gr_vh2pv_body       !! Converts from heliocentric velocity to psudeovelocity for GR calculations
+      procedure, public :: pv2vh          => gr_pv2vh_body       !! Converts from psudeovelocity to heliocentric velocity for GR calculations
       procedure, public :: initialize     => io_read_body_in     !! Read in body initial conditions from a file
       procedure, public :: read_frame     => io_read_frame_body  !! I/O routine for writing out a single frame of time-series data for the central body
       procedure, public :: write_frame    => io_write_frame_body !! I/O routine for writing out a single frame of time-series data for the central body
@@ -407,6 +409,12 @@ module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh)
          real(DP), dimension(:),     intent(out) :: vh    !! Heliocentric velocity vector 
       end subroutine gr_pseudovel2vel
 
+      module pure subroutine gr_pv2vh_body(self, param)
+         implicit none
+         class(swiftest_body),       intent(inout) :: self  !! Swiftest particle object
+         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
+      end subroutine gr_pv2vh_body
+
       module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
          implicit none
          class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
@@ -416,6 +424,12 @@ module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
          real(DP), dimension(:),     intent(out) :: pv    !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
       end subroutine gr_vel2pseudovel
 
+      module pure subroutine gr_vh2pv_body(self, param)
+         implicit none
+         class(swiftest_body),       intent(inout) :: self  !! Swiftest particle object
+         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters
+      end subroutine gr_vh2pv_body
+
       module subroutine io_dump_param(self, param_file_name)
          implicit none
          class(swiftest_parameters),intent(in)    :: self            !! Output collection of parameters

src/modules/whm_classes.f90

@@ -30,21 +30,19 @@ module whm_classes
       !! Note to developers: If you add componenets to this class, be sure to update methods and subroutines that traverse the
       !!    component list, such as whm_setup_pl and whm_spill_pl
    contains
-      procedure, public :: h2j      => whm_coord_h2j_pl        !! Convert position and velcoity vectors from heliocentric to Jacobi coordinates 
-      procedure, public :: j2h      => whm_coord_j2h_pl        !! Convert position and velcoity vectors from Jacobi to helliocentric coordinates 
-      procedure, public :: vh2vj    => whm_coord_vh2vj_pl      !! Convert velocity vectors from heliocentric to Jacobi coordinates 
-      procedure, public :: drift    => whm_drift_pl            !! Loop through massive bodies and call Danby drift routine
-      procedure, public :: fill     => whm_fill_pl             !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic)
-      procedure, public :: accel    => whm_getacch_pl          !! Compute heliocentric accelerations of massive bodies
-      procedure, public :: accel_gr => whm_gr_getacch_pl       !! Acceleration term arising from the post-Newtonian correction
-      procedure, public :: p4       => whm_gr_p4_pl            !! Position kick due to p**4 term in the post-Newtonian correction
-      procedure, public :: vh2pv    => whm_gr_vh2pv_pl         !! Converts from heliocentric velocity to psudeovelocity for GR calculations
-      procedure, public :: pv2vh    => whm_gr_pv2vh_pl         !! Converts from psudeovelocity to heliocentric velocity for GR calculations
-      procedure, public :: setup    => whm_setup_pl            !! Constructor method - Allocates space for number of particles
-      procedure, public :: set_mu   => whm_util_set_mu_eta_pl !! Sets the Jacobi mass value for all massive bodies.
-      procedure, public :: set_ir3  => whm_setup_set_ir3j     !! Sets both the heliocentric and jacobi inverse radius terms (1/rj**3 and 1/rh**3)
-      procedure, public :: step     => whm_step_pl             !! Steps the body forward one stepsize
-      procedure, public :: spill    => whm_spill_pl            !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
+      procedure, public :: h2j         => whm_coord_h2j_pl        !! Convert position and velcoity vectors from heliocentric to Jacobi coordinates 
+      procedure, public :: j2h         => whm_coord_j2h_pl        !! Convert position and velcoity vectors from Jacobi to helliocentric coordinates 
+      procedure, public :: vh2vj       => whm_coord_vh2vj_pl      !! Convert velocity vectors from heliocentric to Jacobi coordinates 
+      procedure, public :: drift       => whm_drift_pl            !! Loop through massive bodies and call Danby drift routine
+      procedure, public :: fill        => whm_fill_pl             !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic)
+      procedure, public :: accel       => whm_getacch_pl          !! Compute heliocentric accelerations of massive bodies
+      procedure, public :: accel_gr    => whm_gr_getacch_pl       !! Acceleration term arising from the post-Newtonian correction
+      procedure, public :: gr_pos_kick => whm_gr_p4_pl            !! Position kick due to p**4 term in the post-Newtonian correction
+      procedure, public :: setup       => whm_setup_pl            !! Constructor method - Allocates space for number of particles
+      procedure, public :: set_mu      => whm_util_set_mu_eta_pl !! Sets the Jacobi mass value for all massive bodies.
+      procedure, public :: set_ir3     => whm_setup_set_ir3j     !! Sets both the heliocentric and jacobi inverse radius terms (1/rj**3 and 1/rh**3)
+      procedure, public :: step        => whm_step_pl             !! Steps the body forward one stepsize
+      procedure, public :: spill       => whm_spill_pl            !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
    end type whm_pl
 
    !********************************************************************************************************************************
@@ -57,14 +55,12 @@ module whm_classes
       !!    component list, such as whm_setup_tp and whm_spill_tp
    contains
       private
-      procedure, public :: drift    => whm_drift_tp        !! Loop through test particles and call Danby drift routine
-      procedure, public :: accel    => whm_getacch_tp      !! Compute heliocentric accelerations of test particles
-      procedure, public :: accel_gr => whm_gr_getacch_tp   !! Acceleration term arising from the post-Newtonian correction
-      procedure, public :: p4       => whm_gr_p4_tp        !! Position kick due to p**4 term in the post-Newtonian correction
-      procedure, public :: vh2pv    => whm_gr_vh2pv_tp     !! Converts from heliocentric velocity to psudeovelocity for GR calculations
-      procedure, public :: pv2vh    => whm_gr_pv2vh_tp     !! Converts from psudeovelocity to heliocentric velocity for GR calculations
-      procedure, public :: setup    => whm_setup_tp        !! Allocates new components of the whm class and recursively calls parent allocations
-      procedure, public :: step     => whm_step_tp         !! Steps the particle forward one stepsize
+      procedure, public :: drift       => whm_drift_tp        !! Loop through test particles and call Danby drift routine
+      procedure, public :: accel       => whm_getacch_tp      !! Compute heliocentric accelerations of test particles
+      procedure, public :: accel_gr    => whm_gr_getacch_tp   !! Acceleration term arising from the post-Newtonian correction
+      procedure, public :: gr_pos_kick => whm_gr_p4_tp        !! Position kick due to p**4 term in the post-Newtonian correction
+      procedure, public :: setup       => whm_setup_tp        !! Allocates new components of the whm class and recursively calls parent allocations
+      procedure, public :: step        => whm_step_tp         !! Steps the particle forward one stepsize
    end type whm_tp
 
    !********************************************************************************************************************************
@@ -179,34 +175,6 @@ module pure subroutine whm_gr_p4_tp(self, param, dt)
          real(DP),                    intent(in)   :: dt    !! Step size
       end subroutine whm_gr_p4_tp
 
-      module pure subroutine whm_gr_pv2vh_pl(self, param)
-         use swiftest_classes, only : swiftest_parameters
-         implicit none
-         class(whm_pl),              intent(inout) :: self   !! Swiftest particle object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
-      end subroutine whm_gr_pv2vh_pl
-
-      module pure subroutine whm_gr_pv2vh_tp(self, param)
-         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 of on parameters 
-      end subroutine whm_gr_pv2vh_tp
-
-      module pure subroutine whm_gr_vh2pv_pl(self, param)
-         use swiftest_classes, only : swiftest_parameters
-         implicit none
-         class(whm_pl),              intent(inout) :: self   !! Swiftest particle object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
-      end subroutine whm_gr_vh2pv_pl
-
-      module pure subroutine whm_gr_vh2pv_tp(self, param)
-         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 of on parameters 
-      end subroutine whm_gr_vh2pv_tp
-
       !> Reads WHM massive body object in from file
       module subroutine whm_setup_pl(self,n)
          implicit none

src/whm/whm_gr.f90

@@ -17,20 +17,19 @@ module subroutine whm_gr_getacch_pl(self, param) !! author: David A. Minton
       real(DP), dimension(:, :), allocatable       :: aj
       real(DP)                                     :: beta, rjmag4
 
-      associate(n => self%nbody, mu => self%muj, c2 => param%inv_c2, &
-         ah => self%ah, xj => self%xj, GMpl => self%Gmass, eta => self%eta)
-         if (n == 0) return
-         allocate(aj, mold = ah)
-         do i = 1, n
-            rjmag4 = (dot_product(xj(:, i), xj(:, i)))**2
-            beta   = - mu(i)**2 * c2 
-            aj(:, i)  = 2 * beta * xj(:, i) / rjmag4
+      associate(pl => self, npl => self%nbody, inv_c2 => param%inv_c2)
+         if (npl == 0) return
+         allocate(aj, mold = pl%ah)
+         do i = 1, npl
+            rjmag4 = (dot_product(pl%xj(:, i), pl%xj(:, i)))**2
+            beta   = -pl%muj(i)**2 * inv_c2 
+            aj(:, i)  = 2 * beta * pl%xj(:, i) / rjmag4
          end do
          suma(:) = 0.0_DP
-         ah(:, 1) = ah(:, 1) + aj(:, 1)
-         do i = 2, n
-            suma(:) = suma(:) + GMpl(i) * aj(:, i) / eta(i)
-            ah(:, i) = ah(:, i) + aj(:, i) + suma(:)
+         pl%ah(:, 1) = pl%ah(:, 1) + aj(:, 1)
+         do i = 2, npl
+            suma(:) = suma(:) + pl%Gmass(i) * aj(:, i) / pl%eta(i)
+            pl%ah(:, i) = pl%ah(:, i) + aj(:, i) + suma(:)
          end do
       end associate
       return
@@ -51,13 +50,12 @@ module subroutine whm_gr_getacch_tp(self, param)
       integer(I4B)                                 :: i
       real(DP)                                     :: rjmag4, beta
 
-      associate(n => self%nbody, mu => self%mu,& 
-         c2 => param%inv_c2, ah => self%ah, xh => self%xh, status => self%status)
-         if (n == 0) return
-         do i = 1, n
-            rjmag4 = (dot_product(xh(:, i), xh(:, i)))**2
-            beta = - mu(i)**2 * c2 
-            ah(:, i) = ah(:, i) + beta * xh(:, i) / rjmag4
+      associate(tp => self, ntp => self%nbody, inv_c2 => param%inv_c2)
+         if (ntp == 0) return
+         do i = 1, ntp
+            rjmag4 = (dot_product(tp%xh(:, i), tp%xh(:, i)))**2
+            beta = - tp%mu(i)**2 * inv_c2 
+            tp%ah(:, i) = tp%ah(:, i) + beta * tp%xh(:, i) / rjmag4
          end do
       end associate
       return
@@ -113,100 +111,4 @@ module pure subroutine whm_gr_p4_tp(self, param, dt)
      return
    end subroutine whm_gr_p4_tp
 
-   module pure subroutine whm_gr_pv2vh_pl(self, param)
-      !! author: David A. Minton
-      !!
-      !! Wrapper function that converts from pseudovelocity to heliocentric velocity for massive bodies
-      !! in a WHM object
-      implicit none
-      ! Arguments
-      class(whm_pl),              intent(inout) :: self  !! Swiftest particle object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
-      ! Internals
-      integer(I4B)                              :: i
-      real(DP), dimension(:,:), allocatable     :: vh    !! Temporary holder of pseudovelocity for in-place conversion
-
-      associate(pl => self, npl => self%nbody)
-         if (npl == 0) return
-         allocate(vh, mold = pl%vh)
-         do i = 1, npl
-            call gr_pseudovel2vel(param, pl%mu(i), pl%xh(:, i), pl%vh(:, i), vh(:, i))
-         end do
-         call move_alloc(vh, pl%vh)
-      end associate
-      return
-   end subroutine whm_gr_pv2vh_pl
-
-   module pure subroutine whm_gr_pv2vh_tp(self, param)
-      !! author: David A. Minton
-      !!
-      !! Wrapper function that converts from pseudovelocity to heliocentric velocity for test particles bodies
-      !! in a WHM object
-      implicit none
-      ! Arguments
-      class(whm_tp),              intent(inout) :: self   !! Swiftest particle object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
-      ! Internals
-      integer(I4B)                                 :: i
-      real(DP), dimension(:,:), allocatable        :: vh !! Temporary holder of pseudovelocity for in-place conversion
-
-      associate(tp => self, ntp => self%nbody)
-         if (ntp == 0) return
-         allocate(vh, mold = tp%vh)
-         do i = 1, ntp
-            call gr_pseudovel2vel(param, tp%mu(i), tp%xh(:, i), tp%vh(:, i), vh(:, i))
-         end do
-         call move_alloc(vh, tp%vh)
-      end associate
-      return
-   end subroutine whm_gr_pv2vh_tp
-
-   module pure subroutine whm_gr_vh2pv_pl(self, param)
-      !! author: David A. Minton
-      !!
-      !! Wrapper function that converts from heliocentric velocity to pseudovelocity for massive bodies
-      !! in a WHM object
-      implicit none
-      ! Arguments
-      class(whm_pl),              intent(inout) :: self   !! Swiftest particle object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
-      ! Internals
-      integer(I4B)                                 :: i
-      real(DP), dimension(:,:), allocatable        :: pv !! Temporary holder of pseudovelocity for in-place conversion
-
-      associate(pl => self, npl => self%nbody)
-         if (npl == 0) return
-         allocate(pv, mold = pl%vh)
-         do i = 1, npl
-            call gr_vel2pseudovel(param, pl%mu(i), pl%xh(:, i), pl%vh(:, i), pv(:, i))
-         end do
-         call move_alloc(pv, pl%vh)
-      end associate
-      return
-   end subroutine whm_gr_vh2pv_pl
-
-   module pure subroutine whm_gr_vh2pv_tp(self, param)
-      !! author: David A. Minton
-      !!
-      !! Wrapper function that converts from heliocentric velocity to pseudovelocity for teset particles
-      !! in a WHM object
-      implicit none
-      ! Arguments
-      class(whm_tp),              intent(inout) :: self   !! Swiftest particle object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
-      ! Internals
-      integer(I4B)                                 :: i
-      real(DP), dimension(:,:), allocatable        :: pv !! Temporary holder of pseudovelocity for in-place conversion
-
-      associate(tp => self, ntp => self%nbody)
-         if (ntp == 0) return
-         allocate(pv, mold = tp%vh)
-         do i = 1, ntp
-            call gr_vel2pseudovel(param, tp%mu(i), tp%xh(:, i), tp%vh(:, i), pv(:, i))
-         end do
-         call move_alloc(pv, tp%vh)
-      end associate
-      return
-   end subroutine whm_gr_vh2pv_tp
-
 end submodule s_whm_gr