From b0d38888190a056d25f02f8fa2c7d741c7282de5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 29 Dec 2022 11:11:34 -0500 Subject: [PATCH] Better angular momentum control on the basic disruption model --- examples/Fragmentation/Fragmentation_Movie.py | 2 +- src/collision/collision_generate.f90 | 66 +++++++++---------- src/fraggle/fraggle_generate.f90 | 2 +- 3 files changed, 34 insertions(+), 36 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index ebad58734..f216fc242 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -232,7 +232,7 @@ def data_stream(self, frame=0): # Set fragmentation parameters minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades - sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(collision_model="disruption", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) sim.run(dt=1e-3, tstop=1.0e-3, istep_out=1, dump_cadence=0) print("Generating animation") diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 1f588b4ad..c3348841d 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -534,10 +534,10 @@ module subroutine collision_generate_simple_vel_vec(collider) ! Internals integer(I4B) :: i,j, loop logical :: lhitandrun, lnoncat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear + 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 + real(DP) :: vmag, vdisp, Lmag, Lscale, rotmag integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale @@ -553,6 +553,30 @@ module subroutine collision_generate_simple_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 if (lhitandrun) then vdisp = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) @@ -560,9 +584,6 @@ module subroutine collision_generate_simple_vel_vec(collider) vdisp = 2 * sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) end if - ! 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) - ! 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 = 1:nfrag) @@ -582,7 +603,7 @@ module subroutine collision_generate_simple_vel_vec(collider) do concurrent(i = 2:nfrag) j = fragments%origin_body(i) vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) - vmag = .mag.impactors%vc(:,j) * vscale(i) * mass_vscale(i) + vmag = .mag.fragments%vc(:,i) * vscale(i) * mass_vscale(i) if (lhitandrun) then fragments%vc(:,i) = vmag * 0.5_DP * impactors%bounce_unit(:) * vsign(i) + vrot(:) else @@ -593,40 +614,17 @@ module subroutine collision_generate_simple_vel_vec(collider) end if end do - ! ! 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 - - ! Check for any residual angular momentum, and if there is any, put it into velocity shear of the fragments + ! 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() Lscale = fragments%mtot * sum(fragments%radius(:)) * sum(vimp(:)) Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) Lmag = .mag.Lresidual(:)/Lscale - do concurrent(i = 3:nfrag) - vshear(:) = (Lresidual(:) / (nfrag-2)) / (fragments%mass(i) * fragments%rmag(i)) - fragments%vc(:,i) = fragments%vc(:,i) + vshear(:) - end do - ! Recompute the collision system coordinates, which will also compute the unit basis vectors for each fragment - call collider%set_coordinate_system() + 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) + call fragments%get_angular_momentum() Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) Lmag = .mag.Lresidual(:)/Lscale diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index f8718e25c..6fe15633b 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -94,7 +94,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur ! Minimize difference between energy/momentum and budgets call fraggle_generate_minimize(self, lfailure_local) - ! call fraggle_generate_tan_vel(self, lfailure_local) + call fraggle_generate_tan_vel(self, lfailure_local) ! call fraggle_generate_rad_vel(self, lfailure_local) call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)