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

Commit

Permalink
More restructuring to put most of the basic infrastructure of generat…
Browse files Browse the repository at this point in the history 
…ing fragments into the simple disruption model
  • Loading branch information
@daminton
daminton committed Dec 23, 2022
1 parent 3350420 commit 9eaa799
Show file tree
Hide file tree
Showing 9 changed files with 326 additions and 379 deletions.
212 changes: 176 additions & 36 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,26 @@
use swiftest
contains

module subroutine collision_generate_basic(self, nbody_system, param, t)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Merge massive bodies no matter the regime
!!
!! 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_basic), 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

call self%merge(nbody_system, param, t)

return
end subroutine collision_generate_basic

module subroutine collision_generate_bounce(self, nbody_system, param, t)
!! author: David A. Minton
!!
Expand Down Expand Up @@ -55,9 +75,9 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t)
end select
end select
case (COLLRESOLVE_REGIME_HIT_AND_RUN)
call collision_generate_hitandrun(self, nbody_system, param, t)
call self%hitandrun(nbody_system, param, t)
case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE)
call self%collision_merge%generate(nbody_system, param, t) ! Use the default collision model, which is merge
call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge
case default
write(*,*) "Error in swiftest_collision, unrecognized collision regime"
call util_exit(FAILURE)
Expand All @@ -70,54 +90,178 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t)
end subroutine collision_generate_bounce


module subroutine collision_generate_hitandrun(collider, nbody_system, param, t)
module subroutine collision_generate_hitandrun(self, nbody_system, param, t)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Create the fragments resulting from a non-catastrophic hit-and-run collision
!!
implicit none
! Arguments
class(collision_merge), intent(inout) :: collider !! Fraggle collision system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
class(collision_basic), intent(inout) :: self !! Collision system object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
! Result
integer(I4B) :: status !! Status flag assigned to this outcome
! Internals
integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj
logical :: lpure
character(len=STRMAX) :: message
real(DP) :: dpe


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

associate(impactors => collider%impactors, fragments => collider%fragments, pl => nbody_system%pl)
associate(impactors => self%impactors)
message = "Hit and run between"
call collision_io_collider_message(nbody_system%pl, impactors%id, message)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message)))

collider%status = HIT_AND_RUN_PURE
pl%status(impactors%id(:)) = ACTIVE
pl%ldiscard(impactors%id(:)) = .false.
pl%lcollision(impactors%id(:)) = .false.

! Be sure to save the pl so that snapshots still work
select type(before => collider%before)
class is (swiftest_nbody_system)
select type(after => collider%after)
class is (swiftest_nbody_system)
allocate(after%pl, source=before%pl)
end select
end select
if (impactors%mass(1) > impactors%mass(2)) then
jtarg = 1
jproj = 2
else
jtarg = 2
jproj = 1
end if

! The simple disruption model (and its extended types allow for non-pure hit and run.
!For the basic merge model, all hit and runs are pure
select type(self)
class is (collision_simple_disruption)
if (impactors%mass_dist(2) > 0.9_DP * impactors%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Pure hit and run. No new fragments generated.")
nfrag = 0
lpure = .true.
else ! Imperfect hit and run, so we'll keep the largest body and destroy the other
lpure = .false.
call self%set_mass_dist(param)

! Generate the position and velocity distributions of the fragments
call self%disrupt(nbody_system, param, t, lpure)

dpe = self%pe(2) - self%pe(1)
nbody_system%Ecollisions = nbody_system%Ecollisions - dpe
nbody_system%Euntracked = nbody_system%Euntracked + dpe

if (lpure) then
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Should have been a pure hit and run instead")
nfrag = 0
else
nfrag = self%fragments%nbody
write(message, *) nfrag
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments")
end if
end if
class default
lpure = .true.
end select

if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them
status = HIT_AND_RUN_PURE
pl%status(impactors%id(:)) = ACTIVE
pl%ldiscard(impactors%id(:)) = .false.
pl%lcollision(impactors%id(:)) = .false.
! Be sure to save the pl so that snapshots still work
select type(before => self%before)
class is (swiftest_nbody_system)
select type(after => self%after)
class is (swiftest_nbody_system)
allocate(after%pl, source=before%pl)
end select
end select
else
ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1))
self%fragments%id(1) = pl%id(ibiggest)
self%fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)]
param%maxid = self%fragments%id(nfrag)
status = HIT_AND_RUN_DISRUPT
call collision_resolve_mergeaddsub(nbody_system, param, t, status)
end if
end associate
end select
end select

return
end subroutine collision_generate_hitandrun


module subroutine collision_generate_simple(self, nbody_system, param, t)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Create the fragments resulting from a non-catastrophic disruption collision
!!
implicit none
! Arguments
class(collision_simple_disruption), intent(inout) :: self
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions
real(DP), intent(in) :: t !! Time of collision
! Internals
integer(I4B) :: i, ibiggest, nfrag, status
logical :: lfailure
character(len=STRMAX) :: message
real(DP) :: dpe

select type(nbody_system)
class is (swiftest_nbody_system)
select type(pl => nbody_system%pl)
class is (swiftest_pl)
associate(impactors => self%impactors, fragments => self%fragments, status => self%status)

select case (impactors%regime)
case (COLLRESOLVE_REGIME_HIT_AND_RUN)
call self%hitandrun(nbody_system, param, t)
return
case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE)
call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge
return
case(COLLRESOLVE_REGIME_DISRUPTION)
message = "Disruption between"
case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
message = "Supercatastrophic disruption between"
case default
write(*,*) "Error in swiftest_collision, unrecognized collision regime"
call util_exit(FAILURE)
end select
call collision_io_collider_message(nbody_system%pl, impactors%id, message)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)

! Collisional fragments will be uniformly distributed around the pre-impact barycenter
call self%set_mass_dist(param)

! Generate the position and velocity distributions of the fragments
call self%disrupt(nbody_system, param, t)

dpe = self%pe(2) - self%pe(1)
nbody_system%Ecollisions = nbody_system%Ecollisions - dpe
nbody_system%Euntracked = nbody_system%Euntracked + dpe

! Populate the list of new bodies
nfrag = fragments%nbody
write(message, *) nfrag
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments")
select case(impactors%regime)
case(COLLRESOLVE_REGIME_DISRUPTION)
status = DISRUPTED
ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1))
fragments%id(1) = pl%id(ibiggest)
fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)]
param%maxid = fragments%id(nfrag)
case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
status = SUPERCATASTROPHIC
fragments%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
param%maxid = fragments%id(nfrag)
end select

call collision_resolve_mergeaddsub(nbody_system, param, t, status)
end associate
end select
end select
return
end subroutine collision_generate_hitandrun
end subroutine collision_generate_simple


module subroutine collision_generate_merge(self, nbody_system, param, t)
Expand All @@ -130,7 +274,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
!! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f
implicit none
! Arguments
class(collision_merge), intent(inout) :: self !! Merge fragment system object
class(collision_basic), 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
Expand Down Expand Up @@ -220,7 +364,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
end subroutine collision_generate_merge


module subroutine collision_generate_simple_disruption(self, nbody_system, param, t)
module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfailure)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Generates a simple fragment position and velocity distribution based on the collision
Expand All @@ -231,18 +375,16 @@ module subroutine collision_generate_simple_disruption(self, nbody_system, param
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
logical, optional, intent(out) :: lfailure
! Internals
real(DP) :: r_max_start

associate(impactors => self%impactors, fragments => self%fragments, nfrag => self%fragments%nbody)
r_max_start = 1.1_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1))
call collision_generate_simple_pos_vec(self, r_max_start)
call self%set_coordinate_system()
call collision_generate_simple_vel_vec(self)
end associate

r_max_start = 1.1_DP * .mag.(self%impactors%rb(:,2) - self%impactors%rb(:,1))
call collision_generate_simple_pos_vec(self, r_max_start)
call self%set_coordinate_system()
call collision_generate_simple_vel_vec(self)
return
end subroutine collision_generate_simple_disruption
end subroutine collision_generate_disrupt


module subroutine collision_generate_simple_pos_vec(collider, r_max_start)
Expand Down Expand Up @@ -371,6 +513,4 @@ module subroutine collision_generate_simple_vel_vec(collider)
end subroutine collision_generate_simple_vel_vec




end submodule s_collision_generate
Loading

0 comments on commit 9eaa799

Please sign in to comment.