Updated the fragmentation model to do a better job of spanning the av…

…ailble mass range when generating collisional fragments.

src/fraggle/fraggle_util.f90

@@ -86,19 +86,20 @@ module subroutine fraggle_util_set_mass_dist(self, param)
       class(collision_fraggle),   intent(inout) :: self  !! Fraggle collision system object
       class(swiftest_parameters), intent(in)    :: param !! Current Swiftest run configuration parameters
       ! Internals
-      integer(I4B)              :: i, j, jproj, jtarg, nfrag, istart, nfragmax
+      integer(I4B)              :: i, j, jproj, jtarg, nfrag, istart, nfragmax, nrem
       real(DP), dimension(2)    :: volume
       real(DP), dimension(NDIM) :: Ip_avg
-      real(DP) :: mfrag, mremaining, mtot, mcumul, G, mass_noise
+      real(DP) :: mfrag, mremaining, mtot, mcumul, G, mass_noise, Mslr, mscale, y, yp, Mrat, m0
       real(DP), dimension(:), allocatable :: mass
-      real(DP), parameter :: BETA = 2.85_DP
+      real(DP)  :: beta 
       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 = 3 !! Minimum number of fragments that can be generated 
       integer(I4B), parameter :: iMlr = 1
       integer(I4B), parameter :: iMslr = 2
       integer(I4B), parameter :: iMrem = 3
+      integer(I4B), parameter :: NFRAGMIN = iMrem + 2 !! Minimum number of fragments that can be generated 
       integer(I4B), dimension(:), allocatable :: ind
+      integer(I4B), parameter :: MAXLOOP = 10000
       logical :: flipper
 
       associate(impactors => self%impactors, min_mfrag => self%min_mfrag)
@@ -107,7 +108,6 @@ module subroutine fraggle_util_set_mass_dist(self, param)
          mtot = sum(impactors%mass(:))
          G = impactors%Gmass(1) / impactors%mass(1)
          Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / mtot
-         nfragmax = ceiling(NFRAGMAX_UNSCALED / param%nfrag_reduction)
 
          if (impactors%mass(1) >= impactors%mass(2)) then
             jtarg = 1
@@ -116,7 +116,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
             jtarg = 2
             jproj = 1
          end if
-  
+
          select case(impactors%regime)
          case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC, COLLRESOLVE_REGIME_HIT_AND_RUN)
             ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of collision_regime.
@@ -127,16 +127,23 @@ module subroutine fraggle_util_set_mass_dist(self, param)
             ! ! Add a small amount of noise to the last digits of the minimum mass value so that multiple fragments don't get generated with identical mass values
             call random_number(mass_noise)
             mass_noise = 1.0_DP + mass_noise * epsilon(1.0_DP) * 10**(MASS_NOISE_FACTOR)
-            min_mfrag = (param%min_GMfrag / param%GU)
-
+            min_mfrag = (param%min_GMfrag / param%GU) * mass_noise
+            
             mremaining = impactors%mass_dist(iMrem)
             nfrag = iMrem - 1 
+            beta = 2.85_DP ! From Leinhardt & Stewart (2012)
+            Mslr = impactors%mass_dist(iMslr)
+
+            mfrag = Mslr
             do while (mremaining > 0.0_DP)
-               nfrag = nfrag + 1
-               mfrag = max((nfrag - 1)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag)
+               mfrag = (nfrag)**(-3._DP / BETA) * impactors%mass_dist(iMslr)
+               mfrag = max(mfrag, min_mfrag)
                mremaining = mremaining - mfrag
+               nfrag = nfrag + 1
             end do
+            nfragmax = ceiling(NFRAGMAX_UNSCALED / param%nfrag_reduction)
             nfrag = max(min(ceiling(nfrag / param%nfrag_reduction), nfragmax), NFRAGMIN)
+
             call self%setup_fragments(nfrag)
 
          case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) 
@@ -162,28 +169,39 @@ module subroutine fraggle_util_set_mass_dist(self, param)
             mass(1) = impactors%mass_dist(iMlr) 
             mass(2) = impactors%mass_dist(iMslr)
 
-            ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA 
+            ! Recompute the slope parameter beta so that we span the complete size range
+            if (Mslr == min_mfrag) Mslr = Mslr + impactors%mass_dist(iMrem) / nfrag
             mremaining = impactors%mass_dist(iMrem)
-            do i = iMrem, nfrag
-               call random_number(mass_noise)
-               mass_noise = 1.0_DP + mass_noise * epsilon(1.0_DP) * 10**(MASS_NOISE_FACTOR)
-               mfrag = max((i - 1)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) * mass_noise
+
+            ! The mass will be distributed evenly in logspace between the second-largest remnant and the minimum mass
+            ! Use Newton's method solver to get the logspace slope of the mass function
+            Mrat = mremaining / min_mfrag
+            nrem = nfrag - 2
+            mscale = 1.0_DP - 1.0_DP/Mrat
+            do j = 1, MAXLOOP
+               y = Mrat - (1.0_DP - mscale**nrem)/(1.0_DP - mscale)
+               yp = (mscale + mscale**nrem * ((nrem - 1) * mscale - nrem)) / (mscale * (mscale - 1.0_DP)**2)
+               m0 = mscale
+               mscale = m0 + y/yp
+               if (abs((mscale - m0)/mscale) < epsilon(1.0_DP)) exit
+            end do
+            Mslr = impactors%mass_dist(iMslr)
+
+            mass(iMslr) = Mslr
+            mfrag = min_mfrag
+            do i = iMrem,nfrag
                mass(i) = mfrag
                mremaining = mremaining - mfrag
+               mfrag = max(mfrag * mscale, min_mfrag)
             end do
 
-            ! If there is any residual mass (either positive or negative) we will distribute remaining mass proportionally among the the fragments
-            if (mremaining < 0.0_DP) then ! If the remainder is negative, this means that that the number of fragments required by the SFD is smaller than our lower limit set by fraggle_generate. 
-               istart = iMrem ! We will reduce the mass of the 3:nfrag bodies to prevent the second-largest fragment from going smaller
-            else ! If the remainder is postiive, this means that the number of fragments required by the SFD is larger than our upper limit set by computational expediency. 
-               istart = iMslr ! We will increase the mass of the 2:nfrag bodies to compensate, which ensures that the second largest fragment either remains the second largest or becomes the largest
-            end if
-            mfrag = 1._DP + mremaining / sum(mass(istart:nfrag))
-            mass(istart:nfrag) = mass(istart:nfrag) * mfrag
-
             ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution.
             mremaining = fragments%mtot - sum(mass(1:nfrag))
-            mass(nfrag) = mass(nfrag) + mremaining
+            if (mremaining < 0.0_DP) then
+               mass(iMlr) = mass(iMlr) + mremaining
+            else
+               mass(iMslr) = mass(iMslr) + mremaining
+            end if
 
             ! Sort the distribution in descending order by mass so that the largest fragment is always the first
             call swiftest_util_sort(-mass, ind)

src/swiftest/swiftest_discard.f90

@@ -43,12 +43,10 @@ module subroutine swiftest_discard_system(self, param)
          if (ltp_discards.or.lpl_discards) then
             call nc%open(param)
             if (lpl_discards) then
-               call pl_discards%write_info(self%system_history%nc, param)
                if (param%lenergy) call self%conservation_report(param, lterminal=.false.)
                call pl_discards%setup(0,param) 
             end if
             if (ltp_discards) then
-               call tp_discards%write_info(self%system_history%nc, param)
                call tp_discards%setup(0,param) 
             end if
             call nc%close()

src/swiftest/swiftest_io.f90

@@ -1421,7 +1421,6 @@ module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param,
          call cb%info%set_value(status="ACTIVE")
 
          if (param%lclose) then
-
             status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid)
             if (status == NF90_NOERR) then
                call netcdf_io_check( nf90_get_var(nc%id, nc%origin_type_varid, ctemp, count=[NAMELEN, idmax]), "netcdf_io_read_particle_info_system nf90_getvar origin_type_varid"  )
@@ -1552,7 +1551,6 @@ module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param)
       real(DP), dimension(NDIM)                 :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs
       real(DP)                                  :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf
 
-
       call self%write_info(nc, param)
 
       call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "netcdf_io_write_frame_body nf90_set_fill"  )
@@ -1749,6 +1747,7 @@ module subroutine swiftest_io_netcdf_write_info_body(self, nc, param)
       integer(I4B)                              :: i, j, idslot, old_mode
       integer(I4B), dimension(:), allocatable   :: ind
       character(len=NAMELEN) :: charstring
+      character(len=NAMELEN), dimension(self%nbody) :: origin_type
 
       call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "netcdf_io_write_info_body nf90_set_fill NF90_NOFILL"  )
 
@@ -1773,6 +1772,7 @@ module subroutine swiftest_io_netcdf_write_info_body(self, nc, param)
 
                if (param%lclose) then
                   charstring = trim(adjustl(self%info(j)%origin_type))
+                  origin_type(i) = charstring
                   call netcdf_io_check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "netcdf_io_write_info_body nf90_put_var origin_type_varid"  )
                   call netcdf_io_check( nf90_put_var(nc%id, nc%origin_time_varid,  self%info(j)%origin_time,  start=[idslot]), "netcdf_io_write_info_body nf90_put_var origin_time_varid"  )
                   call netcdf_io_check( nf90_put_var(nc%id, nc%origin_rh_varid,    self%info(j)%origin_rh(:), start=[1,idslot], count=[NDIM,1]), "netcdf_io_write_info_body nf90_put_var origin_rh_varid"  )