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

Commit

Permalink
Lots more restructuring, refactoring, and cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 22, 2022
1 parent 12f3f53 commit 0896ffb
Show file tree
Hide file tree
Showing 56 changed files with 1,436 additions and 1,410 deletions.
36 changes: 18 additions & 18 deletions src/collision/collision_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ pure elemental subroutine collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmtot,
end subroutine collision_check_one


module subroutine collision_check_plpl(self, system, param, t, dt, irec, lany_collision)
module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, lany_collision)
!! author: David A. Minton
!!
!! Check for merger between massive bodies and test particles in SyMBA
Expand All @@ -68,7 +68,7 @@ module subroutine collision_check_plpl(self, system, param, t, dt, irec, lany_co
implicit none
! Arguments
class(collision_list_plpl), intent(inout) :: self !! SyMBA pl-tp encounter list object
class(base_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(base_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! current time
real(DP), intent(in) :: dt !! step size
Expand All @@ -86,9 +86,9 @@ module subroutine collision_check_plpl(self, system, param, t, dt, irec, lany_co
lany_collision = .false.
if (self%nenc == 0) return

select type(system)
select type(nbody_system)
class is (swiftest_nbody_system)
associate(pl => system%pl)
associate(pl => nbody_system%pl)

nenc = self%nenc
allocate(lmask(nenc))
Expand Down Expand Up @@ -118,18 +118,18 @@ module subroutine collision_check_plpl(self, system, param, t, dt, irec, lany_co


if (lany_collision .or. lany_closest) then
call pl%rh2rb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
do k = 1, nenc
if (.not.lcollision(k) .and. .not. self%lclosest(k)) cycle
i = self%index1(k)
j = self%index2(k)
self%r1(:,k) = pl%rh(:,i) + system%cb%rb(:)
self%r1(:,k) = pl%rh(:,i) + nbody_system%cb%rb(:)
self%v1(:,k) = pl%vb(:,i)
if (lcollision(k)) then
self%status(k) = COLLIDED
self%tcollision(k) = t
end if
self%r2(:,k) = pl%rh(:,j) + system%cb%rb(:)
self%r2(:,k) = pl%rh(:,j) + nbody_system%cb%rb(:)
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 collider pair
Expand All @@ -145,19 +145,19 @@ module subroutine collision_check_plpl(self, system, param, t, dt, irec, lany_co
end do

! Extract the pl-pl encounter list and return the pl-pl collision_list
call self%extract_collisions(system, param)
call self%extract_collisions(nbody_system, param)
end if

! Take snapshots of pairs of bodies at close approach (but not collision) if requested
if (lany_closest) call system%encounter_history%take_snapshot(param, system, t, "closest")
if (lany_closest) call nbody_system%encounter_history%take_snapshot(param, nbody_system, t, "closest")

end associate
end select
return
end subroutine collision_check_plpl


module subroutine collision_check_pltp(self, system, param, t, dt, irec, lany_collision)
module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, lany_collision)
!! author: David A. Minton
!!
!! Check for merger between massive bodies and test particles in SyMBA
Expand All @@ -168,7 +168,7 @@ module subroutine collision_check_pltp(self, system, param, t, dt, irec, lany_co
implicit none
! Arguments
class(collision_list_pltp), intent(inout) :: self !! SyMBA pl-tp encounter list object
class(base_nbody_system), intent(inout) :: system !! SyMBA nbody system object
class(base_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! current time
real(DP), intent(in) :: dt !! step size
Expand All @@ -185,12 +185,12 @@ module subroutine collision_check_pltp(self, system, param, t, dt, irec, lany_co

lany_collision = .false.
if (self%nenc == 0) return
select type(system)
select type(nbody_system)
class is (swiftest_nbody_system)
select type(param)
class is (swiftest_parameters)

associate(pl => system%pl, tp => system%tp)
associate(pl => nbody_system%pl, tp => nbody_system%tp)

nenc = self%nenc
allocate(lmask(nenc))
Expand Down Expand Up @@ -222,19 +222,19 @@ module subroutine collision_check_pltp(self, system, param, t, dt, irec, lany_co


if (lany_collision .or. lany_closest) then
call pl%rh2rb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
do k = 1, nenc
if (.not.lcollision(k) .and. .not. self%lclosest(k)) cycle
i = self%index1(k)
j = self%index2(k)
self%r1(:,k) = pl%rh(:,i) + system%cb%rb(:)
self%r1(:,k) = pl%rh(:,i) + nbody_system%cb%rb(:)
self%v1(:,k) = pl%vb(:,i)
if (lcollision(k)) then
self%status(k) = COLLIDED
self%tcollision(k) = t
end if

self%r2(:,k) = tp%rh(:,j) + system%cb%rb(:)
self%r2(:,k) = tp%rh(:,j) + nbody_system%cb%rb(:)
self%v2(:,k) = tp%vb(:,j)
if (lcollision(k)) then
tp%status(j) = DISCARDED_PLR
Expand All @@ -246,7 +246,7 @@ module subroutine collision_check_pltp(self, system, param, t, dt, irec, lany_co
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 swiftest_io_log_one_message(FRAGGLE_LOG_OUT, message)
!call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
end if
end do

Expand All @@ -256,7 +256,7 @@ module subroutine collision_check_pltp(self, system, param, t, dt, irec, lany_co
end if

! Take snapshots of pairs of bodies at close approach (but not collision) if requested
if (lany_closest) call system%encounter_history%take_snapshot(param, system, t, "closest")
if (lany_closest) call nbody_system%encounter_history%take_snapshot(param, nbody_system, t, "closest")
end associate
end select
end select
Expand Down
139 changes: 116 additions & 23 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,29 +12,122 @@
use swiftest
contains

module subroutine collision_generate_merge_system(self, nbody_system, param, t)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Merge massive bodies in a "MERGE" style collision model (only pure mergers)
implicit none
class(collision_merge), intent(inout) :: self !! Merge fragment system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! The time of the collision

module subroutine collision_generate_merge_system(self, nbody_system, param, t)
implicit none
class(collision_merge), intent(inout) :: self !! Merge fragment system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! The time of the collision
end subroutine collision_generate_merge_system

module subroutine collision_generate_bounce_system(self, nbody_system, param, t)
implicit none
class(collision_bounce), intent(inout) :: self !! Bounce fragment system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! The time of the collision
end subroutine collision_generate_bounce_system

module subroutine collision_generate_simple_system(self, nbody_system, param, t)
implicit none
class(collision_simple), intent(inout) :: self !! Simple fragment system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! The time of the collision
end subroutine collision_generate_simple_system
call collision_generate_merge_any(self, nbody_system, param, t)

return
end subroutine collision_generate_merge_system


module subroutine collision_generate_merge_any(self, nbody_system, param, t)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Merge massive bodies in any collisionals ystem.
!!
!! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90
!!
!! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f
implicit none
! Arguments
class(collision_system), intent(inout) :: self !! Merge fragment system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! The time of the collision
! Internals
integer(I4B) :: i, j, k, ibiggest
real(DP), dimension(NDIM) :: Lspin_new
real(DP) :: dpe
character(len=STRMAX) :: message

select type(nbody_system)
class is (swiftest_nbody_system)
associate(impactors => nbody_system%collider%impactors, fragments => nbody_system%collider%fragments)
message = "Merging"
call collision_io_collider_message(nbody_system%pl, impactors%id, message)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)

select type(pl => nbody_system%pl)
class is (swiftest_pl)

!call self%set_mass_dist(param)

! Calculate the initial energy of the nbody_system without the collisional family
call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)

ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1))
fragments%id(1) = pl%id(ibiggest)
fragments%rb(:,1) = impactors%rbcom(:)
fragments%vb(:,1) = impactors%vbcom(:)

if (param%lrotation) then
! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body
Lspin_new(:) = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2)

! Assume prinicpal axis rotation on 3rd Ip axis
fragments%rot(:,1) = Lspin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2)
else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable
nbody_system%Lescape(:) = nbody_system%Lescape(:) + impactors%Lorbit(:,1) + impactors%Lorbit(:,2)
end if

! Keep track of the component of potential energy due to the pre-impact impactors%id for book-keeping
! Get the energy of the system after the collision
call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
dpe = self%pe(2) - self%pe(1)
nbody_system%Ecollisions = nbody_system%Ecollisions - dpe
nbody_system%Euntracked = nbody_system%Euntracked + dpe


! Update any encounter lists that have the removed bodies in them so that they instead point to the new body
do k = 1, nbody_system%plpl_encounter%nenc
do j = 1, impactors%ncoll
i = impactors%id(j)
if (i == ibiggest) cycle
if (nbody_system%plpl_encounter%id1(k) == pl%id(i)) then
nbody_system%plpl_encounter%id1(k) = pl%id(ibiggest)
nbody_system%plpl_encounter%index1(k) = i
end if
if (nbody_system%plpl_encounter%id2(k) == pl%id(i)) then
nbody_system%plpl_encounter%id2(k) = pl%id(ibiggest)
nbody_system%plpl_encounter%index2(k) = i
end if
if (nbody_system%plpl_encounter%id1(k) == nbody_system%plpl_encounter%id2(k)) nbody_system%plpl_encounter%status(k) = INACTIVE
end do
end do

self%status = MERGED

call collision_resolve_mergeaddsub(nbody_system, param, t, self%status)

end select
end associate
end select
return
end subroutine collision_generate_merge_any


module subroutine collision_generate_bounce_system(self, nbody_system, param, t)
implicit none
class(collision_bounce), intent(inout) :: self !! Bounce fragment system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! The time of the collision
end subroutine collision_generate_bounce_system

module subroutine collision_generate_simple_system(self, nbody_system, param, t)
implicit none
class(collision_simple), intent(inout) :: self !! Simple fragment system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! The time of the collision
end subroutine collision_generate_simple_system

end submodule s_collision_model
Loading

0 comments on commit 0896ffb

Please sign in to comment.