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

Commit

Permalink
Fixed it so particles can only have closest encounters after their or…
Browse files Browse the repository at this point in the history
…igin time. Also ensured that the actual collision time gets communicated down the chain.
  • Loading branch information
daminton committed Dec 14, 2022
1 parent 60277ba commit 00e5173
Show file tree
Hide file tree
Showing 3 changed files with 55 additions and 43 deletions.
1 change: 1 addition & 0 deletions src/encounter/encounter_util.f90
Original file line number Diff line number Diff line change
Expand Up @@ -690,6 +690,7 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg)
vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k)
call orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi)
snapshot%t = t + tperi
if (snapshot%t < maxval(pl_snap%info(:)%origin_time)) cycle

! Computer the center mass of the pair
rcom(:) = (plplenc_list%r1(:,k) * pl_snap%Gmass(1) + plplenc_list%r2(:,k) * pl_snap%Gmass(2)) / Gmtot
Expand Down
31 changes: 18 additions & 13 deletions src/modules/symba_classes.f90
Original file line number Diff line number Diff line change
Expand Up @@ -232,18 +232,20 @@ module subroutine symba_collision_make_colliders_pl(self,idx)
integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision
end subroutine symba_collision_make_colliders_pl

module subroutine symba_resolve_collision_fragmentations(self, system, param)
module subroutine symba_resolve_collision_fragmentations(self, system, param, t)
implicit none
class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
end subroutine symba_resolve_collision_fragmentations

module subroutine symba_resolve_collision_mergers(self, system, param)
module subroutine symba_resolve_collision_mergers(self, system, param, t)
implicit none
class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
end subroutine symba_resolve_collision_mergers

module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, irec)
Expand Down Expand Up @@ -343,25 +345,28 @@ pure module subroutine symba_gr_p4_tp(self, system, param, dt)
real(DP), intent(in) :: dt !! Step size
end subroutine symba_gr_p4_tp

module function symba_collision_casedisruption(system, param) result(status)
module function symba_collision_casedisruption(system, param, t) result(status)
implicit none
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
integer(I4B) :: status !! Status flag assigned to this outcome
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
integer(I4B) :: status !! Status flag assigned to this outcome
end function symba_collision_casedisruption

module function symba_collision_casehitandrun(system, param) result(status)
module function symba_collision_casehitandrun(system, param, t) result(status)
implicit none
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
integer(I4B) :: status !! Status flag assigned to this outcome
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
integer(I4B) :: status !! Status flag assigned to this outcome
end function symba_collision_casehitandrun

module function symba_collision_casemerge(system, param) result(status)
module function symba_collision_casemerge(system, param, t) result(status)
implicit none
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
integer(I4B) :: status !! Status flag assigned to this outcome
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
integer(I4B) :: status !! Status flag assigned to this outcome
end function symba_collision_casemerge

module subroutine symba_util_set_renc(self, scale)
Expand Down
66 changes: 36 additions & 30 deletions src/symba/symba_collision.f90
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,16 @@

contains

module function symba_collision_casedisruption(system, param) result(status)
module function symba_collision_casedisruption(system, param, t) result(status)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Create the fragments resulting from a non-catastrophic disruption collision
!!
implicit none
! Arguments
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
! Result
integer(I4B) :: status !! Status flag assigned to this outcome
! Internals
Expand Down Expand Up @@ -73,23 +74,24 @@ module function symba_collision_casedisruption(system, param) result(status)
param%maxid = fragments%id(nfrag)
end select

call symba_collision_mergeaddsub(system, param, status)
call symba_collision_mergeaddsub(system, param, t, status)
end if
end associate

return
end function symba_collision_casedisruption


module function symba_collision_casehitandrun(system, param) result(status)
module function symba_collision_casehitandrun(system, param, t) result(status)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Create the fragments resulting from a non-catastrophic hit-and-run collision
!!
implicit none
! Arguments
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
! Result
integer(I4B) :: status !! Status flag assigned to this outcom
! Internals
Expand Down Expand Up @@ -144,7 +146,7 @@ module function symba_collision_casehitandrun(system, param) result(status)
fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)]
param%maxid = fragments%id(nfrag)
status = HIT_AND_RUN_DISRUPT
call symba_collision_mergeaddsub(system, param, status)
call symba_collision_mergeaddsub(system, param, t, status)
end if

end associate
Expand All @@ -153,7 +155,7 @@ module function symba_collision_casehitandrun(system, param) result(status)
end function symba_collision_casehitandrun


module function symba_collision_casemerge(system, param) result(status)
module function symba_collision_casemerge(system, param, t) result(status)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Merge massive bodies.
Expand All @@ -163,8 +165,9 @@ module function symba_collision_casemerge(system, param) result(status)
!! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f
implicit none
! Arguments
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
! Result
integer(I4B) :: status !! Status flag assigned to this outcome
! Internals
Expand Down Expand Up @@ -226,7 +229,7 @@ module function symba_collision_casemerge(system, param) result(status)

status = MERGED

call symba_collision_mergeaddsub(system, param, status)
call symba_collision_mergeaddsub(system, param, t, status)

end select
end associate
Expand Down Expand Up @@ -719,16 +722,17 @@ module subroutine symba_collision_make_colliders_pl(self, idx)
end subroutine symba_collision_make_colliders_pl


subroutine symba_collision_mergeaddsub(system, param, status)
subroutine symba_collision_mergeaddsub(system, param, t, status)
!! author: David A. Minton
!!
!! Fills the pl_discards and pl_adds with removed and added bodies
!!
implicit none
! Arguments
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
integer(I4B), intent(in) :: status !! Status flag to assign to adds
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
integer(I4B), intent(in) :: status !! Status flag to assign to adds
! Internals
integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag
logical, dimension(system%pl%nbody) :: lmask
Expand Down Expand Up @@ -775,7 +779,7 @@ subroutine symba_collision_mergeaddsub(system, param, status)
plnew%status(1:nfrag) = NEW_PARTICLE
do i = 1, nfrag
write(newname, FRAGFMT) fragments%id(i)
call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, &
call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=t, name=newname, &
origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), &
collision_id=param%maxid_collision)
end do
Expand All @@ -785,7 +789,7 @@ subroutine symba_collision_mergeaddsub(system, param, status)
else
iother = ibiggest
end if
call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, &
call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=t, &
discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), &
discard_body_id=iother)
end do
Expand All @@ -799,14 +803,14 @@ subroutine symba_collision_mergeaddsub(system, param, status)
plnew%status(1) = OLD_PARTICLE
do i = 2, nfrag
write(newname, FRAGFMT) fragments%id(i)
call plnew%info(i)%set_value(origin_type=origin_type, origin_time=system%t, name=newname, &
call plnew%info(i)%set_value(origin_type=origin_type, origin_time=t, name=newname, &
origin_rh=plnew%rh(:,i), origin_vh=plnew%vh(:,i), &
collision_id=param%maxid_collision)
end do
do i = 1, ncolliders
if (colliders%idx(i) == ibiggest) cycle
iother = ibiggest
call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=system%t, &
call pl%info(colliders%idx(i))%set_value(status=origin_type, discard_time=t, &
discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i), &
discard_body_id=iother)
end do
Expand All @@ -817,7 +821,7 @@ subroutine symba_collision_mergeaddsub(system, param, status)
if (colliders%idx(i) == ibiggest) cycle

iother = ibiggest
call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_rh=pl%rh(:,i), &
call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=t, discard_rh=pl%rh(:,i), &
discard_vh=pl%vh(:,i), discard_body_id=iother)
end do
end select
Expand Down Expand Up @@ -880,7 +884,7 @@ subroutine symba_collision_mergeaddsub(system, param, status)
end subroutine symba_collision_mergeaddsub


module subroutine symba_resolve_collision_fragmentations(self, system, param)
module subroutine symba_resolve_collision_fragmentations(self, system, param, t)
!! author: David A. Minton
!!
!! Process list of collisions, determine the collisional regime, and then create fragments.
Expand All @@ -890,13 +894,14 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)
class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
! Internals
! Internals
integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision
logical :: lgoodcollision
integer(I4B) :: i

associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, t => system%t, collision_history => param%collision_history)
associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, collision_history => param%collision_history)
select type(pl => system%pl)
class is (symba_pl)
select type (cb => system%cb)
Expand All @@ -914,11 +919,11 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)
if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "before")
select case (system%fragments%regime)
case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
plplcollision_list%status(i) = symba_collision_casedisruption(system, param)
plplcollision_list%status(i) = symba_collision_casedisruption(system, param, t)
case (COLLRESOLVE_REGIME_HIT_AND_RUN)
plplcollision_list%status(i) = symba_collision_casehitandrun(system, param)
plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, t)
case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE)
plplcollision_list%status(i) = symba_collision_casemerge(system, param)
plplcollision_list%status(i) = symba_collision_casemerge(system, param, t)
case default
write(*,*) "Error in symba_collision, unrecognized collision regime"
call util_exit(FAILURE)
Expand All @@ -934,7 +939,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)
end subroutine symba_resolve_collision_fragmentations


module subroutine symba_resolve_collision_mergers(self, system, param)
module subroutine symba_resolve_collision_mergers(self, system, param, t)
!! author: David A. Minton
!!
!! Process list of collisions and merge colliding bodies together.
Expand All @@ -944,6 +949,7 @@ module subroutine symba_resolve_collision_mergers(self, system, param)
class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
! Internals
integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision
logical :: lgoodcollision
Expand All @@ -968,7 +974,7 @@ module subroutine symba_resolve_collision_mergers(self, system, param)
fragments%mass_dist(3) = 0.0_DP
fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot
fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot
plplcollision_list%status(i) = symba_collision_casemerge(system, param)
plplcollision_list%status(i) = symba_collision_casemerge(system, param, t)
end do
end select
end select
Expand Down Expand Up @@ -1024,9 +1030,9 @@ module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, ir
"***********************************************************")
allocate(tmp_param, source=param)
if (param%lfragmentation) then
call plplcollision_list%resolve_fragmentations(system, param)
call plplcollision_list%resolve_fragmentations(system, param, t)
else
call plplcollision_list%resolve_mergers(system, param)
call plplcollision_list%resolve_mergers(system, param, t)
end if

! Destroy the collision list now that the collisions are resolved
Expand Down

0 comments on commit 00e5173

Please sign in to comment.