From 8c8c8a0d400cf92b89055a05298bce81815f116b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 4 Aug 2021 13:24:28 -0400 Subject: [PATCH] Added in the family consolidation code --- src/symba/symba_collision.f90 | 118 ++++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 11175fe3a..13c874b0f 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -143,6 +143,124 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt end function symba_collision_check_one + function symba_collision_consolidate_familes(pl, idx_parent, family, x, v, mass, radius, L_spin, Ip) result(lflag) + !! author: David A. Minton + !! + !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all family members, + !! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome. + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl !! SyMBA massive body object + integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + integer(I4B), dimension(:), allocatable, intent(out) :: family !! List of indices of all bodies inovlved in the collision + real(DP), dimension(NDIM,2), intent(out) :: x, v, L_spin, Ip !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + real(DP), dimension(2), intent(out) :: mass, radius !! Output values that represent a 2-body equivalent of a possibly 2+ body collision + ! Result + logical :: lflag !! Logical flag indicating whether a family was successfully created or not + ! Internals + type family_array + integer(I4B), dimension(:), allocatable :: id + integer(I4B), dimension(:), allocatable :: idx + end type family_array + type(family_array), dimension(2) :: parent_child_index_array + integer(I4B), dimension(2) :: nchild + integer(I4B) :: i, j, fam_size, idx_child + real(DP), dimension(2) :: volume, density + real(DP) :: mchild, mtot, volchild + real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv + + nchild(:) = pl%kin(idx_parent(:))%nchild + ! If all of these bodies share a parent, but this is still a unique collision, move the last child + ! out of the parent's position and make it the secondary body + if (idx_parent(1) == idx_parent(2)) then + if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring of the kinship relationships, though it should be rare) + lflag = .false. + return + end if + idx_parent(2) = pl%kin(idx_parent(1))%child(nchild(1)) + nchild(1) = nchild(1) - 1 + nchild(2) = 0 + pl%kin(idx_parent(:))%nchild = nchild(:) + pl%kin(idx_parent(2))%parent = idx_parent(1) + end if + + mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si + radius(:) = pl%radius(idx_parent(:)) + volume(:) = (4.0_DP / 3.0_DP) * PI * radius(:)**3 + + ! Group together the ids and indexes of each collisional parent and its children + do j = 1, 2 + allocate(parent_child_index_array(j)%idx(nchild(j)+ 1)) + allocate(parent_child_index_array(j)%id(nchild(j)+ 1)) + associate(idx_arr => parent_child_index_array(j)%idx, & + id_arr => parent_child_index_array(j)%id, & + ncj => nchild(j), & + pl => pl, & + plkinj => pl%kin(idx_parent(j))) + idx_arr(1) = idx_parent(j) + if (ncj > 0) idx_arr(2:ncj + 1) = plkinj%child(1:ncj) + id_arr(:) = pl%id(idx_arr(:)) + end associate + end do + + ! Consolidate the groups of collsional parents with any children they may have into a single "family" index array + fam_size = 2 + sum(nchild(:)) + allocate(family(fam_size)) + family = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] + fam_size = count(pl%status(family(:)) == ACTIVE) + family = pack(family(:), pl%status(family(:)) == ACTIVE) + + ! Find the barycenter of each body along with its children, if it has any + 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) + L_spin(:, j) = Ip(3, j) * radius(j)**2 * pl%rot(:, idx_parent(j)) + + ! Use conservation of energy and angular momentum to adjust velocities and distances to be those equivalent to where they would be when the mutual radii are touching + !call vector_adjust(mass(j), x(:,j), v(:,j)) + if (nchild(j) > 0) then + do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties + idx_child = parent_child_index_array(j)%idx(i + 1) + if ((idx_child) /= COLLISION) cycle + mchild = pl%Gmass(idx_child) + xchild(:) = pl%xb(:, idx_child) + vchild(:) = pl%vb(:, idx_child) + volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 + volume(j) = volume(j) + volchild + ! Get angular momentum of the child-parent pair and add that to the spin + xcom(:) = (mass(j) * x(:,j) + mchild * xchild(:)) / (mass(j) + mchild) + vcom(:) = (mass(j) * v(:,j) + mchild * vchild(:)) / (mass(j) + mchild) + xc(:) = x(:, j) - xcom(:) + vc(:) = v(:, j) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + L_spin(:, j) = L_spin(:, j) + mass(j) * xcrossv(:) + + xc(:) = xchild(:) - xcom(:) + vc(:) = vchild(:) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + L_spin(:, j) = L_spin(:, j) + mchild * xcrossv(:) + + ! Add the child's spin + L_spin(:, j) = L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child) + + ! Merge the child and parent + mass(j) = mass(j) + mchild + x(:, j) = xcom(:) + v(:, j) = vcom(:) + Ip(:, j) = Ip(:, j) + mchild * pl%Ip(:, idx_child) + end do + end if + density(j) = mass(j) / volume(j) + radius(j) = ((3 * mass(j)) / (density(j) * 4 * pi))**(1.0_DP / 3.0_DP) + Ip(:, j) = Ip(:, j) / mass(j) + end do + lflag = .true. + + return + end function symba_collision_consolidate_familes + module subroutine symba_collision_encounter_scrub(self, system, param) !! author: David A. Minton !!