diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 5afb9484d..3a8f32b50 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -114,6 +114,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() @@ -338,7 +348,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/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)