diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index c30947d68..0fe43e391 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -297,9 +297,10 @@ module swiftest_classes real(DP), dimension(:,:), allocatable :: v1 !! the velocity of body 1 in the encounter real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter contains - procedure :: copy => util_copy_encounter - procedure :: resize => util_resize_encounter - procedure :: setup => setup_encounter + procedure :: setup => setup_encounter !! A constructor that sets the number of encounters and allocates and initializes all arrays + procedure :: copy => util_copy_encounter !! Copies elements from the source encounter list into self. + procedure :: spill => util_spill_encounter !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure :: resize => util_resize_encounter !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. end type swiftest_encounter abstract interface @@ -1152,6 +1153,14 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not end subroutine util_spill_body + module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive) + implicit none + class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list + class(swiftest_encounter), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + end subroutine util_spill_encounter + module subroutine util_spill_pl(self, discards, lspill_list, ldestructive) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index ee032c7f3..d682020ae 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -133,7 +133,7 @@ module symba_classes procedure :: encounter_check => symba_encounter_check_pltpenc !! Checks if massive bodies are going through close encounters with each other procedure :: kick => symba_kick_pltpenc !! Kick barycentric velocities of active test particles within SyMBA recursion procedure :: setup => symba_setup_pltpenc !! A constructor that sets the number of encounters and allocates and initializes all arrays - procedure :: copy => symba_util_copy_pltpenc !! Copies all elements of one pltpenc list to another + procedure :: spill => symba_util_spill_pltpenc !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) end type symba_pltpenc !******************************************************************************************************************************** @@ -447,13 +447,6 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) class(swiftest_body), intent(in) :: source !! Source object to append logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_tp - - module subroutine symba_util_copy_pltpenc(self, source) - use swiftest_classes, only : swiftest_encounter - implicit none - class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list - class(swiftest_encounter), intent(in) :: source !! Source object to copy into - end subroutine symba_util_copy_pltpenc end interface interface util_fill @@ -580,6 +573,15 @@ module subroutine symba_util_spill_pl(self, discards, lspill_list, ldestructive) logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not end subroutine symba_util_spill_pl + module subroutine symba_util_spill_pltpenc(self, discards, lspill_list, ldestructive) + use swiftest_classes, only : swiftest_encounter + implicit none + class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list + class(swiftest_encounter), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + end subroutine symba_util_spill_pltpenc + module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) use swiftest_classes, only : swiftest_body implicit none diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 5cc886485..739ca9778 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -143,6 +143,74 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt end function symba_collision_check_one + module subroutine symba_collision_encounter_scrub(self, 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(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + logical, dimension(self%nenc) :: lplpl_collision + logical, dimension(:), allocatable :: lplpl_unique_parent + integer(I4B), dimension(:), pointer :: plparent + integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx + integer(I4B) :: i, index_coll, ncollisions, nunique_parent + type(symba_plplenc) :: plplenc_noncollision + + + select type (pl => system%pl) + class is (symba_pl) + associate(plplenc_list => self, nplplenc => self%nenc, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) + lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION + if (.not.any(lplpl_collision)) return + + ! Get the subset of pl-pl encounters that lead to a collision + ncollisions = count(lplpl_collision(:)) + allocate(collision_idx(ncollisions)) + collision_idx = pack([(i, i=1, 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 hapepn 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, 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(:) = .true. + lplpl_collision(collision_idx(:)) = .false. + call plplenc_list%spill(plplenc_noncollision, lplpl_collision, ldestructive = .true.) + end associate + end select + + return + end subroutine symba_collision_encounter_scrub + + module subroutine symba_collision_make_family_pl(self, idx) !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton !! @@ -204,20 +272,4 @@ module subroutine symba_collision_make_family_pl(self, idx) return end subroutine symba_collision_make_family_pl - module subroutine symba_collision_encounter_scrub(self, 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(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - - return - end subroutine - - end submodule s_symba_collision \ No newline at end of file diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 7050bf050..56eaacc2c 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -142,25 +142,6 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) return end subroutine symba_util_append_tp - module subroutine symba_util_copy_pltpenc(self, source) - !! author: David A. Minton - !! - !! Copies elements from the source encounter list into self. - implicit none - ! Arguments - class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list - class(swiftest_encounter), intent(in) :: source !! Source object to copy into - - call util_copy_encounter(self, source) - associate(n => source%nenc) - select type(source) - class is (symba_pltpenc) - self%level(1:n) = source%level(1:n) - end select - end associate - - return - end subroutine symba_util_copy_pltpenc module subroutine symba_util_fill_arr_info(keeps, inserts, lfill_list) !! author: David A. Minton @@ -730,6 +711,35 @@ module subroutine symba_util_spill_pl(self, discards, lspill_list, ldestructive) end subroutine symba_util_spill_pl + module subroutine symba_util_spill_pltpenc(self, discards, lspill_list, ldestructive) + !! author: David A. Minton + !! + !! Move spilled (discarded) SyMBA encounter structure from active list to discard list + !! Note: Because the symba_plplenc currently does not contain any additional variable components, this method can recieve it as an input as well. + implicit none + ! Arguments + class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list + class(swiftest_encounter), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + ! Internals + integer(I4B) :: i + + associate(keeps => self) + select type(discards) + class is (symba_pltpenc) + call util_spill(keeps%level, discards%level, lspill_list, ldestructive) + call util_spill_encounter(keeps, discards, lspill_list, ldestructive) + class default + write(*,*) "Invalid object passed to the spill method. Source must be of class symba_pltpenc or its descendents!" + call util_exit(FAILURE) + end select + end associate + + return + end subroutine symba_util_spill_pltpenc + + module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) !! author: David A. Minton !! diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index 3772a8207..c6d5aa34f 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -211,6 +211,9 @@ module subroutine util_resize_encounter(self, nnew) !! author: David A. Minton !! !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. + !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an + !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment + !! Memory usage grows by a factor of 2 each time it fills up, but no more. implicit none ! Arguments class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 5f942854a..fc945765c 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -194,6 +194,39 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) return end subroutine util_spill_body + module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive) + !! author: David A. Minton + !! + !! Move spilled (discarded) Swiftest encounter structure from active list to discard list + implicit none + ! Arguments + class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list + class(swiftest_encounter), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + ! Internals + integer(I4B) :: i + + associate(keeps => self) + + call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive) + call util_spill(keeps%status, discards%status, lspill_list, ldestructive) + call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive) + call util_spill(keeps%index2, discards%index2, lspill_list, ldestructive) + call util_spill(keeps%x1, discards%x1, lspill_list, ldestructive) + call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive) + call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive) + call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive) + + ! This is the base class, so will be the last to be called in the cascade. + ! Therefore we need to set the nenc values for both the keeps and discareds + discards%nenc = count(lspill_list(:)) + keeps%nenc = count(.not.lspill_list(:)) + end associate + + return + end subroutine util_spill_encounter + module subroutine util_spill_pl(self, discards, lspill_list, ldestructive) !! author: David A. Minton