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

Commit

Permalink
More tweaks to find the best fragment distribution for the tricky cases
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 31, 2022
1 parent cf7f9f6 commit 28e5cde
Showing 1 changed file with 64 additions and 57 deletions.
121 changes: 64 additions & 57 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -370,13 +370,14 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object
logical, intent(out) :: lfailure !! Did the velocity computation fail?
! Internals
integer(I4B) :: i, j, loop, istart, n, ndof
integer(I4B) :: i, j, loop, try, istart, n, ndof
logical :: lhitandrun, lsupercat
real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear, vunit
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 = 1000
integer(I4B), parameter :: MAXLOOP = 100
integer(I4B), parameter :: MAXTRY = 100
real(DP), parameter :: TOL = 1e-6
class(collision_fragments(:)), allocatable :: fragments

Expand Down Expand Up @@ -442,68 +443,74 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
end if
call fragments%set_coordinate_system()
ke_residual_min = -huge(1.0_DP)
do loop = 1, MAXLOOP
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)) .or. ((ke_residual >= 0.0_DP) .and. (ke_residual_min < 0.0_DP))) 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
if ((ke_residual > 0.0_DP) .and. (ke_residual < TOL * fragments%ke_budget)) exit
end if
! 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(:)
ke_tot = 0.0_DP
ke_per_dof = -ke_residual
do i = 1, 2*(nfrag - istart + 1)
n = count(ke_avail(istart:nfrag) > -ke_residual/i)
if (ke_residual < 0.0_DP) n = n + count(fragments%ke_spin(istart:nfrag) > -ke_residual/i)
if (abs(n * ke_per_dof) > ke_tot) then
ke_per_dof = -ke_residual/i
ke_tot = n * ke_per_dof
ndof = i
if (abs(ke_tot) > abs(ke_residual)) then
ke_tot = -ke_residual
ke_per_dof = ke_tot/n
exit
end if
outer: do try = 1, MAXTRY
do loop = 1, MAXLOOP
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)) .or. ((ke_residual >= 0.0_DP) .and. (ke_residual_min < 0.0_DP))) 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
if ((ke_residual > 0.0_DP) .and. (ke_residual < TOL * fragments%ke_budget)) exit outer
end if
end do
do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof)
vmag = max(fragments%vmag(i)**2 - 2*ke_per_dof/fragments%mass(i),vesc**2)
fragments%vmag(i) = sqrt(vmag)
fragments%vc(:,i) = fragments%vmag(i) * .unit.fragments%vc(:,i)
end do
do concurrent(i = istart:nfrag, fragments%ke_spin(i) > ke_per_dof)
rotmag = fragments%rotmag(i)**2 - 2*ke_per_dof/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i))
rotmag = max(rotmag, 0.0_DP)
fragments%rotmag(i) = sqrt(rotmag)
fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i)
end do
! 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(:)
ke_tot = 0.0_DP
ke_per_dof = -ke_residual
do i = 1, 2*(nfrag - istart + 1)
n = count(ke_avail(istart:nfrag) > -ke_residual/i)
if (ke_residual < 0.0_DP) n = n + count(fragments%ke_spin(istart:nfrag) > -ke_residual/i)
if (abs(n * ke_per_dof) > ke_tot) then
ke_per_dof = -ke_residual/i
ke_tot = n * ke_per_dof
ndof = i
if (abs(ke_tot) > abs(ke_residual)) then
ke_tot = -ke_residual
ke_per_dof = ke_tot/n
exit
end if
end if
end do
do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof)
vmag = max(fragments%vmag(i)**2 - 2*ke_per_dof/fragments%mass(i),vesc**2)
fragments%vmag(i) = sqrt(vmag)
fragments%vc(:,i) = fragments%vmag(i) * .unit.fragments%vc(:,i)
end do
do concurrent(i = istart:nfrag, fragments%ke_spin(i) > ke_per_dof)
rotmag = fragments%rotmag(i)**2 - 2*ke_per_dof/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i))
rotmag = max(rotmag, 0.0_DP)
fragments%rotmag(i) = sqrt(rotmag)
fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i)
end do

call fragments%set_coordinate_system()
if (lsupercat) then
! Put some of the residual angular momentum into velocity shear. Not too much, or we get some weird trajectories
call fragments%set_coordinate_system()
if (lsupercat) then
! Put some of the residual angular momentum into velocity shear. Not too much, or we get some weird trajectories
call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:))
do concurrent(i = istart:nfrag)
vunit(:) = .unit. (Lresidual(:) .cross. fragments%r_unit(:,i))
vshear(:) = vunit(:) * (.mag.Lresidual(:) / ((nfrag-istart+1)*fragments%mass(i) * fragments%rmag(i)))
fragments%vc(:,i) = fragments%vc(:,i) + vshear(:)
end do
end if
! Check for any residual angular momentum, and if there is any, put it into spin
call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:))
do concurrent(i = istart:nfrag)
vunit(:) = .unit. (Lresidual(:) .cross. fragments%r_unit(:,i))
vshear(:) = vunit(:) * (.mag.Lresidual(:) / ((nfrag-istart+1)*fragments%mass(i) * fragments%rmag(i)))
fragments%vc(:,i) = fragments%vc(:,i) + vshear(:)
do concurrent(i = 1:nfrag)
fragments%Lspin(:,i) = fragments%Lspin(:,i) + Lresidual(:) / nfrag
fragments%rot(:,i) = fragments%Lspin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i))
end do
end if
! Check for any residual angular momentum, and if there is any, put it into spin
call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:))
do concurrent(i = 1:nfrag)
fragments%Lspin(:,i) = fragments%Lspin(:,i) + Lresidual(:) / nfrag
fragments%rot(:,i) = fragments%Lspin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i))
end do

call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:))
call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:))

end do
end do
! We didn't converge. Try another configuration and see if we get a better result
call fraggle_generate_pos_vec(collider)
call fraggle_generate_rot_vec(collider)
collider%fail_scale = collider%fail_scale*1.01_DP
end do outer
lfailure = ke_residual < 0.0_DP

do concurrent(i = 1:nfrag)
Expand Down

0 comments on commit 28e5cde

Please sign in to comment.