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

Commit

Permalink
Fixed bugs related to not properly removing pl-tp encounters from the…
Browse files Browse the repository at this point in the history 
… list after a pl-tp collision takes place.
  • Loading branch information
@daminton
daminton committed Dec 14, 2021
1 parent 72904e2 commit 2fbd367
Show file tree
Hide file tree
Showing 2 changed files with 53 additions and 46 deletions.
90 changes: 47 additions & 43 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec
real(DP) :: rlim, Gmtot
logical :: isplpl
character(len=STRMAX) :: timestr, idstri, idstrj, message

class(symba_encounter), allocatable :: tmp

lany_collision = .false.
if (self%nenc == 0) return
Expand Down Expand Up @@ -324,54 +324,58 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec
end do
end if

if (any(lcollision(1:nenc))) call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
do k = 1, nenc
i = self%index1(k)
j = self%index2(k)
if (lcollision(k)) self%status(k) = COLLISION
self%t(k) = t
self%x1(:,k) = pl%xh(:,i) + system%cb%xb(:)
self%v1(:,k) = pl%vb(:,i)
if (isplpl) then
self%x2(:,k) = pl%xh(:,j) + system%cb%xb(:)
self%v2(:,k) = pl%vb(:,j)
if (lcollision(k)) then
! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx
if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j])

! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step
pl%lcollision([i, j]) = .true.
pl%status([i, j]) = COLLISION
call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i))
call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,j), discard_vh=pl%vh(:,j))
end if
else
self%x2(:,k) = tp%xh(:,j) + system%cb%xb(:)
self%v2(:,k) = tp%vb(:,j)
if (lcollision(k)) then
tp%status(j) = DISCARDED_PLR
tp%ldiscard(j) = .true.
write(idstri, *) pl%id(i)
write(idstrj, *) tp%id(j)
write(timestr, *) t
call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_xh=tp%xh(:,j), discard_vh=tp%vh(:,j))
write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" &
// " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" &
// " at t = " // trim(adjustl(timestr))
call io_log_one_message(FRAGGLE_LOG_OUT, message)
lany_collision = any(lcollision(:))
if (lany_collision) then
call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
do k = 1, nenc
i = self%index1(k)
j = self%index2(k)
if (lcollision(k)) self%status(k) = COLLISION
self%t(k) = t
self%x1(:,k) = pl%xh(:,i) + system%cb%xb(:)
self%v1(:,k) = pl%vb(:,i)
if (isplpl) then
self%x2(:,k) = pl%xh(:,j) + system%cb%xb(:)
self%v2(:,k) = pl%vb(:,j)
if (lcollision(k)) then
! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx
if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j])

! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step
pl%lcollision([i, j]) = .true.
pl%status([i, j]) = COLLISION
call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i))
call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,j), discard_vh=pl%vh(:,j))
end if
else
self%x2(:,k) = tp%xh(:,j) + system%cb%xb(:)
self%v2(:,k) = tp%vb(:,j)
if (lcollision(k)) then
tp%status(j) = DISCARDED_PLR
tp%ldiscard(j) = .true.
write(idstri, *) pl%id(i)
write(idstrj, *) tp%id(j)
write(timestr, *) t
call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_xh=tp%xh(:,j), discard_vh=tp%vh(:,j))
write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" &
// " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" &
// " at t = " // trim(adjustl(timestr))
call io_log_one_message(FRAGGLE_LOG_OUT, message)
end if
end if
end if
end do
end do
end if
end select
end select

lany_collision = any(lcollision(:))

! Extract the pl-pl encounter list and return the plplcollision_list
! Extract the pl-pl or pl-tpencounter list and return the plplcollision_list
if (lany_collision) then
select type(plplenc_list => self)
select type(self)
class is (symba_plplenc)
call plplenc_list%extract_collisions(system, param)
call self%extract_collisions(system, param)
class default
allocate(tmp, mold=self)
call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list
end select
end if

Expand Down
9 changes: 6 additions & 3 deletions src/symba/symba_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,6 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn)

if (self%nenc == 0) return


select type(self)
class is (symba_plplenc)
isplpl = .true.
Expand All @@ -175,8 +174,12 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn)
select type(tp => system%tp)
class is (symba_tp)
associate(npl => pl%nbody, ntp => tp%nbody, nenc => self%nenc)
if (npl > 0) pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE
if (ntp > 0) tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE
if (npl == 0) return
pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE
if (.not. isplpl) then
if (ntp == 0) return
tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE
end if
allocate(lgoodlevel(nenc))

irm1 = irec - 1
Expand Down

0 comments on commit 2fbd367

Please sign in to comment.