From 110d3343282abd59911b2170ffde5b9293dc7c78 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 20:03:33 -0500 Subject: [PATCH] Updated the SFD code to pin the top of the SFD rather than the bottom. This means that bodies smaller than the minimum can be made, but the minimum still sets the total number --- examples/Fragmentation/Fragmentation_Movie.py | 4 ++-- src/fraggle/fraggle_util.f90 | 9 +++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index adb37d371..b60dbc1db 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -166,7 +166,7 @@ def encounter_combiner(sim): enc = enc.sel(time=tgood) # The following will combine the two datasets along the time dimension, sort the time dimension, and then fill in any time gaps with interpolation - ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").chunk(dict(time=-1)).interpolate_na(dim="time", method="akima") + ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time", method="akima") # Rename the merged Target body so that their data can be combined tname=[n for n in ds['name'].data if names[0] in n] @@ -359,7 +359,7 @@ def vec_props(self, c): minimum_fragment_gmass = 0.01 * body_Gmass[style][1] gmtiny = 0.10 * body_Gmass[style][1] sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style]) - sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0, dask=False) + sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0) print("Generating animation") anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index d1327d723..92bbb4aba 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -175,15 +175,16 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! 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 + Mrat = (mremaining + Mslr) / Mslr nrem = nfrag - 2 - x0 = 1.0_DP + 100*epsilon(1.0_DP) + x0 = 1.0_DP - 100*epsilon(1.0_DP) x1 = Mrat**(1.0/nrem) do j = 1, MAXLOOP y0 = Mrat - (1.0_DP - x0**nrem)/(1.0_DP - x0) y1 = Mrat - (1.0_DP - x1**nrem)/(1.0_DP - x1) if (y0*y1 < 0.0_DP) exit x1 = x1 * 1.6_DP + x0 = x0 / 1.6_DP end do ! Find the mass scaling factor with bisection @@ -214,11 +215,11 @@ module subroutine fraggle_util_set_mass_dist(self, param) Mslr = impactors%mass_dist(iMslr) mass(iMslr) = Mslr - mfrag = min_mfrag + mfrag = Mslr do i = iMrem,nfrag + mfrag = mfrag * mscale mass(i) = mfrag mremaining = mremaining - mfrag - mfrag = max(mfrag * mscale, min_mfrag) end do ! 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.