Fixed issues converging on position and fragment number

MintonGroup · Feb 21, 2023 · d3b9a93 · d3b9a93
1 parent 110d334
commit d3b9a93
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 10 deletions.
diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90
@@ -309,10 +309,8 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
       integer(I4B), parameter :: MAXLOOP = 10000
       real(DP), parameter :: cloud_size_scale_factor = 3.0_DP ! Scale factor to apply to the size of the cloud relative to the distance from the impact point. 
                                                               ! A larger value puts more space between fragments initially
-      real(DP) :: rbuffer  ! Body radii are inflated by this scale factor to prevent secondary collisions 
-      real(DP), parameter :: rbuffer_max = 1.02_DP
+      real(DP), parameter :: rbuffer = 1.01_DP ! Body radii are inflated by this scale factor to prevent secondary collisions 
       real(DP), parameter :: pack_density = 0.5236_DP ! packing density of loose spheres
-      rbuffer = 1.01_DP 
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody, &
          pl => nbody_system%pl, tp => nbody_system%tp, npl => nbody_system%pl%nbody, ntp => nbody_system%tp%nbody)
@@ -415,24 +413,24 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
             fragments%rmag(:) = .mag. fragments%rc(:,:)
 
             ! Check for any overlapping bodies.
-            !loverlap(:) = .false.
-            do j = 1, nfrag
+            do j = nfrag, 1, -1
                if (.not.loverlap(j)) cycle
                loverlap(j) = .false.
                ! Check for overlaps between fragments
-               do concurrent(i = 1:nfrag, i/=j)
+               do i = 1,nfrag
+                  if (i == j) cycle
                   dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i))
-                  loverlap(j) = loverlap(j) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) 
+                  loverlap(j) = (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j)))
+                  if (loverlap(j)) exit
                end do
-               ! Check for overlaps with existing bodies that are not involved in the collision 
+               if (loverlap(j)) cycle
                do i = 1, npl
                   if (any(impactors%id(:) == i)) cycle
                   dis = .mag. (fragments%rc(:,j) - (pl%rb(:,i) / collider%dscale - impactors%rbcom(:)))
                   loverlap(j) = loverlap(j) .or. (dis <= rbuffer * (pl%radius(i) / collider%dscale + fragments%radius(j))) 
                end do
             end do
             rdistance = rdistance * collider%fail_scale
-            rbuffer = min(rbuffer * collider%fail_scale, rbuffer_max)
          end do
 
          lfailure = any(loverlap(:))

diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90
@@ -177,7 +177,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
             ! Use Newton's method solver to get the logspace slope of the mass function
             Mrat = (mremaining + Mslr) / Mslr
             nrem = nfrag - 2
-            x0 = 1.0_DP - 100*epsilon(1.0_DP)
+            x0 = 1.0_DP - 1.0_DP/Mrat 
             x1 = Mrat**(1.0/nrem)
             do j = 1, MAXLOOP 
                y0 = Mrat - (1.0_DP - x0**nrem)/(1.0_DP - x0)