diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 75ab1521b..62e182fbb 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -536,7 +536,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu 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 logical :: lhitandrun, lsupercat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new, dL_metric + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, vdisp, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new, dL_metric real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, rn, dL1_mag, dE_conv integer(I4B), dimension(:), allocatable :: vsign real(DP), dimension(:), allocatable :: vscale @@ -569,6 +569,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) E_residual_best = huge(1.0_DP) + L_residual_best(:) = 0.0_DP lfailure = .false. dE_metric = huge(1.0_DP) dE_best = huge(1.0_DP) @@ -581,7 +582,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu if (allocated(dLi_mag)) deallocate(dLi_mag); allocate(dLi_mag(fragments%nbody)) if (lhitandrun) then - vesc = sqrt(2 * sum(fragments%Gmass(istart:fragments%nbody)) / sum(fragments%radius(istart:fragments%nbody))) + vesc = sqrt(2 * sum(fragments%Gmass(istart:fragments%nbody)) / impactors%radius(2)) vmin_guess = .mag.impactors%vc(:,2) * (1.0_DP - hitandrun_vscale) vmax_guess = .mag.impactors%vc(:,2) * (1.0_DP + hitandrun_vscale) else @@ -615,7 +616,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu j = fragments%origin_body(i) vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) if (lhitandrun) then - fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) + vrot(:) + vdisp(:) = .unit.(fragments%rc(:,i) - impactors%rc(:,2)) * vesc + fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) + vrot(:) + vdisp(:) else vmag = vscale(i) rimp(:) = fragments%rc(:,i) - impactors%rcimp(:) @@ -638,6 +640,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) L_residual_unit(:) = .unit. L_residual(:) + if (nsteps == 1) L_residual_best(:) = L_residual(:) ! Use equipartition of spin kinetic energy to distribution spin angular momentum do concurrent(i = istart:fragments%nbody) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 98ebe0654..2f498a0cf 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -94,7 +94,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) real(DP), parameter :: BETA = 2.85_DP integer(I4B), parameter :: MASS_NOISE_FACTOR = 5 !! The number of digits of random noise that get added to the minimum mass value to prevent identical masses from being generated in a single run integer(I4B), parameter :: NFRAGMAX_UNSCALED = 3000 !! Maximum number of fragments that can be generated - integer(I4B), parameter :: NFRAGMIN = 1 !! Minimum number of fragments that can be generated (set by the fraggle_generate algorithm for constraining momentum and energy) + integer(I4B), parameter :: NFRAGMIN = 3 !! Minimum number of fragments that can be generated integer(I4B), parameter :: iMlr = 1 integer(I4B), parameter :: iMslr = 2 integer(I4B), parameter :: iMrem = 3 @@ -136,7 +136,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) mfrag = max((nfrag - 1)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) mremaining = mremaining - mfrag end do - nfrag = min(ceiling(nfrag / param%nfrag_reduction), nfragmax) + nfrag = max(min(ceiling(nfrag / param%nfrag_reduction), nfragmax), NFRAGMIN) call self%setup_fragments(nfrag) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE)