diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index fb5f6ec06..2b131ef76 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -86,13 +86,14 @@ module symba_classes type(symba_kinship), dimension(:), allocatable :: kin !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step type(symba_particle_info), dimension(:), allocatable :: info contains - procedure :: discard => symba_discard_pl !! Process massive body discards - procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level - procedure :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other - procedure :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies - procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for number of particle - procedure :: sort => symba_util_sort_pl !! Sorts body arrays by a sortable componen - procedure :: rearrange => symba_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods + procedure :: make_family => symba_collision_make_family_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family + procedure :: discard => symba_discard_pl !! Process massive body discards + procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level + procedure :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other + procedure :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies + procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for number of particle + procedure :: sort => symba_util_sort_pl !! Sorts body arrays by a sortable componen + procedure :: rearrange => symba_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods end type symba_pl !******************************************************************************************************************************** @@ -189,6 +190,12 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec integer(I4B), intent(in) :: irec !! Current recursion level end subroutine symba_collision_check_plplenc + module subroutine symba_collision_make_family_pl(self,idx) + implicit none + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision + end subroutine symba_collision_make_family_pl + module subroutine symba_discard_pl(self, system, param) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index bdf83a595..c15e5df81 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -29,41 +29,38 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec select type(pl => system%pl) class is (symba_pl) - select type(tp => system%tp) - class is (symba_tp) - associate(plplenc_list => self, nplplenc => self%nenc, ind1 => self%index1, ind2 => self%index2) - allocate(lmask(nplplenc)) - lmask(:) = ((plplenc_list%status(1:nplplenc) == ACTIVE) & - .and. (pl%levelg(ind1(1:nplplenc)) >= irec) & - .and. (pl%levelg(ind2(1:nplplenc)) >= irec)) - if (.not.any(lmask(:))) return - - allocate(lcollision(nplplenc)) - lcollision(:) = .false. - - do concurrent(k = 1:nplplenc, lmask(k)) - xr(:) = pl%xh(:, ind1(k)) - pl%xh(:, ind2(k)) - vr(:) = pl%vb(:, ind1(k)) - pl%vb(:, ind2(k)) - rlim = pl%radius(ind1(k)) + pl%radius(ind2(k)) - mtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), mtot, rlim, dt, plplenc_list%lvdotr(k)) + associate(plplenc_list => self, nplplenc => self%nenc, ind1 => self%index1, ind2 => self%index2) + allocate(lmask(nplplenc)) + lmask(:) = ((plplenc_list%status(1:nplplenc) == ACTIVE) & + .and. (pl%levelg(ind1(1:nplplenc)) >= irec) & + .and. (pl%levelg(ind2(1:nplplenc)) >= irec)) + if (.not.any(lmask(:))) return + + allocate(lcollision(nplplenc)) + lcollision(:) = .false. + + do concurrent(k = 1:nplplenc, lmask(k)) + xr(:) = pl%xh(:, ind1(k)) - pl%xh(:, ind2(k)) + vr(:) = pl%vb(:, ind1(k)) - pl%vb(:, ind2(k)) + rlim = pl%radius(ind1(k)) + pl%radius(ind2(k)) + mtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), mtot, rlim, dt, plplenc_list%lvdotr(k)) + end do + + if (any(lcollision(:))) then + do k = 1, nplplenc + if (plplenc_list%status(k) /= COLLISION) cycle + plplenc_list%status(k) = COLLISION + plplenc_list%xh1(:,k) = pl%xh(:,ind1(k)) + plplenc_list%vb1(:,k) = pl%vb(:,ind1(k)) + plplenc_list%xh2(:,k) = pl%xh(:,ind2(k)) + plplenc_list%vb2(:,k) = pl%vb(:,ind2(k)) + if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call pl%make_family([ind1(k),ind2(k)]) + pl%lcollision(ind1(k)) = .true. + pl%lcollision(ind2(k)) = .true. end do - - if (any(lcollision(:))) then - do k = 1, nplplenc - if (plplenc_list%status(k) /= COLLISION) cycle - plplenc_list%status(k) = COLLISION - plplenc_list%xh1(:,k) = pl%xh(:,ind1(k)) - plplenc_list%vb1(:,k) = pl%vb(:,ind1(k)) - plplenc_list%xh2(:,k) = pl%xh(:,ind2(k)) - plplenc_list%vb2(:,k) = pl%vb(:,ind2(k)) - if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call plplenc_list%make_family(system) - pl%lcollision(ind1(k)) = .true. - pl%lcollision(ind2(k)) = .true. - end do - end if - end associate - end select + end if + end associate end select return @@ -177,4 +174,23 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt return end function symba_collision_check_one + module subroutine symba_collision_make_family_pl(self, 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 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 symba_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routine symba5_merge.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision + ! Internals + + return + end subroutine symba_collision_make_family_pl + end submodule s_symba_collision \ No newline at end of file