Changes in a number of places to improve the robustness of the code

src/discard/discard.f90

@@ -236,13 +236,12 @@ subroutine discard_pl_tp(tp, system, param)
                      tp%status(i) = DISCARDED_PLR
                      tp%lmask(i) = .false.
                      pl%ldiscard(j) = .true.
-                     write(*, *) "Particle ", tp%id(i), " too close to massive body ", pl%id(j), " at t = ", t
                      write(idstri, *) tp%id(i)
                      write(idstrj, *) pl%id(j)
                      write(timestr, *) param%t
-                     write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstri)) // ")" &
-                                            // "  too close to massive body " // trim(adjustl(pl%info(i)%name)) // " ("  // trim(adjustl(idstrj)) &
-                                            // " at t = " // trim(adjustl(timestr))
+                     write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstri)) // ")" &
+                                                  // "  too close to massive body " // trim(adjustl(pl%info(i)%name)) // " ("  // trim(adjustl(idstrj)) &
+                                                  // " at t = " // trim(adjustl(timestr))
                      tp%ldiscard(i) = .true.
                      exit
                   end if

src/fragmentation/fragmentation.f90

@@ -152,12 +152,12 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       ! write(*,        "(' -------------------------------------------------------------------------------------')")
       ! call calculate_system_energy(linclude_fragments=.true.)
       if (lfailure) then
-         ! write(*,*) "symba_frag_pos failed after: ",try," tries"
+          write(*,*) "symba_frag_pos failed after: ",try," tries"
          do ii = 1, nfrag
             vb_frag(:, ii) = vcom(:)
          end do
-      ! else
-      !    write(*,*) "symba_frag_pos succeeded after: ",try," tries"
+       else
+          write(*,*) "symba_frag_pos succeeded after: ",try," tries"
       !    write(*,        "(' dL_tot should be very small' )")
       !    write(*,fmtlabel) ' dL_tot      |', dLmag / Lmag_before
       !    write(*,        "(' dE_tot should be negative and equal to Qloss' )")

src/rmvs/rmvs_discard.f90

@@ -16,6 +16,7 @@ module subroutine rmvs_discard_tp(self, system, param)
       class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
       ! Internals
       integer(I4B)                                 :: i
+      character(len=STRMAX) :: timestr, idstri, idstrj
 
       if (self%nbody == 0) return
 
@@ -25,7 +26,12 @@ module subroutine rmvs_discard_tp(self, system, param)
                if ((tp%status(i) == ACTIVE) .and. (tp%lperi(i))) then 
                   if ((tp%peri(i) < pl%radius(iplperP))) then
                      tp%status(i) = DISCARDED_PLQ
-                     write(*, *) "Particle ",tp%id(i)," q with respect to Planet ",pl%id(iplperP)," is too small at t = ",t
+                     write(idstri, *) tp%id(i)
+                     write(idstrj, *) pl%id(iplperP)
+                     write(timestr, *) t
+                     write(*, *) "Particle "  // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) &
+                              // ") q with respect to massive body " // trim(adjustl(pl%info(iplperP)%name)) // " (" // trim(adjustl(idstrj)) &
+                              // ") is too small at t = " // trim(adjustl(timestr))
                      tp%ldiscard(i) = .true.
                      tp%lmask(i) = .false.
                   end if

src/rmvs/rmvs_util.f90

@@ -189,16 +189,19 @@ module subroutine rmvs_util_sort_pl(self, sortby, ascending)
       character(*),   intent(in)     :: sortby    !! Sorting attribute
       logical,        intent(in)     :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order
       ! Internals
-      integer(I4B), dimension(self%nbody) :: ind
+      integer(I4B), dimension(:), allocatable :: ind
       integer(I4B) :: direction
 
+      if (self%nbody == 0) return
+
       if (ascending) then
          direction = 1
       else
          direction = -1
       end if
 
       associate(pl => self, npl => self%nbody)
+         allocate(ind(npl))
          select case(sortby)
          case("nenc")
             call util_sort(direction * pl%nenc(1:npl), ind(1:npl))
@@ -231,8 +234,10 @@ module subroutine rmvs_util_sort_tp(self, sortby, ascending)
       character(*),   intent(in)    :: sortby    !! Sorting attribute
       logical,        intent(in)    :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order
       ! Internals
-      integer(I4B), dimension(self%nbody) :: ind
-      integer(I4B)                        :: direction
+      integer(I4B), dimension(:), allocatable :: ind
+      integer(I4B)                            :: direction
+
+      if (self%nbody == 0) return
 
       if (ascending) then
          direction = 1
@@ -241,6 +246,7 @@ module subroutine rmvs_util_sort_tp(self, sortby, ascending)
       end if
 
       associate(tp => self, ntp => self%nbody)
+         allocate(ind(ntp))
          select case(sortby)
          case("plperP")
             call util_sort(direction * tp%plperP(1:ntp), ind(1:ntp))

src/symba/symba_collision.f90

@@ -740,7 +740,7 @@ module subroutine symba_collision_encounter_extract_collisions(self, system, par
       class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
       class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters
       ! Internals
-      logical,      dimension(self%nenc)      :: lplpl_collision
+      logical,      dimension(:), allocatable :: lplpl_collision
       logical,      dimension(:), allocatable :: lplpl_unique_parent
       integer(I4B), dimension(:), pointer     :: plparent
       integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx
@@ -749,6 +749,7 @@ module subroutine symba_collision_encounter_extract_collisions(self, system, par
       select type (pl => system%pl)
       class is (symba_pl)
          associate(plplenc_list => self, nplplenc => self%nenc, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent)
+            allocate(lplpl_collision(nplplenc))
             lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION
             if (.not.any(lplpl_collision)) return 
             ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies.
@@ -1139,6 +1140,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir
       ! Internals
       real(DP) :: Eorbit_before, Eorbit_after
       logical :: lplpl_collision
+      character(len=STRMAX) :: timestr
 
       associate(plplenc_list => self, plplcollision_list => system%plplcollision_list)
          select type(pl => system%pl)
@@ -1157,7 +1159,8 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir
                end if
 
                do
-                  write(*, *) "Collision between massive bodies detected at time t = ", t
+                  write(timestr,*) t
+                  write(*, *) "Collision between massive bodies detected at time t = " // trim(adjustl(timestr))
                   if (param%lfragmentation) then
                      call plplcollision_list%resolve_fragmentations(system, param)
                   else

src/symba/symba_discard.f90

@@ -36,14 +36,14 @@ subroutine symba_discard_cb_pl(pl, system, param)
                   pl%status(i) = DISCARDED_RMAX
                   write(idstr, *) pl%id(i)
                   write(timestr, *) param%t
-                  write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr))
+                  write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr))
                else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then
                   pl%ldiscard(i) = .true.
                   pl%lcollision(i) = .false. 
                   pl%status(i) = DISCARDED_RMIN
                   write(idstr, *) pl%id(i)
                   write(timestr, *) param%t
-                  write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr))
+                  write(*, *) trim(adjustl(pl%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr))
                else if (param%rmaxu >= 0.0_DP) then
                   rb2 = dot_product(pl%xb(:,i), pl%xb(:,i))
                   vb2 = dot_product(pl%vb(:,i), pl%vb(:,i))
@@ -54,7 +54,7 @@ subroutine symba_discard_cb_pl(pl, system, param)
                      pl%status(i) = DISCARDED_RMAXU
                      write(idstr, *) pl%id(i)
                      write(timestr, *) param%t
-                     write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
+                     write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
                   end if
                end if
             end if
@@ -255,6 +255,7 @@ subroutine symba_discard_peri_pl(pl, system, param)
       logical, save      :: lfirst = .true.
       logical            :: lfirst_orig
       integer(I4B)       :: i
+      character(len=STRMAX) :: timestr, idstr
 
 
       lfirst_orig = pl%lfirst
@@ -271,7 +272,9 @@ subroutine symba_discard_peri_pl(pl, system, param)
                      pl%ldiscard(i) = .true.
                      pl%lcollision(i) = .false.
                      pl%status(i) = DISCARDED_PERI
-                     write(*, *) "Particle ", pl%id(i), " perihelion distance too small at t = ", param%t
+                     write(timestr, *) param%t
+                     write(idstr, *) pl%id(i)
+                     write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ") perihelion distance too small at t = " // trim(adjustl(timestr)) 
                   end if
                end if
             end if

src/symba/symba_encounter_check.f90

@@ -99,7 +99,7 @@ module function symba_encounter_check(self, system, dt, irec) result(lany_encoun
       integer(I4B),              intent(in)    :: irec       !! Current recursion level 
       logical                                  :: lany_encounter !! Returns true if there is at least one close encounter  
       ! Internals
-      integer(I4B)              :: k
+      integer(I4B)              :: i, j,k
       real(DP), dimension(NDIM) :: xr, vr
       logical                   :: lencounter, isplpl
       real(DP)                  :: rlim2, rji2
@@ -122,40 +122,40 @@ module function symba_encounter_check(self, system, dt, irec) result(lany_encoun
             allocate(lencmask(self%nenc))
             lencmask(:) = (self%status(1:self%nenc) == ACTIVE) .and. (self%level(1:self%nenc) == irec - 1)
             if (.not.any(lencmask(:))) return
-            associate(ind1 => self%index1, ind2 => self%index2) 
-               do concurrent(k = 1:self%nenc, lencmask(k))
+            do concurrent(k = 1:self%nenc, lencmask(k))
+               i = self%index1(k)
+               j = self%index2(k)
+               if (isplpl) then
+                  xr(:) = pl%xh(:,j) - pl%xh(:,i)
+                  vr(:) = pl%vb(:,j) - pl%vb(:,i)
+                  call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), pl%rhill(j), dt, irec, lencounter, self%lvdotr(k))
+               else
+                  xr(:) = tp%xh(:,j) - pl%xh(:,i)
+                  vr(:) = tp%vb(:,j) - pl%vb(:,i)
+                  call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), 0.0_DP, dt, irec, lencounter, self%lvdotr(k))
+               end if
+               if (lencounter) then
                   if (isplpl) then
-                     xr(:) = pl%xh(:,ind2(k)) - pl%xh(:,ind1(k))
-                     vr(:) = pl%vb(:,ind2(k)) - pl%vb(:,ind1(k))
-                     call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(ind1(k)), pl%rhill(ind2(k)), dt, irec, lencounter, self%lvdotr(k))
+                     rlim2 = (pl%radius(i) + pl%radius(j))**2
                   else
-                     xr(:) = tp%xh(:,ind2(k)) - pl%xh(:,ind1(k))
-                     vr(:) = tp%vb(:,ind2(k)) - pl%vb(:,ind1(k))
-                     call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(ind1(k)), 0.0_DP, dt, irec, lencounter, self%lvdotr(k))
+                     rlim2 = (pl%radius(i))**2
                   end if
-                  if (lencounter) then
+                  rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore
+                  if (rji2 > rlim2) then
+                     lany_encounter = .true.
+                     pl%levelg(i) = irec
+                     pl%levelm(i) = MAX(irec, pl%levelm(i))
                      if (isplpl) then
-                        rlim2 = (pl%radius(ind1(k)) + pl%radius(ind2(k)))**2
+                        pl%levelg(j) = irec
+                        pl%levelm(j) = MAX(irec, pl%levelm(j))
                      else
-                        rlim2 = (pl%radius(ind1(k)))**2
+                        tp%levelg(j) = irec
+                        tp%levelm(j) = MAX(irec, tp%levelm(j))
                      end if
-                     rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore
-                     if (rji2 > rlim2) then
-                        lany_encounter = .true.
-                        pl%levelg(ind1(k)) = irec
-                        pl%levelm(ind1(k)) = MAX(irec, pl%levelm(ind1(k)))
-                        if (isplpl) then
-                           pl%levelg(ind2(k)) = irec
-                           pl%levelm(ind2(k)) = MAX(irec, pl%levelm(ind2(k)))
-                        else
-                           tp%levelg(ind2(k)) = irec
-                           tp%levelm(ind2(k)) = MAX(irec, tp%levelm(ind2(k)))
-                        end if
-                        self%level(k) = irec
-                     end if
-                  end if   
-               end do
-            end associate
+                     self%level(k) = irec
+                  end if
+               end if   
+            end do
          end select
       end select