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

Commit

Permalink
Started work on having pl-tp collisions recorded in the collision.nc …
Browse files Browse the repository at this point in the history 
…file. Still need to get tp snapshots saved
  • Loading branch information
@daminton
daminton committed Feb 22, 2024
1 parent f198460 commit e85df57
Show file tree
Hide file tree
Showing 4 changed files with 112 additions and 52 deletions.
2 changes: 2 additions & 0 deletions src/collision/collision_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -389,8 +389,10 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param)
if (allocated(pl)) deallocate(pl)
select case(stage)
case(1)
if (.not. allocated(before%pl)) cycle
allocate(pl, source=before%pl)
case(2)
if (.not. allocated(after%pl)) cycle
allocate(pl, source=after%pl)
end select
npl = pl%nbody
Expand Down
51 changes: 51 additions & 0 deletions src/collision/collision_resolve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -707,6 +707,14 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec)
real(DP), intent(in) :: t !! Current simulation tim
real(DP), intent(in) :: dt !! Current simulation step size
integer(I4B), intent(in) :: irec !! Current recursion level
! Internals
class(swiftest_pl), allocatable :: plsub
logical :: lpltp_collision
character(len=STRMAX) :: timestr, idstr
integer(I4B) :: i, j, nnew, loop
integer(I8B) :: k, ncollisions
integer(I4B), dimension(:), allocatable :: idnew
logical, dimension(:), allocatable :: lmask

! Make sure coordinate systems are all synced up due to being inside the recursion at this point
select type(nbody_system)
Expand All @@ -726,6 +734,49 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec)

! Discard the collider
call nbody_system%tp%discard(nbody_system, param)

associate(idx1 => pltp_collision%index1, idx2 => pltp_collision%index2)
ncollisions = pltp_collision%nenc
write(timestr,*) t
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "")
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // &
"***********************************************************")
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision between test particle and massive body detected at time t = " // &
trim(adjustl(timestr)))
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // &
"***********************************************************")

do k = 1_I8B, ncollisions
! Advance the collision id number and save it
collider%maxid_collision = max(collider%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id))
collider%maxid_collision = collider%maxid_collision + 1
collider%collision_id = collider%maxid_collision
write(idstr,*) collider%collision_id
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "collision_id " // trim(adjustl(idstr)))
collider%impactors%regime = COLLRESOLVE_REGIME_MERGE
allocate(lmask, mold=pl%lmask)
lmask(:) = .false.
lmask(idx1(k)) = .true.

allocate(plsub, mold=pl)
call pl%spill(plsub, lmask, ldestructive=.false.)

! Save the before snapshots
select type(before => collider%before)
class is (swiftest_nbody_system)
call move_alloc(plsub, before%pl)
end select

call collision_history%take_snapshot(param,nbody_system, t, "particle")

call impactors%dealloc()
end do

! Destroy the collision list now that the collisions are resolved
call pltp_collision%setup(0_I8B)

end associate

end associate
end select
end select
Expand Down
109 changes: 58 additions & 51 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -756,8 +756,8 @@ end subroutine collision_util_shift_vector_to_origin
module subroutine collision_util_snapshot(self, param, nbody_system, t, arg)
!! author: David A. Minton
!!
!! Takes a minimal snapshot of the state of the nbody_system during an encounter so that the trajectories
!! can be played back through the encounter
!! Takes a minimal snapshot of the state of the nbody_system during a collision to record the before and after states of the
!! system through the collision.
implicit none
! Internals
class(collision_storage), intent(inout) :: self !! Swiftest storage object
Expand Down Expand Up @@ -791,63 +791,70 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg)
snapshot%t = t
case ("after")
phase_val = 2
case ("particle")
phase_val = -1
allocate(collision_snapshot :: snapshot)
allocate(snapshot%collider, source=nbody_system%collider)
case default
write(*,*) "collision_util_snapshot requies either 'before' or 'after' passed to 'arg'"
write(*,*) "collision_util_snapshot requies either 'before', 'after', or 'particle' passed to 'arg'"
return
end select

! Get and record the energy of the system before the collision
call nbody_system%get_energy_and_momentum(param)
snapshot%collider%L_orbit(:,phase_val) = nbody_system%L_orbit(:)
snapshot%collider%L_spin(:,phase_val) = nbody_system%L_spin(:)
snapshot%collider%L_total(:,phase_val) = nbody_system%L_total(:)
snapshot%collider%ke_orbit(phase_val) = nbody_system%ke_orbit
snapshot%collider%ke_spin(phase_val) = nbody_system%ke_spin
snapshot%collider%pe(phase_val) = nbody_system%pe
snapshot%collider%be(phase_val) = nbody_system%be
snapshot%collider%te(phase_val) = nbody_system%te

if (stage == "after") then
select type(before_snap => snapshot%collider%before )
class is (swiftest_nbody_system)
select type(before_orig => nbody_system%collider%before)
if (stage /= "particle" ) then
! Get and record the energy of the system before the collision
call nbody_system%get_energy_and_momentum(param)
snapshot%collider%L_orbit(:,phase_val) = nbody_system%L_orbit(:)
snapshot%collider%L_spin(:,phase_val) = nbody_system%L_spin(:)
snapshot%collider%L_total(:,phase_val) = nbody_system%L_total(:)
snapshot%collider%ke_orbit(phase_val) = nbody_system%ke_orbit
snapshot%collider%ke_spin(phase_val) = nbody_system%ke_spin
snapshot%collider%pe(phase_val) = nbody_system%pe
snapshot%collider%be(phase_val) = nbody_system%be
snapshot%collider%te(phase_val) = nbody_system%te

if (stage == "after") then
select type(before_snap => snapshot%collider%before )
class is (swiftest_nbody_system)
select type(plsub => before_orig%pl)
class is (swiftest_pl)
! Log the properties of the old and new bodies
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Removing bodies:")
do i = 1, plsub%nbody
write(message,*) trim(adjustl(plsub%info(i)%name)), " (", trim(adjustl(plsub%info(i)%particle_type)),")"
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end do

allocate(before_snap%pl, source=plsub)
select type(before_orig => nbody_system%collider%before)
class is (swiftest_nbody_system)
select type(plsub => before_orig%pl)
class is (swiftest_pl)
! Log the properties of the old and new bodies
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Removing bodies:")
do i = 1, plsub%nbody
write(message,*) trim(adjustl(plsub%info(i)%name)), " (", trim(adjustl(plsub%info(i)%particle_type)),")"
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end do

allocate(before_snap%pl, source=plsub)
end select
deallocate(before_orig%pl)
end select
deallocate(before_orig%pl)
end select
end select


select type(after_snap => snapshot%collider%after )
class is (swiftest_nbody_system)
select type(after_orig => nbody_system%collider%after)
class is (swiftest_nbody_system)
select type(plnew => after_orig%pl)
class is (swiftest_pl)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Adding bodies:")
do i = 1, plnew%nbody
write(message,*) trim(adjustl(plnew%info(i)%name)), " (", trim(adjustl(plnew%info(i)%particle_type)),")"
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end do
call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
"***********************************************************" // &
"***********************************************************")
allocate(after_snap%pl, source=plnew)
end select
deallocate(after_orig%pl)
end select
end select

select type(after_snap => snapshot%collider%after )
class is (swiftest_nbody_system)
select type(after_orig => nbody_system%collider%after)
class is (swiftest_nbody_system)
select type(plnew => after_orig%pl)
class is (swiftest_pl)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Adding bodies:")
do i = 1, plnew%nbody
write(message,*) trim(adjustl(plnew%info(i)%name)), " (", trim(adjustl(plnew%info(i)%particle_type)),")"
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end do
call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
"***********************************************************" // &
"***********************************************************")
allocate(after_snap%pl, source=plnew)
end select
deallocate(after_orig%pl)
end select
end select
end if
end if

if ((stage == "after") .or. (stage == "particle")) then
! Save the snapshot for posterity
call self%save(snapshot)
deallocate(snapshot)
Expand Down
2 changes: 1 addition & 1 deletion src/swiftest/swiftest_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,7 @@ subroutine swiftest_discard_pl_tp(tp, nbody_system, param)
write(idstri, *) tp%id(i)
write(idstrj, *) pl%id(j)
write(timestr, *) nbody_system%t
write(message, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" &
write(message, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" &
// " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" &
// " at t = " // trim(adjustl(timestr))
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
Expand Down

0 comments on commit e85df57

Please sign in to comment.