Merge remote-tracking branch 'origin/master'

MintonGroup · Feb 13, 2023 · 2f14424 · 2f14424
2 parents 85bcfd7 + 7b220e3
commit 2f14424
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 5 deletions.
diff --git 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
@@ -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)