Improved robustness of spill methods when keeps and lspill_list are d…

…ifferent sizes. Created new plplcollision_list structure in order to prevent resizing of the plplenc_list every time there is a collision

src/io/io.f90

@@ -574,7 +574,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
                iostat = -1
                return
             end if
-
          end if
          if (param%qmin > 0.0_DP) then
             if ((param%qmin_coord /= "HELIO") .and. (param%qmin_coord /= "BARY")) then
@@ -622,6 +621,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
          write(*,*) "MU2KG          = ",param%MU2KG       
          write(*,*) "TU2S           = ",param%TU2S        
          write(*,*) "DU2M           = ",param%DU2M        
+         if (trim(adjustl(param%outfile)) == "") then
+            param%outfile = BIN_OUTFILE
+         end if
          if (trim(adjustl(param%enc_out)) /= "") then
             write(*,*) "ENC_OUT        = ",trim(adjustl(param%enc_out))
          else
@@ -632,6 +634,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
             write(*,*) "BIG_DISCARD    = ",param%lbig_discard
          else
             write(*,*) "! DISCARD_OUT not set: Discards will not be recorded to file"
+            param%lbig_discard = .false.
             write(*,*) "! BIG_DISCARD    = ",param%lbig_discard
          end if
 
@@ -643,10 +646,11 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
 
          ! Calculate the G for the system units
          param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2))
-
-         ! Calculate the inverse speed of light in the system units
-         param%inv_c2 = einsteinC * param%TU2S / param%DU2M
-         param%inv_c2 = (param%inv_c2)**(-2)
+         if (param%lgr) then
+            ! Calculate the inverse speed of light in the system units
+            param%inv_c2 = einsteinC * param%TU2S / param%DU2M
+            param%inv_c2 = (param%inv_c2)**(-2)
+         end if
 
          associate(integrator => v_list(1))
             if (integrator == RMVS) then

src/modules/swiftest_globals.f90

@@ -107,9 +107,6 @@ module swiftest_globals
    character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.dat', 'dump_param2.dat']
 
    !> Default file names that can be changed by the user in the parameters file
-   character(*), parameter :: ENC_OUTFILE  = 'encounter.out'
-   character(*), parameter :: DISCARD_FILE = 'discard.out'
-   character(*), parameter :: ENERGY_FILE  = 'energy.out'
    character(*), parameter :: CB_INFILE    = 'cb.in'
    character(*), parameter :: PL_INFILE    = 'pl.in'
    character(*), parameter :: TP_INFILE    = 'tp.in'

src/modules/symba_classes.f90

@@ -149,7 +149,7 @@ module symba_classes
    !> SyMBA class for tracking pl-pl close encounters in a step
    type, extends(symba_pltpenc) :: symba_plplenc
    contains
-      procedure :: scrub_non_collision    => symba_collision_encounter_scrub        !! Processes the pl-pl encounter list remove only those encounters that led to a collision
+      procedure :: extract_collisions     => symba_collision_encounter_extract_collisions        !! Processes the pl-pl encounter list remove only those encounters that led to a collision
       procedure :: resolve_fragmentations => symba_collision_resolve_fragmentations !! Process list of collisions, determine the collisional regime, and then create fragments
       procedure :: resolve_mergers        => symba_collision_resolve_mergers        !! Process list of collisions and merge colliding bodies together
    end type symba_plplenc
@@ -161,6 +161,7 @@ module symba_classes
       class(symba_merger),  allocatable :: pl_adds !! List of added bodies in mergers or collisions
       class(symba_pltpenc), allocatable :: pltpenc_list  !! List of massive body-test particle encounters in a single step 
       class(symba_plplenc), allocatable :: plplenc_list  !! List of massive body-massive body encounters in a single step
+      class(symba_plplenc), allocatable :: plplcollision_list  !! List of massive body-massive body collisions in a single step
       integer(I4B)                      :: irec          !! System recursion level
    contains
       procedure :: write_discard    => symba_io_write_discard             !! Write out information about discarded and merged planets and test particles in SyMBA
@@ -184,11 +185,11 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec
          integer(I4B),               intent(in)    :: irec   !! Current recursion level
       end subroutine symba_collision_check_pltpenc
 
-      module subroutine symba_collision_encounter_scrub(self, system, param)
+      module subroutine symba_collision_encounter_extract_collisions(self, system, param)
          implicit none
          class(symba_plplenc),       intent(inout) :: self   !! SyMBA pl-pl encounter list
          class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
-         class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameterss
+         class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters
       end subroutine
 
       module subroutine symba_collision_make_family_pl(self,idx)

src/setup/setup.f90

@@ -56,8 +56,9 @@ module subroutine setup_construct_system(system, param)
             allocate(symba_tp :: system%tp_discards)
             allocate(symba_merger :: system%pl_adds)
             allocate(symba_merger :: system%pl_discards)
-            allocate(symba_plplenc :: system%plplenc_list)
             allocate(symba_pltpenc :: system%pltpenc_list)
+            allocate(symba_plplenc :: system%plplenc_list)
+            allocate(symba_plplenc :: system%plplcollision_list)
          end select
       case (RINGMOONS)
          write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled'

src/symba/symba_collision.f90

@@ -267,7 +267,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v
    end function symba_collision_consolidate_familes
 
 
-   module subroutine symba_collision_encounter_scrub(self, system, param)
+   module subroutine symba_collision_encounter_extract_collisions(self, system, param)
       !! author: David A. Minton
       !! 
       !! Processes the pl-pl encounter list remove only those encounters that led to a collision
@@ -283,56 +283,55 @@ module subroutine symba_collision_encounter_scrub(self, system, param)
       integer(I4B), dimension(:), pointer     :: plparent
       integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx
       integer(I4B)                            :: i, index_coll, ncollisions, nunique_parent
-      type(symba_plplenc)                     :: plplenc_noncollision
 
       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)
             lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION
-            if (any(lplpl_collision)) then ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies.
-
-               ! Get the subset of pl-pl encounters that lead to a collision
-               ncollisions = count(lplpl_collision(:))
-               allocate(collision_idx(ncollisions))
-               collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision)
-
-               ! Get the subset of collisions that involve a unique pair of parents
-               allocate(lplpl_unique_parent(ncollisions))
-
-               lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:)))
-               nunique_parent = count(lplpl_unique_parent(:))
-               allocate(unique_parent_idx(nunique_parent))
-               unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:))
-
-               ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind
-               ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases
-               ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single
-               ! step
-               lplpl_unique_parent(:) = .true.
-               do index_coll = 1, ncollisions
-                  associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll))))
-                     lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. &
-                                                      any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. &
-                                                      any(plparent(idx1(unique_parent_idx(:))) == ip2) .or. &
-                                                      any(plparent(idx2(unique_parent_idx(:))) == ip2) )
-                  end associate
-               end do
-
-               ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't
-               ! contain a parent body on the unique parent list.
-               ncollisions = nunique_parent + count(lplpl_unique_parent)
-               collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))]
-
-               ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them
-               lplpl_collision(:) = .false.
-               lplpl_collision(collision_idx(:)) = .true.
-            end if
-            call plplenc_list%spill(plplenc_noncollision, .not.lplpl_collision, ldestructive=.true.) ! Remove any encounters that are not collisions from the list.
+            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.
+
+            ! Get the subset of pl-pl encounters that lead to a collision
+            ncollisions = count(lplpl_collision(:))
+            allocate(collision_idx(ncollisions))
+            collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision)
+
+            ! Get the subset of collisions that involve a unique pair of parents
+            allocate(lplpl_unique_parent(ncollisions))
+
+            lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:)))
+            nunique_parent = count(lplpl_unique_parent(:))
+            allocate(unique_parent_idx(nunique_parent))
+            unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:))
+
+            ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind
+            ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases
+            ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single
+            ! step
+            lplpl_unique_parent(:) = .true.
+            do index_coll = 1, ncollisions
+               associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll))))
+                  lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. &
+                                                   any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. &
+                                                   any(plparent(idx1(unique_parent_idx(:))) == ip2) .or. &
+                                                   any(plparent(idx2(unique_parent_idx(:))) == ip2) )
+               end associate
+            end do
+
+            ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't
+            ! contain a parent body on the unique parent list.
+            ncollisions = nunique_parent + count(lplpl_unique_parent)
+            collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))]
+
+            ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them
+            lplpl_collision(:) = .false.
+            lplpl_collision(collision_idx(:)) = .true.
+            call plplenc_list%spill(system%plplcollision_list, lplpl_collision, ldestructive=.false.) ! Extract any encounters that are not collisions from the list.
          end associate
       end select
 
       return
-   end subroutine symba_collision_encounter_scrub
+   end subroutine symba_collision_encounter_extract_collisions
 
 
    module subroutine symba_collision_make_family_pl(self, idx)

src/symba/symba_discard.f90

@@ -291,21 +291,21 @@ module subroutine symba_discard_pl(self, system, param)
       class is (symba_nbody_system)
          select type(param)
          class is (symba_parameters)
-            associate(pl => self, plplenc_list => system%plplenc_list)
+            associate(pl => self, plplenc_list => system%plplenc_list, plplcollision_list => system%plplcollision_list)
                call pl%h2b(system%cb) 
 
                ! First deal with the non pl-pl collisions
                call symba_discard_nonplpl(self, system, param)
 
-               ! Scrub the pl-pl encounter list of any encounters that did not lead to a collision
-               call plplenc_list%scrub_non_collision(system, param)
+               ! Extract the pl-pl encounter list and return the plplcollision_list
+               call plplenc_list%extract_collisions(system, param)
 
-               if ((plplenc_list%nenc > 0) .and. any(pl%lcollision(:))) then 
+               if ((plplcollision_list%nenc > 0) .and. any(pl%lcollision(:))) then 
                   write(*, *) "Collision between massive bodies detected at time t = ",param%t
                   if (param%lfragmentation) then
-                     call plplenc_list%resolve_fragmentations(system, param)
+                     call plplcollision_list%resolve_fragmentations(system, param)
                   else
-                     call plplenc_list%resolve_mergers(system, param)
+                     call plplcollision_list%resolve_mergers(system, param)
                   end if
                end if
 

src/symba/symba_setup.f90

@@ -76,6 +76,7 @@ module subroutine symba_setup_initialize_system(self, param)
          call whm_setup_initialize_system(system, param)
          call system%pltpenc_list%setup(0)
          call system%plplenc_list%setup(0)
+         call system%plplcollision_list%setup(0)
          select type(pl => system%pl)
          class is (symba_pl)
             call pl%sort("mass", ascending=.false.)

src/symba/symba_step.f90

@@ -237,7 +237,7 @@ module subroutine symba_step_reset_system(self)
       ! Internals
       integer(I4B) :: i
 
-      associate(system => self, pltpenc_list => self%pltpenc_list, plplenc_list => self%plplenc_list, pl_adds => self%pl_adds, pl_discards => self%pl_discards)
+      associate(system => self)
          select type(pl => system%pl)
          class is (symba_pl)
             select type(tp => system%tp)
@@ -255,18 +255,19 @@ module subroutine symba_step_reset_system(self)
                   pl%levelm(:) = 0
                   pl%lencounter = .false.
                   pl%lcollision = .false.
-                  plplenc_list%nenc = 0
+                  system%plplenc_list%nenc = 0
+                  system%plplcollision_list%nenc = 0
                end if
 
                if (tp%nbody > 0) then
                   tp%nplenc(:) = 0 
                   tp%levelg(:) = 0
                   tp%levelm(:) = 0
-               pltpenc_list%nenc = 0
+               system%pltpenc_list%nenc = 0
                end if
 
-               call pl_adds%resize(0)
-               call pl_discards%resize(0)
+               call system%pl_adds%resize(0)
+               call system%pl_discards%resize(0)
             end select
          end select
       end associate

src/symba/symba_util.f90

@@ -704,14 +704,25 @@ module subroutine symba_util_spill_arr_info(keeps, discards, lspill_list, ldestr
       type(symba_particle_info), dimension(:), allocatable, intent(inout) :: discards     !! Array of discards
       logical,                   dimension(:),              intent(in)    :: lspill_list  !! Logical array of bodies to spill into the discardss
       logical,                                              intent(in)    :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not
+      ! Internals
+      integer(I4B) :: nspill, nkeep, nlist
+
+      nkeep = count(.not.lspill_list(:))
+      nspill = count(lspill_list(:))
+      nlist = size(lspill_list(:))
+
+      if (.not.allocated(keeps) .or. nspill == 0) return
+      if (.not.allocated(discards)) then
+         allocate(discards(nspill))
+      else if (size(discards) /= nspill) then
+         deallocate(discards)
+         allocate(discards(nspill))
+      end if
 
-      if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return
-      if (.not.allocated(discards)) allocate(discards(count(lspill_list(:))))
-
-      discards(:) = pack(keeps(:), lspill_list(:))
+      discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist))
       if (ldestructive) then
-         if (count(.not.lspill_list(:)) > 0) then
-            keeps(:) = pack(keeps(:), .not. lspill_list(:))
+         if (nkeep > 0) then
+            keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist))
          else
             deallocate(keeps)
          end if
@@ -732,14 +743,25 @@ module subroutine symba_util_spill_arr_kin(keeps, discards, lspill_list, ldestru
       type(symba_kinship), dimension(:), allocatable, intent(inout) :: discards     !! Array of discards
       logical,             dimension(:),              intent(in)    :: lspill_list  !! Logical array of bodies to spill into the discardss
       logical,                                        intent(in)    :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not
+      ! Internals
+      integer(I4B) :: nspill, nkeep, nlist
+
+      nkeep = count(.not.lspill_list(:))
+      nspill = count(lspill_list(:))
+      nlist = size(lspill_list(:))
+
+      if (.not.allocated(keeps) .or. nspill == 0) return
+      if (.not.allocated(discards)) then
+         allocate(discards(nspill))
+      else if (size(discards) /= nspill) then
+         deallocate(discards)
+         allocate(discards(nspill))
+      end if
 
-      if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return
-      if (.not.allocated(discards)) allocate(discards(count(lspill_list(:))))
-
-      discards(:) = pack(keeps(:), lspill_list(:))
+      discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist))
       if (ldestructive) then
-         if (count(.not.lspill_list(:)) > 0) then
-            keeps(:) = pack(keeps(:), .not. lspill_list(:))
+         if (nkeep > 0) then
+            keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist))
          else
             deallocate(keeps)
          end if
@@ -806,7 +828,7 @@ module subroutine symba_util_spill_pltpenc(self, discards, lspill_list, ldestruc
       ! Internals
       integer(I4B) :: i
 
-      associate(keeps => self)
+      associate(keeps => self, nenc => self%nenc)
          select type(discards)
          class is (symba_pltpenc)
             call util_spill(keeps%level, discards%level, lspill_list, ldestructive)

src/util/util_rescale.f90

@@ -20,9 +20,11 @@ module subroutine util_rescale_system(self, param, mscale, dscale, tscale)
       ! Calculate the G for the system units
       param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2))
 
-      ! Calculate the inverse speed of light in the system units
-      param%inv_c2 = einsteinC * param%TU2S / param%DU2M
-      param%inv_c2 = (param%inv_c2)**(-2)
+      if (param%lgr) then
+         ! Calculate the inverse speed of light in the system units
+         param%inv_c2 = einsteinC * param%TU2S / param%DU2M
+         param%inv_c2 = (param%inv_c2)**(-2)
+      end if
 
       vscale = dscale / tscale