diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 2f498a0cf..1f1e51b72 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/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) diff --git a/src/swiftest/swiftest_discard.f90 b/src/swiftest/swiftest_discard.f90 index bc3d9b1f0..707ed19ca 100644 --- a/src/swiftest/swiftest_discard.f90 +++ b/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() diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 00e882ecd..744c204ed 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/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" )