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

Commit

Permalink
Consolidated pltp/plpl encounter lists into a single encounter object
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 3, 2021
1 parent 947c4d4 commit 9efb88c
Show file tree
Hide file tree
Showing 8 changed files with 159 additions and 240 deletions.
20 changes: 19 additions & 1 deletion src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -292,8 +292,14 @@ module swiftest_classes
integer(I4B), dimension(:), allocatable :: status !! status of the interaction
integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter
integer(I4B), dimension(:), allocatable :: index2 !! position of the second body in the encounter
real(DP), dimension(:,:), allocatable :: x1 !! the position of body 1 in the encounter
real(DP), dimension(:,:), allocatable :: x2 !! the position of body 2 in the encounter
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 :: copy => util_copy_encounter
procedure :: resize => util_resize_encounter
procedure :: setup => setup_encounter
end type swiftest_encounter

abstract interface
Expand Down Expand Up @@ -707,6 +713,12 @@ module subroutine setup_construct_system(system, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine setup_construct_system

module subroutine setup_encounter(self, n)
implicit none
class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter structure
integer(I4B), intent(in) :: n !! Number of encounters to allocate space for
end subroutine setup_encounter

module subroutine setup_initialize_system(self, param)
implicit none
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object
Expand Down Expand Up @@ -952,6 +964,12 @@ module subroutine util_resize_body(self, nnew)
integer(I4B), intent(in) :: nnew !! New size neded
end subroutine util_resize_body

module subroutine util_resize_encounter(self, nnew)
implicit none
class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list
integer(I4B), intent(in) :: nnew !! New size of list needed
end subroutine util_resize_encounter

module subroutine util_resize_pl(self, nnew)
implicit none
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
Expand Down
33 changes: 0 additions & 33 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,22 +133,14 @@ module symba_classes
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 :: resize => symba_util_resize_pltpenc !! Checks the current size of the pltpenc_list against the required size and extends it by a factor of 2 more than requested if it is too small
end type symba_pltpenc

!********************************************************************************************************************************
! symba_plplenc class definitions and method interfaces
!*******************************************************************************************************************************
!> SyMBA class for tracking pl-pl close encounters in a step
type, extends(symba_pltpenc) :: symba_plplenc
real(DP), dimension(:,:), allocatable :: xh1 !! the heliocentric position of parent 1 in encounter
real(DP), dimension(:,:), allocatable :: xh2 !! the heliocentric position of parent 2 in encounter
real(DP), dimension(:,:), allocatable :: vb1 !! the barycentric velocity of parent 1 in encounter
real(DP), dimension(:,:), allocatable :: vb2 !! the barycentric velocity of parent 2 in encounter
contains
procedure :: collision_check => symba_collision_check_plplenc !! Checks if two massive bodies are going to collide
procedure :: setup => symba_setup_plplenc !! A constructor that sets the number of encounters and allocates and initializes all arrays
procedure :: copy => symba_util_copy_plplenc !! Copies all elements of one plplenc list to another
end type symba_plplenc

!********************************************************************************************************************************
Expand Down Expand Up @@ -182,17 +174,6 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec
integer(I4B), intent(in) :: irec !! Current recursion level
end subroutine symba_collision_check_pltpenc

module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec)
use swiftest_classes, only : swiftest_parameters
implicit none
class(symba_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! current time
real(DP), intent(in) :: dt !! step size
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
Expand Down Expand Up @@ -343,7 +324,6 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
end subroutine symba_kick_pltpenc


module subroutine symba_io_write_frame_info(self, iu, param)
use swiftest_classes, only : swiftest_parameters
implicit none
Expand Down Expand Up @@ -465,13 +445,6 @@ module subroutine symba_util_copy_pltpenc(self, source)
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

module subroutine symba_util_copy_plplenc(self, source)
use swiftest_classes, only : swifest_encounter
implicit none
class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list
class(swiftest_encounter), intent(in) :: source !! Source object to copy into
end subroutine symba_util_copy_plplenc
end interface

interface util_fill
Expand Down Expand Up @@ -529,12 +502,6 @@ module subroutine symba_util_resize_pl(self, nnew)
integer(I4B), intent(in) :: nnew !! New size neded
end subroutine symba_util_resize_pl

module subroutine symba_util_resize_pltpenc(self, nnew)
implicit none
class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list
integer(I4B), intent(in) :: nnew !! New size of list needed
end subroutine symba_util_resize_pltpenc

module subroutine symba_util_resize_tp(self, nnew)
implicit none
class(symba_tp), intent(inout) :: self !! SyMBA massive body object
Expand Down
44 changes: 44 additions & 0 deletions src/setup/setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,50 @@ module subroutine setup_construct_system(system, param)
end subroutine setup_construct_system


module subroutine setup_encounter(self, n)
!! author: David A. Minton
!!
!! A constructor that sets the number of encounters and allocates and initializes all arrays
!!
implicit none
! Arguments
class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter structure
integer(I4B), intent(in) :: n !! Number of encounters to allocate space for

self%nenc = n
if (n == 0) return

if (allocated(self%lvdotr)) deallocate(self%lvdotr)
if (allocated(self%status)) deallocate(self%status)
if (allocated(self%index1)) deallocate(self%index1)
if (allocated(self%index2)) deallocate(self%index2)
if (allocated(self%x1)) deallocate(self%x1)
if (allocated(self%x2)) deallocate(self%x2)
if (allocated(self%v1)) deallocate(self%v1)
if (allocated(self%v2)) deallocate(self%v2)

allocate(self%lvdotr(n))
allocate(self%status(n))
allocate(self%index1(n))
allocate(self%index2(n))
allocate(self%x1(NDIM,n))
allocate(self%x2(NDIM,n))
allocate(self%v1(NDIM,n))
allocate(self%v2(NDIM,n))

self%lvdotr(:) = .false.
self%status(:) = INACTIVE
self%index1(:) = 0
self%index2(:) = 0
self%x1(:,:) = 0.0_DP
self%x2(:,:) = 0.0_DP
self%v1(:,:) = 0.0_DP
self%v2(:,:) = 0.0_DP

return
end subroutine setup_encounter


module subroutine setup_initialize_system(self, param)
!! author: David A. Minton
!!
Expand Down
154 changes: 56 additions & 98 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,89 +2,12 @@
use swiftest
contains

module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec)
!! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton
!!
!! Check for merger between massive bodies in SyMBA. If the user has turned on the FRAGMENTATION feature, it will call the
!! symba_regime subroutine to determine what kind of collision will occur.
!!
!! 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_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! current time
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
! Internals
logical, dimension(:), allocatable :: lcollision, lmask
real(DP), dimension(NDIM) :: xr, vr
integer(I4B) :: k
real(DP) :: rlim, mtot

if (self%nenc == 0) return

select type(pl => system%pl)
class is (symba_pl)
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
deallocate(lmask)

if (any(lcollision(:))) then
allocate(lmask(pl%nbody))
do k = 1, nplplenc
if (.not.lcollision(k)) cycle

! Set this encounter as a collision and save the position and velocity vectors at the time of the collision
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))

! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional family
if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call pl%make_family([ind1(k),ind2(k)])

! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step
pl%lcollision([ind1(k), ind2(k)]) = .true.
pl%lcollision([ind1(k), ind2(k)]) = .true.
pl%status([ind1(k), ind2(k)]) = COLLISION
end do

end if
end associate
end select

return

return
end subroutine symba_collision_check_plplenc


module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec)
!! author: David A. Minton
!!
!! Check for merger between massive bodies and test particles in SyMBA
!!
!! Adapted from David E. Kaufmann's Swifter routine symba_merge_tp.f90
!! Adapted from David E. Kaufmann's Swifter routine symba_merge.f90 and symba_merge_tp.f90
!!
!! Adapted from Hal Levison's Swift routine symba5_merge.f
implicit none
Expand All @@ -99,39 +22,74 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec
logical, dimension(:), allocatable :: lcollision, lmask
real(DP), dimension(NDIM) :: xr, vr
integer(I4B) :: k
real(DP) :: rlim, mtot
logical :: isplpl

if (self%nenc == 0) return
select type(self)
class is (symba_plplenc)
isplpl = .true.
class default
isplpl = .false.
end select

select type(pl => system%pl)
class is (symba_pl)
select type(tp => system%tp)
class is (symba_tp)
associate(pltpenc_list => self, npltpenc => self%nenc, plind => self%index1, tpind => self%index2)
allocate(lmask(npltpenc))
lmask(:) = ((pltpenc_list%status(1:npltpenc) == ACTIVE) &
.and. (pl%levelg(plind(1:npltpenc)) >= irec) &
.and. (tp%levelg(tpind(1:npltpenc)) >= irec))
associate(nenc => self%nenc, ind1 => self%index1, ind2 => self%index2)
allocate(lmask(nenc))
lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(ind1(1:nenc)) >= irec))
if (isplpl) then
lmask(:) = lmask(:) .and. (pl%levelg(ind2(1:nenc)) >= irec)
else
lmask(:) = lmask(:) .and. (tp%levelg(ind2(1:nenc)) >= irec)
end if
if (.not.any(lmask(:))) return

allocate(lcollision(npltpenc))
allocate(lcollision(nenc))
lcollision(:) = .false.

do concurrent(k = 1:npltpenc, lmask(k))
xr(:) = pl%xh(:, plind(k)) - tp%xh(:, tpind(k))
vr(:) = pl%vb(:, plind(k)) - tp%vb(:, tpind(k))
lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(plind(k)), pl%radius(plind(k)), dt, pltpenc_list%lvdotr(k))
end do
if (isplpl) then
do concurrent(k = 1:nenc, 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, self%lvdotr(k))
end do
else
do concurrent(k = 1:nenc, lmask(k))
xr(:) = pl%xh(:, ind1(k)) - tp%xh(:, ind2(k))
vr(:) = pl%vb(:, ind1(k)) - tp%vb(:, ind2(k))
lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(ind1(k)), pl%radius(ind1(k)), dt, self%lvdotr(k))
end do
end if

if (any(lcollision(:))) then
where(lcollision(1:npltpenc))
pltpenc_list%status(1:npltpenc) = COLLISION
tp%status(tpind(1:npltpenc)) = DISCARDED_PLR
tp%ldiscard(tpind(1:npltpenc)) = .true.
end where

do k = 1, npltpenc
if (pltpenc_list%status(k) /= COLLISION) cycle
write(*,*) 'Test particle ',tp%id(tpind(k)), ' collided with massive body ',pl%id(plind(k)), ' at time ',t
do k = 1, nenc
if (.not.lcollision(k)) cycle
self%status(k) = COLLISION
self%x1(:,k) = pl%xh(:,ind1(k))
self%v1(:,k) = pl%vb(:,ind1(k))
if (isplpl) then
self%x2(:,k) = pl%xh(:,ind2(k))
self%v2(:,k) = pl%vb(:,ind2(k))

! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional family
if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call pl%make_family([ind1(k),ind2(k)])

! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step
pl%lcollision([ind1(k), ind2(k)]) = .true.
pl%ldiscard([ind1(k), ind2(k)]) = .true.
pl%status([ind1(k), ind2(k)]) = COLLISION
else
self%x2(:,k) = tp%xh(:,ind2(k))
self%v2(:,k) = tp%vb(:,ind2(k))
tp%status(ind2(k)) = DISCARDED_PLR
tp%ldiscard(ind2(k)) = .true.
write(*,*) 'Test particle ',tp%id(ind2(k)), ' collided with massive body ',pl%id(ind1(k)), ' at time ',t
end if
end do
end if
end associate
Expand Down
Loading

0 comments on commit 9efb88c

Please sign in to comment.