Added specific subroutines for evaluating pl-plm encounters in SyMBA

src/io/io.f90

@@ -650,9 +650,16 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
                case ("INTERACTION_LOOPS")
                   call io_toupper(param_value)
                   param%interaction_loops = param_value
+               case ("ENCOUNTER_CHECK_PLPL")
+                  call io_toupper(param_value)
+                  param%encounter_check_plpl = param_value
+               case ("ENCOUNTER_CHECK_PLTP")
+                  call io_toupper(param_value)
+                  param%encounter_check_pltp = param_value
                case ("ENCOUNTER_CHECK")
                   call io_toupper(param_value)
-                  param%encounter_check = param_value
+                  param%encounter_check_plpl = param_value
+                  param%encounter_check_pltp = param_value
                case ("FIRSTKICK")
                   call io_toupper(param_value)
                   if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false. 
@@ -836,28 +843,50 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
             call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric")
          end select
 
+         select case(trim(adjustl(param%encounter_check_plpl)))
+         case("ADAPTIVE")
+            param%ladaptive_encounters_plpl = .true.
+            param%lencounter_sas_plpl = .true.
+            call io_log_start(param, ENCOUNTER_PLPL_TIMER_LOG_OUT, "Encounter check loop timer logfile")
+            call io_log_one_message(ENCOUNTER_PLPL_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric")
+         case("TRIANGULAR")
+            param%ladaptive_encounters_plpl = .false.
+            param%lencounter_sas_plpl = .false.
+         case("SORTSWEEP")
+            param%ladaptive_encounters_plpl = .false.
+            param%lencounter_sas_plpl = .true.
+         case default
+            write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLPL: -> ",trim(adjustl(param%encounter_check_plpl))
+            write(*,*) "Must be one of the following: TRIANGULAR, SORTSWEEP, or ADAPTIVE"
+            write(*,*) "Using default value of ADAPTIVE"
+            param%encounter_check_plpl = "ADAPTIVE"
+            param%ladaptive_encounters_plpl = .true.
+            param%lencounter_sas_plpl = .true.
+            call io_log_start(param, ENCOUNTER_PLPL_TIMER_LOG_OUT, "Encounter check loop timer logfile")
+            call io_log_one_message(ENCOUNTER_PLPL_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric")
+         end select
 
-         select case(trim(adjustl(param%encounter_check)))
+         select case(trim(adjustl(param%encounter_check_pltp)))
          case("ADAPTIVE")
-            param%ladaptive_encounters = .true.
-            param%lencounter_sas = .true.
-            call io_log_start(param, ENCOUNTER_TIMER_LOG_OUT, "Encounter check loop timer logfile")
-            call io_log_one_message(ENCOUNTER_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric")
+            param%ladaptive_encounters_pltp = .true.
+            param%lencounter_sas_pltp = .true.
+            call io_log_start(param, ENCOUNTER_PLTP_TIMER_LOG_OUT, "Encounter check loop timer logfile")
+            call io_log_one_message(ENCOUNTER_PLTP_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, npltp, metric")
          case("TRIANGULAR")
-            param%ladaptive_encounters = .false.
-            param%lencounter_sas = .false.
+            param%ladaptive_encounters_pltp = .false.
+            param%lencounter_sas_pltp = .false.
          case("SORTSWEEP")
-            param%ladaptive_encounters = .false.
-            param%lencounter_sas = .true.
+            param%ladaptive_encounters_pltp = .false.
+            param%lencounter_sas_pltp = .true.
          case default
-            write(*,*) "Unknown value for parameter ENCOUNTER_CHECK: -> ",trim(adjustl(param%encounter_check))
+            write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLTP: -> ",trim(adjustl(param%encounter_check_pltp))
             write(*,*) "Must be one of the following: TRIANGULAR, SORTSWEEP, or ADAPTIVE"
             write(*,*) "Using default value of ADAPTIVE"
-            param%encounter_check = "ADAPTIVE"
-            param%ladaptive_encounters = .true.
-            param%lencounter_sas = .true.
-            call io_log_start(param, ENCOUNTER_TIMER_LOG_OUT, "Encounter check loop timer logfile")
-            call io_log_one_message(ENCOUNTER_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric")
+            param%encounter_check_pltp = "ADAPTIVE"
+            param%ladaptive_encounters_pltp = .true.
+            param%lencounter_sas_pltp = .true.
+            call io_log_start(param, ENCOUNTER_PLTP_TIMER_LOG_OUT, "Encounter check loop timer logfile")
+            call io_log_one_message(ENCOUNTER_PLTP_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, npltp, metric")
          end select
 
          iostat = 0
@@ -950,7 +979,8 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
          call io_param_writer_one("ROTATION", param%lrotation, unit)
          call io_param_writer_one("TIDES", param%ltides, unit)
          call io_param_writer_one("INTERACTION_LOOPS", param%interaction_loops, unit)
-         call io_param_writer_one("ENCOUNTER_CHECK", param%encounter_check, unit)
+         call io_param_writer_one("ENCOUNTER_CHECK_PLPL", param%encounter_check_plpl, unit)
+         call io_param_writer_one("ENCOUNTER_CHECK_PLTP", param%encounter_check_pltp, unit)
 
          if (param%lenergy) then
             call io_param_writer_one("FIRSTENERGY", param%lfirstenergy, unit)

src/modules/encounter_classes.f90

@@ -64,12 +64,11 @@ module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc
          logical,      dimension(:),   intent(out) :: lvdotr     !! Logical array indicating which pairs are approaching
       end subroutine encounter_check_all
 
-      module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc)
+      module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, index1, index2, nenc)
          import swiftest_parameters
          implicit none
          class(swiftest_parameters),              intent(inout) :: param  !! Current Swiftest run configuration parameter5s
          integer(I4B),                            intent(in)    :: npl    !! Total number of massive bodies
-         integer(I4B),                            intent(in)    :: nplm   !! Number of fully interacting massive bodies
          real(DP),     dimension(:,:),            intent(in)    :: x      !! 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 
@@ -80,6 +79,25 @@ module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvd
          integer(I4B),                            intent(out)   :: nenc   !! Total number of encounters
       end subroutine encounter_check_all_plpl
 
+      module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc)
+         import swiftest_parameters
+         implicit none
+         class(swiftest_parameters),              intent(inout) :: param  !! Current Swiftest run configuration parameter5s
+         integer(I4B),                            intent(in)    :: nplm   !! Total number of fully interacting massive bodies 
+         integer(I4B),                            intent(in)    :: nplt   !! Total number of partially interacting masive bodies (GM < GMTINY) 
+         real(DP),     dimension(:,:),            intent(in)    :: xplm   !! Position vectors of fully interacting massive bodies
+         real(DP),     dimension(:,:),            intent(in)    :: vplm   !! Velocity vectors of fully interacting massive bodies
+         real(DP),     dimension(:,:),            intent(in)    :: xplt   !! Position vectors of partially interacting massive bodies
+         real(DP),     dimension(:,:),            intent(in)    :: vplt   !! Velocity vectors of partially interacting massive bodies
+         real(DP),     dimension(:),              intent(in)    :: rencm  !! Critical radii of fully interacting massive bodies that defines an encounter
+         real(DP),     dimension(:),              intent(in)    :: renct  !! Critical radii of partially interacting massive bodies that defines an encounter
+         real(DP),                                intent(in)    :: dt     !! Step size
+         logical,      dimension(:), allocatable, intent(out)   :: lvdotr !! Logical flag indicating the sign of v .dot. x
+         integer(I4B), dimension(:), allocatable, intent(out)   :: index1 !! List of indices for body 1 in each encounter
+         integer(I4B), dimension(:), allocatable, intent(out)   :: index2 !! List of indices for body 2 in each encounter
+         integer(I4B),                            intent(out)   :: nenc   !! Total number of encounters
+      end subroutine encounter_check_all_plplm
+
       module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc)
          import swiftest_parameters
          implicit none

src/modules/swiftest_classes.f90

@@ -126,12 +126,15 @@ module swiftest_classes
       real(DP)             :: inv_c2         = -1.0_DP            !! Inverse speed of light squared in the system units
       character(STRMAX)    :: energy_out     = ""                 !! Name of output energy and momentum report file
       character(NAMELEN)   :: interaction_loops = "ADAPTIVE"      !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" 
-      character(NAMELEN)   :: encounter_check = "ADAPTIVE"        !! Method used to compute encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" 
+      character(NAMELEN)   :: encounter_check_plpl = "ADAPTIVE"   !! Method used to compute pl-pl encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" 
+      character(NAMELEN)   :: encounter_check_pltp = "ADAPTIVE"   !! Method used to compute pl-tp encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" 
       ! The following are used internally, and are not set by the user, but instead are determined by the input value of INTERACTION_LOOPS
       logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interaction loops
       logical :: ladaptive_interactions = .false. !! Adaptive interaction loop is turned on (choose between TRIANGULAR and FLAT based on periodic timing tests)
-      logical :: lencounter_sas = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters
-      logical :: ladaptive_encounters = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests)
+      logical :: lencounter_sas_plpl = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters
+      logical :: lencounter_sas_pltp = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters
+      logical :: ladaptive_encounters_plpl = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests)
+      logical :: ladaptive_encounters_pltp = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests)
 
       ! Logical flags to turn on or off various features of the code
       logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually)

src/modules/walltime_classes.f90

@@ -9,7 +9,8 @@ module walltime_classes
 
    integer(I4B) :: INTERACTION_TIMER_CADENCE = 1000 !! Minimum number of steps to wait before timing an interaction loop in ADAPTIVE mode
    character(len=*), parameter :: INTERACTION_TIMER_LOG_OUT  = "interaction_timer.log" !! Name of log file for recording results of interaction loop timing
-   character(len=*), parameter :: ENCOUNTER_TIMER_LOG_OUT  = "encounter_check_timer.log" !! Name of log file for recording results of encounter check method timing
+   character(len=*), parameter :: ENCOUNTER_PLPL_TIMER_LOG_OUT  = "encounter_check_plpl_timer.log" !! Name of log file for recording results of encounter check method timing
+   character(len=*), parameter :: ENCOUNTER_PLTP_TIMER_LOG_OUT  = "encounter_check_pltp_timer.log" !! Name of log file for recording results of encounter check method timing
 
    type :: walltimer
       integer(I8B) :: count_rate                 !! Rate at wich the clock ticks

src/symba/symba_encounter_check.f90

@@ -18,7 +18,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
       logical                                   :: lany_encounter !! Returns true if there is at least one close encounter      
       ! Internals
       integer(I8B) :: k, nplplm, kenc
-      integer(I4B) :: i, j, nenc, npl, nplm
+      integer(I4B) :: i, j, nenc, npl, nplm, nplt
       logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr
       integer(I4B), dimension(:), allocatable :: index1, index2
       integer(I4B), dimension(:,:), allocatable :: k_plpl_enc 
@@ -31,7 +31,12 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
 
          npl = pl%nbody
          nplm = pl%nplm
-         call encounter_check_all_plpl(param, npl, nplm, pl%xh, pl%vh, pl%renc, dt, lvdotr, index1, index2, nenc)
+         nplt = npl - nplm
+         if (nplt == 0) then
+            call encounter_check_all_plpl(param, npl, pl%xh, pl%vh, pl%renc, dt, lvdotr, index1, index2, nenc)
+         else
+            call encounter_check_all_plplm(param, nplm, nplt, pl%xh(:,1:nplm), pl%vh(:,1:nplm), pl%xh(:,nplm+1:npl), pl%vh(:,nplm+1:npl), pl%renc(1:nplm), pl%renc(nplm+1:npl), dt, lvdotr, index1, index2, nenc)
+         end if
          lany_encounter = nenc > 0
          if (lany_encounter) then
             call plplenc_list%resize(nenc)

src/symba/symba_util.f90

@@ -290,6 +290,10 @@ module subroutine symba_util_flatten_eucl_plpl(self, param)
 
       associate(pl => self, nplpl => self%nplpl, nplplm => self%nplplm)
          npl = int(self%nbody, kind=I8B)
+         select type(param)
+         class is (symba_parameters)
+            pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY 
+         end select
          nplm = count(.not. pl%lmtiny(1:npl))
          pl%nplm = int(nplm, kind=I4B)
          nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column

src/walltime/walltime.f90

@@ -155,7 +155,7 @@ module subroutine walltime_interaction_adapt(self, param, ninteractions, pl)
       case("ENCOUNTER")
          write(advancedstyle, *) "SORTSWEEP "
          write(standardstyle, *) "TRIANGULAR"
-         write(logfile,*) ENCOUNTER_TIMER_LOG_OUT
+         write(logfile,*) ENCOUNTER_PLPL_TIMER_LOG_OUT
       case default
          write(logfile,*) "unknown_looptimer.log"
       end select
@@ -251,8 +251,10 @@ module subroutine walltime_interaction_flip_loop_style(self, param, pl)
       select case(trim(adjustl(self%looptype)))
       case("INTERACTION")
          param%lflatten_interactions = .not. param%lflatten_interactions
-      case("ENCOUNTER")
-         param%lencounter_sas = .not. param%lencounter_sas
+      case("ENCOUNTER_PLPL")
+         param%lencounter_sas_plpl= .not. param%lencounter_sas_plpl
+      case("ENCOUNTER_PLTP")
+         param%lencounter_sas_pltp= .not. param%lencounter_sas_pltp
       end select
 
       if (present(pl)) then
@@ -286,7 +288,7 @@ module subroutine walltime_interaction_time_this_loop(self, param, ninteractions
       case("INTERACTION")
          write(logfile,*) INTERACTION_TIMER_LOG_OUT
       case("ENCOUNTER")
-         write(logfile,*) ENCOUNTER_TIMER_LOG_OUT
+         write(logfile,*) ENCOUNTER_PLPL_TIMER_LOG_OUT
       case default
          write(logfile,*) "unknown_looptimer.log"
       end select
@@ -299,16 +301,20 @@ module subroutine walltime_interaction_time_this_loop(self, param, ninteractions
          select case(trim(adjustl(self%looptype)))
          case("INTERACTION")
             self%stage1_is_advanced = param%lflatten_interactions
-         case("ENCOUNTER")
-            self%stage1_is_advanced = param%lencounter_sas
+         case("ENCOUNTER_PLPL")
+            self%stage1_is_advanced = param%lencounter_sas_plpl
+         case("ENCOUNTER_PLTP")
+            self%stage1_is_advanced = param%lencounter_sas_pltp
          end select
          call io_log_one_message(logfile, trim(adjustl(self%loopname)) // ": loop timer turned on at t = " // trim(adjustl(tstr)))
       case(2)
          select case(trim(adjustl(self%looptype)))
          case("INTERACTION")
             param%lflatten_interactions = self%stage1_is_advanced 
-         case("ENCOUNTER")
-            param%lencounter_sas = self%stage1_is_advanced 
+         case("ENCOUNTER_PLPL")
+            param%lencounter_sas_plpl= self%stage1_is_advanced 
+         case("ENCOUNTER_PLTP")
+            param%lencounter_sas_pltp= self%stage1_is_advanced 
          end select
          call self%flip(param, pl) 
       case default