diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 3eaa65ee4..9fdb9e8ce 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -239,6 +239,11 @@ def __init__(self,read_param: bool = False, read_data: bool = False, simdir: os. fragmentation model is enabled. Ignored otherwise *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass Parameter input file equivalent: `MIN_GMFRAG` + nfrag_reduction : float, optional + If fragmentation is turne don, this is a reduction factor used to limit the number of fragments generated in a collision. + For instance, if the SFD of the collision would generated 300 fragments above the `minimum_fragment_mass`, then a value + of `nfrag_reduction = 30.0` would reduce it to 10. + *Note.* Currently only used by the Fraggle collision model. rotation : bool, default False If set to True, this turns on rotation tracking and radius, rotation vector, and moments of inertia values must be included in the initial conditions. @@ -764,6 +769,7 @@ def set_parameter(self, verbose: bool = True, **kwargs): "qmin_coord": "HELIO", "gmtiny": 0.0, "mtiny": None, + "nfrag_reduction": 30.0, "close_encounter_check": True, "general_relativity": True, "collision_model": "FRAGGLE", @@ -1009,6 +1015,7 @@ def set_feature(self, collision_model: Literal["MERGE","BOUNCE","FRAGGLE"] | None = None, minimum_fragment_gmass: float | None = None, minimum_fragment_mass: float | None = None, + nfrag_reduction: float | None = None, rotation: bool | None = None, compute_conservation_values: bool | None = None, extra_force: bool | None = None, @@ -1054,6 +1061,11 @@ def set_feature(self, fragmentation model is enabled. Ignored otherwise *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass Parameter input file equivalent: `MIN_GMFRAG` + nfrag_reduction : float, optional + If fragmentation is turne don, this is a reduction factor used to limit the number of fragments generated in a collision. + For instance, if the SFD of the collision would generated 300 fragments above the `minimum_fragment_mass`, then a value + of `nfrag_reduction = 30.0` would reduce it to 10. + *Note.* Currently only used by the Fraggle collision model. rotation : bool, optional If set to True, this turns on rotation tracking and radius, rotation vector, and moments of inertia values must be included in the initial conditions. @@ -1151,6 +1163,10 @@ def set_feature(self, self.param["MIN_GMFRAG"] = minimum_fragment_mass * self.GU if "minimum_fragment_gmass" not in update_list: update_list.append("minimum_fragment_gmass") + + if nfrag_reduction is not None: + self.param["NFRAG_REDUCTION"] = nfrag_reduction + update_list.append("nfrag_reduction") if rotation is not None: self.param['ROTATION'] = rotation @@ -1167,6 +1183,7 @@ def set_feature(self, if extra_force is not None: self.param["EXTRA_FORCE"] = extra_force update_list.append("extra_force") + if big_discard is not None: if self.codename != "Swifter": @@ -1264,6 +1281,7 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N "collision_model": "COLLISION_MODEL", "encounter_save": "ENCOUNTER_SAVE", "minimum_fragment_gmass": "MIN_GMFRAG", + "nfrag_reduction": "NFRAG_REDUCTION", "rotation": "ROTATION", "general_relativity": "GR", "compute_conservation_values": "ENERGY", diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index 0669f62f3..51d632625 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -54,6 +54,7 @@ module base real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units real(DP) :: GMTINY = -1.0_DP !! Smallest G*mass that is fully gravitating real(DP) :: min_GMfrag = -1.0_DP !! Smallest G*mass that can be produced in a fragmentation event + real(DP) :: nfrag_reduction = 30.0_DP !! Reduction factor for limiting the number of fragments in a collision integer(I4B), dimension(:), allocatable :: seed !! Random seeds for fragmentation modeling logical :: lmtiny_pl = .false. !! Include semi-interacting massive bodies character(STRMAX) :: collision_model = "MERGE" !! The Coll diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index fd68287c2..150eb06c9 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -34,6 +34,8 @@ module subroutine fraggle_generate(self, nbody_system, param, t) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) class is (swiftest_pl) + select type(param) + class is (swiftest_parameters) associate(impactors => self%impactors, status => self%status, maxid => nbody_system%maxid) ! Set the coordinate system of the impactors call impactors%set_coordinate_system() @@ -99,6 +101,7 @@ module subroutine fraggle_generate(self, nbody_system, param, t) end associate end select end select + end select return end subroutine fraggle_generate @@ -228,6 +231,8 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) class is (swiftest_pl) + select type(param) + class is (swiftest_parameters) associate(impactors => self%impactors, maxid => nbody_system%maxid) call collision_io_collider_message(nbody_system%pl, impactors%id, message) if (impactors%mass(1) > impactors%mass(2)) then @@ -273,6 +278,7 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) end associate end select end select + end select return end subroutine fraggle_generate_hitandrun diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 37905a188..7cdfbd925 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -85,7 +85,7 @@ end subroutine fraggle_util_restructure module subroutine fraggle_util_set_mass_dist(self, param) implicit none class(collision_fraggle), intent(inout) :: self !! Fraggle collision object - class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine fraggle_util_set_mass_dist end interface diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 82f026621..f5f84ed8f 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -83,8 +83,8 @@ module subroutine fraggle_util_set_mass_dist(self, param) !! implicit none ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + 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 real(DP), dimension(2) :: volume @@ -124,35 +124,19 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! the limits bracketed above and the model size distribution of fragments. ! Check to see if our size distribution would give us a smaller number of fragments than the maximum number - select type(param) - class is (swiftest_parameters) - ! 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) * mass_noise - ! The number of fragments we generate is bracked by the minimum required by fraggle_generate nd the - ! maximum set by the NFRAG_SIZE_MULTIPLIER which limits the total number of fragments to prevent the nbody - ! code from getting an overwhelmingly large number of fragments - nfrag = ceiling(NFRAG_SIZE_MULTIPLIER * log(mtot / min_mfrag)) - nfrag = max(min(nfrag, NFRAGMAX), NFRAGMIN) - class default - min_mfrag = 0.0_DP - nfrag = NFRAGMAX - end select + ! ! 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) mremaining = impactors%mass_dist(iMrem) - if (mremaining > 0.0_DP) then - i = iMrem - do while (i <= nfrag) - mfrag = max((i - 1)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) - if (mremaining - mfrag <= 0.0_DP) exit - mremaining = mremaining - mfrag - i = i + 1 - end do - else - i = iMrem - 1 - end if - nfrag = i + nfrag = iMrem - 1 + do while (mremaining > 0.0_DP) + nfrag = nfrag + 1 + mfrag = max((nfrag - 1)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) + mremaining = mremaining - mfrag + end do + nfrag = ceiling(nfrag / param%nfrag_reduction) call self%setup_fragments(nfrag) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) @@ -167,7 +151,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) end associate return case default - write(*,*) "collision_util_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime + write(*,*) "fraggle_util_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime end select associate(fragments => self%fragments) @@ -176,14 +160,14 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! Make the first two bins the same as the Mlr and Mslr values that came from collision_regime mass(1) = impactors%mass_dist(iMlr) - mass(2) = impactors%mass_dist(iMslr) + mass(2) = impactors%mass_dist(iMslr) ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA 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((1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) * mass_noise + mfrag = max((i - 1)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) * mass_noise mass(i) = mfrag mremaining = mremaining - mfrag end do @@ -192,7 +176,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) 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 remains the second largest + 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 diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 79daa0220..6f4845065 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -2054,6 +2054,8 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i read(param_value, *) param%GMTINY case ("MIN_GMFRAG") read(param_value, *) param%min_GMfrag + case ("NFRAG_REDUCTION") + read(param_value, *) param%nfrag_reduction case ("ENCOUNTER_SAVE") call swiftest_io_toupper(param_value) read(param_value, *) param%encounter_save @@ -2229,6 +2231,10 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i call random_seed(get = param%seed) end if if (param%min_GMfrag < 0.0_DP) param%min_GMfrag = param%GMTINY + if (param%nfrag_reduction < 1.0_DP) then + write(iomsg,*) "Warning: NFRAG_REDUCTION value invalid. Setting to 1.0" + param%nfrag_reduction = 1.0_DP + end if end if ! Determine if the GR flag is set correctly for this integrator @@ -2406,6 +2412,7 @@ module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, i if (param%min_GMfrag >= 0.0_DP) call io_param_writer_one("MIN_GMFRAG",param%min_GMfrag, unit) call io_param_writer_one("COLLISION_MODEL",param%collision_model, unit) if (param%collision_model == "FRAGGLE" ) then + call io_param_writer_one("NFRAG_REDUCTION",param%nfrag_reduction, unit) nseeds = size(param%seed) call random_seed(get = param%seed) call io_param_writer_one("SEED", [nseeds, param%seed(:)], unit)