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

Commit

Permalink
Fixed issue with rotation velocity calculation causing floating point…
Browse files Browse the repository at this point in the history 
… errors. Made other improvements, including setting collisions to merge if they gain energy.
  • Loading branch information
@daminton
daminton committed Jan 12, 2023
1 parent d4753a8 commit 9079568
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 71 deletions.
2 changes: 1 addition & 1 deletion src/base/base_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -346,7 +346,7 @@ subroutine base_util_resize_storage(self, nnew)
integer(I4B), intent(in) :: nnew !! New size
! Internals
class(base_storage_frame), dimension(:), allocatable :: tmp
integer(I4B) :: i, n, nold, nbig
integer(I4B) :: i, nold, nbig

nold = self%nframes
if (nnew <= nold) return
Expand Down
38 changes: 16 additions & 22 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,12 @@ module collision
integer(I4B), parameter :: COLLRESOLVE_REGIME_SUPERCATASTROPHIC = 3
integer(I4B), parameter :: COLLRESOLVE_REGIME_GRAZE_AND_MERGE = 4
integer(I4B), parameter :: COLLRESOLVE_REGIME_HIT_AND_RUN = 5
character(len=*),dimension(5), parameter :: REGIME_NAMES = ["Merge", "Disruption", "Supercatastrophic", "Graze and Merge", "Hit and Run"]
character(len=NAMELEN),parameter :: REGIME_NAME_MERGE = "Merge"
character(len=NAMELEN),parameter :: REGIME_NAME_DISRUPTION = "Disruption"
character(len=NAMELEN),parameter :: REGIME_NAME_SUPERCATASTROPHIC = "Supercatastrophic"
character(len=NAMELEN),parameter :: REGIME_NAME_GRAZE_AND_MERGE = "Graze and Merge"
character(len=NAMELEN),parameter :: REGIME_NAME_HIT_AND_RUN = "Hit and Run"
character(len=NAMELEN),dimension(5), parameter :: REGIME_NAMES = [REGIME_NAME_MERGE, REGIME_NAME_DISRUPTION, REGIME_NAME_SUPERCATASTROPHIC, REGIME_NAME_GRAZE_AND_MERGE, REGIME_NAME_HIT_AND_RUN]

!> Swiftest class for tracking pl-pl close encounters in a step when collisions are possible
type, extends(encounter_list) :: collision_list_plpl
Expand All @@ -52,14 +57,16 @@ module collision
!> Class definition for the variables that describe the bodies involved in the collision
type, extends(base_object) :: collision_impactors
integer(I4B) :: ncoll !! Number of bodies involved in the collision
integer(I4B), dimension(:), allocatable :: id !! Index of bodies involved in the collision
integer(I4B), dimension(:), allocatable :: id !! Index of bodies involved in the collision
real(DP), dimension(NDIM,2) :: rb !! Two-body equivalent position vectors of the collider bodies prior to collision in system barycentric coordinates
real(DP), dimension(NDIM,2) :: vb !! Two-body equivalent velocity vectors of the collider bodies prior to collision in system barycentric coordinate
real(DP), dimension(NDIM,2) :: rc !! Two-body equivalent position vectors of the collider bodies prior to collision in collision center of mass coordinates
real(DP), dimension(NDIM,2) :: vc !! Two-body equivalent velocity vectors of the collider bodies prior to collision in collision center of mass coordinates
real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision
real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision
real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision
real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision
real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision
real(DP), dimension(2) :: ke_orbit !! Orbital kinetic energy of each individual impactor
real(DP), dimension(2) :: ke_spin !! Spin kinetic energy of each individual impactors
real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision
real(DP), dimension(2) :: Gmass !! Two-body equivalent G*mass of the collider bodies prior to the collision
real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision
Expand Down Expand Up @@ -113,10 +120,10 @@ module collision
real(DP), dimension(NDIM,nbody) :: t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame
real(DP), dimension(NDIM,nbody) :: n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame
integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from
real(DP), dimension(NDIM) :: L_orbit_tot !! Orbital angular momentum vector of all fragments
real(DP), dimension(NDIM) :: L_spin_tot !! Spin angular momentum vector of all fragments
real(DP), dimension(NDIM,nbody) :: L_orbit !! Orbital angular momentum vector of each individual fragment
real(DP), dimension(NDIM,nbody) :: L_spin !! Spin angular momentum vector of each individual fragment
real(DP), dimension(NDIM) :: L_orbit_tot !! Orbital angular momentum vector of all fragments
real(DP), dimension(NDIM) :: L_spin_tot !! Spin angular momentum vector of all fragments
real(DP), dimension(NDIM,nbody) :: L_orbit !! Orbital angular momentum vector of each individual fragment
real(DP), dimension(NDIM,nbody) :: L_spin !! Spin angular momentum vector of each individual fragment
real(DP) :: ke_orbit_tot !! Orbital kinetic energy of all fragments
real(DP) :: ke_spin_tot !! Spin kinetic energy of all fragments
real(DP) :: pe !! Potential energy of all fragments
Expand All @@ -127,7 +134,6 @@ module collision
procedure :: dealloc => collision_util_dealloc_fragments !! Deallocates all allocatable arrays and sets everything else to 0
procedure :: reset => collision_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate)
procedure :: set_coordinate_system => collision_util_set_coordinate_fragments !! Sets the coordinate system of the fragments
final :: collision_final_fragments !! Finalizer deallocates all allocatables
end type collision_fragments


Expand Down Expand Up @@ -189,7 +195,7 @@ module collision
integer(I4B) :: stage_dimid !! ID for the stage dimension
integer(I4B) :: stage_varid !! ID for the stage variable
character(NAMELEN) :: stage_dimname = "stage" !! name of the stage dimension (before/after)
character(len=6), dimension(2) :: stage_coords = ["before", "after"] !! The stage coordinate labels
character(len=6), dimension(2) :: stage_coords = ["before", "after "] !! The stage coordinate labels

integer(I4B) :: collision_id_dimid !! ID for the collision event dimension

Expand Down Expand Up @@ -500,18 +506,6 @@ end subroutine collision_util_set_original_scale_factors

contains

subroutine collision_final_fragments(self)
!! author: David A. Minton
!!
!! Finalizer will deallocate all allocatables
implicit none
! Arguments
type(collision_fragments(*)), intent(inout) :: self

if (allocated(self%info)) deallocate(self%info)
return
end subroutine collision_final_fragments

subroutine collision_final_impactors(self)
!! author: David A. Minton
!!
Expand Down
4 changes: 3 additions & 1 deletion src/collision/collision_resolve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -509,7 +509,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
real(DP) :: E_before, E_after, dLmag
real(DP), dimension(NDIM) :: L_orbit_before, L_orbit_after, L_spin_before, L_spin_after, L_before, L_after, dL_orbit, dL_spin, dL
logical :: lplpl_collision
character(len=STRMAX) :: timestr
character(len=STRMAX) :: timestr, idstr
integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision
logical :: lgoodcollision
integer(I4B) :: i, loop, ncollisions
Expand Down Expand Up @@ -562,6 +562,8 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
param%maxid_collision = max(param%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id))
param%maxid_collision = param%maxid_collision + 1
collider%collision_id = param%maxid_collision
write(idstr,*) collider%collision_id
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Collision id " // trim(adjustl(idstr)))

! Get the collision regime
call impactors%get_regime(nbody_system, param)
Expand Down
15 changes: 14 additions & 1 deletion src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,7 +233,14 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par
self%pe(phase_val) = constraint_system%pe / self%Escale
self%be(phase_val) = constraint_system%be / self%Escale
self%te(phase_val) = constraint_system%te / self%Escale
if (phase_val == 2) then
if (phase_val == 1) then
do concurrent(i = 1:2)
impactors%ke_orbit(i) = 0.5_DP * impactors%mass(i) * dot_product(impactors%vc(:,i), impactors%vc(:,i))
impactors%ke_spin(i) = 0.5_DP * impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * dot_product(impactors%rot(:,i), impactors%rot(:,i))
impactors%L_orbit(:,i) = impactors%mass(i) * impactors%rc(:,i) .cross. impactors%vc(:,i)
impactors%L_spin(:,i) = impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * impactors%rot(:,i)
end do
else if (phase_val == 2) then
do concurrent(i = 1:nfrag)
fragments%ke_orbit(i) = 0.5_DP * fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i))
fragments%ke_spin(i) = 0.5_DP * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i), fragments%rot(:,i))
Expand Down Expand Up @@ -315,6 +322,8 @@ module subroutine collision_util_dealloc_impactors(self)
self%rot(:,:) = 0.0_DP
self%L_spin(:,:) = 0.0_DP
self%L_orbit(:,:) = 0.0_DP
self%ke_spin(:) = 0.0_DP
self%ke_orbit(:) = 0.0_DP
self%Ip(:,:) = 0.0_DP
self%mass(:) = 0.0_DP
self%radius(:) = 0.0_DP
Expand Down Expand Up @@ -810,6 +819,8 @@ module subroutine collision_util_set_natural_scale_factors(self)
impactors%radius(:) = impactors%radius(:) / collider%dscale
impactors%L_spin(:,:) = impactors%L_spin(:,:) / collider%Lscale
impactors%L_orbit(:,:) = impactors%L_orbit(:,:) / collider%Lscale
impactors%ke_orbit(:) = impactors%ke_orbit(:) / collider%Escale
impactors%ke_spin(:) = impactors%ke_spin(:) / collider%Escale

do concurrent(i = 1:2)
impactors%rot(:,i) = impactors%L_spin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i))
Expand Down Expand Up @@ -859,6 +870,8 @@ module subroutine collision_util_set_original_scale_factors(self)
impactors%vc = impactors%vc * collider%vscale
impactors%L_spin = impactors%L_spin * collider%Lscale
impactors%L_orbit = impactors%L_orbit * collider%Lscale
impactors%ke_orbit = impactors%ke_orbit * collider%Escale
impactors%ke_spin = impactors%ke_spin * collider%Escale
do concurrent(i = 1:2)
impactors%rot(:,i) = impactors%L_spin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i))
end do
Expand Down
Loading

0 comments on commit 9079568

Please sign in to comment.