From d3b9a93d094e35660d1ab70492c6adc7d854161a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 20:44:10 -0500 Subject: [PATCH] Fixed issues converging on position and fragment number --- src/fraggle/fraggle_generate.f90 | 16 +++++++--------- src/fraggle/fraggle_util.f90 | 2 +- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 5ef782c0a..554da1d31 100644 --- 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,16 +413,17 @@ 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(:))) @@ -432,7 +431,6 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu 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 index 92bbb4aba..5863a3073 100644 --- 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)