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

Commit

Permalink
Made a number of changes to improve the robustness of Fraggle
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 31, 2022
1 parent af0a15e commit f1fa891
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 31 deletions.
16 changes: 10 additions & 6 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,10 +122,11 @@ module collision
real(DP), dimension(nbody) :: ke_orbit !! Orbital kinetic energy of each individual fragment
real(DP), dimension(nbody) :: ke_spin !! Spin kinetic energy of each individual fragment
contains
procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0
procedure :: get_angular_momentum => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments
procedure :: get_kinetic_energy => collision_util_get_kinetic_energy !! Calcualtes the current kinetic energy of the fragments
final :: collision_final_fragments !! Finalizer deallocates all allocatables
procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0
procedure :: get_angular_momentum => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments
procedure :: get_kinetic_energy => collision_util_get_kinetic_energy !! Calcualtes the current kinetic energy of the fragments
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 @@ -231,7 +232,6 @@ 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_hitandrun(self, nbody_system, param, t)
implicit none
class(collision_basic), intent(inout) :: self
Expand All @@ -247,7 +247,6 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
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


module subroutine collision_io_collider_message(pl, collidx, collider_message)
implicit none
Expand Down Expand Up @@ -406,6 +405,11 @@ module subroutine collision_util_set_coordinate_collider(self)
class(collision_basic), intent(inout) :: self !! collisional system
end subroutine collision_util_set_coordinate_collider

module subroutine collision_util_set_coordinate_fragments(self)
implicit none
class(collision_fragments(*)), intent(inout) :: self !! Collisional nbody_system
end subroutine collision_util_set_coordinate_fragments

module subroutine collision_util_set_coordinate_impactors(self)
implicit none
class(collision_impactors), intent(inout) :: self !! collisional system
Expand Down
23 changes: 20 additions & 3 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -449,22 +449,39 @@ module subroutine collision_util_set_coordinate_collider(self)
call impactors%set_coordinate_system()

if (.not.allocated(self%fragments)) return
call fragments%set_coordinate_system()


end associate

return
end subroutine collision_util_set_coordinate_collider


module subroutine collision_util_set_coordinate_fragments(self)
!! author: David A. Minton
!!
!! Defines the collisional coordinate nbody_system, including the unit vectors of both the nbody_system and individual fragments.
implicit none
! Arguments
class(collision_fragments(*)), intent(inout) :: self !! Collisional nbody_system

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

fragments%rmag(:) = .mag. fragments%rc(:,:)
fragments%vmag(:) = .mag. fragments%vc(:,:)
fragments%rotmag(:) = .mag. fragments%rot(:,:)

! Define the radial, normal, and tangential unit vectors for each individual fragment
fragments%r_unit(:,:) = .unit. fragments%rc(:,:)
fragments%v_unit(:,:) = .unit. fragments%vc(:,:)
fragments%n_unit(:,:) = .unit. (fragments%rc(:,:) .cross. fragments%vc(:,:))
fragments%t_unit(:,:) = -.unit. (fragments%r_unit(:,:) .cross. fragments%n_unit(:,:))

end associate

return
end subroutine collision_util_set_coordinate_collider
end subroutine collision_util_set_coordinate_fragments


module subroutine collision_util_set_coordinate_impactors(self)
Expand Down
49 changes: 27 additions & 22 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
logical :: lk_plpl, lfailure_local
logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes
character(len=STRMAX) :: message
real(DP), parameter :: fail_scale_initial = 1.001_DP
real(DP), parameter :: fail_scale_initial = 1.01_DP

! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we
! can simply fail the attempt and try again. So we need to turn off any floating point exception halting modes temporarily
Expand Down Expand Up @@ -136,23 +136,16 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
! Compute the "before" energy/momentum and the budgets
call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)
call self%set_budgets()
do try = 1, MAXTRY
self%fail_scale = (fail_scale_initial)**try
write(message,*) try
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message)))
call fraggle_generate_pos_vec(self)
call fraggle_generate_rot_vec(self)
call fraggle_generate_vel_vec(self,lfailure_local)
if (.not.lfailure_local) exit
end do
self%fail_scale = fail_scale_initial
call fraggle_generate_pos_vec(self)
call fraggle_generate_rot_vec(self)
call fraggle_generate_vel_vec(self,lfailure_local)
call self%set_original_scale()

if (lfailure_local) then
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle fragment generation failed after " // &
trim(adjustl(message)) // " tries. This collision may have gained energy.")
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "This collision may have gained energy.")
else
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle fragment generation succeeded after " // &
trim(adjustl(message)) // " tries")
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle fragment generation succeeded.")
end if


Expand Down Expand Up @@ -253,8 +246,8 @@ module subroutine fraggle_generate_pos_vec(collider)
logical, dimension(collider%fragments%nbody) :: loverlap
integer(I4B) :: i, j, loop
logical :: lcat, lhitandrun
integer(I4B), parameter :: MAXLOOP = 10000
real(DP), parameter :: rdistance_scale_factor = 0.01_DP ! Scale factor to apply to distance scaling in the event of a fail
integer(I4B), parameter :: MAXLOOP = 20000
real(DP), parameter :: rdistance_scale_factor = 0.1_DP ! Scale factor to apply to distance scaling in the event of a fail

associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
Expand All @@ -269,6 +262,8 @@ module subroutine fraggle_generate_pos_vec(collider)
if (.not.any(loverlap(:))) exit
fragment_cloud_center(:,1) = impactors%rc(:,1) - rdistance(1) * impactors%bounce_unit(:)
fragment_cloud_center(:,2) = impactors%rc(:,2) + rdistance(2) * impactors%bounce_unit(:)
fragment_cloud_radius(1) = .mag.(fragment_cloud_center(:,1) - impactors%rbimp(:))
fragment_cloud_radius(2) = .mag.(fragment_cloud_center(:,2) - impactors%rbimp(:))
do concurrent(i = 1:nfrag, loverlap(i))
if (i < 3) then
fragments%rc(:,i) = fragment_cloud_center(:,i)
Expand Down Expand Up @@ -298,7 +293,6 @@ module subroutine fraggle_generate_pos_vec(collider)
end do
end do
rdistance(:) = rdistance(:) * collider%fail_scale
fragment_cloud_radius(:) = fragment_cloud_radius(:) * collider%fail_scale
end do

call collision_util_shift_vector_to_origin(fragments%mass, fragments%rc)
Expand Down Expand Up @@ -375,16 +369,19 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
integer(I4B) :: i, j, loop, istart, n, ndof
logical :: lhitandrun, lnoncat
real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear, vunit
real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof, ke_tot
real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof, ke_tot, ke_residual_min
integer(I4B), dimension(collider%fragments%nbody) :: vsign
real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail
integer(I4B), parameter :: MAXLOOP = 100
integer(I4B), parameter :: MAXLOOP = 2000
real(DP), parameter :: TOL = 1e-2
class(collision_fragments(:)), allocatable :: fragments

associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
associate(impactors => collider%impactors, nfrag => collider%fragments%nbody)
lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN)
lnoncat = (impactors%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point

allocate(fragments, source=collider%fragments)

! The fragments will be divided into two "clouds" based on identified origin body.
! These clouds will collectively travel like two impactors bouncing off of each other.
where(fragments%origin_body(:) == 1)
Expand Down Expand Up @@ -414,6 +411,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence
mass_vscale(:) = mass_vscale(:) / minval(mass_vscale(:))


! Set the velocities of all fragments using all of the scale factors determined above
do concurrent(i = 1:nfrag)
j = fragments%origin_body(i)
Expand All @@ -439,10 +437,17 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
else
istart = 1
end if
ke_residual_min = huge(1.0_DP)
do loop = 1, MAXLOOP
call collider%set_coordinate_system()
call fragments%set_coordinate_system()
call fragments%get_kinetic_energy()
ke_residual = fragments%ke_budget - (fragments%ke_orbit_tot + fragments%ke_spin_tot)
if (abs(ke_residual) < abs(ke_residual_min)) then ! This is our best case so far. Save it for posterity
if (allocated(collider%fragments)) deallocate(collider%fragments)
allocate(collider%fragments, source=fragments)
ke_residual_min = ke_residual
end if
ke_residual_min = max(ke_residual_min,ke_residual)
if (ke_residual > 0.0_DP) exit
! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape
ke_avail(:) = fragments%ke_orbit(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%rmag(:)
Expand Down Expand Up @@ -475,7 +480,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
end do

! Check for any residual angular momentum, and if there is any, put it into spin
call collider%set_coordinate_system()
call fragments%set_coordinate_system()
call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:))
do concurrent(i = istart:nfrag)
Expand Down

0 comments on commit f1fa891

Please sign in to comment.