Restructured kick methods and refactored getacch to get_accel method.…

… Created basic helio setup method

src/drift/drift.f90

@@ -18,10 +18,10 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag
       !! Adapted from Hal Levison and Martin Duncan's Swift routine drift_one.f
       implicit none
       ! Arguments
-      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
-      real(DP), intent(in)                   :: dt    !! Step size
-      integer(I4B), intent(out)              :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
+      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
+      real(DP), intent(in)      :: dt                     !! Step size
+      integer(I4B), intent(out) :: iflag                  !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
       ! Internals
       integer(I4B) :: i
       real(DP)   :: dttmp

src/gr/gr.f90

@@ -11,11 +11,11 @@ subroutine gr_getaccb_ns_body(self, cb, param, agr, agr0)
       !! 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_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
+      real(DP), dimension(:, :),  intent(inout) :: agr
+      real(DP), dimension(NDIM),  intent(out)   :: agr0
       ! Internals
       real(DP), dimension(NDIM) :: xh, vh
       real(DP)                  :: rmag, rdotv, vmag2

src/helio/helio_getacch.f90

@@ -27,7 +27,7 @@ module subroutine helio_getacch_pl(self, system, param, t)
          pl%ah(:,:) = pl%ahi(:,:)
          if (param%loblatecb) call pl%obl_acc(cb)
          if (param%lextra_force) call pl%user_getacch(system, param, t)
-         if (param%lgr) call pl%gr_getacch(param)
+         !if (param%lgr) call pl%gr_get_accel(param)
       end associate
 
       return
@@ -61,7 +61,7 @@ module subroutine helio_getacch_tp(self, system, param, t, xhp)
                tp%ah(:,:) = tp%ahi(:,:)
                if (param%loblatecb) call tp%obl_acc(cb) 
                if (param%lextra_force) call tp%user_getacch(system, param, t)
-               if (param%lgr) call tp%gr_getacch(param)
+               !if (param%lgr) call tp%gr_get_accel(param)
             end select
          end associate
          return

src/helio/helio_kick.f90

@@ -0,0 +1,53 @@
+submodule(helio_classes) s_helio_kick
+   use swiftest
+contains
+   module subroutine helio_kickvb_pl(self, dt)
+      !! 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
+      implicit none
+      ! Arguments
+      class(helio_pl), intent(inout) :: self !! Swiftest generic body object
+      real(DP),        intent(in)    :: dt   !! Stepsize
+      ! Internals
+      integer(I4B) :: i
+
+      associate(pl => self, npl => self%nbody)
+         if (npl ==0) return
+         do concurrent(i = 1:npl, pl%status(i) == ACTIVE) 
+            pl%vb(:, i) = pl%vb(:, i) + pl%ah(:, i) * dt
+         end do
+      end associate
+
+      return
+
+   end subroutine helio_kickvb_pl
+
+   module subroutine helio_kickvb_tp(self, dt)
+      !! 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
+      implicit none
+      ! Arguments
+      class(helio_tp), intent(inout) :: self !! Swiftest generic body object
+      real(DP),        intent(in)    :: dt   !! Stepsize
+      ! Internals
+      integer(I4B) :: i
+
+      associate(tp => self, ntp => self%nbody)
+         if (ntp ==0) return
+         do concurrent(i = 1:ntp, tp%status(i) == ACTIVE) 
+            tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt
+         end do
+      end associate
+
+      return
+
+   end subroutine helio_kickvb_tp
+end submodule s_helio_kick

src/helio/helio_setup.f90

@@ -1,6 +1,20 @@
 submodule(helio_classes) s_helio_setup
    use swiftest
 contains
+   module subroutine helio_setup_system(self, param)
+      !! author: David A. Minton
+      !!
+      !! Initialize a Helio nbody system from files 
+      implicit none
+      ! Arguments
+      class(helio_nbody_system),  intent(inout) :: self   !! Helio system object
+      class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
+
+      call io_read_initialize_system(self, param)
+      ! Make sure that the discard list gets allocated initially
+      call self%tp_discards%setup(self%tp%nbody)
+   end subroutine helio_setup_system
+
    module procedure helio_setup_pl
       !! author: David A. Minton & Carlisle A. Wishard
       !!

src/helio/helio_step.f90

@@ -64,12 +64,12 @@ module subroutine helio_step_pl(self, system, param, t, dt)
             lfirst = .false.
          end if
          call pl%lindrift(cb, dth, ptbeg)
-         call pl%getacch(system, param, t)
-         call pl%kickvb(dth)
+         call pl%get_accel(system, param, t)
+         call pl%kick(dth)
 
          call pl%drift(system, param, dt)
-         call pl%getacch(system, param, t + dt)
-         call pl%kickvb(dth)
+         call pl%get_accel(system, param, t + dt)
+         call pl%kick(dth)
          call pl%lindrift(cb, dth, ptend)
          call pl%vb2vh(cb)
       end associate
@@ -106,11 +106,11 @@ module subroutine helio_step_tp(self, system, param, t, dt)
                lfirst = .false.
             end if
             call tp%lindrift(dth, tp%ptbeg)
-            call tp%getacch(system, param, t, xbeg)
-            call tp%kickvb(dth)
+            call tp%get_accel(system, param, t, xbeg)
+            call tp%kick(dth)
             call tp%drift(system, param, dt)
-            call tp%getacch(system, param, t + dt, xend)
-            call tp%kickvb(dth)
+            call tp%get_accel(system, param, t + dt, xend)
+            call tp%kick(dth)
             call tp%lindrift(dth, tp%ptend)
             call tp%vb2vh(vbcb = -tp%ptend)
          end associate

src/kick/kick.f90

@@ -1,7 +1,7 @@
 submodule(swiftest_classes) s_kick
    use swiftest
 contains
-   module subroutine kick_vh_body(self, dt)
+   module subroutine kickvh_body(self, dt)
       !! author: David A. Minton
       !!
       !! Kick heliocentric velocities of bodies
@@ -23,30 +23,6 @@ module subroutine kick_vh_body(self, dt)
       end associate
 
       return
-   end subroutine kick_vh_body
+   end subroutine kickvh_body
 
-   module subroutine kick_vb_body(self, dt)
-      !! author: David A. Minton
-      !!
-      !! Kick barycentric velocities of bodies
-      !!
-      !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f and kickvh_tp.f
-      !! Adapted from David E. Kaufmann's Swifter routine helio_kickvb.f90 and helio_kickvb_tp.f90
-      implicit none
-      ! Arguments
-      class(swiftest_body),         intent(inout) :: self !! Swiftest generic body object
-      real(DP),                     intent(in)    :: dt   !! Stepsize
-      ! Internals
-      integer(I4B) :: i
-
-      associate(n => self%nbody, vb => self%vb, ah => self%ah, status => self%status)
-         if (n ==0) return
-         do concurrent(i = 1:n, status(i) == ACTIVE) 
-            vb(:, i) = vb(:, i) + ah(:, i) * dt
-         end do
-      end associate
-
-      return
-
-   end subroutine kick_vb_body
 end submodule s_kick

src/modules/helio_classes.f90

@@ -4,24 +4,26 @@ module helio_classes
    !! Definition of classes and methods specific to the Democratic Heliocentric Method
    !! Adapted from David E. Kaufmann's Swifter routine: helio.f90
    use swiftest_globals
-   use rmvs_classes, only : rmvs_cb, rmvs_pl, rmvs_tp, rmvs_nbody_system
+   use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_tp
+   use whm_classes, only : whm_nbody_system
    implicit none
 
 
    !********************************************************************************************************************************
    !  helio_nbody_system class definitions and method interfaces
    !********************************************************************************************************************************
-   type, public, extends(rmvs_nbody_system) :: helio_nbody_system
+   type, public, extends(whm_nbody_system) :: helio_nbody_system
    contains
       private
+      procedure, public :: initialize    => helio_setup_system  !! Performs Helio-specific initilization steps,
       procedure, public :: step          => helio_step_system
    end type helio_nbody_system
 
    !********************************************************************************************************************************
    ! helio_cb class definitions and method interfaces
    !*******************************************************************************************************************************
    !> Helio central body particle class
-   type, public, extends(rmvs_cb) :: helio_cb
+   type, public, extends(swiftest_cb) :: helio_cb
    contains
    end type helio_cb
 
@@ -30,14 +32,15 @@ module helio_classes
    !*******************************************************************************************************************************
 
    !! Helio massive body particle class
-   type, public, extends(rmvs_pl) :: helio_pl
+   type, public, extends(swiftest_pl) :: helio_pl
       real(DP),     dimension(:,:), allocatable :: ahi         !!  heliocentric acceleration due to interactions
    contains
       procedure, public :: vh2vb       => helio_coord_vh2vb_pl  !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
       procedure, public :: vb2vh       => helio_coord_vb2vh_pl  !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
       procedure, public :: drift       => helio_drift_pl        !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure, public :: lindrift    => helio_drift_linear_pl !! Method for linear drift of massive bodies due to barycentric momentum of Sun
-      procedure, public :: getacch     => helio_getacch_pl      !! Compute heliocentric accelerations of massive bodies
+      procedure, public :: get_accel   => helio_getacch_pl      !! Compute heliocentric accelerations of massive bodies
+      procedure, public :: kick        => helio_kickvb_pl       !! Kicks the barycentric velocities
       procedure, public :: setup       => helio_setup_pl        !! Constructor method - Allocates space for number of particles
       procedure, public :: step        => helio_step_pl         !! Steps the body forward one stepsize
    end type helio_pl
@@ -47,18 +50,19 @@ module helio_classes
    !*******************************************************************************************************************************
 
    !! Helio test particle class
-   type, public, extends(rmvs_tp) :: helio_tp
+   type, public, extends(swiftest_tp) :: helio_tp
       real(DP),     dimension(:,:), allocatable :: ahi     !!  heliocentric acceleration due to interactions
       real(DP),     dimension(NDIM) :: ptbeg   !! negative barycentric velocity of the Sun at beginning of time step
       real(DP),     dimension(NDIM) :: ptend   !! negative barycentric velocity of the Sun at beginning of time step
    contains
-      procedure, public :: vh2vb       => helio_coord_vh2vb_tp        !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
-      procedure, public :: vb2vh       => helio_coord_vb2vh_tp        !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
+      procedure, public :: vh2vb       => helio_coord_vh2vb_tp  !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
+      procedure, public :: vb2vh       => helio_coord_vb2vh_tp  !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
       procedure, public :: drift       => helio_drift_tp        !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure, public :: lindrift    => helio_drift_linear_tp !! Method for linear drift of massive bodies due to barycentric momentum of Sun
-      procedure, public :: getacch     => helio_getacch_tp      !! Compute heliocentric accelerations of massive bodies
+      procedure, public :: get_accel   => helio_getacch_tp      !! Compute heliocentric accelerations of massive bodies
+      procedure, public :: kick        => helio_kickvb_tp       !! Kicks the barycentric velocities
       procedure, public :: setup       => helio_setup_tp        !! Constructor method - Allocates space for number of particles
-      procedure, public :: step        => helio_step_tp             !! Steps the body forward one stepsize
+      procedure, public :: step        => helio_step_tp         !! Steps the body forward one stepsize
    end type helio_tp
 
    interface
@@ -109,23 +113,23 @@ end subroutine helio_drift_tp
       module subroutine helio_drift_linear_pl(self, cb, dt, pt)
          use swiftest_classes, only : swiftest_cb
          implicit none
-         class(helio_pl),             intent(inout) :: self !! Helio test particle data structure
-         class(swiftest_cb),          intent(in)    :: cb   !! Helio central body data structure
+         class(helio_pl),             intent(inout) :: self !! Helio test particle object
+         class(swiftest_cb),          intent(in)    :: cb   !! Helio central body object
          real(DP),                    intent(in)    :: dt   !! Stepsize
          real(DP), dimension(:),      intent(out)   :: pt   !! negative barycentric velocity of the central body
       end subroutine helio_drift_linear_pl 
 
       module subroutine helio_drift_linear_tp(self, dt, pt)
          implicit none
-         class(helio_tp),             intent(inout) :: self   !! Helio test particle data structure
+         class(helio_tp),             intent(inout) :: self   !! Helio test particle object
          real(DP),                    intent(in)    :: dt     !! Stepsize
          real(DP), dimension(:),      intent(in)    :: pt     !! negative barycentric velocity of the Sun
       end subroutine helio_drift_linear_tp
 
       module subroutine helio_getacch_pl(self, system, param, t)
          use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
          implicit none
-         class(helio_pl),              intent(inout) :: self   !! Helio massive body particle data structure
+         class(helio_pl),              intent(inout) :: self   !! Helio massive body object
          class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
          class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters of 
          real(DP),                     intent(in)    :: t      !! Current simulation time
@@ -134,19 +138,38 @@ end subroutine helio_getacch_pl
       module subroutine helio_getacch_tp(self, system, param, t, xhp)
          use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
          implicit none
-         class(helio_tp),              intent(inout) :: self   !! Helio test particle data structure
+         class(helio_tp),              intent(inout) :: self   !! Helio test particle object
          class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
          class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters
          real(DP),                     intent(in)    :: t      !! Current time
          real(DP), dimension(:,:),     intent(in)    :: xhp    !! Heliocentric positions of planets
       end subroutine helio_getacch_tp
 
+      module subroutine helio_kickvb_pl(self, dt)
+         implicit none
+         class(helio_pl), intent(inout) :: self !! Helio massive body object
+         real(DP),        intent(in)    :: dt   !! Stepsize
+      end subroutine helio_kickvb_pl
+
+      module subroutine helio_kickvb_tp(self, dt)
+         implicit none
+         class(helio_tp), intent(inout) :: self !! Helio test particle object
+         real(DP),        intent(in)    :: dt   !! Stepsize
+      end subroutine helio_kickvb_tp
+
       module subroutine helio_setup_pl(self, n)
          implicit none
          class(helio_pl), intent(inout) :: self !! Helio massive body object
          integer,         intent(in)    :: n    !! Number of test particles to allocate
       end subroutine helio_setup_pl
 
+      module subroutine helio_setup_system(self, param)
+         use swiftest_classes, only : swiftest_parameters
+         implicit none
+         class(helio_nbody_system),  intent(inout) :: self   !! Helio system object
+         class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
+      end subroutine helio_setup_system
+
       module subroutine helio_setup_tp(self,n)
          implicit none
          class(helio_tp), intent(inout) :: self !! Helio test particle object
@@ -165,7 +188,7 @@ end subroutine helio_step_system
       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 data structure
+         class(helio_pl),              intent(inout) :: self   !! Helio massive body particle object
          class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system
          class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
          real(DP),                     intent(in)    :: t      !! Current simulation time
@@ -175,7 +198,7 @@ end subroutine helio_step_pl
       module subroutine helio_step_tp(self, system, param, t, dt)
          use swiftest_classes, only : swiftest_cb, swiftest_parameters
          implicit none
-         class(helio_tp),              intent(inout) :: self   !! Helio test particle data structure
+         class(helio_tp),              intent(inout) :: self   !! Helio test particle object
          class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
          class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
          real(DP),                     intent(in)    :: t      !! Current simulation time