! Copyright 2024 - The Minton Group at Purdue University ! This file is part of Swiftest. ! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. ! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty ! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! You should have received a copy of the GNU General Public License along with Swiftest. ! If not, see: https://www.gnu.org/licenses. submodule (collision) s_collision_resolve use swiftest contains module subroutine collision_resolve_consolidate_impactors(self, nbody_system, param, idx_parent, lflag) !! author: David A. Minton !! !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all impactors%id !! members, and pairs of quantities (x and v vectors, mass, radius, L_rot, and Ip) that can be used to resolve the !! collisional outcome. implicit none ! Arguments class(collision_impactors),intent(out) :: self !! Collision impactors object class(base_nbody_system),intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(in) :: param !! Current run configuration parameters with Swiftest additions integer(I4B), dimension(:), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical, intent(out) :: lflag !! Logical flag indicating whether a impactors%id was successfully created or not ! Internals type collidx_array integer(I4B), dimension(:), allocatable :: id integer(I4B), dimension(:), allocatable :: idx end type collidx_array type(collidx_array), dimension(2) :: parent_child_index_array integer(I4B), dimension(2) :: nchild integer(I4B) :: i, j, nimpactors, idx_child real(DP), dimension(2) :: volume, density real(DP) :: mchild, volchild, rrel_mag, rlim, mtot, vdotr real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv, rrel, vrel, rrel_unit, vrel_unit, dr real(DP), dimension(NDIM,2) :: mxc, vcc select type(nbody_system) class is (swiftest_nbody_system) associate(impactors => self, pl => nbody_system%pl, cb => nbody_system%cb) 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. call pl%reset_kinship([idx_parent(1)]) 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 impactors%Gmass(:) = pl%Gmass(idx_parent(:)) impactors%mass(:) = pl%mass(idx_parent(:)) impactors%radius(:) = pl%radius(idx_parent(:)) volume(:) = (4.0_DP / 3.0_DP) * PI * impactors%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), & 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 "impactors%id" index array nimpactors = 2 + sum(nchild(:)) allocate(impactors%id(nimpactors)) impactors%id = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] impactors%ncoll = count(pl%lcollision(impactors%id(:))) impactors%id = pack(impactors%id(:), pl%lcollision(impactors%id(:))) impactors%L_rot(:,:) = 0.0_DP impactors%Ip(:,:) = 0.0_DP ! Find the barycenter of each body along with its children, if it has any do j = 1, 2 impactors%rb(:, j) = pl%rh(:, idx_parent(j)) + cb%rb(:) impactors%vb(:, j) = pl%vb(:, idx_parent(j)) ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) if (param%lrotation) then impactors%Ip(:, j) = impactors%mass(j) * pl%Ip(:, idx_parent(j)) impactors%L_rot(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) * DEG2RAD 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 idx_child = parent_child_index_array(j)%idx(i + 1) if (.not. pl%lcollision(idx_child)) cycle mchild = pl%mass(idx_child) xchild(:) = pl%rh(:, idx_child) + cb%rb(:) 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 rotation ! Add the child's rotation if (param%lrotation) then xcom(:) = (impactors%mass(j) * impactors%rb(:,j) + mchild * xchild(:)) / (impactors%mass(j) + mchild) vcom(:) = (impactors%mass(j) * impactors%vb(:,j) + mchild * vchild(:)) / (impactors%mass(j) + mchild) xc(:) = impactors%rb(:, j) - xcom(:) vc(:) = impactors%vb(:, j) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) impactors%L_rot(:, j) = impactors%L_rot(:, j) + impactors%mass(j) * xcrossv(:) xc(:) = xchild(:) - xcom(:) vc(:) = vchild(:) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) impactors%L_rot(:, j) = impactors%L_rot(:, j) + mchild * xcrossv(:) impactors%L_rot(:, j) = impactors%L_rot(:, j) + mchild * pl%Ip(3, idx_child) & * pl%radius(idx_child)**2 & * pl%rot(:, idx_child) * DEG2RAD impactors%Ip(:, j) = impactors%Ip(:, j) + mchild * pl%Ip(:, idx_child) end if ! Merge the child and parent impactors%mass(j) = impactors%mass(j) + mchild impactors%rb(:, j) = xcom(:) impactors%vb(:, j) = vcom(:) end do end if density(j) = impactors%mass(j) / volume(j) impactors%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) if (param%lrotation) then impactors%Ip(:, j) = impactors%Ip(:, j) / impactors%mass(j) impactors%rot(:,j) = impactors%L_rot(:, j) / (impactors%Ip(3,j) * impactors%mass(j) * impactors%radius(j)**2) & * RAD2DEG end if end do lflag = .true. mtot = sum(impactors%mass(:)) xcom(:) = (impactors%mass(1) * impactors%rb(:, 1) + impactors%mass(2) * impactors%rb(:, 2)) / mtot vcom(:) = (impactors%mass(1) * impactors%vb(:, 1) + impactors%mass(2) * impactors%vb(:, 2)) / mtot ! Shift the impactors so that they are not overlapping rlim = sum(impactors%radius(1:2)) rrel = impactors%rb(:,2) - impactors%rb(:,1) rrel_mag = .mag. rrel if (rrel_mag < rlim) then rrel_unit = .unit.rrel vrel = impactors%vb(:,2) - impactors%vb(:,1) vrel_unit = .unit.vrel vdotr = dot_product(vrel_unit, rrel) dr(:) = -(vdotr - sign(1.0_DP, vdotr) * sqrt(rlim**2 - rrel_mag**2 + vdotr**2)) * vrel_unit(:) dr(:) = (1.0_DP + 2*epsilon(1.0_DP)) * dr(:) impactors%rb(:,1) = impactors%rb(:,1) - dr(:) * impactors%mass(2) / mtot impactors%rb(:,2) = impactors%rb(:,2) + dr(:) * impactors%mass(1) / mtot rrel = impactors%rb(:,2) - impactors%rb(:,1) rrel_mag = .mag. rrel end if mxc(:, 1) = impactors%mass(1) * (impactors%rb(:, 1) - xcom(:)) mxc(:, 2) = impactors%mass(2) * (impactors%rb(:, 2) - xcom(:)) vcc(:, 1) = impactors%vb(:, 1) - vcom(:) vcc(:, 2) = impactors%vb(:, 2) - vcom(:) impactors%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) ! Destroy the kinship relationships for all members of this impactors%id call pl%reset_kinship(impactors%id(:)) end associate end select return end subroutine collision_resolve_consolidate_impactors module subroutine collision_resolve_extract_plpl(self, nbody_system, param) !! author: David A. Minton !! !! Processes the pl-pl encounter list remove only those encounters that led to a collision !! implicit none ! Arguments class(collision_list_plpl), intent(inout) :: self !! pl-pl encounter list class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(in) :: param !! Current run configuration parameters ! Internals logical, dimension(:), allocatable :: lplpl_collision logical, dimension(:), allocatable :: lplpl_unique_parent integer(I4B), dimension(:), pointer :: plparent integer(I4B) :: nunique_parent integer(I8B) :: ncollisions, index_coll, k, nplplenc integer(I8B), dimension(:), allocatable :: unique_parent_idx, collision_idx select type(nbody_system) class is (swiftest_nbody_system) select type (pl => nbody_system%pl) class is (swiftest_pl) associate(idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) nplplenc = self%nenc allocate(lplpl_collision(nplplenc)) lplpl_collision(:) = self%status(1_I8B:nplplenc) == COLLIDED 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([(k, k=1_I8B, 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 happen 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_I8B, 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. ! Extract any encounters that are not collisions from the list. call self%spill(nbody_system%plpl_collision, lplpl_collision, ldestructive=.true.) end associate end select end select return end subroutine collision_resolve_extract_plpl module subroutine collision_resolve_extract_pltp(self, nbody_system, param) implicit none class(collision_list_pltp), intent(inout) :: self !! pl-tp encounter list class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(in) :: param !! Current run configuration parameters ! Internals logical, dimension(:), allocatable :: lpltp_collision integer(I8B) :: ncollisions, k, npltpenc integer(I8B), dimension(:), allocatable :: collision_idx select type(nbody_system) class is (swiftest_nbody_system) select type (pl => nbody_system%pl) class is (swiftest_pl) select type (tp => nbody_system%tp) class is (swiftest_tp) associate(idx1 => self%index1, idx2 => self%index2) npltpenc = self%nenc allocate(lpltp_collision(npltpenc)) lpltp_collision(:) = self%status(1_I8B:npltpenc) == COLLIDED if (.not.any(lpltp_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-tp encounters that lead to a collision ncollisions = count(lpltp_collision(:)) allocate(collision_idx(ncollisions)) collision_idx = pack([(k, k=1_I8B, npltpenc)], lpltp_collision) ! Create a mask that contains only the pl-tp encounters that did not result in a collision, and then discard them lpltp_collision(:) = .false. lpltp_collision(collision_idx(:)) = .true. ! Extract any encounters that are not collisions from the list. call self%spill(nbody_system%pltp_collision, lpltp_collision, ldestructive=.true.) end associate end select end select end select return end subroutine collision_resolve_extract_pltp module subroutine collision_resolve_make_impactors_pl(pl, idx) !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton !! !! When a single body is involved in more than one collision in a single step, it becomes part of a collision family !! The largest body involved in a multi-body collision is the "parent" and all bodies that collide with it are its "children," !! including those that collide with the children. !! !! Adapted from David E. Kaufmann's Swifter routine swiftest_merge_pl.f90 !! !! Adapted from Hal Levison's Swift routine symba5_merge.f implicit none ! Arguments class(base_object), intent(inout) :: pl !! Swiftest massive body object integer(I4B), dimension(:), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision ! Internals integer(I4B) :: i, j, index_parent, index_child, p1, p2 integer(I4B) :: nchild_inherit, nchild_orig, nchild_new integer(I4B), dimension(:), allocatable :: temp select type(pl) class is (swiftest_pl) p1 = pl%kin(idx(1))%parent p2 = pl%kin(idx(2))%parent if (p1 == p2) return ! This is a collision between to children of a shared parent. We will ignore it. if (pl%mass(p1) > pl%mass(p2)) then index_parent = p1 index_child = p2 else index_parent = p2 index_child = p1 end if ! Expand the child array (or create it if necessary) and copy over the previous lists of children nchild_orig = pl%kin(index_parent)%nchild nchild_inherit = pl%kin(index_child)%nchild nchild_new = nchild_orig + nchild_inherit + 1 allocate(temp(nchild_new)) if (nchild_orig > 0) temp(1:nchild_orig) = pl%kin(index_parent)%child(1:nchild_orig) ! Find out if the child body has any children of its own. The new parent wil inherit these children if (nchild_inherit > 0) then temp(nchild_orig+1:nchild_orig+nchild_inherit) = pl%kin(index_child)%child(1:nchild_inherit) do i = 1, nchild_inherit j = pl%kin(index_child)%child(i) ! Set the childrens' parent to the new parent pl%kin(j)%parent = index_parent end do end if call pl%reset_kinship([index_child]) ! Add the new child to its parent pl%kin(index_child)%parent = index_parent temp(nchild_new) = index_child ! Save the new child array to the parent pl%kin(index_parent)%nchild = nchild_new call move_alloc(from=temp, to=pl%kin(index_parent)%child) end select return end subroutine collision_resolve_make_impactors_pl module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) !! author: David A. Minton !! !! Fills the pl_discards and pl_adds with removed and added bodies !! use symba, only : symba_pl implicit none ! Arguments class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions real(DP), intent(in) :: t !! Time of collision integer(I4B), intent(in) :: status !! Status flag to assign to adds ! Internals integer(I4B) :: i, ibiggest, ismallest, iother, nimpactors, nfrag logical, dimension(:), allocatable :: lmask class(swiftest_pl), allocatable :: plnew, plsub character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' character(len=NAMELEN) :: newname, origin_type real(DP) :: volume select type(nbody_system) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, & pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, collider => nbody_system%collider, & impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments) ! Add the impactors%id bodies to the subtraction list nimpactors = impactors%ncoll nfrag = fragments%nbody ! Setup new bodies allocate(plnew, mold=pl) call plnew%setup(nfrag, param) ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) ismallest = impactors%id(minloc(pl%Gmass(impactors%id(:)), dim=1)) ! Copy over identification, information, and physical properties of the new bodies from the fragment list plnew%id(1:nfrag) = fragments%id(1:nfrag) plnew%rb(:, 1:nfrag) = fragments%rb(:, 1:nfrag) plnew%vb(:, 1:nfrag) = fragments%vb(:, 1:nfrag) plnew%status(1:nfrag) = ACTIVE call plnew%b2h(cb) plnew%mass(1:nfrag) = fragments%mass(1:nfrag) plnew%Gmass(1:nfrag) = param%GU * fragments%mass(1:nfrag) plnew%radius(1:nfrag) = fragments%radius(1:nfrag) if (allocated(pl%a)) call plnew%xv2el(cb) if (param%lrotation) then plnew%Ip(:, 1:nfrag) = fragments%Ip(:, 1:nfrag) plnew%rot(:, 1:nfrag) = fragments%rot(:, 1:nfrag) end if call plnew%set_rhill(cb) ! if (param%ltides) then ! plnew%Q = pl%Q(ibiggest) ! plnew%k2 = pl%k2(ibiggest) ! plnew%tlag = pl%tlag(ibiggest) ! end if #ifdef DOCONLOC do concurrent(i = 1:nfrag) shared(plnew,fragments) local(volume) #else do concurrent(i = 1:nfrag) #endif volume = 4.0_DP/3.0_DP * PI * plnew%radius(i)**3 plnew%density(i) = fragments%mass(i) / volume end do select case(status) case(SUPERCATASTROPHIC) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) fragments%id(i) call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=t, name=newname, & origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=collider%maxid_collision) end do do i = 1, nimpactors if (impactors%id(i) == ibiggest) then iother = ismallest else iother = ibiggest end if call pl%info(impactors%id(i))%set_value(status="Supercatastrophic", discard_time=t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do case(DISRUPTED,HIT_AND_RUN_DISRUPT) if (status == DISRUPTED) then write(origin_type,*) "Disruption" else if (status == HIT_AND_RUN_DISRUPT) then write(origin_type,*) "Hit and run fragmentation" end if call plnew%info(1)%copy(pl%info(ibiggest)) do i = 2, nfrag write(newname, FRAGFMT) fragments%id(i) call plnew%info(i)%set_value(origin_type=origin_type, origin_time=t, name=newname, & origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), & collision_id=collider%maxid_collision) end do do i = 1, nimpactors if (impactors%id(i) == ibiggest) cycle iother = ibiggest call pl%info(impactors%id(i))%set_value(status=origin_type, discard_time=t, & discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do case(MERGED) call plnew%info(1)%copy(pl%info(ibiggest)) do i = 1, nimpactors if (impactors%id(i) == ibiggest) cycle iother = ibiggest call pl%info(impactors%id(i))%set_value(status="MERGED", discard_time=t, discard_rh=pl%rh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select !Copy over or set integration parameters for new bodies plnew%lcollision(1:nfrag) = .false. plnew%ldiscard(1:nfrag) = .false. select type(pl) class is (symba_pl) select type(plnew) class is (symba_pl) plnew%levelg(1:nfrag) = pl%levelg(ibiggest) plnew%levelm(1:nfrag) = pl%levelm(ibiggest) end select end select plnew%lmtiny(1:nfrag) = plnew%Gmass(1:nfrag) < param%GMTINY where(plnew%lmtiny(1:nfrag)) plnew%info(1:nfrag)%particle_type = PL_TINY_TYPE_NAME elsewhere plnew%info(1:nfrag)%particle_type = PL_TYPE_NAME end where ! Append the new merged body to the list call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)]) ! Add the discarded bodies to the discard list pl%status(impactors%id(:)) = MERGED pl%ldiscard(impactors%id(:)) = .true. pl%lcollision(impactors%id(:)) = .true. allocate(lmask, mold=pl%lmask) lmask(:) = .false. lmask(impactors%id(:)) = .true. allocate(plsub, mold=pl) call pl%spill(plsub, lmask, ldestructive=.false.) ! Save the before/after snapshots select type(before => collider%before) class is (swiftest_nbody_system) call move_alloc(plsub, before%pl) end select select type(after => collider%after) class is (swiftest_nbody_system) call move_alloc(plnew, after%pl) end select end associate end select end select return end subroutine collision_resolve_mergeaddsub module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) !! author: David A. Minton !! !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the collision !! implicit none ! Arguments class(collision_list_plpl), intent(inout) :: self !! Swiftest pl-pl encounter list class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Current simulation step size integer(I4B), intent(in) :: irec !! Current recursion level ! Internals real(DP) :: E_before, E_after, mnew logical :: lplpl_collision character(len=STRMAX) :: timestr integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision integer(I4B) :: nnew, loop, iframe_start, iframe_end integer(I8B) :: k, ncollisions integer(I4B), dimension(:), allocatable :: idnew integer(I4B), parameter :: MAXCASCADE = 1000 select type (nbody_system) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) class is (swiftest_pl) select type(param) class is (swiftest_parameters) associate(plpl_collision => nbody_system%plpl_collision, & collision_history => nbody_system%collision_history, & pl => nbody_system%pl, cb => nbody_system%cb, & collider => nbody_system%collider, & fragments => nbody_system%collider%fragments, & impactors => nbody_system%collider%impactors) if (plpl_collision%nenc == 0) return ! No collisions to resolve ! Make sure that the heliocentric and barycentric coordinates are consistent with each other call pl%vb2vh(nbody_system%cb) call pl%rh2rb(nbody_system%cb) ! Get the energy before the collision is resolved if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) E_before = nbody_system%te end if do loop = 1, MAXCASCADE associate( idx1 => plpl_collision%index1, idx2 => plpl_collision%index2) ncollisions = plpl_collision%nenc if (ncollisions == 0) exit write(timestr,*) t call swiftest_io_log_one_message(COLLISION_LOG_OUT, "") call swiftest_io_log_one_message(COLLISION_LOG_OUT,& "***********************************************************" // & "***********************************************************") call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision between massive bodies detected at time t = " // & trim(adjustl(timestr))) call swiftest_io_log_one_message(COLLISION_LOG_OUT, & "***********************************************************" // & "***********************************************************") iframe_start = collision_history%iframe + 1 do k = 1_I8B, ncollisions idx_parent(1) = pl%kin(idx1(k))%parent idx_parent(2) = pl%kin(idx2(k))%parent call impactors%consolidate(nbody_system, param, idx_parent, lgoodcollision) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle ! Get the collision regime call collider%get_regime(nbody_system, param) ! Temporarily store the original status of the parent bodies so we can compute the energy change correctly call collision_history%take_snapshot(param,nbody_system, t, "before") ! Generate the new bodies resulting from the collision call collider%generate(nbody_system, param, t) call collision_history%take_snapshot(param,nbody_system, t, "after") plpl_collision%status(k) = collider%status call impactors%dealloc() end do iframe_end = collision_history%iframe if (param%lenergy) call collision_history%save_energy_snapshot("before", nbody_system, iframe_start, iframe_end) ! Destroy the collision list now that the collisions are resolved call plpl_collision%setup(0_I8B) if ((nbody_system%pl_adds%nbody == 0) .and. (.not.any(pl%ldiscard(:)))) exit if (allocated(idnew)) deallocate(idnew) nnew = nbody_system%pl_adds%nbody allocate(idnew, source=nbody_system%pl_adds%id) mnew = sum(nbody_system%pl_adds%mass(:)) ! Rearrange the arrays: Remove discarded bodies, add any new bodies, re-sort, and recompute all indices and ! encounter lists call pl%rearray(nbody_system, param) ! Destroy the add/discard list so that we don't append the same body multiple times if another collision ! is detected call nbody_system%pl_discards%setup(0, param) call nbody_system%pl_adds%setup(0, param) if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) ! call nbody_system%get_energy_and_momentum(param) E_after = nbody_system%te nbody_system%E_collisions = nbody_system%E_collisions + (E_after - E_before) call collision_history%save_energy_snapshot("after", nbody_system, iframe_start, iframe_end) end if ! Check whether or not any of the particles that were just added are themselves in a collision state. This will ! generate a new plpl_collision call self%collision_check(nbody_system, param, t, dt, irec, lplpl_collision) if (.not.lplpl_collision) exit if (loop == MAXCASCADE) then call swiftest_io_log_one_message(COLLISION_LOG_OUT,"A runaway collisional cascade has been detected in " // & "collision_resolve_plpl.") call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Consider reducing the step size or changing the " // & "parameters in the collisional model to reduce the " // & "number of fragments.") call base_util_exit(FAILURE,unit=param%display_unit) end if end associate end do end associate end select end select end select return end subroutine collision_resolve_plpl module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) !! author: David A. Minton !! !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the collision implicit none ! Arguments class(collision_list_pltp), intent(inout) :: self !! Swiftest pl-pl encounter list class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Current simulation step size integer(I4B), intent(in) :: irec !! Current recursion level ! Internals character(len=STRMAX) :: timestr integer(I8B) :: k, ncollisions logical, dimension(:), allocatable :: ldiscard_pl, ldiscard_tp ! Make sure coordinate systems are all synced up due to being inside the recursion at this point select type(nbody_system) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) associate(pltp_collision => nbody_system%pltp_collision, & collision_history => nbody_system%collision_history, & pl => nbody_system%pl, & cb => nbody_system%cb, & tp => nbody_system%tp, & collider => nbody_system%collider, & impactors => nbody_system%collider%impactors) call pl%vb2vh(nbody_system%cb) call tp%vb2vh(nbody_system%cb%vb) call pl%b2h(nbody_system%cb) call tp%b2h(nbody_system%cb) ! Restructure the massive bodies based on the outcome of the collision call tp%rearray(nbody_system, param) ! Check for discards call nbody_system%tp%discard(nbody_system, param) associate(idx1 => pltp_collision%index1, idx2 => pltp_collision%index2) ncollisions = pltp_collision%nenc write(timestr,*) t do k = 1_I8B, ncollisions call swiftest_io_log_one_message(COLLISION_LOG_OUT, "") call swiftest_io_log_one_message(COLLISION_LOG_OUT,"***********************************************************" & // "***********************************************************") call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision between test particle and massive body detected " & // "at time t = " // trim(adjustl(timestr))) call swiftest_io_log_one_message(COLLISION_LOG_OUT,"***********************************************************" & // "***********************************************************") impactors%regime = COLLRESOLVE_REGIME_MERGE allocate(ldiscard_pl, mold=pl%ldiscard) ldiscard_pl(:) = .false. ldiscard_pl(idx1(k)) = .true. allocate(ldiscard_tp, mold=tp%ldiscard) ldiscard_tp(:) = .false. ldiscard_tp(idx2(k)) = .true. ! Save the system snapshot call tp%save_discard(ldiscard_tp,nbody_system,collider%before) call pl%save_discard(ldiscard_pl,nbody_system,collider%before) call collision_history%take_snapshot(param,nbody_system, t, "before") call pl%save_discard(ldiscard_pl,nbody_system,collider%after) call collision_history%take_snapshot(param,nbody_system, t, "after") deallocate(ldiscard_pl,ldiscard_tp) end do ! Destroy the collision list now that the collisions are resolved call pltp_collision%setup(0_I8B) end associate end associate end select end select return end subroutine collision_resolve_pltp end submodule s_collision_resolve