Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed some issues with masks in the encounter to collision scrub
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 4, 2021
1 parent 96e3170 commit 8ea46e1
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 47 deletions.
84 changes: 43 additions & 41 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
2 changes: 2 additions & 0 deletions src/symba/symba_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
9 changes: 3 additions & 6 deletions src/util/util_spill.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 8ea46e1

Please sign in to comment.