From 9eb5246a58a61af4a46248b2cfda1b2e40f618b2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 21 Feb 2024 13:36:43 -0500 Subject: [PATCH] Added missing pltp collision routines --- src/collision/collision_check.f90 | 8 ++------ src/collision/collision_resolve.f90 | 32 +++++++++++++++++++++++++++++ src/swiftest/swiftest_module.f90 | 4 ++-- src/swiftest/swiftest_util.f90 | 1 + src/symba/symba_step.f90 | 2 ++ src/symba/symba_util.f90 | 3 +-- 6 files changed, 40 insertions(+), 10 deletions(-) diff --git a/src/collision/collision_check.f90 b/src/collision/collision_check.f90 index be11ee9a8..c95547a6e 100644 --- a/src/collision/collision_check.f90 +++ b/src/collision/collision_check.f90 @@ -187,11 +187,9 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l lany_collision = .false. if (self%nenc == 0) return + select type(nbody_system) class is (swiftest_nbody_system) - select type(param) - class is (swiftest_parameters) - associate(pl => nbody_system%pl, tp => nbody_system%tp) nenc = self%nenc @@ -256,15 +254,13 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l end do ! Extract the pl-tp encounter list and return the pl-tp collision_list - allocate(tmp, mold=self) - call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list + call self%extract_collisions(nbody_system, param) end if ! Take snapshots of pairs of bodies at close approach (but not collision) if requested if (lany_closest) call nbody_system%encounter_history%take_snapshot(param, nbody_system, t, "closest") end associate end select - end select return end subroutine collision_check_pltp diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index a3daf7122..ab27b72e4 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -255,6 +255,38 @@ module subroutine collision_resolve_extract_pltp(self, nbody_system, param) 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, index_coll, 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. + call self%spill(nbody_system%pltp_collision, lpltp_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. + end associate + end select + end select + end select + return end subroutine collision_resolve_extract_pltp diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index fdefc6ec0..70c49a619 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -346,10 +346,10 @@ module swiftest class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure class(swiftest_pl), allocatable :: pl_adds !! List of added bodies in mergers or collisions class(swiftest_tp), allocatable :: tp_adds !! List of added bodies in mergers or collisions - class(encounter_list), allocatable :: pltp_encounter !! List of massive body-test particle encounters in a single step + class(encounter_list), allocatable :: pltp_encounter !! List of massive body-test particle encounters in a single step class(encounter_list), allocatable :: plpl_encounter !! List of massive body-massive body encounters in a single step class(collision_list_plpl), allocatable :: plpl_collision !! List of massive body-massive body collisions in a single step - class(collision_list_plpl), allocatable :: pltp_collision !! List of massive body-massive body collisions in a single step + class(collision_list_pltp), allocatable :: pltp_collision !! List of massive body-test particle collisions in a single step class(collision_basic), allocatable :: collider !! Collision system object class(encounter_storage), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file class(collision_storage), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 6b6b00b65..78bcb7b01 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -2412,6 +2412,7 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) allocate(symba_list_pltp :: nbody_system%pltp_encounter) allocate(symba_list_plpl :: nbody_system%plpl_encounter) allocate(collision_list_plpl :: nbody_system%plpl_collision) + allocate(collision_list_pltp :: nbody_system%pltp_collision) end select case (INT_RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index e647e4d37..948eeaacc 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -305,7 +305,9 @@ module subroutine symba_step_reset_system(self, param) nenc_old = nbody_system%pltp_encounter%nenc call nbody_system%pltp_encounter%setup(0_I8B) + call nbody_system%pltp_collision%setup(0_I8B) if (ntp > 0) then + tp%lcollision(1:ntp) = .false. tp%nplenc(1:ntp) = 0 tp%levelg(1:ntp) = -1 tp%levelm(1:ntp) = -1 diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 8ed213ba8..8e4b45e37 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -157,7 +157,6 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) associate(keeps => self) select type(inserts) class is (symba_tp) - call util_fill(keeps%nplenc, inserts%nplenc, lfill_list) call util_fill(keeps%levelg, inserts%levelg, lfill_list) call util_fill(keeps%levelm, inserts%levelm, lfill_list) @@ -290,6 +289,7 @@ module subroutine symba_util_setup_initialize_system(self, system_history, param call nbody_system%pltp_encounter%setup(0_I8B) call nbody_system%plpl_encounter%setup(0_I8B) call nbody_system%plpl_collision%setup(0_I8B) + call nbody_system%pltp_collision%setup(0_I8B) if (param%lenc_save_trajectory .or. param%lenc_save_closest) then allocate(encounter_netcdf_parameters :: encounter_history%nc) @@ -563,7 +563,6 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) associate(keeps => self) select type(discards) class is (symba_tp) - call util_spill(keeps%nplenc, discards%nplenc, lspill_list, ldestructive) call util_spill(keeps%levelg, discards%levelg, lspill_list, ldestructive) call util_spill(keeps%levelm, discards%levelm, lspill_list, ldestructive)