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

Commit

Permalink
Improved reporting of collisions with body names and ids
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 28, 2021
1 parent 24ad0d6 commit 3a79624
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 11 deletions.
18 changes: 14 additions & 4 deletions src/discard/discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,14 +125,14 @@ subroutine discard_cb_tp(tp, system, param)
tp%status(i) = DISCARDED_RMAX
write(idstr, *) tp%id(i)
write(timestr, *) param%t
write(*, *) "Particle " // trim(adjustl(idstr)) // " too far from the central body at t = " // trim(adjustl(timestr))
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
tp%lmask(i) = .false.
else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then
tp%status(i) = DISCARDED_RMIN
write(idstr, *) tp%id(i)
write(timestr, *) param%t
write(*, *) "Particle " // trim(adjustl(idstr)) // " too close to the central body at t = " // trim(adjustl(timestr))
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
tp%lmask(i) = .false.
else if (param%rmaxu >= 0.0_DP) then
Expand All @@ -143,7 +143,7 @@ subroutine discard_cb_tp(tp, system, param)
tp%status(i) = DISCARDED_RMAXU
write(idstr, *) tp%id(i)
write(timestr, *) param%t
write(*, *) "Particle " // trim(adjustl(idstr)) // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
tp%lmask(i) = .false.
end if
Expand Down Expand Up @@ -173,6 +173,7 @@ subroutine discard_peri_tp(tp, system, param)
integer(I4B) :: i, j, ih
real(DP) :: r2
real(DP), dimension(NDIM) :: dx
character(len=STRMAX) :: idstr, timestr

associate(cb => system%cb, ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t)
call tp%get_peri(system, param)
Expand All @@ -190,7 +191,9 @@ subroutine discard_peri_tp(tp, system, param)
(tp%atp(i) <= param%qmin_ahi) .and. &
(tp%peri(i) <= param%qmin)) then
tp%status(i) = DISCARDED_PERI
write(*, *) "Particle ", tp%id(i), " perihelion distance too small at t = ", t
write(idstr, *) tp%id(i)
write(timestr, *) param%t
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " perihelion distance too small at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
end if
end if
Expand Down Expand Up @@ -219,6 +222,7 @@ subroutine discard_pl_tp(tp, system, param)
integer(I4B) :: i, j, isp
real(DP) :: r2min, radius
real(DP), dimension(NDIM) :: dx, dv
character(len=STRMAX) :: idstri, idstrj, timestr

associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t, dt => param%dt)
do i = 1, ntp
Expand All @@ -233,6 +237,12 @@ subroutine discard_pl_tp(tp, system, param)
tp%lmask(i) = .false.
pl%ldiscard(j) = .true.
write(*, *) "Particle ", tp%id(i), " too close to massive body ", pl%id(j), " at t = ", t
write(idstri, *) tp%id(i)
write(idstrj, *) pl%id(j)
write(timestr, *) param%t
write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" &
// " too close to massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstrj)) &
// " at t = " // trim(adjustl(timestr))
tp%ldiscard(i) = .true.
exit
end if
Expand Down
52 changes: 48 additions & 4 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,11 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass
real(DP), dimension(:), allocatable :: m_frag, rad_frag
integer(I4B), dimension(:), allocatable :: id_frag
logical :: lfailure
character(len=STRMAX) :: collider_message

write(*, '("Disruption between bodies ",I8,99(:,",",I8))') system%pl%id(family(:))
collider_message = "Disruption between"
call symba_collision_collider_message(system%pl, family, collider_message)
write(*,*) trim(adjustl(collider_message))

! Collisional fragments will be uniformly distributed around the pre-impact barycenter
nfrag = NFRAG_DISRUPT
Expand Down Expand Up @@ -125,8 +128,12 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass,
integer(I4B), dimension(:), allocatable :: id_frag
logical :: lpure
logical, dimension(system%pl%nbody) :: lmask
character(len=STRMAX) :: collider_message
character(len=NAMELEN) :: idstr

write(*, '("Hit and run between bodies ",I8,99(:,",",I8))') system%pl%id(family(:))
collider_message = "Hit and run between"
call symba_collision_collider_message(system%pl, family, collider_message)
write(*,*) trim(adjustl(collider_message))

mtot = sum(mass(:))
xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot
Expand Down Expand Up @@ -232,10 +239,15 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad
real(DP), dimension(NDIM, 1) :: vb_frag, xb_frag, rot_frag, Ip_frag
real(DP), dimension(1) :: m_frag, rad_frag
integer(I4B), dimension(1) :: id_frag
character(len=STRMAX) :: collider_message
character(len=NAMELEN) :: idstr

collider_message = "Merging"
call symba_collision_collider_message(system%pl, family, collider_message)
write(*,*) trim(adjustl(collider_message))

select type(pl => system%pl)
class is (symba_pl)
write(*, '("Merging bodies ",I8,99(:,",",I8))') pl%id(family(:))

ibiggest = family(maxloc(pl%Gmass(family(:)), dim=1))
id_frag(1) = pl%id(ibiggest)
Expand Down Expand Up @@ -338,8 +350,12 @@ module function symba_collision_casesupercatastrophic(system, param, family, x,
integer(I4B), dimension(:), allocatable :: id_frag
logical :: lfailure
logical, dimension(system%pl%nbody) :: lmask
character(len=STRMAX) :: collider_message
character(len=NAMELEN) :: idstr

write(*, '("Supercatastrophic disruption between bodies ",I8,99(:,",",I8))') system%pl%id(family(:))
collider_message = "Supercatastrophic disruption between"
call symba_collision_collider_message(system%pl, family, collider_message)
write(*,*) trim(adjustl(collider_message))

! Collisional fragments will be uniformly distributed around the pre-impact barycenter
nfrag = NFRAG_SUPERCAT
Expand Down Expand Up @@ -403,6 +419,34 @@ module function symba_collision_casesupercatastrophic(system, param, family, x,
end function symba_collision_casesupercatastrophic


subroutine symba_collision_collider_message(pl, family, collider_message)
!! author: David A. Minton
!!
!! Prints a nicely formatted message about which bodies collided, including their names and ids.
!! This subroutine appends the body names and ids to an input message.
implicit none
! Arguments
class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object
integer(I4B), dimension(:), intent(in) :: family !! Index of collisional family members
character(*), intent(inout) :: collider_message !! The message to print to the screen.
! Internals
integer(I4B) :: i, n
character(len=STRMAX) :: idstr

n = size(family)
if (n == 0) return

do i = 1, n
if (i > 1) collider_message = trim(adjustl(collider_message)) // " and "
collider_message = " " // trim(adjustl(collider_message)) // " " // trim(adjustl(pl%info(family(i))%name))
write(idstr, '(I10)') pl%id(family(i))
collider_message = trim(adjustl(collider_message)) // " (" // trim(adjustl(idstr)) // ") "
end do

return
end subroutine symba_collision_collider_message


module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision)
!! author: David A. Minton
!!
Expand Down
6 changes: 3 additions & 3 deletions src/symba/symba_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,14 +36,14 @@ subroutine symba_discard_cb_pl(pl, system, param)
pl%status(i) = DISCARDED_RMAX
write(idstr, *) pl%id(i)
write(timestr, *) param%t
write(*, *) "Massive body " // trim(adjustl(idstr)) // " too far from the central body at t = " // trim(adjustl(timestr))
write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr))
else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then
pl%ldiscard(i) = .true.
pl%lcollision(i) = .false.
pl%status(i) = DISCARDED_RMIN
write(idstr, *) pl%id(i)
write(timestr, *) param%t
write(*, *) "Massive body " // trim(adjustl(idstr)) // " too close to the central body at t = " // trim(adjustl(timestr))
write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr))
else if (param%rmaxu >= 0.0_DP) then
rb2 = dot_product(pl%xb(:,i), pl%xb(:,i))
vb2 = dot_product(pl%vb(:,i), pl%vb(:,i))
Expand All @@ -54,7 +54,7 @@ subroutine symba_discard_cb_pl(pl, system, param)
pl%status(i) = DISCARDED_RMAXU
write(idstr, *) pl%id(i)
write(timestr, *) param%t
write(*, *) "Massive body " // trim(adjustl(idstr)) // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
end if
end if
end if
Expand Down

0 comments on commit 3a79624

Please sign in to comment.