Refactored util_index to util_flatten as it is more descriptive and c…

…onsistent with the rest of the terminology. Put in safety check in case allocating k_plpl fails

src/fraggle/fraggle_generate.f90

@@ -18,7 +18,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa
       class(fraggle_fragments),     intent(inout) :: self      !! Fraggle system object the outputs will be the fragmentation 
       class(fraggle_colliders),     intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies 
       class(swiftest_nbody_system), intent(inout) :: system    !! Swiftest nbody system object
-      class(swiftest_parameters),   intent(in)    :: param     !! Current run configuration parameters 
+      class(swiftest_parameters),   intent(inout) :: param     !! Current run configuration parameters 
       logical,                      intent(out)   :: lfailure  !! Answers the question: Should this have been a merger instead?
       ! Internals
       integer(I4B)                         :: i
@@ -140,7 +140,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa
          call frag%set_original_scale(colliders)
 
          ! Restore the big array
-         if (lk_plpl) call pl%index(param)
+         if (lk_plpl) call pl%flatten(param)
       end associate
       call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes)  ! Save the current halting modes so we can turn them off temporarily
 

src/fraggle/fraggle_util.f90

@@ -135,7 +135,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para
       class(fraggle_fragments),     intent(inout) :: self      !! Fraggle fragment system object
       class(fraggle_colliders),     intent(inout) :: colliders !! Fraggle collider system object
       class(swiftest_nbody_system), intent(inout) :: system    !! Swiftest nbody system object
-      class(swiftest_parameters),   intent(in)    :: param     !! Current swiftest run configuration parameters
+      class(swiftest_parameters),   intent(inout) :: param     !! Current swiftest run configuration parameters
       logical,                      intent(in)    :: lbefore   !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa
       ! Internals
       integer(I4B) :: i, nplm
@@ -172,7 +172,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para
             call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam)
          end if 
 
-         call tmpsys%pl%index(param)
+         call tmpsys%pl%flatten(param)
 
          call tmpsys%get_energy_and_momentum(param) 
 

src/helio/helio_setup.f90

@@ -16,7 +16,7 @@ module subroutine helio_setup_initialize_system(self, param)
       call self%pl%h2b(self%cb)
       call self%tp%h2b(self%cb)
       call self%pl%sort("mass", ascending=.false.)
-      call self%pl%index(param)
+      call self%pl%flatten(param)
 
       return
    end subroutine helio_setup_initialize_system

src/modules/fraggle_classes.f90

@@ -104,7 +104,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa
          class(fraggle_fragments),     intent(inout) :: self      !! Fraggle fragment system object 
          class(fraggle_colliders),     intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies 
          class(swiftest_nbody_system), intent(inout) :: system    !! Swiftest nbody system object
-         class(swiftest_parameters),   intent(in)    :: param     !! Current run configuration parameters 
+         class(swiftest_parameters),   intent(inout) :: param     !! Current run configuration parameters 
          logical,                      intent(out)   :: lfailure  !! Answers the question: Should this have been a merger instead?
       end subroutine fraggle_generate_fragments
 
@@ -249,7 +249,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para
          class(fraggle_fragments),     intent(inout) :: self      !! Fraggle fragment system object
          class(fraggle_colliders),     intent(inout) :: colliders !! Fraggle collider system object
          class(swiftest_nbody_system), intent(inout) :: system    !! Swiftest nbody system object
-         class(swiftest_parameters),   intent(in)    :: param     !! Current swiftest run configuration parameters
+         class(swiftest_parameters),   intent(inout) :: param     !! Current swiftest run configuration parameters
          logical,                      intent(in)    :: lbefore   !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa
       end subroutine fraggle_util_get_energy_momentum
 

src/modules/swiftest_classes.f90

@@ -133,7 +133,7 @@ module swiftest_classes
       logical :: loblatecb      = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2 is input)
       logical :: lrotation      = .false. !! Include rotation states of big bodies
       logical :: ltides         = .false. !! Include tidal dissipation 
-      logical :: lflatten_interactions = .true. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory)
+      logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory)
 
       ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values)
       real(DP)                  :: Eorbit_orig = 0.0_DP   !! Initial orbital energy
@@ -324,7 +324,7 @@ module swiftest_classes
       ! Massive body-specific concrete methods 
       ! These are concrete because they are the same implemenation for all integrators
       procedure :: discard      => discard_pl             !! Placeholder method for discarding massive bodies 
-      procedure :: index        => util_index_eucl_plpl   !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
+      procedure :: flatten      => util_flatten_eucl_plpl   !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
       procedure :: accel_int    => kick_getacch_int_pl    !! Compute direct cross (third) term heliocentric accelerations of massive bodies
       procedure :: accel_obl    => obl_acc_pl             !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
       procedure :: setup        => setup_pl               !! A base constructor that sets the number of bodies and allocates and initializes all arrays  
@@ -561,34 +561,34 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag
          integer(I4B), intent(out)      :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
       end subroutine drift_one
 
-      module pure subroutine util_index_eucl_ij_to_k(n, i, j, k)
-         !$omp declare simd(util_index_eucl_ij_to_k)
+      module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
+         !$omp declare simd(util_flatten_eucl_ij_to_k)
          implicit none
          integer(I4B), intent(in)  :: n !! Number of bodies
          integer(I4B), intent(in)  :: i !! Index of the ith body
          integer(I4B), intent(in)  :: j !! Index of the jth body
          integer(I8B), intent(out) :: k !! Index of the flattened matrix
-      end subroutine util_index_eucl_ij_to_k
+      end subroutine util_flatten_eucl_ij_to_k
 
-      module pure subroutine util_index_eucl_k_to_ij(n, k, i, j)
+      module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j)
          implicit none
          integer(I4B), intent(in)  :: n !! Number of bodies
          integer(I8B), intent(in)  :: k !! Index of the flattened matrix
          integer(I4B), intent(out) :: i !! Index of the ith body
          integer(I4B), intent(out) :: j !! Index of the jth body
-      end subroutine util_index_eucl_k_to_ij
+      end subroutine util_flatten_eucl_k_to_ij
 
-      module subroutine util_index_eucl_plpl(self, param)
+      module subroutine util_flatten_eucl_plpl(self, param)
          implicit none
          class(swiftest_pl),         intent(inout) :: self  !! Swiftest massive body object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
+         class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
       end subroutine
 
-      module subroutine util_index_eucl_pltp(self, pl, param)
+      module subroutine util_flatten_eucl_pltp(self, pl, param)
          implicit none
          class(swiftest_tp),         intent(inout) :: self  !! Swiftest test particle object
          class(swiftest_pl),         intent(in)    :: pl    !! Swiftest massive body object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
+         class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
       end subroutine
 
       module pure subroutine gr_kick_getaccb_ns_body(self, system, param)

src/modules/symba_classes.f90

@@ -67,7 +67,7 @@ module symba_classes
       type(symba_kinship),       dimension(:), allocatable :: kin        !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step
    contains
       procedure :: make_colliders  => symba_collision_make_colliders_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family
-      procedure :: index           => symba_util_index_eucl_plpl     !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
+      procedure :: flatten         => symba_util_flatten_eucl_plpl     !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
       procedure :: discard         => symba_discard_pl               !! Process massive body discards
       procedure :: drift           => symba_drift_pl                 !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level
       procedure :: encounter_check => symba_encounter_check_pl       !! Checks if massive bodies are going through close encounters with each other
@@ -316,12 +316,12 @@ module function symba_collision_casemerge(system, param, colliders, frag) result
          integer(I4B)                             :: status    !! Status flag assigned to this outcome
       end function symba_collision_casemerge
 
-      module subroutine symba_util_index_eucl_plpl(self, param)
+      module subroutine symba_util_flatten_eucl_plpl(self, param)
          use swiftest_classes, only : swiftest_parameters
          implicit none
          class(symba_pl),            intent(inout) :: self  !! SyMBA massive body object
-         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
-      end subroutine symba_util_index_eucl_plpl
+         class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
+      end subroutine symba_util_flatten_eucl_plpl
 
       module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) 
          implicit none

src/modules/walltime_classes.f90

@@ -6,6 +6,8 @@ module walltime_classes
    implicit none
    public
 
+   integer(I4B) :: INTERACTION_TIMER_CADENCE = 1000 !! Minimum number of steps to wait before timing an interaction loop in ADAPTIVE mode
+
    type :: walltimer
       integer(I8B) :: count_rate                 !! Rate at wich the clock ticks
       integer(I8B) :: count_max                  !! Maximum value of the clock ticker

src/modules/whm_classes.f90

@@ -108,13 +108,6 @@ module subroutine whm_drift_pl(self, system, param, dt)
          real(DP),                     intent(in)    :: dt     !! Stepsize
       end subroutine whm_drift_pl
 
-      module subroutine whm_util_index_eucl_plpl(self, param)
-         use swiftest_classes, only : swiftest_parameters
-         implicit none
-         class(whm_pl),            intent(inout) :: self  !! WHM massive body object
-         class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters
-      end subroutine whm_util_index_eucl_plpl
-
       !> Get heliocentric accelration of massive bodies
       module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg)
          use swiftest_classes, only : swiftest_cb, swiftest_parameters

src/setup/setup.f90

@@ -175,7 +175,7 @@ module subroutine setup_initialize_system(self, param)
          call self%pl%el2xv(self%cb)
          call self%tp%el2xv(self%cb)
       end if
-      call self%pl%index(param)
+      call self%pl%flatten(param)
       if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb)
       self%pl%lfirst = param%lfirstkick
       self%tp%lfirst = param%lfirstkick

src/symba/symba_util.f90

@@ -269,7 +269,7 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list)
    end subroutine symba_util_fill_tp
 
 
-   module subroutine symba_util_index_eucl_plpl(self, param)
+   module subroutine symba_util_flatten_eucl_plpl(self, param)
       !! author: Jacob R. Elliott and David A. Minton
       !!
       !! Turns i,j indices into k index for use in the Euclidean distance matrix. This also sets the lmtiny flag and computes the
@@ -283,7 +283,7 @@ module subroutine symba_util_index_eucl_plpl(self, param)
       implicit none
       ! Arguments
       class(symba_pl),            intent(inout) :: self  !! SyMBA massive body object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
+      class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
       ! Internals
       integer(I8B) :: k, nplpl, nplplm
       integer(I4B) :: i, j, npl, nplm, ip, jp
@@ -299,7 +299,7 @@ module subroutine symba_util_index_eucl_plpl(self, param)
             allocate(self%k_plpl(2, pl%nplpl))
             do concurrent (i = 1:npl)
                do concurrent (j = i+1:npl)
-                  call util_index_eucl_ij_to_k(npl, i, j, k)
+                  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
@@ -308,7 +308,7 @@ module subroutine symba_util_index_eucl_plpl(self, param)
       end associate
 
       return
-   end subroutine symba_util_index_eucl_plpl
+   end subroutine symba_util_flatten_eucl_plpl
 
 
    module subroutine symba_util_peri_pl(self, system, param)
@@ -478,7 +478,7 @@ module subroutine symba_util_rearray_pl(self, system, param)
 
          ! Reindex the new list of bodies 
          call pl%sort("mass", ascending=.false.)
-         call pl%index(param)
+         call pl%flatten(param)
 
          ! Reset the kinship trackers
          call pl%reset_kinship([(i, i=1, npl)])

src/util/util_index.f90 → src/util/util_flatten.f90

@@ -2,8 +2,8 @@
    use swiftest
 contains
 
-   module pure subroutine util_index_eucl_ij_to_k(n, i, j, k)
-      !$omp declare simd(util_index_eucl_ij_to_k)
+   module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
+      !$omp declare simd(util_flatten_eucl_ij_to_k)
       !! author: Jacob R. Elliott and David A. Minton
       !!
       !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions.
@@ -27,10 +27,10 @@ module pure subroutine util_index_eucl_ij_to_k(n, i, j, k)
       k = (i8 - 1_I8B) * n8 - i8 * (i8 - 1_I8B) / 2_I8B + (j8 - i8)
 
       return
-   end subroutine util_index_eucl_ij_to_k
+   end subroutine util_flatten_eucl_ij_to_k
 
 
-   module pure subroutine util_index_eucl_k_to_ij(n, k, i, j)
+   module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j)
       !! author: Jacob R. Elliott and David A. Minton
       !!
       !! Turns k index into i,j indices for use in the Euclidean distance matrix for pl-pl interactions.
@@ -59,10 +59,10 @@ module pure subroutine util_index_eucl_k_to_ij(n, k, i, j)
       j = int(j8, kind=I4B)
 
       return
-   end subroutine util_index_eucl_k_to_ij
+   end subroutine util_flatten_eucl_k_to_ij
 
 
-   module subroutine util_index_eucl_plpl(self, param)
+   module subroutine util_flatten_eucl_plpl(self, param)
       !! author: Jacob R. Elliott and David A. Minton
       !!
       !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions for a Swiftest massive body object
@@ -74,30 +74,34 @@ module subroutine util_index_eucl_plpl(self, param)
       implicit none
       ! Arguments
       class(swiftest_pl),         intent(inout) :: self  !! Swiftest massive body object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
+      class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
       ! Internals
-      integer(I4B) :: i, j, npl
+      integer(I4B) :: i, j, npl, err
       integer(I8B) :: k
 
       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 (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously
-            allocate(self%k_plpl(2, nplpl))
-            do concurrent (i=1:npl, j=1:npl, j>i)
-               call util_index_eucl_ij_to_k(npl, i, j, k)
-               self%k_plpl(1, k) = i
-               self%k_plpl(2, k) = j
-            end do
+            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.
+            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 if
          end if
       end associate
 
       return
-   end subroutine util_index_eucl_plpl
+   end subroutine util_flatten_eucl_plpl
 
 
-   module subroutine util_index_eucl_pltp(self, pl, param)
+   module subroutine util_flatten_eucl_pltp(self, pl, param)
       !! author: Jacob R. Elliott and David A. Minton
       !!
       !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-tp interactions
@@ -110,7 +114,7 @@ module subroutine util_index_eucl_pltp(self, pl, param)
       ! Arguments
       class(swiftest_tp),         intent(inout) :: self  !! Swiftest test particle object
       class(swiftest_pl),         intent(in)    :: pl    !! Swiftest massive body object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
+      class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
       ! Internals
       integer(I8B) :: i, j, counter, npl, ntp
 
@@ -131,6 +135,6 @@ module subroutine util_index_eucl_pltp(self, pl, param)
       end associate
 
       return
-   end subroutine util_index_eucl_pltp
+   end subroutine util_flatten_eucl_pltp
 
 end submodule s_util_index

src/whm/whm_setup.f90

@@ -82,7 +82,7 @@ module subroutine whm_setup_initialize_system(self, param)
       ! First we need to make sure that the massive bodies are sorted by heliocentric distance before computing jacobies
       call util_set_ir3h(self%pl)
       call self%pl%sort("ir3h", ascending=.false.)
-      call self%pl%index(param)
+      call self%pl%flatten(param)
 
       ! Make sure that the discard list gets allocated initially
       call self%tp_discards%setup(0, param)