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

Commit

Permalink
Cleanup and fixed up some of the coordinate system vector calcs
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 23, 2022
1 parent 13ae3bb commit bed6564
Show file tree
Hide file tree
Showing 4 changed files with 81 additions and 77 deletions.
120 changes: 60 additions & 60 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,64 @@
use swiftest
contains

module subroutine collision_generate_bounce(self, nbody_system, param, t)
!! author: David A. Minton
!!
!! In this collision model, if the collision would result in a disruption, the bodies are instead "bounced" off
!! of the center of mass. This is done as a reflection in the 2-body equivalent distance vector direction.
implicit none
! Arguments
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
! Internals
integer(I4B) :: i,j,nfrag
real(DP), dimension(NDIM) :: vcom, rnorm

select type(nbody_system)
class is (swiftest_nbody_system)
select type (pl => nbody_system%pl)
class is (swiftest_pl)
associate(impactors => nbody_system%collider%impactors, fragments => nbody_system%collider%fragments)
select case (impactors%regime)
case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
nfrag = size(impactors%id(:))
do i = 1, nfrag
j = impactors%id(i)
vcom(:) = pl%vb(:,j) - impactors%vbcom(:)
rnorm(:) = .unit. (impactors%rb(:,2) - impactors%rb(:,1))
! Do the reflection
vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:)
pl%vb(:,j) = impactors%vbcom(:) + vcom(:)
self%status = DISRUPTED
pl%status(j) = ACTIVE
pl%ldiscard(j) = .false.
pl%lcollision(j) = .false.
end do
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) ! Be sure to save the pl so that snapshots still work
end select
end select
case (COLLRESOLVE_REGIME_HIT_AND_RUN)
call collision_generate_hitandrun(self, 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
case default
write(*,*) "Error in swiftest_collision, unrecognized collision regime"
call util_exit(FAILURE)
end select
end associate
end select
end select

return
end subroutine collision_generate_bounce


module subroutine collision_generate_hitandrun(collider, nbody_system, param, t)
!! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
Expand Down Expand Up @@ -159,70 +217,12 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
end subroutine collision_generate_merge


module subroutine collision_generate_bounce(self, nbody_system, param, t)
!! author: David A. Minton
!!
!! In this collision model, if the collision would result in a disruption, the bodies are instead "bounced" off
!! of the center of mass. This is done as a reflection in the 2-body equivalent distance vector direction.
implicit none
! Arguments
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
! Internals
integer(I4B) :: i,j,nfrag
real(DP), dimension(NDIM) :: vcom, rnorm

select type(nbody_system)
class is (swiftest_nbody_system)
select type (pl => nbody_system%pl)
class is (swiftest_pl)
associate(impactors => nbody_system%collider%impactors, fragments => nbody_system%collider%fragments)
select case (impactors%regime)
case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
nfrag = size(impactors%id(:))
do i = 1, nfrag
j = impactors%id(i)
vcom(:) = pl%vb(:,j) - impactors%vbcom(:)
rnorm(:) = .unit. (impactors%rb(:,2) - impactors%rb(:,1))
! Do the reflection
vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:)
pl%vb(:,j) = impactors%vbcom(:) + vcom(:)
self%status = DISRUPTED
pl%status(j) = ACTIVE
pl%ldiscard(j) = .false.
pl%lcollision(j) = .false.
end do
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) ! Be sure to save the pl so that snapshots still work
end select
end select
case (COLLRESOLVE_REGIME_HIT_AND_RUN)
call collision_generate_hitandrun(self, 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
case default
write(*,*) "Error in swiftest_collision, unrecognized collision regime"
call util_exit(FAILURE)
end select
end associate
end select
end select

return
end subroutine collision_generate_bounce


module subroutine collision_generate_simple(self, nbody_system, param, t)
module subroutine collision_generate_simple_disruption(self, nbody_system, param, t)
implicit none
class(collision_simple_disruption), 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
end subroutine collision_generate_simple_disruption

end submodule s_collision_generate
21 changes: 11 additions & 10 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,13 +66,14 @@ module collision
real(DP) :: Mcb !! Mass of central body (used to compute potential energy in regime determination)

! Values in a coordinate frame centered on the collider barycenter and collisional nbody_system unit vectors
real(DP), dimension(NDIM) :: x_unit !! x-direction unit vector of collisional nbody_system
real(DP), dimension(NDIM) :: y_unit !! y-direction unit vector of collisional nbody_system
real(DP), dimension(NDIM) :: z_unit !! z-direction unit vector of collisional nbody_system
real(DP), dimension(NDIM) :: v_unit !! z-direction unit vector of collisional nbody_system
real(DP), dimension(NDIM) :: rbcom !! Center of mass position vector of the collider nbody_system in nbody_system barycentric coordinates
real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider nbody_system in nbody_system barycentric coordinates
real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider nbody_system in nbody_system barycentric coordinates
real(DP), dimension(NDIM) :: x_unit !! x-direction unit vector of collisional nbody_system
real(DP), dimension(NDIM) :: y_unit !! y-direction unit vector of collisional nbody_system
real(DP), dimension(NDIM) :: z_unit !! z-direction unit vector of collisional nbody_system
real(DP), dimension(NDIM) :: v_unit !! velocity direction unit vector of collisional nbody_system
real(DP), dimension(NDIM) :: rbcom !! Center of mass position vector of the collider nbody_system in nbody_system barycentric coordinates
real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider nbody_system in nbody_system barycentric coordinates
real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider nbody_system in nbody_system barycentric coordinates
real(DP), dimension(NDIM) :: vbimp !! The impact point velocity vector is the component of the velocity in the distance vector direction

contains
procedure :: consolidate => collision_resolve_consolidate_impactors !! Consolidates a multi-body collision into an equivalent 2-body collision
Expand Down Expand Up @@ -149,7 +150,7 @@ module collision

type, extends(collision_merge) :: collision_simple_disruption
contains
procedure :: generate => collision_generate_simple !! If a collision would result in a disruption [TODO: SOMETHING LIKE CHAMBERS 2012]
procedure :: generate => collision_generate_simple_disruption !! If a collision would result in a disruption [TODO: SOMETHING LIKE CHAMBERS 2012]
procedure :: set_mass_dist => collision_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type
final :: collision_final_simple !! Finalizer will deallocate all allocatables
end type collision_simple_disruption
Expand Down Expand Up @@ -223,13 +224,13 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t)
real(DP), intent(in) :: t !! The time of the collision
end subroutine collision_generate_bounce

module subroutine collision_generate_simple(self, nbody_system, param, t)
module subroutine collision_generate_simple_disruption(self, nbody_system, param, t)
implicit none
class(collision_simple_disruption), intent(inout) :: self !! Simple 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_simple
end subroutine collision_generate_simple_disruption

module subroutine collision_io_collider_message(pl, collidx, collider_message)
implicit none
Expand Down
11 changes: 7 additions & 4 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ module subroutine collision_util_construct_temporary_system(self, nbody_system,
! Set up a new system based on the original
if (allocated(tmpparam)) deallocate(tmpparam)
if (allocated(tmpsys)) deallocate(tmpsys)
allocate(tmpparam, source=param)
allocate(tmpparam_local, source=param)
call swiftest_util_setup_construct_system(tmpsys_local, tmpparam_local)

! No test particles necessary for energy/momentum calcs
Expand Down Expand Up @@ -411,10 +411,13 @@ module subroutine collision_util_set_coordinate_system(self)

! The cross product of the y- by z-axis will give us the x-axis
impactors%x_unit(:) = impactors%y_unit(:) .cross. impactors%z_unit(:)

impactors%v_unit(:) = .unit.delta_v(:)

if (.not.any(fragments%rc(:,:) > 0.0_DP)) return
! The "bounce" unit vector is the projection of
impactors%vbimp(:) = dot_product(delta_v(:),impactors%x_unit(:)) * impactors%x_unit(:)

if ((.not.allocated(self%fragments)) .or. (nfrag == 0) .or. (.not.any(fragments%rc(:,:) > 0.0_DP))) return

fragments%rmag(:) = .mag. fragments%rc(:,:)

! Randomize the tangential velocity direction.
Expand Down Expand Up @@ -623,7 +626,7 @@ module subroutine collision_util_setup_fragments_system(self, nfrag)
return
end subroutine collision_util_setup_fragments_system


subroutine collision_util_save_snapshot(collision_history, snapshot)
!! author: David A. Minton
!!
Expand Down
6 changes: 3 additions & 3 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,9 +252,6 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai
return
end if

! This is a factor that will "distort" the shape of the frgment cloud in the direction of the impact velocity
f_spin= .mag. (runit(:) .cross. vunit(:))

if (param%lflatten_interactions) then
lk_plpl = allocated(pl%k_plpl)
if (lk_plpl) deallocate(pl%k_plpl)
Expand Down Expand Up @@ -288,6 +285,9 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai
call fraggle_generate_pos_vec(collider, r_max_start)
call collider%set_coordinate_system()

! This is a factor that will "distort" the shape of the fragment cloud in the direction of the impact velocity
f_spin= .mag. (impactors%y_unit(:) .cross. impactors%v_unit(:))

! Initial velocity guess will be the barycentric velocity of the colliding nbody_system so that the budgets are based on the much smaller collisional-frame velocities
do concurrent (i = 1:nfrag)
fragments%vb(:, i) = impactors%vbcom(:)
Expand Down

0 comments on commit bed6564

Please sign in to comment.