Moved set_beg_end method to swiftest_pl structure for better readabil…

…ity. Fixed h2b method.

src/helio/helio_step.f90

@@ -53,9 +53,9 @@ module subroutine helio_step_pl(self, system, param, t, dt)
             call pl%lindrift(system, dth, ptb)
             call pl%get_accel(system, param, t)
             call pl%kick(dth)
-            call system%set_beg_end(xbeg = pl%xh)
+            call pl%set_beg_end(xbeg = pl%xh)
             call pl%drift(system, param, dt)
-            call system%set_beg_end(xend = pl%xh)
+            call pl%set_beg_end(xend = pl%xh)
             call pl%get_accel(system, param, t + dt)
             call pl%kick(dth)
             call pl%lindrift(system, dth, pte)
@@ -87,17 +87,17 @@ module subroutine helio_step_tp(self, system, param, t, dt)
 
       select type(system)
       class is (helio_nbody_system)
-         associate(tp => self, cb => system%cb, pl => system%pl, xbeg => system%xbeg, xend => system%xend, ptb => system%ptb, pte => system%pte)
+         associate(tp => self, cb => system%cb, pl => system%pl, ptb => system%ptb, pte => system%pte)
             dth = 0.5_DP * dt
             if (tp%lfirst) then
                call tp%vh2vb(vbcb = -ptb)
                tp%lfirst = .false.
             end if
             call tp%lindrift(system, dth, ptb)
-            call tp%get_accel(system, param, t, xbeg)
+            call tp%get_accel(system, param, t, pl%xbeg)
             call tp%kick(dth)
             call tp%drift(system, param, dt)
-            call tp%get_accel(system, param, t + dt, xend)
+            call tp%get_accel(system, param, t + dt, pl%xend)
             call tp%kick(dth)
             call tp%lindrift(system, dth, pte)
             call tp%vb2vh(vbcb = -pte)

src/modules/rmvs_classes.f90

@@ -29,7 +29,6 @@ module rmvs_classes
       !> Replace the abstract procedures with concrete ones
       procedure, public :: initialize    => rmvs_setup_system      !! Performs RMVS-specific initilization steps, including generating the close encounter planetocentric structures
       procedure, public :: step          => rmvs_step_system       !! Advance the RMVS nbody system forward in time by one step
-      procedure, public :: set_beg_end   => rmvs_setup_set_beg_end !! Sets the beginning and ending values of planet positions. Also adds the end velocity for RMVS.
    end type rmvs_nbody_system
 
    type, private :: rmvs_interp
@@ -145,12 +144,6 @@ module subroutine rmvs_setup_pl(self,n)
          integer,             intent(in)    :: n    !! Number of test particles to allocate
       end subroutine rmvs_setup_pl
 
-      module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg)
-         implicit none
-         class(rmvs_nbody_system), intent(inout)          :: self !! RMVS nbody system object
-         real(DP), dimension(:,:), intent(in),   optional :: xbeg, xend, vbeg
-      end subroutine rmvs_setup_set_beg_end
-
       module subroutine rmvs_setup_system(self, param)
          use swiftest_classes, only : swiftest_parameters
          implicit none

src/modules/swiftest_classes.f90

@@ -19,8 +19,8 @@ module swiftest_classes
              setup_set_rhill, setup_tp
    public :: user_getacch_body
    public :: util_coord_b2h_pl, util_coord_b2h_tp, util_coord_h2b_pl, util_coord_h2b_tp, util_copy_body, util_copy_cb, util_copy_pl, &
-             util_copy_tp, util_copy_system, util_fill_body, util_fill_pl, util_fill_tp, util_reverse_status, util_spill_body, &
-             util_spill_pl, util_spill_tp
+             util_copy_tp, util_copy_system, util_fill_body, util_fill_pl, util_fill_tp, util_reverse_status, util_set_beg_end, & 
+             util_spill_body, util_spill_pl, util_spill_tp
 
    !********************************************************************************************************************************
    ! swiftest_parameters class definitions 
@@ -198,23 +198,27 @@ module swiftest_classes
       integer(I4B), dimension(:,:), allocatable :: k_eucl          !! Index array that converts i, j array indices into k index for use in 
                                                                    !!  the Euclidean distance matrix
       real(DP),     dimension(:),   allocatable :: irij3           !! 1.0_DP / (rji2 * sqrt(rji2)) where rji2 is the square of the Euclidean distance
+      real(DP),     dimension(:,:), allocatable :: xbeg            !! Position at beginning of step
+      real(DP),     dimension(:,:), allocatable :: xend            !! Position at end of step
+      real(DP),     dimension(:,:), allocatable :: vbeg            !! Velocity at beginning of step
       !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the
       !!    component list, such as setup_pl and util_spill_pl
    contains
       private
       ! Massive body-specific concrete methods 
       ! These are concrete because they are the same implemenation for all integrators
-      procedure, public :: discard    => discard_pl           !! Placeholder method for discarding massive bodies 
-      procedure, public :: eucl_index => eucl_dist_index_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
-      procedure, public :: eucl_irij3 => eucl_irij3_plpl      !! Parallelized single loop blocking for Euclidean distance matrix calcualtion
-      procedure, public :: setup      => setup_pl             !! A base constructor that sets the number of bodies and allocates and initializes all arrays  
-      procedure, public :: set_mu     => setup_set_mu_pl      !! Method used to construct the vectorized form of the central body mass
-      procedure, public :: set_rhill  => setup_set_rhill      !! Calculates the Hill's radii for each body
-      procedure, public :: h2b        => util_coord_h2b_pl    !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity)
-      procedure, public :: b2h        => util_coord_b2h_pl    !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity)
-      procedure, public :: copy       => util_copy_pl         !! Copies elements of one object to another.
-      procedure, public :: fill       => util_fill_pl         !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic)
-      procedure, public :: spill      => util_spill_pl        !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
+      procedure, public :: discard     => discard_pl           !! Placeholder method for discarding massive bodies 
+      procedure, public :: eucl_index  => eucl_dist_index_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
+      procedure, public :: eucl_irij3  => eucl_irij3_plpl      !! Parallelized single loop blocking for Euclidean distance matrix calcualtion
+      procedure, public :: setup       => setup_pl             !! A base constructor that sets the number of bodies and allocates and initializes all arrays  
+      procedure, public :: set_mu      => setup_set_mu_pl      !! Method used to construct the vectorized form of the central body mass
+      procedure, public :: set_rhill   => setup_set_rhill      !! Calculates the Hill's radii for each body
+      procedure, public :: h2b         => util_coord_h2b_pl    !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity)
+      procedure, public :: b2h         => util_coord_b2h_pl    !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity)
+      procedure, public :: copy        => util_copy_pl         !! Copies elements of one object to another.
+      procedure, public :: fill        => util_fill_pl         !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic)
+      procedure, public :: set_beg_end => util_set_beg_end     !! Sets the beginning and ending positions of planets.
+      procedure, public :: spill       => util_spill_pl        !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
    end type swiftest_pl
 
    !********************************************************************************************************************************
@@ -263,7 +267,6 @@ module swiftest_classes
       private
       !> Each integrator will have its own version of the step
       procedure(abstract_step_system), public, deferred :: step
-      procedure(abstract_setup_set_beg_end), public, deferred :: set_beg_end !! Sets the beginning and ending positions of planets.
 
       ! Concrete classes that are common to the basic integrator (only test particles considered for discard)
       procedure, public :: discard        => discard_system               !! Perform a discard step on the system
@@ -312,20 +315,13 @@ subroutine abstract_set_mu(self, cb)
          class(swiftest_cb),           intent(inout) :: cb   !! Swiftest central body object
       end subroutine abstract_set_mu
 
-      subroutine abstract_setup_set_beg_end(self, xbeg, xend, vbeg)
-         import swiftest_nbody_system, DP
-         class(swiftest_nbody_system),  intent(inout)     :: self !! WHM nbody system object
-         real(DP), dimension(:,:), intent(in),   optional :: xbeg, xend
-         real(DP), dimension(:,:), intent(in),   optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS
-      end subroutine abstract_setup_set_beg_end
-
       subroutine abstract_step_body(self, system, param, t, dt)
          import DP, swiftest_body, swiftest_nbody_system, swiftest_parameters
          implicit none
          class(swiftest_body),         intent(inout) :: self   !! Swiftest particle object
          class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object
          class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
-         real(DP),                     intent(in)    :: t       !! Simulation time
+         real(DP),                     intent(in)    :: t      !! Simulation time
          real(DP),                     intent(in)    :: dt     !! Current stepsize
       end subroutine abstract_step_body
 
@@ -700,7 +696,7 @@ end subroutine util_fill_body
 
       module subroutine util_fill_pl(self, inserts, lfill_list)
          implicit none
-         class(swiftest_pl),    intent(inout) :: self        !! Swiftest massive body object
+         class(swiftest_pl),    intent(inout) :: self       !! Swiftest massive body object
          class(swiftest_body),  intent(inout) :: inserts    !! Inserted object 
          logical, dimension(:), intent(in)    :: lfill_list !! Logical array of bodies to merge into the keeps
       end subroutine util_fill_pl
@@ -717,6 +713,13 @@ module subroutine util_reverse_status(self)
          class(swiftest_body),         intent(inout) :: self
       end subroutine util_reverse_status
 
+      module subroutine util_set_beg_end(self, xbeg, xend, vbeg)
+         implicit none
+         class(swiftest_pl),       intent(inout) :: self       !! Swiftest massive body object
+         real(DP), dimension(:,:), intent(in),   optional :: xbeg, xend !! Positions at beginning and end of step
+         real(DP), dimension(:,:), intent(in),   optional :: vbeg       !! vbeg is an unused variable to keep this method forward compatible with RMVS
+      end subroutine util_set_beg_end
+
       module subroutine util_spill_body(self, discards, lspill_list)
          implicit none
          class(swiftest_body), intent(inout) :: self       !! Swiftest generic body object

src/modules/whm_classes.f90

@@ -72,14 +72,11 @@ module whm_classes
    !********************************************************************************************************************************
    !> An abstract class for the WHM integrator nbody system 
    type, public, extends(swiftest_nbody_system) :: whm_nbody_system
-      !> In the WHM integrator, only test particles are discarded
-      real(DP), dimension(:,:), allocatable :: xbeg, xend    !! Positions of massive bodies at beginning and end of a step. Required in order to separate the test particle step from the massive body step
    contains
       private
       !> Replace the abstract procedures with concrete ones
       procedure, public :: initialize   => whm_setup_system      !! Performs WHM-specific initilization steps, like calculating the Jacobi masses
       procedure, public :: step         => whm_step_system       !! Advance the WHM nbody system forward in time by one step
-      procedure, public :: set_beg_end  => whm_setup_set_beg_end !! Sets the beginning and ending positions of planets.
    end type whm_nbody_system
 
    interface
@@ -216,13 +213,6 @@ module subroutine whm_setup_pl(self,n)
          integer(I4B),  intent(in)       :: n    !! Number of test particles to allocate
       end subroutine whm_setup_pl
 
-      module subroutine whm_setup_set_beg_end(self, xbeg, xend, vbeg)
-         implicit none
-         class(whm_nbody_system),  intent(inout)          :: self !! WHM nbody system object
-         real(DP), dimension(:,:), intent(in),   optional :: xbeg, xend
-         real(DP), dimension(:,:), intent(in),   optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS
-      end subroutine whm_setup_set_beg_end
-
       module subroutine whm_setup_set_ir3j(self)
          implicit none
          class(whm_pl),                intent(inout) :: self    !! WHM massive body object

src/rmvs/rmvs_encounter_check.f90

@@ -24,14 +24,14 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter)
 
       select type(pl => system%pl)
       class is (rmvs_pl)
-         associate(tp => self, ntp => self%nbody, npl => pl%nbody, xbeg => system%xbeg, vbeg => system%vbeg, rts => system%rts)
+         associate(tp => self, ntp => self%nbody, npl => pl%nbody, rts => system%rts)
             r2crit(:) = (rts * pl%rhill(:))**2
             tp%plencP(:) = 0
             do j = 1, npl
                do i = 1, ntp
                   if ((tp%status(i) /= ACTIVE).or.(tp%plencP(i) /= 0)) cycle
-                  xr(:) = tp%xh(:, i) - xbeg(:, j)
-                  vr(:) = tp%vh(:, i) - vbeg(:, j)
+                  xr(:) = tp%xh(:, i) - pl%xbeg(:, j)
+                  vr(:) = tp%vh(:, i) - pl%vbeg(:, j)
                   r2 = dot_product(xr(:), xr(:))
                   v2 = dot_product(vr(:), vr(:))
                   vdotr = dot_product(vr(:), xr(:))

src/rmvs/rmvs_setup.f90

@@ -145,32 +145,4 @@ module subroutine rmvs_setup_tp(self,n)
       return
    end subroutine rmvs_setup_tp
 
-
-
-   module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg)
-      !! author: David A. Minton
-      !! 
-      !! Sets one or more of the values of xbeg, xend, and vbeg
-      implicit none
-      ! Arguments
-      class(rmvs_nbody_system), intent(inout)          :: self !! RMVS test particle object
-      real(DP), dimension(:,:), intent(in),   optional :: xbeg, xend, vbeg
-
-      if (present(xbeg)) then
-         if (allocated(self%xbeg)) deallocate(self%xbeg)
-         allocate(self%xbeg, source=xbeg)
-      end if
-      if (present(xend)) then
-         if (allocated(self%xend)) deallocate(self%xend)
-         allocate(self%xend, source=xend)
-      end if
-      if (present(vbeg)) then
-         if (allocated(self%vbeg)) deallocate(self%vbeg)
-         allocate(self%vbeg, source=vbeg)
-      end if
-
-      return
-
-   end subroutine rmvs_setup_set_beg_end
-
 end submodule s_rmvs_setup

src/rmvs/rmvs_step.f90

@@ -30,7 +30,7 @@ module subroutine rmvs_step_system(self, param, t, dt)
                   allocate(xbeg, source=pl%xh)
                   allocate(vbeg, source=pl%vh)
                   call pl%set_rhill(cb)
-                  call system%set_beg_end(xbeg = xbeg, vbeg = vbeg)
+                  call pl%set_beg_end(xbeg = xbeg, vbeg = vbeg)
                   ! ****** Check for close encounters ***** !
                   system%rts = RHSCALE
                   lencounter = tp%encounter_check(system, dt)
@@ -42,11 +42,11 @@ module subroutine rmvs_step_system(self, param, t, dt)
                      call pl%step(system, param, t, dt) 
                      pl%outer(NTENC)%x(:,:) = pl%xh(:,:)
                      pl%outer(NTENC)%v(:,:) = pl%vh(:,:)
-                     call system%set_beg_end(xend = pl%xh)
+                     call pl%set_beg_end(xend = pl%xh)
                      call rmvs_interp_out(system, param, dt)
                      call rmvs_step_out(system, param, t, dt) 
                      call tp%reverse_status()
-                     call system%set_beg_end(xbeg = xbeg, xend = xend)
+                     call pl%set_beg_end(xbeg = xbeg, xend = xend)
                      tp%lfirst = .true.
                      call tp%step(system, param, t, dt)
                      where (tp%status(:) == INACTIVE) tp%status(:) = ACTIVE
@@ -95,7 +95,7 @@ subroutine rmvs_step_out(system, param, t, dt)
                end where
                do outer_index = 1, NTENC
                   outer_time = t + (outer_index - 1) * dto
-                  call system%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, :), &
+                  call pl%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, :), &
                                           vbeg = pl%outer(outer_index - 1)%v(:, :), &
                                           xend = pl%outer(outer_index    )%x(:, :))
                   system%rts = RHPSCALE
@@ -167,8 +167,8 @@ subroutine rmvs_step_in(system, param, outer_time, dto)
                               lfirsttp = .true.
                               do inner_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level
                                  plenci%xh(:,:) = plenci%inner(inner_index - 1)%x(:,:)
-                                 call planetocen_system%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, &
-                                                                    xend = plenci%inner(inner_index)%x)
+                                 call plenci%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, &
+                                                         xend = plenci%inner(inner_index)%x)
                                  call tpenci%step(planetocen_system, param, inner_time, dti)
                                  do j = 1, pl%nenc(i)
                                     tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(inner_index)%x(:,i)

src/util/util_coord.f90

@@ -15,16 +15,22 @@ module subroutine util_coord_h2b_pl(self, cb)
       ! Internals
       integer(I4B)  :: i
       real(DP)      :: msys
+      real(DP), dimension(NDIM) :: xtmp, vtmp
 
-      associate(npl => self%nbody, xbcb => cb%xb, vbcb => cb%vb, Mcb => cb%Gmass, &
-         xb => self%xb, xh => self%xh, vb => self%vb, vh => self%vh, Mpl => self%Gmass)
-
-         msys = Mcb + sum(Mpl(1:npl))
-         do i = 1, NDIM
-            xbcb(i) = -sum(Mpl(1:npl) * xh(i, 1:npl)) / Mcb
-            vbcb(i) = -sum(Mpl(1:npl) * vh(i, 1:npl)) / Mcb
-            xb(i, 1:npl) = xh(i, 1:npl) + xbcb(i)
-            vb(i, 1:npl) = vh(i, 1:npl) + vbcb(i)
+      associate(pl => self, npl => self%nbody)
+         msys = cb%Gmass
+         xtmp(:) = 0.0_DP
+         vtmp(:) = 0.0_DP
+         do i = 1, npl
+            msys = msys + pl%Gmass(i)
+            xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i)
+            vtmp(:) = vtmp(:) + pl%Gmass(i) * pl%vh(:,i)
+         end do
+         cb%xb(:) = -xtmp(:) / msys
+         cb%vb(:) = -vtmp(:) / msys
+         do i = 1, npl
+            pl%xb(:,i) = pl%xh(:,i) + cb%xb(:)
+            pl%vb(:,i) = pl%vh(:,i) + cb%vb(:)
          end do
       end associate