Merge branch 'debug'

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",

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

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

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
 

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
@@ -95,7 +95,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
       integer(I4B), parameter :: MASS_NOISE_FACTOR = 4  !! 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 = 100  !! Maximum number of fragments that can be generated
       integer(I4B), parameter :: NFRAGMIN = 1 !! Minimum number of fragments that can be generated (set by the fraggle_generate algorithm for constraining momentum and energy)
-      integer(I4B), parameter :: NFRAG_SIZE_MULTIPLIER = 3 !! Log-space scale factor that scales the number of fragments by the collisional system mass
+      integer(I4B), parameter :: NFRAG_SIZE_MULTIPLIER = 10 !! Log-space scale factor that scales the number of fragments by the collisional system mass
       integer(I4B), parameter :: iMlr = 1
       integer(I4B), parameter :: iMslr = 2
       integer(I4B), parameter :: iMrem = 3
@@ -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 (7) and 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((1 + i - iMslr)**(-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

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)