Restructured so that kicks and encounters can be timed and adapted se…

…parately

src/io/io.f90

@@ -815,21 +815,25 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
          case("ADAPTIVE")
             param%ladaptive_interactions = .true.
             param%lflatten_interactions = .true.
+            param%lflatten_encounters = .true.
             call io_log_start(param, INTERACTION_TIMER_LOG_OUT, "Interaction loop timer logfile")
             call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric")
          case("TRIANGULAR")
             param%ladaptive_interactions = .false.
             param%lflatten_interactions = .false.
+            param%lflatten_encounters = .false.
          case("FLAT")
             param%ladaptive_interactions = .false.
             param%lflatten_interactions = .true.
+            param%lflatten_encounters = .true.
          case default
             write(*,*) "Unknown value for parameter INTERACTION_LOOPS: -> ",trim(adjustl(param%interaction_loops))
             write(*,*) "Must be one of the following: TRIANGULAR, FLAT, or ADAPTIVE"
             write(*,*) "Using default value of ADAPTIVE"
             param%interaction_loops = "ADAPTIVE"
             param%ladaptive_interactions = .true.
             param%lflatten_interactions = .true.
+            param%lflatten_encounters = .true.
          end select
 
          iostat = 0

src/kick/kick.f90

@@ -22,8 +22,9 @@ module subroutine kick_getacch_int_pl(self, param)
 
       if (param%ladaptive_interactions) then
          if (lfirst) then
+            write(itimer%loopname, *) "kick_getacch_int_pl"
+            write(itimer%looptype, *) "INTERACTION"
             call itimer%time_this_loop(param, self, self%nplpl)
-            write(itimer%loopname, *)  "kick_getacch_int_pl"
             lfirst = .false.
          else
             if (itimer%check(param, self%nplpl)) call itimer%time_this_loop(param, self, self%nplpl)

src/modules/symba_classes.f90

@@ -259,7 +259,7 @@ end subroutine symba_drift_tp
       module function symba_encounter_check_pl(self, param, system, dt, irec) result(lany_encounter)
          implicit none
          class(symba_pl),            intent(inout) :: self           !! SyMBA test particle object  
-         class(swiftest_parameters), intent(in)    :: param          !! Current swiftest run configuration parameters
+         class(swiftest_parameters), intent(inout) :: param          !! Current swiftest run configuration parameters
          class(symba_nbody_system),  intent(inout) :: system         !! SyMBA nbody system object
          real(DP),                   intent(in)    :: dt             !! step size
          integer(I4B),               intent(in)    :: irec           !! Current recursion level

src/modules/walltime_classes.f90

@@ -24,7 +24,8 @@ module walltime_classes
    end type walltimer
 
    type, extends(walltimer) :: interaction_timer
-      character(len=STRMAX) :: loopname             !! Stores the name of the loop being timed for logging purposes
+      character(len=STRMAX) :: loopname              !! Stores the name of the loop being timed for logging purposes
+      character(len=NAMELEN) :: looptype             !! Stores the type of loop (e.g. INTERACTION or ENCOUNTER)
       integer(I8B) :: max_interactions = huge(1_I8B) !! Stores the number of pl-pl interactions that failed when attempting to flatten (e.g. out of memory). Adapting won't occur if ninteractions > max_interactions
       integer(I8B) :: last_interactions = 0          !! Number of interactions that were computed last time. The timer is only run if there has been a change to the number of interactions
       integer(I4B) :: step_counter = 0               !! Number of steps that have elapsed since the last timed loop

src/symba/symba_encounter_check.f90

@@ -130,7 +130,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
       implicit none
       ! Arguments
       class(symba_pl),            intent(inout)  :: self   !! SyMBA test particle object  
-      class(swiftest_parameters), intent(in)     :: param  !! Current swiftest run configuration parameters
+      class(swiftest_parameters), intent(inout)  :: param  !! Current swiftest run configuration parameters
       class(symba_nbody_system),  intent(inout)  :: system !! SyMBA nbody system object
       real(DP),                   intent(in)     :: dt     !! step size
       integer(I4B),               intent(in)     :: irec   !! Current recursion level
@@ -142,13 +142,28 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
       logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr
       integer(I4B), dimension(:), allocatable :: index1, index2
       integer(I4B), dimension(:,:), allocatable :: k_plpl_enc 
+      type(interaction_timer), save :: itimer
+      logical, save :: lfirst = .true.
 
       if (self%nbody == 0) return
 
       associate(pl => self, plplenc_list => system%plplenc_list)
+
+         if (param%ladaptive_interactions) then
+            if (lfirst) then
+               write(itimer%loopname, *)  "symba_encounter_check_pl"
+               write(itimer%looptype, *)  "ENCOUNTERS"
+               call itimer%time_this_loop(param, pl, pl%nplplm)
+               lfirst = .false.
+            else
+               if (itimer%check(param, pl%nplplm)) call itimer%time_this_loop(param, pl, pl%nplplm)
+            end if
+         end if
+
          npl = pl%nbody
          if (param%lflatten_interactions) then
             nplplm = pl%nplplm
+
             allocate(lencounter(nplplm))
             allocate(loc_lvdotr(nplplm))
 
@@ -200,6 +215,11 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
                pl%nplenc(j) = pl%nplenc(j) + 1
             end do
          end if
+
+         if (param%ladaptive_interactions) then 
+            if (itimer%is_on) call itimer%adapt(param, pl, pl%nplplm)
+         end if
+
       end associate
 
       return

src/symba/symba_kick.f90

@@ -19,11 +19,12 @@ module subroutine symba_kick_getacch_int_pl(self, param)
 
       if (param%ladaptive_interactions) then
          if (lfirst) then
-            call itimer%time_this_loop(param, self, self%nplpl)
             write(itimer%loopname, *)  "symba_kick_getacch_int_pl"
+            write(itimer%looptype, *)  "INTERACTION"
+            call itimer%time_this_loop(param, self, self%nplplm)
             lfirst = .false.
          else
-            if (itimer%check(param, self%nplpl)) call itimer%time_this_loop(param, self, self%nplpl)
+            if (itimer%check(param, self%nplplm)) call itimer%time_this_loop(param, self, self%nplplm)
          end if
       end if
 
@@ -34,7 +35,7 @@ module subroutine symba_kick_getacch_int_pl(self, param)
       end if
 
       if (param%ladaptive_interactions) then 
-         if (itimer%is_on) call itimer%adapt(param, self, self%nplpl)
+         if (itimer%is_on) call itimer%adapt(param, self, self%nplplm)
       end if
 
       return

src/symba/symba_util.f90

@@ -285,25 +285,28 @@ module subroutine symba_util_flatten_eucl_plpl(self, param)
       class(symba_pl),            intent(inout) :: self  !! SyMBA massive body object
       class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
       ! Internals
-      integer(I8B) :: k, nplpl, nplplm
-      integer(I4B) :: i, j, npl, nplm, ip, jp
+      integer(I8B) :: k
+      integer(I4B) :: i, j, npl, nplm, err
 
-      associate(pl => self)
+      associate(pl => self, nplpl => self%nplpl, nplplm => self%nplplm)
          npl = int(self%nbody, kind=I8B)
          nplm = count(.not. pl%lmtiny(1:npl))
          pl%nplm = int(nplm, kind=I4B)
-         pl%nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column
-         pl%nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies
-         if (param%lflatten_interactions) then
+         nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column
+         nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies
+         if ((param%lflatten_interactions) .or. (param%lflatten_encounters)) then
             if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously
-            allocate(self%k_plpl(2, pl%nplpl))
-            do concurrent (i = 1:npl)
-               do concurrent (j = i+1:npl)
+            allocate(self%k_plpl(2, nplpl), stat=err)
+            if (err /=0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode
+               param%lflatten_interactions = .false.
+               param%lflatten_encounters = .false.
+            else
+               do concurrent (i=1:npl, j=1:npl, j>i)
                   call util_flatten_eucl_ij_to_k(npl, i, j, k)
                   self%k_plpl(1, k) = i
                   self%k_plpl(2, k) = j
                end do
-            end do
+            end if
          end if
       end associate
 

src/util/util_flatten.f90

@@ -82,11 +82,12 @@ module subroutine util_flatten_eucl_plpl(self, param)
       npl = int(self%nbody, kind=I8B)
       associate(nplpl => self%nplpl)
          nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl
-         if (param%lflatten_interactions) then
+         if ((param%lflatten_interactions) .or. (param%lflatten_encounters)) then
             if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously
             allocate(self%k_plpl(2, nplpl), stat=err)
             if (err /=0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode
                param%lflatten_interactions = .false.
+               param%lflatten_encounters = .false.
             else
                do concurrent (i=1:npl, j=1:npl, j>i)
                   call util_flatten_eucl_ij_to_k(npl, i, j, k)

src/walltime/walltime.f90

@@ -93,41 +93,53 @@ module subroutine walltime_interaction_adapt(self, param, pl, ninteractions)
       character(len=STRMAX) :: tstr, nstr, cstr, mstr
       character(len=11) :: lstyle
       character(len=1) :: schar
+      logical :: lflatten_final
 
       ! Record the elapsed time 
       call system_clock(self%count_finish_step)
 
-      if (param%lflatten_interactions) then
-         write(lstyle,*) "FLAT      "
-      else
-         write(lstyle,*) "TRIANGULAR"
-      end if 
       write(schar,'(I1)') self%stage
       write(nstr,*) ninteractions
 
       select case(self%stage)
       case(1)
+         if (self%stage1_is_flattened) then
+            write(lstyle,*) "FLAT      "
+         else
+            write(lstyle,*) "TRIANGULAR"
+         end if 
          self%stage1_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP)
-         write(mstr,*) self%stage2_metric
+         write(mstr,*) self%stage1_metric
       case(2)
+         if (.not.self%stage1_is_flattened) then
+            write(lstyle,*) "FLAT      "
+         else
+            write(lstyle,*) "TRIANGULAR"
+         end if 
+
          self%stage2_metric = (self%count_finish_step - self%count_start_step) / real(ninteractions, kind=DP)
          self%is_on = .false.
          self%step_counter = 0
-         if (self%stage1_metric < self%stage2_metric) call self%flip(param, pl)  ! Go back to the original style, otherwise keep the stage2 style
-         write(mstr,*) self%stage1_metric
+         if (self%stage1_metric < self%stage2_metric) then
+            lflatten_final = self%stage1_is_flattened
+            call self%flip(param, pl)  ! Go back to the original style, otherwise keep the stage2 style
+         else
+            lflatten_final = .not.self%stage1_is_flattened
+         end if
+         write(mstr,*) self%stage2_metric
       end select
 
       write(cstr,*) self%count_finish_step - self%count_start_step
 
       call io_log_one_message(INTERACTION_TIMER_LOG_OUT, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr)))
 
       if (self%stage == 2) then
-         if (param%lflatten_interactions) then
+         if (lflatten_final) then
             write(lstyle,*) "FLAT      "
          else
             write(lstyle,*) "TRIANGULAR"
          end if 
-         call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "The fastest loop method tested is " // trim(adjustl(lstyle)))
+         call io_log_one_message(INTERACTION_TIMER_LOG_OUT, trim(adjustl(self%loopname)) // ": the fastest loop method tested is " // trim(adjustl(lstyle)))
       end if
 
       return
@@ -174,8 +186,14 @@ module subroutine walltime_interaction_flip_loop_style(self, param, pl)
       class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
       class(swiftest_pl),         intent(inout) :: pl    !! Swiftest massive body object 
 
-      param%lflatten_interactions = .not. param%lflatten_interactions
-      if (param%lflatten_interactions) then
+      select case(trim(adjustl(self%looptype)))
+      case("INTERACTIONS")
+         param%lflatten_interactions = .not. param%lflatten_interactions
+      case("ENCOUNTERS")
+         param%lflatten_encounters = .not. param%lflatten_encounters
+      end select
+
+      if ((param%lflatten_interactions) .or. (param%lflatten_encounters)) then
          call pl%flatten(param)
       else
          if (allocated(pl%k_plpl)) deallocate(pl%k_plpl)
@@ -203,11 +221,21 @@ module subroutine walltime_interaction_time_this_loop(self, param, pl, ninteract
       write(tstr,*) param%t
       select case(self%stage)
       case(1)
-         self%stage1_ninteractions = ninteractions      
-         self%stage1_is_flattened = param%lflatten_interactions
+         self%stage1_ninteractions = ninteractions 
+         select case(trim(adjustl(self%looptype)))
+         case("INTERACTIONS")
+            self%stage1_is_flattened = param%lflatten_interactions
+         case("ENCOUNTERS")
+            self%stage1_is_flattened = param%lflatten_encounters
+         end select
          call io_log_one_message(INTERACTION_TIMER_LOG_OUT, trim(adjustl(self%loopname)) // ": loop timer turned on at t = " // trim(adjustl(tstr)))
       case(2)
-         param%lflatten_interactions = self%stage1_is_flattened
+         select case(trim(adjustl(self%looptype)))
+         case("INTERACTIONS")
+            param%lflatten_interactions = self%stage1_is_flattened 
+         case("ENCOUNTERS")
+            param%lflatten_encounters = self%stage1_is_flattened 
+         end select
          call self%flip(param, pl) 
       case default
          self%stage = 1