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

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 14, 2021
2 parents 5eded5f + 2fbd367 commit d4b22d1
Show file tree
Hide file tree
Showing 4 changed files with 56 additions and 49 deletions.
2 changes: 1 addition & 1 deletion examples/symba_gr_test/init_cond.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@
GMcb = sim.ds.isel(id=0)['Gmass'].values
GU = swiftest.GC / (sim.param['DU2M']**3 / (sim.param['MU2KG'] * sim.param['TU2S']**2))
dens = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M']**3) # Assume a bulk density of 3 g/cm^3
GM_pl = 1e-11
GM_pl = sim.ds.sel(id=1)['Gmass'].values * 0.01
M_pl = GM_pl / GU
R_pl = (3 * M_pl / (4 * np.pi * dens))**(1.0 / 3.0)
Rh_pl = Me_a * (GM_pl / (3 * GMcb))**(1.0/3.0)
Expand Down
4 changes: 2 additions & 2 deletions examples/symba_gr_test/pl.swiftest.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ Neptune 0.0020336100526728302882 0.7749718408665500834
0.000164587904124493665
30.038959912561129073 0.008955570159296157365 1.7711193542961420899
131.82211597488270627 284.47484279411258967 308.45137222693909962
Planetesimal 9.999999999999999395e-12 1.698243864025464137e-05
2.2876635862715332271e-07
Planetesimal 6.553709809565314101e-08 0.00031780626665496970825
4.281092009908029071e-06
0.38709858430953941744 0.20562340108970009189 7.0033025080013837638
48.296118373786072198 29.20442403952453958 339.33948746828792764
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 d4b22d1

Please sign in to comment.