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

Commit

Permalink
Finally worked out a scheme that might yield good constraints without…
Browse files Browse the repository at this point in the history 
… needing the BFGS minimizer
  • Loading branch information
@daminton
daminton committed Dec 29, 2022
1 parent aeef4fe commit 963bd11
Show file tree
Hide file tree
Showing 6 changed files with 95 additions and 107 deletions.
90 changes: 36 additions & 54 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -535,11 +535,10 @@ module subroutine collision_generate_disruption_vel_vec(collider)
integer(I4B) :: i,j, loop
logical :: lhitandrun, lnoncat
real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual
real(DP), dimension(NDIM,2) :: vbounce, vcloud
real(DP), dimension(2) :: vimp, mcloud
real(DP) :: vmag, vdisp, Lmag, Lscale, rotmag
real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof
integer(I4B), dimension(collider%fragments%nbody) :: vsign
real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale
real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail
integer(I4B), parameter :: MAXLOOP = 100

associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN)
Expand All @@ -553,35 +552,11 @@ module subroutine collision_generate_disruption_vel_vec(collider)
vsign(:) = 1
end where

! The dominant velocities of the two clouds will be scaled by the CoM velocities of the two bodies
vimp(1:2) = .mag.impactors%vc(:,1:2)

! Now shift the CoM of each fragment cloud to what the origin body would have been in a bounce
vbounce(:,1) = -.mag.impactors%vc(:,1) * impactors%bounce_unit(:)
vbounce(:,2) = .mag.impactors%vc(:,2) * impactors%bounce_unit(:)

! Compute the current CoM of the fragment clouds
vcloud(:,:) = 0.0_DP
do concurrent(j = 1:2)
mcloud(j) = sum(fragments%mass(:), fragments%origin_body(:) == j)
do concurrent(i = 1:nfrag, fragments%origin_body(i) == j)
vcloud(:,j) = vcloud(:,j) + fragments%mass(i) * fragments%vc(:,i)
end do
vcloud(:,j) = vcloud(:,j) / mcloud(j)
end do

! Subtract off the difference between the cloud CoM velocity and the expected CoM velocity from bouncing
vcloud(:,:) = vcloud(:,:) - vbounce(:,:)
do concurrent(i = 1:nfrag)
j = fragments%origin_body(i)
fragments%vc(:,i) = fragments%vc(:,i) - vcloud(:,j)
end do

! Compute the velocity dispersion based on the escape velocity
! The minimum fragment velocity will be set by the escape velocity
if (lhitandrun) then
vdisp = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2))
vesc = sqrt(2 * impactors%Gmass(2) / impactors%radius(2))
else
vdisp = 2 * sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:)))
vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:)))
end if

! Scale the magnitude of the velocity by the distance from the impact point
Expand All @@ -590,44 +565,54 @@ module subroutine collision_generate_disruption_vel_vec(collider)
rimp(:) = fragments%rc(:,i) - impactors%rbimp(:)
vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1)))
end do
vscale(:) = vscale(:)/maxval(vscale(:))
vscale(:) = vscale(:)/minval(vscale(:))

! Give the fragment velocities a random value that is scaled with fragment mass
call random_number(mass_vscale)
mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2
mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The 1/8 power is arbitrary. It just gives the velocity a small mass dependence
mass_vscale(:) = 2*mass_vscale(:) / maxval(mass_vscale(:))
mass_vscale(:) = mass_vscale(:) / minval(mass_vscale(:))

! Set the velocities of all fragments using all of the scale factors determined above
fragments%vc(:,1) = .mag.impactors%vc(:,1) * impactors%bounce_unit(:)
do concurrent(i = 2:nfrag)
do concurrent(i = 1:nfrag)
j = fragments%origin_body(i)
vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j))
vmag = .mag.fragments%vc(:,i) * vscale(i) * mass_vscale(i)
vmag = vesc * vscale(i) * mass_vscale(i)
if (lhitandrun) then
fragments%vc(:,i) = vmag * 0.5_DP * impactors%bounce_unit(:) * vsign(i) + vrot(:)
fragments%vc(:,i) = vmag * impactors%bounce_unit(:) * vsign(i) + vrot(:)
else
! Add more velocity dispersion to disruptions vs hit and runs.
rimp(:) = fragments%rc(:,i) - impactors%rbimp(:)
vimp_unit(:) = .unit. rimp(:)
fragments%vc(:,i) = vmag * 0.5_DP * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:)
fragments%vc(:,i) = vmag * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:)
end if
end do

! Check for any residual angular momentum, and if there is any, put it into spin of the biggest body
call collider%set_coordinate_system()
call fragments%get_angular_momentum()
do loop = 1, MAXLOOP
call collider%set_coordinate_system()
call fragments%get_kinetic_energy()
ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)
! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape
ke_avail(:) = fragments%ke_orbit_frag(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%vmag(:)
ke_per_dof = -ke_residual/nfrag
do concurrent(i = 1:nfrag, ke_avail(i) > ke_per_dof)
fragments%vmag(i) = sqrt(2 * (fragments%ke_orbit_frag(i) - ke_per_dof)/fragments%mass(i))
fragments%vc(:,i) = fragments%vmag(i) * fragments%v_unit(:,i)
end do
call fragments%get_kinetic_energy()
ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)

! Check for any residual angular momentum, and if there is any, put it into spin of the biggest body
call collider%set_coordinate_system()
call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:))
rotmag = .mag. fragments%rot(:,1)
fragments%rot(:,1) = fragments%rot(:,1) + Lresidual(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1))
rotmag = .mag. fragments%rot(:,1)

Lscale = fragments%mtot * sum(fragments%radius(:)) * sum(vimp(:))
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:))
Lmag = .mag.Lresidual(:)/Lscale
rotmag = .mag. fragments%rot(:,1)
fragments%rot(:,1) = fragments%rot(:,1) + Lresidual(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1))
rotmag = .mag. fragments%rot(:,1)
if (ke_residual >= 0.0_DP) exit

call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:))
Lmag = .mag.Lresidual(:)/Lscale
end do

do concurrent(i = 1:nfrag)
fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:)
Expand All @@ -639,12 +624,9 @@ module subroutine collision_generate_disruption_vel_vec(collider)
end do
impactors%vbcom(:) = impactors%vbcom(:) / fragments%mtot

! Distribute any remaining angular momentum into fragments pin
call fragments%set_spins()

end associate
return
end subroutine collision_generate_disruption_vel_vec


end submodule s_collision_generate
end submodule s_collision_generate
21 changes: 12 additions & 9 deletions src/collision/collision_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,16 +106,19 @@ module collision
real(DP), dimension(nbody) :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame
real(DP), dimension(nbody) :: vmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame
real(DP), dimension(nbody) :: rotmag !! Array of rotation magnitudes of individual fragments
real(DP), dimension(NDIM,nbody) :: v_r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame
real(DP), dimension(NDIM,nbody) :: v_t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame
real(DP), dimension(NDIM,nbody) :: v_n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame
real(DP), dimension(NDIM,nbody) :: r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame
real(DP), dimension(NDIM,nbody) :: v_unit !! Array of velocity direction unit vectors of individual fragments in the collisional coordinate frame
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) :: Lorbit !! Orbital angular momentum vector of all fragments
real(DP), dimension(NDIM) :: Lspin !! Spin angular momentum vector of all fragments
real(DP) :: ke_orbit !! Orbital kinetic energy of all fragments
real(DP) :: ke_spin !! Spin kinetic energy of all fragments
real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories
real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories
real(DP), dimension(NDIM) :: Lorbit !! Orbital angular momentum vector of all fragments
real(DP), dimension(NDIM) :: Lspin !! Spin angular momentum vector of all fragments
real(DP) :: ke_orbit !! Orbital kinetic energy of all fragments
real(DP) :: ke_spin !! Spin kinetic energy of all fragments
real(DP) :: ke_budget !! Kinetic energy budget for computing fragment trajectories
real(DP), dimension(NDIM) :: L_budget !! Angular momentum budget for computing fragment trajectories
real(DP), dimension(nbody) :: ke_orbit_frag !! Orbital kinetic energy of each individual fragment
real(DP), dimension(nbody) :: ke_spin_frag !! 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
Expand Down
37 changes: 20 additions & 17 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,21 +205,23 @@ module subroutine collision_util_get_kinetic_energy(self)
!! Calculates the current kinetic energy of the fragments
implicit none
! Argument
class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object
class(collision_fragments(*)), intent(inout) :: self !! Fragment system object
! Internals
integer(I4B) :: i

associate(fragments => self, nfrag => self%nbody)
fragments%ke_orbit = 0.0_DP
fragments%ke_spin = 0.0_DP
fragments%ke_orbit_frag(:) = 0.0_DP
fragments%ke_spin_frag(:) = 0.0_DP

do i = 1, nfrag
fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i))
fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) )
do concurrent(i = 1:nfrag)
fragments%ke_orbit_frag(i) = fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i))
fragments%ke_spin_frag(i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) )
end do

fragments%ke_orbit = fragments%ke_orbit / 2
fragments%ke_spin = fragments%ke_spin / 2
fragments%ke_orbit_frag(:) = fragments%ke_orbit_frag(:) / 2
fragments%ke_spin_frag(:) = fragments%ke_spin_frag(:) / 2
fragments%ke_orbit = sum(fragments%ke_orbit_frag(:))
fragments%ke_spin = sum(fragments%ke_spin_frag(:))

end associate

Expand Down Expand Up @@ -374,9 +376,9 @@ module subroutine collision_util_reset_fragments(self)
self%density(:) = 0.0_DP
self%rc(:,:) = 0.0_DP
self%vc(:,:) = 0.0_DP
self%v_r_unit(:,:) = 0.0_DP
self%v_t_unit(:,:) = 0.0_DP
self%v_n_unit(:,:) = 0.0_DP
self%r_unit(:,:) = 0.0_DP
self%t_unit(:,:) = 0.0_DP
self%n_unit(:,:) = 0.0_DP

return
end subroutine collision_util_reset_fragments
Expand Down Expand Up @@ -453,9 +455,10 @@ module subroutine collision_util_set_coordinate_collider(self)
fragments%vmag(:) = .mag. fragments%vc(:,:)

! Define the radial, normal, and tangential unit vectors for each individual fragment
fragments%v_r_unit(:,:) = .unit. fragments%rc(:,:)
fragments%v_t_unit(:,:) = .unit. fragments%vc(:,:)
fragments%v_n_unit(:,:) = .unit. (fragments%v_r_unit(:,:) .cross. fragments%v_t_unit(:,:))
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%r_unit(:,:))

end associate

Expand Down Expand Up @@ -770,9 +773,9 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag)
self%fragments%vc(:,:) = 0.0_DP
self%fragments%rot(:,:) = 0.0_DP
self%fragments%Ip(:,:) = 0.0_DP
self%fragments%v_r_unit(:,:) = 0.0_DP
self%fragments%v_t_unit(:,:) = 0.0_DP
self%fragments%v_n_unit(:,:) = 0.0_DP
self%fragments%r_unit(:,:) = 0.0_DP
self%fragments%t_unit(:,:) = 0.0_DP
self%fragments%n_unit(:,:) = 0.0_DP
self%fragments%mass(:) = 0.0_DP
self%fragments%radius(:) = 0.0_DP
self%fragments%density(:) = 0.0_DP
Expand Down
Loading

0 comments on commit 963bd11

Please sign in to comment.