Simplified sort and sweep to use radial distance and set up for testing

src/encounter/encounter_module.f90

@@ -17,8 +17,8 @@ module encounter
    implicit none
    public
 
-   integer(I4B), parameter :: SWEEPDIM = 3
    character(len=*), parameter :: ENCOUNTER_OUTFILE = 'encounters.nc'  !! Name of NetCDF output file for encounter information
+   real(DP), parameter :: RSWEEP_FACTOR = 1.1_DP
 
    type, abstract :: encounter_list
       integer(I8B)                              :: nenc = 0   !! Total number of encounters
@@ -95,7 +95,7 @@ module encounter
 
 
    type encounter_bounding_box
-      type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb
+      type(encounter_bounding_box_1D) :: aabb
    contains
       procedure :: dealloc      => encounter_util_dealloc_bounding_box    !! Deallocates all allocatables
       procedure :: setup        => encounter_util_setup_aabb              !! Setup a new axis-aligned bounding box structure
@@ -107,12 +107,12 @@ module encounter
 
 
    interface
-      module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr)
+      module subroutine encounter_check_all_plpl(param, npl, r, v, renc, dt, nenc, index1, index2, lvdotr)
          use base, only: base_parameters
          implicit none
          class(base_parameters),                  intent(inout) :: param  !! Current Swiftest run configuration parameter5s
          integer(I4B),                            intent(in)    :: npl    !! Total number of massive bodies
-         real(DP),     dimension(:,:),            intent(in)    :: x      !! Position vectors of massive bodies
+         real(DP),     dimension(:,:),            intent(in)    :: r      !! Position vectors of massive bodies
          real(DP),     dimension(:,:),            intent(in)    :: v      !! Velocity vectors of massive bodies
          real(DP),     dimension(:),              intent(in)    :: renc   !! Critical radii of massive bodies that defines an encounter 
          real(DP),                                intent(in)    :: dt     !! Step size
@@ -142,15 +142,15 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt,
          logical,      dimension(:), allocatable, intent(out)   :: lvdotr !! Logical flag indicating the sign of v .dot. x
       end subroutine encounter_check_all_plplm
 
-      module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr)
+      module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr)
          use base, only: base_parameters
          implicit none
          class(base_parameters),                  intent(inout) :: param  !! Current Swiftest run configuration parameter5s
          integer(I4B),                            intent(in)    :: npl    !! Total number of massive bodies 
          integer(I4B),                            intent(in)    :: ntp    !! Total number of test particles 
          real(DP),     dimension(:,:),            intent(in)    :: rpl    !! Position vectors of massive bodies
          real(DP),     dimension(:,:),            intent(in)    :: vpl    !! Velocity vectors of massive bodies
-         real(DP),     dimension(:,:),            intent(in)    :: xtp    !! Position vectors of massive bodies
+         real(DP),     dimension(:,:),            intent(in)    :: rtp    !! Position vectors of massive bodies
          real(DP),     dimension(:,:),            intent(in)    :: vtp    !! Velocity vectors of massive bodies
          real(DP),     dimension(:),              intent(in)    :: renc   !! Critical radii of massive bodies that defines an encounter
          real(DP),                                intent(in)    :: dt     !! Step size
@@ -205,11 +205,11 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r
          logical,      dimension(:), allocatable, intent(out)   :: lvdotr     !! Logical array indicating which pairs are approaching
       end subroutine encounter_check_sweep_aabb_double_list
 
-      module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt, nenc, index1, index2, lvdotr)
+      module subroutine encounter_check_sweep_aabb_single_list(self, n, r, v, renc, dt, nenc, index1, index2, lvdotr)
          implicit none
          class(encounter_bounding_box),           intent(inout) :: self       !! Multi-dimensional bounding box structure
          integer(I4B),                            intent(in)    :: n          !! Number of bodies
-         real(DP),     dimension(:,:),            intent(in)    :: x, v       !! Array of position and velocity vectors 
+         real(DP),     dimension(:,:),            intent(in)    :: r, v       !! Array of position and velocity vectors 
          real(DP),     dimension(:),              intent(in)    :: renc       !! Radius of encounter regions of bodies 1
          real(DP),                                intent(in)    :: dt         !! Step size
          integer(I8B),                            intent(out)   :: nenc       !! Total number of encounter candidates

src/encounter/encounter_util.f90

@@ -104,9 +104,7 @@ module subroutine encounter_util_dealloc_bounding_box(self)
       ! Internals
       integer(I4B) :: i
 
-      do i = 1,NDIM
-         call self%aabb(i)%dealloc()
-      end do
+      call self%aabb%dealloc()
 
       return
    end subroutine encounter_util_dealloc_bounding_box
@@ -351,26 +349,20 @@ module subroutine encounter_util_setup_aabb(self, n, n_last)
       next_last = 2 * n_last
 
       if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies
-         do dim = 1, SWEEPDIM
-            allocate(itmp(next))
-            if (n_last > 0) itmp(1:next_last) = self%aabb(dim)%ind(1:next_last)
-            call move_alloc(itmp, self%aabb(dim)%ind)
-            self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)]
-         end do
+         allocate(itmp(next))
+         if (n_last > 0) itmp(1:next_last) = self%aabb%ind(1:next_last)
+         call move_alloc(itmp, self%aabb%ind)
+         self%aabb%ind(next_last+1:next) = [(k, k = next_last+1, next)]
       else ! The number of bodies has gone down. Resize and chop of the old indices
-         do dim = 1, SWEEPDIM
-            allocate(itmp(next))
-            itmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next)
-            call move_alloc(itmp, self%aabb(dim)%ind)
-         end do
+         allocate(itmp(next))
+         itmp(1:next) = pack(self%aabb%ind(1:next_last), self%aabb%ind(1:next_last) <= next)
+         call move_alloc(itmp, self%aabb%ind)
       end if
 
-      do dim = 1, SWEEPDIM
-         if (allocated(self%aabb(dim)%ibeg)) deallocate(self%aabb(dim)%ibeg)
-         allocate(self%aabb(dim)%ibeg(n))
-         if (allocated(self%aabb(dim)%iend)) deallocate(self%aabb(dim)%iend)
-         allocate(self%aabb(dim)%iend(n))
-      end do
+      if (allocated(self%aabb%ibeg)) deallocate(self%aabb%ibeg)
+      allocate(self%aabb%ibeg(n))
+      if (allocated(self%aabb%iend)) deallocate(self%aabb%iend)
+      allocate(self%aabb%iend(n))
 
       return
    end subroutine encounter_util_setup_aabb

src/rmvs/rmvs_encounter_check.f90

@@ -40,7 +40,7 @@ module function rmvs_encounter_check_tp(self, param, nbody_system, dt) result(le
 
       select type(pl => nbody_system%pl)
       class is (rmvs_pl)
-         associate(tp => self, ntp => self%nbody, npl => pl%nbody)
+         associate(tp => self, ntp => self%nbody, npl => pl%nbody, cb => nbody_system%cb)
             tp%plencP(1:ntp) = 0
             call encounter_check_all_pltp(param, npl, ntp, pl%rbeg, pl%vbeg, tp%rh, tp%vh, pl%renc, dt, & 
                                           nenc, index1, index2, lvdotr)

src/swiftest/swiftest_module.f90

@@ -140,8 +140,8 @@ module swiftest
       procedure :: write_info      => swiftest_io_netcdf_write_info_body    !! Dump contents of particle information metadata to file
       procedure :: accel_obl       => swiftest_obl_acc_body                 !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
       procedure :: el2xv           => swiftest_orbel_el2xv_vec              !! Convert orbital elements to position and velocity vectors
-      procedure :: xv2el           => swiftest_orbel_xv2el_vec              !! Convert position and velocity vectors to orbital  elements 
-      procedure :: setup           => swiftest_util_setup_body                   !! A constructor that sets the number of bodies and allocates all allocatable arrays
+      procedure :: xv2el           => swiftest_orbel_xv2el_vec              !! Convert position and velocity vectors to orbital elements 
+      procedure :: setup           => swiftest_util_setup_body              !! A constructor that sets the number of bodies and allocates all allocatable arrays
       procedure :: accel_user      => swiftest_user_kick_getacch_body       !! Add user-supplied heliocentric accelerations to planets
       procedure :: append          => swiftest_util_append_body             !! Appends elements from one structure to another
       procedure :: dealloc         => swiftest_util_dealloc_body            !! Deallocates all allocatable arrays

src/swiftest/swiftest_orbel.f90

@@ -680,7 +680,7 @@ real(DP) pure function swiftest_orbel_fhybrid(e,n)
       abn = n
       if(n < 0._DP) abn = -abn
 
-      if(abn < 0.636_DP * e -0.6_DP) then
+      if(abn < 0.636_DP * e - 0.6_DP) then
          swiftest_orbel_fhybrid = swiftest_orbel_flon(e,n)
       else
          swiftest_orbel_fhybrid = swiftest_orbel_fget(e,n)
@@ -715,7 +715,7 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v
       q = 0.0_DP
       x = [px, py, pz]
       v = [vx, vy, vz]
-      r = sqrt(dot_product(x(:), x(:))) 
+      r = .mag.x(:)
       v2 = dot_product(v(:), v(:))
       hvec(:) = x(:) .cross. v(:)
       h2 = dot_product(hvec(:), hvec(:))

src/symba/symba_encounter_check.f90

@@ -249,7 +249,7 @@ module function symba_encounter_check_tp(self, param, nbody_system, dt, irec) re
       lany_encounter = .false.
       if (self%nbody == 0) return
 
-      associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody)
+      associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb)
          call pl%set_renc(irec)
          call encounter_check_all_pltp(param, npl, ntp, pl%rh, pl%vb, tp%rh, tp%vb, pl%renc, dt, nenc, index1, index2, lvdotr) 
 

src/symba/symba_step.f90

@@ -27,6 +27,7 @@ module subroutine symba_step_system(self, param, t, dt)
       real(DP),                   intent(in)    :: dt     !! Current stepsize
       ! Internals
       logical :: lencounter
+      type(walltimer)        :: timer1,timer2,timer3 !! Object used for computing elapsed wall time
 
       select type(pl => self%pl)
       class is (symba_pl)
@@ -48,7 +49,6 @@ module subroutine symba_step_system(self, param, t, dt)
       end select
       end select
       end select
-
       return
    end subroutine symba_step_system