From a5ac1c485ef7683436acc8f330172e47018b6885 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 22 Dec 2022 10:30:39 -0500 Subject: [PATCH] Restructured to try to simplify the code used to resolve collisions --- src/collision/collision_generate.f90 | 53 ++-- src/collision/collision_module.f90 | 73 ++--- src/collision/collision_resolve.f90 | 388 +++++++++++++-------------- src/swiftest/swiftest_setup.f90 | 2 +- 4 files changed, 233 insertions(+), 283 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 281145797..d89353e44 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -13,22 +13,6 @@ contains module subroutine collision_generate_merge_system(self, nbody_system, param, t) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Merge massive bodies in a "MERGE" style collision model (only pure mergers) - implicit none - class(collision_merge), intent(inout) :: self !! Merge fragment system object - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! The time of the collision - - call collision_generate_merge_any(self, nbody_system, param, t) - - return - end subroutine collision_generate_merge_system - - - module subroutine collision_generate_merge_any(self, nbody_system, param, t) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Merge massive bodies in any collisionals ystem. @@ -45,7 +29,7 @@ module subroutine collision_generate_merge_any(self, nbody_system, param, t) ! Internals integer(I4B) :: i, j, k, ibiggest real(DP), dimension(NDIM) :: Lspin_new - real(DP) :: dpe + real(DP) :: volume, dpe character(len=STRMAX) :: message select type(nbody_system) @@ -57,35 +41,46 @@ module subroutine collision_generate_merge_any(self, nbody_system, param, t) select type(pl => nbody_system%pl) class is (swiftest_pl) - - !call self%set_mass_dist(param) + ! Generate the merged body as a single fragment + call self%setup_fragments(1) ! Calculate the initial energy of the nbody_system without the collisional family call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) + ! The new body's metadata will be taken from the largest of the two impactor bodies, so we need + ! its index in the main pl structure ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) fragments%id(1) = pl%id(ibiggest) - fragments%rb(:,1) = impactors%rbcom(:) - fragments%vb(:,1) = impactors%vbcom(:) + allocate(fragments%info, source=pl%info(ibiggest:ibiggest)) + ! Compute the physical properties of the new body after the merge. + volume = 4._DP / 3._DP * PI * sum(impactors%radius(:)**3) + fragments%mass(1) = impactors%mass_dist(1) + fragments%density(1) = fragments%mass(1) / volume + fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD) if (param%lrotation) then - ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - Lspin_new(:) = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2) - - ! Assume prinicpal axis rotation on 3rd Ip axis + do concurrent(i = 1:NDIM) + fragments%Ip(i,1) = sum(impactors%mass(:) * impactors%Ip(i,:)) + Lspin_new(i) = sum(impactors%Lorbit(i,:) + impactors%Lorbit(i,:)) + end do + fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1) fragments%rot(:,1) = Lspin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2) else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable nbody_system%Lescape(:) = nbody_system%Lescape(:) + impactors%Lorbit(:,1) + impactors%Lorbit(:,2) end if - ! Keep track of the component of potential energy due to the pre-impact impactors%id for book-keeping + ! The fragment trajectory will be the barycentric trajectory + fragments%rb(:,1) = impactors%rbcom(:) + fragments%vb(:,1) = impactors%vbcom(:) + ! Get the energy of the system after the collision call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - dpe = self%pe(2) - self%pe(1) + + ! Keep track of the component of potential energy that is now not considered because two bodies became one + dpe = self%pe(2) - self%pe(1) nbody_system%Ecollisions = nbody_system%Ecollisions - dpe nbody_system%Euntracked = nbody_system%Euntracked + dpe - ! Update any encounter lists that have the removed bodies in them so that they instead point to the new body do k = 1, nbody_system%plpl_encounter%nenc do j = 1, impactors%ncoll @@ -111,7 +106,7 @@ module subroutine collision_generate_merge_any(self, nbody_system, param, t) end associate end select return - end subroutine collision_generate_merge_any + end subroutine collision_generate_merge_system module subroutine collision_generate_bounce_system(self, nbody_system, param, t) diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 1eee878b0..5e1166ea3 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -75,9 +75,10 @@ module collision real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider nbody_system in nbody_system barycentric coordinates contains - procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be - procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions - final :: collision_final_impactors !! Finalizer will deallocate all allocatables + procedure :: consolidate => collision_resolve_consolidate_impactors !! Consolidates a multi-body collision into an equivalent 2-body collision + procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be + procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions + final :: collision_final_impactors !! Finalizer will deallocate all allocatables end type collision_impactors @@ -109,9 +110,11 @@ module collision end type collision_fragments - type, abstract :: collision_system + type :: collision_system !! This class defines a collisional nbody_system that stores impactors and fragments. This is written so that various collision models (i.e. Fraggle) could potentially be used !! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments + !! + !! The generate method for this class is the merge model. This allows any extended type to have access to the merge procedure by selecting the collision_system parent class class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision @@ -135,15 +138,9 @@ module collision procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: reset => collision_util_reset_system !! Deallocates all allocatables procedure :: set_coordinate_system => collision_util_set_coordinate_system !! Sets the coordinate nbody_system of the collisional nbody_system - procedure(abstract_generate_system), deferred :: generate !! Generates a nbody_system of fragments depending on collision model + procedure :: generate => collision_generate_merge_system !! Merges the impactors to make a single final body end type collision_system - type, extends(collision_system) :: collision_merge - contains - procedure :: generate => collision_generate_merge_system !! Merges the impactors to make a single final body - final :: collision_final_merge_system !! Finalizer will deallocate all allocatables - end type collision_merge - type, extends(collision_system) :: collision_bounce contains procedure :: generate => collision_generate_bounce_system !! If a collision would result in a disruption, "bounce" the bodies instead. @@ -199,37 +196,10 @@ module collision end type collision_storage - abstract interface - subroutine abstract_generate_system(self, nbody_system, param, t) - import collision_system, base_nbody_system, base_parameters, DP - implicit none - class(collision_system), intent(inout) :: self !! Collision system object - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! The time of the collision - end subroutine abstract_generate_system - - subroutine abstract_set_mass_dist(self, param) - import collision_system, base_parameters - implicit none - class(collision_system), intent(inout) :: self !! Collision system object - class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - end subroutine abstract_set_mass_dist - end interface - - interface - module subroutine collision_generate_merge_any(self, nbody_system, param, t) - implicit none - class(collision_system), intent(inout) :: self !! Merge fragment nbody_system object - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! The time of the collision - end subroutine collision_generate_merge_any - module subroutine collision_generate_merge_system(self, nbody_system, param, t) implicit none - class(collision_merge), intent(inout) :: self !! Merge fragment nbody_system object + class(collision_system), intent(inout) :: self !! Merge fragment nbody_system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! The time of the collision @@ -310,6 +280,15 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l integer(I4B), intent(in) :: irec !! Current recursion level logical, intent(out) :: lany_collision !! Returns true if any pair of encounters resulted in a collision end subroutine collision_check_pltp + + module subroutine collision_resolve_consolidate_impactors(self, nbody_system, param, idx_parent, lflag) + implicit none + 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 + end subroutine collision_resolve_consolidate_impactors module subroutine collision_resolve_extract_plpl(self, nbody_system, param) implicit none @@ -544,22 +523,6 @@ subroutine collision_final_storage(self) end subroutine collision_final_storage - subroutine collision_final_merge_system(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(collision_merge), intent(inout) :: self !! Collision system object - - call self%reset() - if (allocated(self%impactors)) deallocate(self%impactors) - if (allocated(self%fragments)) deallocate(self%fragments) - - return - end subroutine collision_final_merge_system - - subroutine collision_final_bounce_system(self) !! author: David A. Minton !! diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 6cb08a49a..c04da8a62 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -12,20 +12,18 @@ contains - function collision_resolve_consolidate_impactors(pl, cb, param, idx_parent, impactors) result(lflag) + 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, Lspin, and Ip) that can be used to resolve the collisional outcome. implicit none ! Arguments - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - class(base_parameters), intent(in) :: param !! Current run configuration parameters with Swiftest additions - integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - class(collision_impactors), intent(out) :: impactors - ! Result - logical :: lflag !! Logical flag indicating whether a impactors%id was successfully created or not + 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 @@ -39,115 +37,121 @@ function collision_resolve_consolidate_impactors(pl, cb, param, idx_parent, impa real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv real(DP), dimension(NDIM,2) :: mxc, vcc - 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%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 - 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%Lspin(:,:) = 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%Lspin(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) - end if + 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%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 + 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 - 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 spin - ! Add the child's spin + ! 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%Lspin(:,:) = 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 - 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%Lspin(:, j) = impactors%Lspin(:, j) + impactors%mass(j) * xcrossv(:) - - xc(:) = xchild(:) - xcom(:) - vc(:) = vchild(:) - vcom(:) - xcrossv(:) = xc(:) .cross. vc(:) - impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * xcrossv(:) - - impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * pl%Ip(3, idx_child) & - * pl%radius(idx_child)**2 & - * pl%rot(:, idx_child) - impactors%Ip(:, j) = impactors%Ip(:, j) + mchild * pl%Ip(:, idx_child) + impactors%Ip(:, j) = impactors%mass(j) * pl%Ip(:, idx_parent(j)) + impactors%Lspin(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) end if - ! Merge the child and parent - impactors%mass(j) = impactors%mass(j) + mchild - impactors%rb(:, j) = xcom(:) - impactors%vb(:, j) = vcom(:) + 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 spin + ! Add the child's spin + 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%Lspin(:, j) = impactors%Lspin(:, j) + impactors%mass(j) * xcrossv(:) + + xc(:) = xchild(:) - xcom(:) + vc(:) = vchild(:) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * xcrossv(:) + + impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * pl%Ip(3, idx_child) & + * pl%radius(idx_child)**2 & + * pl%rot(:, idx_child) + 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) impactors%Ip(:, j) = impactors%Ip(:, j) / impactors%mass(j) 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) impactors%Ip(:, j) = impactors%Ip(:, j) / impactors%mass(j) - end do - lflag = .true. - - xcom(:) = (impactors%mass(1) * impactors%rb(:, 1) + impactors%mass(2) * impactors%rb(:, 2)) / sum(impactors%mass(:)) - vcom(:) = (impactors%mass(1) * impactors%vb(:, 1) + impactors%mass(2) * impactors%vb(:, 2)) / sum(impactors%mass(:)) - 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%Lorbit(:,:) = mxc(:,:) .cross. vcc(:,:) - - ! Destroy the kinship relationships for all members of this impactors%id - call pl%reset_kinship(impactors%id(:)) + lflag = .true. + + xcom(:) = (impactors%mass(1) * impactors%rb(:, 1) + impactors%mass(2) * impactors%rb(:, 2)) / sum(impactors%mass(:)) + vcom(:) = (impactors%mass(1) * impactors%vb(:, 1) + impactors%mass(2) * impactors%vb(:, 2)) / sum(impactors%mass(:)) + 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%Lorbit(:,:) = 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 function collision_resolve_consolidate_impactors + end subroutine collision_resolve_consolidate_impactors module subroutine collision_resolve_extract_plpl(self, nbody_system, param) @@ -232,48 +236,6 @@ module subroutine collision_resolve_extract_pltp(self, nbody_system, param) end subroutine collision_resolve_extract_pltp - subroutine collision_resolve_list(plpl_collision, nbody_system, param, t) - !! author: David A. Minton - !! - !! Process list of collisions, determine the collisional regime, and then create fragments. - !! - implicit none - ! Arguments - class(collision_list_plpl), intent(inout) :: plpl_collision !! 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 !! Time of collision - ! Internals - ! Internals - integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - logical :: lgoodcollision - integer(I4B) :: i - - select type(nbody_system) - class is (swiftest_nbody_system) - associate(ncollisions => plpl_collision%nenc, idx1 => plpl_collision%index1, idx2 => plpl_collision%index2, & - collision_history => nbody_system%collision_history, impactors => nbody_system%collider%impactors, & - fragments => nbody_system%collider%fragments, & - pl => nbody_system%pl, cb => nbody_system%cb) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = collision_resolve_consolidate_impactors(pl, cb, param, idx_parent, impactors) - if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle - - call impactors%get_regime(nbody_system, param) - call collision_history%take_snapshot(param,nbody_system, t, "before") - call nbody_system%collider%generate(nbody_system, param, t) - call collision_history%take_snapshot(param,nbody_system, t, "after") - call impactors%reset() - - end do - end associate - end select - return - end subroutine collision_resolve_list - - module subroutine collision_resolve_make_impactors_pl(pl, idx) !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton !! @@ -532,6 +494,10 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) logical :: lplpl_collision character(len=STRMAX) :: timestr class(swiftest_parameters), allocatable :: tmp_param + integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + logical :: lgoodcollision + integer(I4B) :: i, loop, ncollisions + integer(I4B), parameter :: MAXCASCADE = 1000 select type (nbody_system) class is (swiftest_nbody_system) @@ -539,59 +505,85 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) class is (swiftest_pl) select type(param) class is (swiftest_parameters) - associate(plpl_encounter => self, plpl_collision => nbody_system%plpl_collision) - 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) - Eorbit_before = nbody_system%te - end if - - do - 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, "***********************************************************" // & - "***********************************************************") - allocate(tmp_param, source=param) - - call collision_resolve_list(plpl_collision, nbody_system, param, t) + associate(plpl_encounter => self, 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) + associate( idx1 => plpl_collision%index1, idx2 => plpl_collision%index2) - ! Destroy the collision list now that the collisions are resolved - call plpl_collision%setup(0_I8B) + if (plpl_collision%nenc == 0) return ! No collisions to resolve - if ((nbody_system%pl_adds%nbody == 0) .and. (nbody_system%pl_discards%nbody == 0)) exit - ! Save the add/discard information to file - call nbody_system%write_discard(tmp_param) + ! 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) - ! Rearrange the arrays: Remove discarded bodies, add any new bodies, resort, and recompute all indices and encounter lists - call pl%rearray(nbody_system, tmp_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) - deallocate(tmp_param) - - ! 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 plpl_encounter%collision_check(nbody_system, param, t, dt, irec, lplpl_collision) + ! Get the energy before the collision is resolved + if (param%lenergy) then + call nbody_system%get_energy_and_momentum(param) + Eorbit_before = nbody_system%te + end if - if (.not.lplpl_collision) exit - end do + do loop = 1, MAXCASCADE + ncollisions = plpl_collision%nenc + 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, "***********************************************************" // & + "***********************************************************") + allocate(tmp_param, source=param) + + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + call impactors%consolidate(nbody_system, param, idx_parent, lgoodcollision) + if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle + + call impactors%get_regime(nbody_system, param) + call collision_history%take_snapshot(param,nbody_system, t, "before") + call collider%generate(nbody_system, param, t) + call collision_history%take_snapshot(param,nbody_system, t, "after") + plpl_collision%status(i) = collider%status + call impactors%reset() + end do + + ! Destroy the collision list now that the collisions are resolved + call plpl_collision%setup(0_I8B) + + if ((nbody_system%pl_adds%nbody == 0) .and. (nbody_system%pl_discards%nbody == 0)) exit + + ! Save the add/discard information to file + call nbody_system%write_discard(tmp_param) + + ! Rearrange the arrays: Remove discarded bodies, add any new bodies, resort, and recompute all indices and encounter lists + call pl%rearray(nbody_system, tmp_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) + deallocate(tmp_param) + + ! 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 plpl_encounter%collision_check(nbody_system, param, t, dt, irec, lplpl_collision) + + if (.not.lplpl_collision) exit + if (loop == MAXCASCADE) then + write(*,*) + write(*,*) "An runaway collisional cascade has been detected in collision_resolve_plpl." + write(*,*) "Consider reducing the step size or changing the parameters in the collisional model to reduce the number of fragments." + call util_exit(FAILURE) + end if + end do - if (param%lenergy) then - call nbody_system%get_energy_and_momentum(param) - Eorbit_after = nbody_system%te - nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before) - end if + if (param%lenergy) then + call nbody_system%get_energy_and_momentum(param) + Eorbit_after = nbody_system%te + nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before) + end if + end associate end associate end select end select diff --git a/src/swiftest/swiftest_setup.f90 b/src/swiftest/swiftest_setup.f90 index a08daf29d..ffa455ac0 100644 --- a/src/swiftest/swiftest_setup.f90 +++ b/src/swiftest/swiftest_setup.f90 @@ -117,7 +117,7 @@ module subroutine swiftest_setup_construct_system(nbody_system, param) select case(param%collision_model) case("MERGE") - allocate(collision_merge :: nbody_system%collider) + allocate(collision_system :: nbody_system%collider) case("BOUNCE") allocate(collision_bounce :: nbody_system%collider) case("SIMPLE")