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

Commit

Permalink
Restructured to try to simplify the code used to resolve collisions
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 22, 2022
1 parent f63621d commit a5ac1c4
Show file tree
Hide file tree
Showing 4 changed files with 233 additions and 283 deletions.
53 changes: 24 additions & 29 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,22 +13,6 @@
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

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.
Expand All @@ -45,7 +29,7 @@ module subroutine collision_generate_merge_any(self, nbody_system, param, t)
! Internals
integer(I4B) :: i, j, k, ibiggest
real(DP), dimension(NDIM) :: Lspin_new
real(DP) :: dpe
real(DP) :: volume, dpe
character(len=STRMAX) :: message

select type(nbody_system)
Expand All @@ -57,35 +41,46 @@ module subroutine collision_generate_merge_any(self, nbody_system, param, t)

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

!call self%set_mass_dist(param)
! Generate the merged body as a single fragment
call self%setup_fragments(1)

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

! The new body's metadata will be taken from the largest of the two impactor bodies, so we need
! its index in the main pl structure
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(:)
allocate(fragments%info, source=pl%info(ibiggest:ibiggest))

! Compute the physical properties of the new body after the merge.
volume = 4._DP / 3._DP * PI * sum(impactors%radius(:)**3)
fragments%mass(1) = impactors%mass_dist(1)
fragments%density(1) = fragments%mass(1) / volume
fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD)
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
do concurrent(i = 1:NDIM)
fragments%Ip(i,1) = sum(impactors%mass(:) * impactors%Ip(i,:))
Lspin_new(i) = sum(impactors%Lorbit(i,:) + impactors%Lorbit(i,:))
end do
fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1)
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
! The fragment trajectory will be the barycentric trajectory
fragments%rb(:,1) = impactors%rbcom(:)
fragments%vb(:,1) = impactors%vbcom(:)

! 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)

! Keep track of the component of potential energy that is now not considered because two bodies became one
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
Expand All @@ -111,7 +106,7 @@ module subroutine collision_generate_merge_any(self, nbody_system, param, t)
end associate
end select
return
end subroutine collision_generate_merge_any
end subroutine collision_generate_merge_system


module subroutine collision_generate_bounce_system(self, nbody_system, param, t)
Expand Down
73 changes: 18 additions & 55 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,9 +75,10 @@ module collision
real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider nbody_system in nbody_system barycentric coordinates

contains
procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be
procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions
final :: collision_final_impactors !! Finalizer will deallocate all allocatables
procedure :: consolidate => collision_resolve_consolidate_impactors !! Consolidates a multi-body collision into an equivalent 2-body collision
procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be
procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions
final :: collision_final_impactors !! Finalizer will deallocate all allocatables
end type collision_impactors


Expand Down Expand Up @@ -109,9 +110,11 @@ module collision
end type collision_fragments


type, abstract :: collision_system
type :: collision_system
!! This class defines a collisional nbody_system that stores impactors and fragments. This is written so that various collision models (i.e. Fraggle) could potentially be used
!! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments
!!
!! The generate method for this class is the merge model. This allows any extended type to have access to the merge procedure by selecting the collision_system parent class
class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system
class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system
class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision
Expand All @@ -135,15 +138,9 @@ module collision
procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.)
procedure :: reset => collision_util_reset_system !! Deallocates all allocatables
procedure :: set_coordinate_system => collision_util_set_coordinate_system !! Sets the coordinate nbody_system of the collisional nbody_system
procedure(abstract_generate_system), deferred :: generate !! Generates a nbody_system of fragments depending on collision model
procedure :: generate => collision_generate_merge_system !! Merges the impactors to make a single final body
end type collision_system

type, extends(collision_system) :: collision_merge
contains
procedure :: generate => collision_generate_merge_system !! Merges the impactors to make a single final body
final :: collision_final_merge_system !! Finalizer will deallocate all allocatables
end type collision_merge

type, extends(collision_system) :: collision_bounce
contains
procedure :: generate => collision_generate_bounce_system !! If a collision would result in a disruption, "bounce" the bodies instead.
Expand Down Expand Up @@ -199,37 +196,10 @@ module collision
end type collision_storage


abstract interface
subroutine abstract_generate_system(self, nbody_system, param, t)
import collision_system, base_nbody_system, base_parameters, DP
implicit none
class(collision_system), 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
real(DP), intent(in) :: t !! The time of the collision
end subroutine abstract_generate_system

subroutine abstract_set_mass_dist(self, param)
import collision_system, base_parameters
implicit none
class(collision_system), intent(inout) :: self !! Collision system object
class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters
end subroutine abstract_set_mass_dist
end interface


interface
module subroutine collision_generate_merge_any(self, nbody_system, param, t)
implicit none
class(collision_system), intent(inout) :: self !! Merge fragment nbody_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_any

module subroutine collision_generate_merge_system(self, nbody_system, param, t)
implicit none
class(collision_merge), intent(inout) :: self !! Merge fragment nbody_system object
class(collision_system), intent(inout) :: self !! Merge fragment nbody_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 @@ -310,6 +280,15 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l
integer(I4B), intent(in) :: irec !! Current recursion level
logical, intent(out) :: lany_collision !! Returns true if any pair of encounters resulted in a collision
end subroutine collision_check_pltp

module subroutine collision_resolve_consolidate_impactors(self, nbody_system, param, idx_parent, lflag)
implicit none
class(collision_impactors), intent(out) :: self !! Collision impactors object
class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(base_parameters), intent(in) :: param !! Current run configuration parameters with Swiftest additions
integer(I4B), dimension(:), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision
logical, intent(out) :: lflag !! Logical flag indicating whether a impactors%id was successfully created or not
end subroutine collision_resolve_consolidate_impactors

module subroutine collision_resolve_extract_plpl(self, nbody_system, param)
implicit none
Expand Down Expand Up @@ -544,22 +523,6 @@ subroutine collision_final_storage(self)
end subroutine collision_final_storage


subroutine collision_final_merge_system(self)
!! author: David A. Minton
!!
!! Finalizer will deallocate all allocatables
implicit none
! Arguments
type(collision_merge), intent(inout) :: self !! Collision system object

call self%reset()
if (allocated(self%impactors)) deallocate(self%impactors)
if (allocated(self%fragments)) deallocate(self%fragments)

return
end subroutine collision_final_merge_system


subroutine collision_final_bounce_system(self)
!! author: David A. Minton
!!
Expand Down
Loading

0 comments on commit a5ac1c4

Please sign in to comment.