diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 7d96c231e..d79615591 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -124,6 +124,16 @@ "merge_spinner" : 5.0e-3, } +nfrag_reduction = {"disruption_headon" : 10.0, + "disruption_off_axis" : 10.0, + "supercatastrophic_headon" : 10.0, + "supercatastrophic_off_axis" : 10.0, + "hitandrun_disrupt" : 1.0, + "hitandrun_pure" : 1.0, + "merge" : 1.0, + "merge_spinner" : 1.0, + } + density = 3000 * swiftest.AU2M**3 / swiftest.MSun GU = swiftest.GMSun * swiftest.YR2S**2 / swiftest.AU2M**3 body_radius = body_Gmass.copy() @@ -348,7 +358,7 @@ def vec_props(self, c): # Set fragmentation parameters minimum_fragment_gmass = 0.01 * body_Gmass[style][1] gmtiny = 0.10 * body_Gmass[style][1] - sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=10.0) + sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style]) sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0) print("Generating animation") diff --git a/src/collision/collision_regime.f90 b/src/collision/collision_regime.f90 index 29f4f60c4..f41f9490a 100644 --- a/src/collision/collision_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -187,7 +187,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, real(DP) :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat real(DP) :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat real(DP) :: Mint, Mtot, Mtmp, Mbig, Msmall - real(DP) :: Rp, rhill + real(DP) :: Rc1, rhill real(DP) :: Mresidual real(DP) :: U_binding @@ -209,14 +209,14 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, alpha = 1.0_DP Mint = m2 end if - Rp = (3 * (m1 / den1 + alpha * m2 / den2) / (4 * PI))**(1.0_DP/3.0_DP) ! (Mustill et al. 2018) - c_star = calc_c_star(Rp) + Rc1 = (3 * Mtot / (4 * PI * DENSITY1))**THIRD ! Stewart and Leinhardt (2009) + c_star = calc_c_star(Rc1) !calculate Vescp - Vescp = sqrt(2 * GC * Mtot / Rp) !Mustill et al. 2018 eq 6 + Vescp = sqrt(2 * GC * Mtot / Rc1) ! Steward and Leinhardt (2012) e.g. Fig. 7 caption. !calculate rhill - rhill = a1 * (m1 / 3.0_DP / (Mcb + m1))**(1.0_DP/3.0_DP) + rhill = a1 * (m1 / (3 *(Mcb + m1)))**(1.0_DP/3.0_DP) !calculate Vhill if ((rad2 + rad1) < rhill) then @@ -245,7 +245,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3) bcrit = rad1 / (rad1 + rad2) ! Specific binding energy - U_binding = (3 * GC * Mtot) / (5 * Rp) ! LS12 eq. 27 + U_binding = (3 * GC * Mtot) / (5 * Rc1) ! LS12 eq. 27 ke = 0.5_DP * Vimp**2 pe = - GC * m1 * m2 / (Mtot * norm2(rh2 - rh1)) @@ -325,17 +325,18 @@ function calc_Qrd_pstar(Mtarg, Mp, alpha, c_star) result(Qrd_pstar) ! Result real(DP) :: Qrd_pstar ! Internals - real(DP) :: Qrd_star1, mu_alpha, mu, Qrd_star + real(DP) :: Qrd_star1, mu_alpha, mu, Qrd_star, gamma ! calc mu, mu_alpha mu = (Mtarg * Mp) / (Mtarg + Mp) ! [kg] mu_alpha = (Mtarg * alpha * Mp) / (Mtarg + alpha * Mp) ! [kg] + gamma = Mp / Mtarg ! calc Qrd_star1 - Qrd_star1 = (c_star * 4 * PI * DENSITY1 * GC * Rp**2) / 5.0_DP + Qrd_star1 = (c_star * 4 * PI * DENSITY1 * GC * Rc1**2) / 5.0_DP ! LS12 eq. 28 ! calc Qrd_star - Qrd_star = Qrd_star1 * (((Mp / Mtarg + 1.0_DP)**2) / (4 * Mp / Mtarg))**(2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 23) + Qrd_star = Qrd_star1 * (((gamma + 1.0_DP)**2) / (4 * gamma))**(2.0_DP / (3 * MU_BAR) - 1.0_DP) !(eq 23) ! calc Qrd_pstar, v_pstar - Qrd_pstar = ((mu / mu_alpha)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) * Qrd_star ! (eq 15) + Qrd_pstar = ((mu / mu_alpha)**(2.0_DP - 3 * MU_BAR / 2.0_DP)) * Qrd_star ! (eq 15) return end function calc_Qrd_pstar @@ -365,8 +366,8 @@ function calc_Qrd_rev(Mp, Mtarg, Mint, den1, den2, Vimp, c_star) result(Mslr) ! calc Qrd_star1, v_star1 Qrd_star1 = (c_star * 4 * PI * mtot_rev * GC ) / Rc1 / 5.0_DP ! calc Qrd_pstar_rev - Qrd_star = Qrd_star1 * (((gamma_rev + 1.0_DP)**2) / (4 * gamma_rev)) ** (2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 52) - Qrd_pstar = Qrd_star * ((mu_rev / mu_alpha_rev)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) + Qrd_star = Qrd_star1 * (((gamma_rev + 1.0_DP)**2) / (4 * gamma_rev)) ** (2.0_DP / (3 * MU_BAR) - 1.0_DP) !(eq 52) + Qrd_pstar = Qrd_star * ((mu_rev / mu_alpha_rev)**(2.0_DP - 3 * MU_BAR / 2.0_DP)) Qrd_pstar_rev = Qrd_pstar * (Vhill / Vescp)**CRUFU !Rufu and Aharaonson eq (3) !calc Qr_supercat_rev Qr_supercat_rev = 1.8_DP * Qrd_pstar_rev @@ -417,8 +418,8 @@ function calc_c_star(Rc1) result(c_star) !! Result real(DP) :: c_star !! Internals - real(DP), parameter :: loR = 1.0e3_DP ! Lower bound of inteRpolation size (m) - real(DP), parameter :: hiR = 1.0e7_DP ! Upper bound of inteRpolation size (m) + real(DP), parameter :: loR = 1.0e3_DP ! Lower bound of interpolation size (m) + real(DP), parameter :: hiR = 1.0e7_DP ! Upper bound of interpolation size (m) real(DP), parameter :: loval = 5.0_DP ! Value of C* at lower bound real(DP), parameter :: hival = 1.9_DP ! Value of C* at upper bound diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 824a69b16..b9b38fc4c 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -533,7 +533,7 @@ 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-3_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) :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_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, nsteps_best, posloop logical :: lhitandrun, lsupercat @@ -542,13 +542,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu integer(I4B), dimension(:), allocatable :: vsign real(DP), dimension(:), allocatable :: vscale real(DP), dimension(:), allocatable :: dLi_mag - real(DP), parameter :: L_ROT_VEL_RATIO = 0.5_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments ! 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), parameter :: hitandrun_vscale = 0.25_DP real(DP) :: vmin_guess real(DP) :: vmax_guess - integer(I4B), parameter :: MAXLOOP = 20 - integer(I4B), parameter :: MAXTRY = 10 + integer(I4B), parameter :: MAXLOOP = 100 + integer(I4B), parameter :: MAXTRY = 100 integer(I4B), parameter :: MAXANGMTM = 1000 class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -653,12 +652,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu do i = istart,fragments%nbody dL(:) = -dL1_mag * dLi_mag(i) * L_residual_unit(:) - drot(:) = L_ROT_VEL_RATIO * dL(:) / (fragments%mass(i) * fragments%Ip(3,i) * fragments%radius(i)**2) + drot(:) = dL(:) / (fragments%mass(i) * fragments%Ip(3,i) * fragments%radius(i)**2) rot_new(:) = fragments%rot(:,i) + drot(:) if (.mag.rot_new(:) < collider_local%max_rot) then fragments%rot(:,i) = rot_new(:) fragments%rotmag(i) = .mag.fragments%rot(:,i) - else ! We would break the spin barrier here. Put less into spin and more into velocity shear. + else ! We would break the spin barrier here. Add a random component of rotation that is less than what would break the limit. The rest will go in velocity shear call random_number(drot) call random_number(rn) drot(:) = (rn * collider_local%max_rot - fragments%rotmag(i)) * 2 * (drot(:) - 0.5_DP) @@ -669,7 +668,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fragments%rot(:,i) = fragments%rotmag(i) * .unit. fragments%rot(:,i) end if end if - L_residual(:) = L_residual(:) - drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 + L_residual(:) = L_residual(:) + drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 end do ! Put any remaining residual into velocity shear @@ -682,7 +681,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu if (all(dL_metric(:) <= 1.0_DP)) exit angmtm do i = istart, fragments%nbody - dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody)) + dL(:) = -L_residual(:) / (fragments%nbody - istart + 1) call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) fragments%vmag(i) = .mag.fragments%vc(:,i)