Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Browse files Browse the repository at this point in the history 
Added make_family method to the pl type in symba
  • Loading branch information
@daminton
daminton committed Jul 31, 2021
1 parent 2e06f93 commit 2389b26
Show file tree
Hide file tree
Showing 2 changed files with 64 additions and 41 deletions.
21 changes: 14 additions & 7 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

!********************************************************************************************************************************
Expand Down Expand Up @@ -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
Expand Down
84 changes: 50 additions & 34 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

0 comments on commit 2389b26

Please sign in to comment.