From 1f9d3c2400a6a145573e3c82de322bef9746c665 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 11 Aug 2021 10:51:00 -0400 Subject: [PATCH] Improved robustness of spill methods when keeps and lspill_list are different 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 | 14 ++-- src/modules/swiftest_globals.f90 | 3 - src/modules/symba_classes.f90 | 7 +- src/setup/setup.f90 | 3 +- src/symba/symba_collision.f90 | 83 ++++++++++---------- src/symba/symba_discard.f90 | 12 +-- src/symba/symba_setup.f90 | 1 + src/symba/symba_step.f90 | 11 +-- src/symba/symba_util.f90 | 48 ++++++++---- src/util/util_rescale.f90 | 8 +- src/util/util_spill.f90 | 125 +++++++++++++++++++++---------- 11 files changed, 195 insertions(+), 120 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 2fa42dcc1..10267c140 100644 --- a/src/io/io.f90 +++ b/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 diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 5ec55f6c6..6cacac789 100644 --- a/src/modules/swiftest_globals.f90 +++ b/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' diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 4628202f8..6a69adcc7 100644 --- a/src/modules/symba_classes.f90 +++ b/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) diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 8f96c48a1..843a48698 100644 --- a/src/setup/setup.f90 +++ b/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' diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 1910411b9..f1b023a53 100644 --- a/src/symba/symba_collision.f90 +++ b/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) diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index c154b4290..f003ab9cb 100644 --- a/src/symba/symba_discard.f90 +++ b/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 diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index e06fb20b5..e4644f25e 100644 --- a/src/symba/symba_setup.f90 +++ b/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.) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 7065625b4..847d6e4ee 100644 --- a/src/symba/symba_step.f90 +++ b/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 diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 4c0f256e3..1335272fe 100644 --- a/src/symba/symba_util.f90 +++ b/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) diff --git a/src/util/util_rescale.f90 b/src/util/util_rescale.f90 index 62a9409ec..50ebc088c 100644 --- a/src/util/util_rescale.f90 +++ b/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 diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 9c76ff5e9..8039d14b4 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -13,14 +13,25 @@ module subroutine util_spill_arr_char_string(keeps, discards, lspill_list, ldest character(len=STRMAX), dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards 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 @@ -40,14 +51,25 @@ module subroutine util_spill_arr_DP(keeps, discards, lspill_list, ldestructive) real(DP), 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 @@ -68,19 +90,30 @@ module subroutine util_spill_arr_DPvec(keeps, discards, lspill_list, ldestructiv logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not ! Internals - integer(I4B) :: i - - if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return - if (.not.allocated(discards)) allocate(discards(NDIM, count(lspill_list(:)))) + integer(I4B) :: i, 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(NDIM, nspill)) + else if (size(discards, dim=2) /= nspill) then + deallocate(discards) + allocate(discards(NDIM, nspill)) + end if do i = 1, NDIM - discards(i,:) = pack(keeps(i,:), lspill_list(:)) + discards(i,:) = pack(keeps(i,1:nlist), lspill_list(1:nlist)) end do if (ldestructive) then - if (count(.not.lspill_list(:)) > 0) then + if (nkeep > 0) then do i = 1, NDIM - keeps(i,:) = pack(keeps(i,:), .not. lspill_list(:)) + keeps(i,:) = pack(keeps(i,1:nlist), .not. lspill_list(1:nlist)) end do + else + deallocate(keeps) end if end if @@ -98,14 +131,25 @@ module subroutine util_spill_arr_I4B(keeps, discards, lspill_list, ldestructive) integer(I4B), dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards 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 @@ -125,14 +169,25 @@ module subroutine util_spill_arr_logical(keeps, discards, lspill_list, ldestruct logical, dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or no + ! 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 @@ -154,7 +209,6 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list ! Internals - integer(I4B) :: i ! For each component, pack the discarded bodies into the discard object and do the inverse with the keeps !> Spill all the common components @@ -185,7 +239,6 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) discards%nbody = count(lspill_list(:)) keeps%nbody = keeps%nbody - discards%nbody if (keeps%nbody > size(keeps%status)) keeps%status(keeps%nbody+1:size(keeps%status)) = INACTIVE - end associate return @@ -201,11 +254,8 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive class(swiftest_encounter), intent(inout) :: discards !! Discarded object logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list - ! Internals - integer(I4B) :: i associate(keeps => self) - call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive) call util_spill(keeps%status, discards%status, lspill_list, ldestructive) call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive) @@ -237,11 +287,8 @@ module subroutine util_spill_pl(self, discards, lspill_list, ldestructive) class(swiftest_body), intent(inout) :: discards !! Discarded object logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list - ! Internals - integer(I4B) :: i associate(keeps => self) - select type (discards) ! The standard requires us to select the type of both arguments in order to access all the components class is (swiftest_pl) !> Spill components specific to the massive body class