diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index edb02dee5..efae0ba39 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -217,9 +217,11 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v do j = 1, 2 x(:, j) = pl%xb(:, idx_parent(j)) v(:, j) = pl%vb(:, idx_parent(j)) - Ip(:, j) = mass(j) * pl%Ip(:, idx_parent(j)) ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) - if (param%lrotation) L_spin(:, j) = Ip(3, j) * radius(j)**2 * pl%rot(:, idx_parent(j)) + if (param%lrotation) then + Ip(:, j) = mass(j) * pl%Ip(:, idx_parent(j)) + L_spin(:, j) = Ip(3, j) * radius(j)**2 * pl%rot(:, idx_parent(j)) + end if if (nchild(j) > 0) then do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties @@ -282,49 +284,49 @@ module subroutine symba_collision_encounter_scrub(self, system, param) 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 (.not.any(lplpl_collision)) return - - ! 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(:) = .true. - lplpl_collision(collision_idx(:)) = .false. - call plplenc_list%spill(plplenc_noncollision, lplpl_collision, ldestructive = .true.) + 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. end associate end select diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 57b26e8ea..2fa789d46 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -152,6 +152,8 @@ module subroutine symba_discard_pl(self, system, param) associate(pl => self, plplenc_list => system%plplenc_list) call symba_discard_nonplpl(self, system, param) call plplenc_list%scrub_non_collision(system, param) + if (plplenc_list%nenc == 0) return ! No collisions to resolve + call pl%h2b(system%cb) if (param%lfragmentation) then call plplenc_list%resolve_fragmentations(system, param) diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index fc945765c..0cef49194 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -183,12 +183,8 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) ! Therefore we need to set the nbody values for both the keeps and discareds discards%nbody = count(lspill_list(:)) keeps%nbody = count(.not.lspill_list(:)) - if (allocated(keeps%ldiscard)) deallocate(keeps%ldiscard) - if (allocated(discards%ldiscard)) deallocate(discards%ldiscard) - allocate(keeps%ldiscard(keeps%nbody)) - allocate(discards%ldiscard(discards%nbody)) - keeps%ldiscard = .false. - discards%ldiscard = .true. + if (keeps%nbody > size(keeps%status)) keeps%status(keeps%nbody+1:size(keeps%status)) = INACTIVE + end associate return @@ -222,6 +218,7 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive ! Therefore we need to set the nenc values for both the keeps and discareds discards%nenc = count(lspill_list(:)) keeps%nenc = count(.not.lspill_list(:)) + if (keeps%nenc > size(keeps%status)) keeps%status(keeps%nenc+1:size(keeps%status)) = INACTIVE end associate return