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

Commit

Permalink
Updated Fraggle to help with convergence and deal with overlaps with …
Browse files Browse the repository at this point in the history 
…bodies 1 and 2.
  • Loading branch information
@daminton
daminton committed Feb 2, 2023
1 parent 7e67f22 commit 512b083
Show file tree
Hide file tree
Showing 4 changed files with 62 additions and 33 deletions.
2 changes: 1 addition & 1 deletion src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ module collision
real(DP), dimension(NDIM) :: v_unit !! velocity direction unit vector of collisional 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) :: rcimp !! Impact point position vector of the collider nbody_system in nbody_system barycentric coordinates
real(DP), dimension(NDIM) :: bounce_unit !! The impact point velocity vector is the component of the velocity in the distance vector direction

contains
Expand Down
10 changes: 5 additions & 5 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -359,7 +359,7 @@ module subroutine collision_util_dealloc_impactors(self)
self%v_unit(:) = 0.0_DP
self%rbcom(:) = 0.0_DP
self%vbcom(:) = 0.0_DP
self%rbimp(:) = 0.0_DP
self%rcimp(:) = 0.0_DP

return
end subroutine collision_util_dealloc_impactors
Expand Down Expand Up @@ -690,8 +690,8 @@ module subroutine collision_util_set_coordinate_impactors(self)
impactors%vc(:,1) = impactors%vb(:,1) - impactors%vbcom(:)
impactors%vc(:,2) = impactors%vb(:,2) - impactors%vbcom(:)

! Find the point of impact between the two bodies
impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * impactors%y_unit(:) - impactors%rbcom(:)
! Find the point of impact between the two bodies, defined as the location (in the collisional coordinate system) at the surface of body 1 along the line connecting the two bodies.
impactors%rcimp(:) = impactors%rb(:,1) + impactors%radius(1) * impactors%y_unit(:) - impactors%rbcom(:)

! Set the velocity direction as the "bounce" direction" for disruptions, and body 2's direction for hit and runs
if (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) then
Expand Down Expand Up @@ -921,7 +921,7 @@ module subroutine collision_util_set_natural_scale_factors(self)
! Scale all dimensioned quantities of impactors and fragments
impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale
impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale
impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale
impactors%rcimp(:) = impactors%rcimp(:) / collider%dscale
impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale
impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale
impactors%rc(:,:) = impactors%rc(:,:) / collider%dscale
Expand Down Expand Up @@ -973,7 +973,7 @@ module subroutine collision_util_set_original_scale_factors(self)
! Restore scale factors
impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale
impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale
impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale
impactors%rcimp(:) = impactors%rcimp(:) * collider%dscale

impactors%mass = impactors%mass * collider%mscale
impactors%Gmass(:) = impactors%Gmass(:) * (collider%dscale**3/collider%tscale**2)
Expand Down
80 changes: 54 additions & 26 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,9 +159,11 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
lfailure = .false.
call self%get_energy_and_momentum(nbody_system, param, phase="before")
self%fail_scale = fail_scale_initial
call fraggle_generate_pos_vec(self, nbody_system, param)
call fraggle_generate_rot_vec(self, nbody_system, param)
call fraggle_generate_vel_vec(self, nbody_system, param, lfailure)
call fraggle_generate_pos_vec(self, nbody_system, param, lfailure)
if (.not.lfailure) then
call fraggle_generate_rot_vec(self, nbody_system, param)
call fraggle_generate_vel_vec(self, nbody_system, param, lfailure)
end if

if (.not.lfailure) then
if (self%fragments%nbody /= nfrag_start) then
Expand Down Expand Up @@ -280,7 +282,7 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t)
end subroutine fraggle_generate_hitandrun


module subroutine fraggle_generate_pos_vec(collider, nbody_system, param)
module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailure)
!! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
!! Initializes the position vectors of the fragments around the center of mass based on the collision style.
Expand All @@ -292,6 +294,7 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param)
class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
logical, intent(out) :: lfailure !! Did the velocity computation fail?
! Internals
real(DP) :: dis, direction, rdistance
real(DP), dimension(NDIM,2) :: fragment_cloud_center
Expand All @@ -300,12 +303,13 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param)
real(DP), dimension(collider%fragments%nbody) :: mass_rscale, phi, theta, u
integer(I4B) :: i, j, loop, istart
logical :: lsupercat, lhitandrun
integer(I4B), parameter :: MAXLOOP = 20000
integer(I4B), parameter :: MAXLOOP = 100
real(DP), parameter :: rdistance_scale_factor = 1.0_DP ! Scale factor to apply to distance scaling of cloud centers in the event of overlap
! The distance is chosen to be close to the original locations of the impactors
! but far enough apart to prevent a collisional cascade between fragments
real(DP), parameter :: cloud_size_scale_factor = 3.0_DP ! Scale factor to apply to the size of the cloud relative to the distance from the impact point.
! A larger value puts more space between fragments initially


associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody, &
pl => nbody_system%pl, tp => nbody_system%tp, npl => nbody_system%pl%nbody, ntp => nbody_system%tp%nbody)
Expand All @@ -328,30 +332,31 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param)
mass_rscale(:) = (mass_rscale(:) + 1.0_DP) / 2
mass_rscale(:) = mass_rscale(:) * (fragments%mtot / fragments%mass(1:nfrag))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence
mass_rscale(:) = mass_rscale(:) / maxval(mass_rscale(:))
istart = 3

do loop = 1, MAXLOOP
if (.not.any(loverlap(:))) exit
if (lhitandrun) then
if (lhitandrun) then ! Keep the target unchanged and place the largest fragment at rdistance away from the projectile along its trajectory
fragment_cloud_radius(:) = impactors%radius(:)
fragment_cloud_center(:,1) = impactors%rc(:,1)
fragment_cloud_center(:,2) = impactors%rc(:,2) + rdistance * impactors%bounce_unit(:)
else
fragment_cloud_center(:,1) = impactors%rbimp(:) - 1.001_DP * max(fragments%radius(1),impactors%radius(1)) * impactors%y_unit(:)
fragment_cloud_center(:,2) = impactors%rbimp(:) + 1.001_DP * max(fragments%radius(2),impactors%radius(2)) * impactors%y_unit(:)
else ! Keep the first and second bodies at approximately their original location, but so as not to be overlapping
fragment_cloud_center(:,1) = impactors%rc(:,1)
fragment_cloud_center(:,2) = impactors%rcimp(:) + 1.001_DP * max(fragments%radius(2),impactors%radius(2)) * impactors%y_unit(:)
fragment_cloud_radius(:) = cloud_size_scale_factor * rdistance
end if
fragments%rc(:,1) = fragment_cloud_center(:,1)
fragments%rc(:,2) = fragment_cloud_center(:,2)
istart = 3

do i = istart, nfrag
do i = 1, nfrag
if (loverlap(i)) then
call random_number(phi(i))
call random_number(theta(i))
call random_number(u(i))
end if
end do

! Randomly place the n>2 fragments inside their cloud until none are overlapping
do concurrent(i = istart:nfrag, loverlap(i))
j = fragments%origin_body(i)

Expand All @@ -369,13 +374,25 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param)
fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j)

! Make sure that the fragments are positioned away from the impact point
direction = dot_product(fragments%rc(:,i) - impactors%rbimp(:), fragment_cloud_center(:,j) - impactors%rbimp(:))
direction = dot_product(fragments%rc(:,i) - impactors%rcimp(:), fragment_cloud_center(:,j) - impactors%rcimp(:))
if (direction < 0.0_DP) then
fragments%rc(:,i) = fragments%rc(:,i) - fragment_cloud_center(:,j)
fragments%rc(:,i) = -fragments%rc(:,i) + fragment_cloud_center(:,j)
end if

end do

! Because body 1 and 2 are initialized near the original impactor positions, then if these bodies are still overlapping
! when the rest are not, we will randomly walk their position in space so as not to move them too far from their starting position
if (all(.not.loverlap(istart:nfrag)) .and. any(loverlap(1:istart-1))) then
do concurrent(i = 1:istart-1,loverlap(i))
dis = 0.1_DP * fragments%radius(i)
fragments%rmag(i) = dis * u(i)**(THIRD)
fragments%rc(1,i) = fragments%rmag(i) * sin(theta(i)) * cos(phi(i))
fragments%rc(2,i) = fragments%rmag(i) * sin(theta(i)) * sin(phi(i))
fragments%rc(3,i) = fragments%rmag(i) * cos(theta(i))
end do
end if

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

! Check for any overlapping bodies.
Expand All @@ -397,6 +414,8 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param)
rdistance = rdistance * collider%fail_scale
end do

lfailure = any(loverlap(:))

call collision_util_shift_vector_to_origin(fragments%mass, fragments%rc)
call collider%set_coordinate_system()

Expand All @@ -421,8 +440,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param)
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
! Internals
real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir
real(DP) :: v_init, v_final, mass_init, mass_final, rotmag, dKE, KE_init, KE_final
real(DP) :: v_init, mass_init, mass_final, KE_init
real(DP), parameter :: random_scale_factor = 0.01_DP !! The relative scale factor to apply to the random component of the rotation vector
integer(I4B) :: i
logical :: lhitandrun
Expand Down Expand Up @@ -476,20 +494,20 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
logical, intent(out) :: lfailure !! Did the velocity computation fail?
! Internals
real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-2_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount)
real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-4_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount)
real(DP), parameter :: MOMENTUM_SUCCESS_METRIC = 1.0e-12_DP !! Relative angular momentum error to accept as a success (should be *much* stricter than energy)
integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps
integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best
logical :: lhitandrun, lsupercat
real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, dL, drot, rot_new
real(DP) :: vimp, vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, dE_metric, dM, mfrag, dL_metric, dL_best
real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_min, ke_avail, ke_remove, dE_best, fscale, dE_metric, dM, mfrag, dL_metric, dL_best, rn
integer(I4B), dimension(collider%fragments%nbody) :: vsign
real(DP), dimension(collider%fragments%nbody) :: vscale, volume
! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have
real(DP) :: vmin_guess = 1.01_DP
real(DP) :: vmax_guess
real(DP) :: delta_v, GC
integer(I4B), parameter :: MAXLOOP = 50
integer(I4B), parameter :: MAXTRY = 50
integer(I4B), parameter :: MAXTRY = 100
real(DP), parameter :: MAX_REDUCTION_RATIO = 0.1_DP ! Ratio of difference between first and second fragment mass to remove from the largest fragment in case of a failure
real(DP), parameter :: ROT_MAX_FRAC = 0.001_DP !! Fraction of difference between current rotation and maximum to add when angular momentum budget gets too high
class(collision_fraggle), allocatable :: collider_local
Expand Down Expand Up @@ -538,12 +556,13 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
dE_metric = huge(1.0_DP)
dE_best = huge(1.0_DP)
dL_best = huge(1.0_DP)

nsteps_best = 0
nsteps = 0
outer: do try = 1, MAXTRY
! Scale the magnitude of the velocity by the distance from the impact point
! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct
do concurrent(i = 2:nfrag)
rimp(:) = fragments%rc(:,i) - impactors%rbimp(:)
rimp(:) = fragments%rc(:,i) - impactors%rcimp(:)
vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1)))
end do
vscale(1) = 0.9_DP * minval(vscale(2:nfrag))
Expand All @@ -565,7 +584,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
end if
else
vmag = vscale(i)
rimp(:) = fragments%rc(:,i) - impactors%rbimp(:)
rimp(:) = fragments%rc(:,i) - impactors%rcimp(:)
vimp_unit(:) = .unit. (rimp(:) + vsign(i) * impactors%bounce_unit(:))
fragments%vc(:,i) = vmag * vimp_unit(:) + vrot(:)
end if
Expand All @@ -577,8 +596,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
call fragments%set_coordinate_system()
ke_min = 0.5_DP * fragments%mtot * vesc**2

E_residual = huge(1.0_DP)
inner: do loop = 1, MAXLOOP
nsteps = loop * try
nsteps = nsteps + 1

! Try to put residual angular momentum into the spin, but if this would go past the spin barrier, then put it into velocity shear instead
angmtm: do j = 1, MAXTRY
Expand All @@ -598,7 +618,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
else ! We would break the spin barrier here. Put less into spin and more into velocity shear.
if (i >= istart) then
call random_number(drot)
drot(:) = (0.5_DP * collider_local%max_rot - fragments%rotmag(i)) * 2 * (drot(:) - 0.5_DP)
call random_number(rn)
drot(:) = (rn * collider_local%max_rot - fragments%rotmag(i)) * 2 * (drot(:) - 0.5_DP)
fragments%rot(:,i) = fragments%rot(:,i) + drot(:)
fragments%rotmag(i) = .mag.fragments%rot(:,i)
if (fragments%rotmag(i) > collider%max_rot) then
Expand Down Expand Up @@ -627,6 +648,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
end do

dE = collider_local%te(2) - collider_local%te(1)
E_residual_last = E_residual
E_residual = dE + impactors%Qloss

L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1))
Expand All @@ -638,6 +660,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
E_residual_best = E_residual
dE_best = dE
dL_best = dL_metric
nsteps_best = nsteps

do concurrent(i = 1:fragments%nbody)
fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:)
Expand All @@ -646,6 +669,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
if (allocated(collider%fragments)) deallocate(collider%fragments)
allocate(collider%fragments, source=fragments)
dE_metric = abs(E_residual) / impactors%Qloss
else if (abs(E_residual) >= abs(E_residual_last)) then
exit inner ! We are no longer converging on a solution. At this point, it's best to give up
end if
if ((dE_best < 0.0_DP) .and. (dE_metric <= ENERGY_SUCCESS_METRIC * try)) exit outer ! As the tries increase, we relax the success metric. What was once a failure might become a success
end if
Expand Down Expand Up @@ -683,11 +708,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
fragments%radius(:) = (3._DP * volume(:) / (4._DP * PI))**(THIRD)

call fragments%reset()
call fraggle_generate_pos_vec(collider_local, nbody_system, param)
call fraggle_generate_pos_vec(collider_local, nbody_system, param, lfailure)
if (lfailure) exit
call fraggle_generate_rot_vec(collider_local, nbody_system, param)

! Increase the spatial size factor to get a less dense cloud
collider_local%fail_scale = collider_local%fail_scale*1.001_DP
collider_local%fail_scale = collider_local%fail_scale * 1.001_DP

! Bring the minimum and maximum velocities closer together
delta_v = 0.125_DP * (vmax_guess - vmin_guess)
Expand All @@ -707,6 +733,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
else
call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Fraggle velocity calculation converged after " // trim(adjustl(message)) // " steps.")
end if
write(message,*) nsteps_best
call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Best solution came after " // trim(adjustl(message)) // " steps.")

end associate
end associate
Expand Down
Loading

0 comments on commit 512b083

Please sign in to comment.