From 2f9696d362fc1e6ec0b5d1b180c6fe1937e1b276 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 21 Feb 2024 14:20:29 -0500 Subject: [PATCH] Added a rearray_tp method to the swiftest_tp class in order to better resolve pl-tp collisions --- src/collision/collision_resolve.f90 | 22 +++++--- src/swiftest/swiftest_module.f90 | 14 +++-- src/swiftest/swiftest_util.f90 | 81 ++++++++++++++++++++++++++++- 3 files changed, 105 insertions(+), 12 deletions(-) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index ab27b72e4..1586073cb 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -220,7 +220,7 @@ module subroutine collision_resolve_extract_plpl(self, nbody_system, param) 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 + ! 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. @@ -713,13 +713,19 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - call nbody_system%pl%vb2vh(nbody_system%cb) - call nbody_system%tp%vb2vh(nbody_system%cb%vb) - call nbody_system%pl%b2h(nbody_system%cb) - call nbody_system%tp%b2h(nbody_system%cb) - - ! Discard the collider - call nbody_system%tp%discard(nbody_system, param) + 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) + + call tp%rearray(nbody_system, param) + + ! Discard the collider + !call nbody_system%tp%discard(nbody_system, param) + end associate end select end select diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 70c49a619..6350a42dd 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -321,6 +321,7 @@ module swiftest procedure :: rh2rb => swiftest_util_coord_rh2rb_tp !! Convert test particles from heliocentric to barycentric coordinates (position only) procedure :: dealloc => swiftest_util_dealloc_tp !! Deallocates all allocatable arrays procedure :: fill => swiftest_util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) + procedure :: rearray => swiftest_util_rearray_tp !! Clean up the test particle structures to remove discarded bodies procedure :: resize => swiftest_util_resize_tp !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. procedure :: set_mu => swiftest_util_set_mu_tp !! Method used to construct the vectorized form of the central body mass procedure :: sort => swiftest_util_sort_tp !! Sorts body arrays by a sortable component @@ -1485,11 +1486,18 @@ end subroutine swiftest_util_peri_tp module subroutine swiftest_util_rearray_pl(self, nbody_system, param) implicit none - class(swiftest_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine swiftest_util_rearray_pl + module subroutine swiftest_util_rearray_tp(self, nbody_system, param) + implicit none + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine swiftest_util_rearray_tp + module subroutine swiftest_util_rescale_system(self, param, mscale, dscale, tscale) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 78bcb7b01..5723db82f 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1704,10 +1704,13 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) if (npl == 0) then if (param%lmtiny_pl) pl%nplm = 0 + ! There are no more massive bodies. Reset the encounter lists and move on + if (allocated(nbody_system%plpl_encounter)) call nbody_system%plpl_encounter%setup(0_I8B) + if (allocated(nbody_system%pltp_encounter)) call nbody_system%pltp_encounter%setup(0_I8B) return end if - ! Reset all of the status flags for this body + ! Reset all of the status flags for the remaining bodies pl%status(1:npl) = ACTIVE do i = 1, npl call pl%info(i)%set_value(status="ACTIVE") @@ -1842,6 +1845,82 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) end subroutine swiftest_util_rearray_pl + + module subroutine swiftest_util_rearray_tp(self, nbody_system, param) + !! Author: David A. Minton + !! + !! Clean up the test particle structures to remove discarded bodies + use symba + implicit none + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + class(swiftest_tp), allocatable :: tmp !! The discarded body list. + integer(I4B) :: i, ntp + integer(I8B) :: k, nenc_old, nencmin + logical, dimension(:), allocatable :: lmask + class(encounter_list), allocatable :: pltpenc_old + logical :: lencounter + + associate(pl => self, tp => nbody_system%tp, cb => nbody_system%cb, pl_adds => nbody_system%pl_adds) + + ntp = tp%nbody + if (ntp == 0) return + + ! Remove the discards and destroy the list, as the nbody_system already tracks tp_discards elsewhere + allocate(lmask(ntp)) + lmask(1:ntp) = tp%ldiscard(1:ntp) + if (count(lmask(:)) > 0) then + allocate(tmp, mold=self) + call tp%spill(tmp, lspill_list=lmask, ldestructive=.true.) + ntp = tp%nbody + call tmp%setup(0,param) + deallocate(tmp) + deallocate(lmask) + end if + ntp = tp%nbody + if (ntp == 0) then + ! There are no more test particles. Reset the encounter list and move on + if (allocated(nbody_system%pltp_encounter)) call nbody_system%pltp_encounter%setup(0_I8B) + return + end if + + ! Reset all of thes tatus flags for the remaining bodies + tp%status(1:ntp) = ACTIVE + do i = 1, ntp + call tp%info(i)%set_value(status="ACTIVE") + end do + tp%ldiscard(1:ntp) = .false. + tp%lcollision(1:ntp) = .false. + tp%lmask(1:ntp) = .true. + + if (allocated(nbody_system%pltp_encounter)) then + ! Store the original pltpenc list so we don't remove any of the original encounters + nenc_old = nbody_system%pltp_encounter%nenc + if (nenc_old > 0_I8B) then + allocate(pltpenc_old, source=nbody_system%pltp_encounter) + call pltpenc_old%copy(nbody_system%pltp_encounter) + end if + + ! Re-build the encounter list + ! Be sure to get the level info if this is a SyMBA nbody_system + select type(nbody_system) + class is (symba_nbody_system) + select type(tp) + class is (symba_tp) + lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec) + end select + end select + end if + + end associate + + return + end subroutine swiftest_util_rearray_tp + + module subroutine swiftest_util_rescale_system(self, param, mscale, dscale, tscale) !! author: David A. Minton !!