Improved the SFD of fragments. Fragment number is now variable, and w…

…eakly depends on the collisional system mass.

examples/symba_energy_momentum/param.disruption_headon.in

@@ -21,7 +21,8 @@ ENC_OUT        enc.disruption_headon.dat
 EXTRA_FORCE    no                ! no extra user-defined forces
 BIG_DISCARD    no                ! output all planets if anything discarded
 RHILL_PRESENT  yes               ! Hill's sphere radii in input file
-GMTINY          1.0e-16
+GMTINY         1.0e-16
+MIN_GMFRAG     1.0e-10
 FRAGMENTATION  yes
 MU2KG          1.98908e30
 TU2S           3.1556925e7

examples/symba_energy_momentum/param.disruption_off_axis.in

@@ -22,7 +22,8 @@ DISCARD_OUT    discard.disruption_off_axis.out
 EXTRA_FORCE    no                ! no extra user-defined forces
 BIG_DISCARD    no                ! output all planets if anything discarded
 RHILL_PRESENT  yes               ! Hill's sphere radii in input file
-GMTINY          1.0e-16
+GMTINY         1.0e-16
+MIN_GMFRAG     1.0e-12
 FRAGMENTATION  yes
 MU2KG          1.98908e30
 TU2S           3.1556925e7

examples/symba_energy_momentum/param.supercatastrophic_headon.in

@@ -17,12 +17,11 @@ CHK_RMIN       0.005
 CHK_RMAX       1e6
 CHK_EJECT      -1.0              ! ignore this check
 CHK_QMIN       -1.0              ! ignore this check
-DISCARD_OUT    discard.supercatastrophic_headon.out
-ENC_OUT        enc.supercatastrophic_headon.dat
 EXTRA_FORCE    no                ! no extra user-defined forces
 BIG_DISCARD    no                ! output all planets if anything discarded
 RHILL_PRESENT  yes               ! Hill's sphere radii in input file
-GMTINY          1.0e-16
+GMTINY         1.0e-16
+MIN_GMFRAG     1.0e-11
 FRAGMENTATION  yes
 MU2KG          1.98908e30
 TU2S           3.1556925e7

examples/symba_energy_momentum/param.supercatastrophic_off_axis.in

@@ -21,7 +21,8 @@ ENC_OUT        enc.supercatastrophic_off_axis.dat
 EXTRA_FORCE    no                ! no extra user-defined forces
 BIG_DISCARD    no                ! output all planets if anything discarded
 RHILL_PRESENT  yes               ! Hill's sphere radii in input file
-GMTINY          1.0e-16
+GMTINY         1.0e-16
+MIN_GMFRAG     1.0e-11
 FRAGMENTATION  yes
 MU2KG          1.98908e30
 TU2S           3.1556925e7

src/fraggle/fraggle_regime.f90

@@ -56,6 +56,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param)
          frag%mass_dist(1) = min(max(mlr, 0.0_DP), mtot)
          frag%mass_dist(2) = min(max(mslr, 0.0_DP), mtot)
          frag%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot)
+         frag%mtot = sum(frag%mass_dist(1:3))
 
          ! Convert quantities back to the system units and save them into the fragment system
          frag%mass_dist(:) = (frag%mass_dist(:) / param%MU2KG) 
@@ -98,6 +99,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb
       real(DP), parameter   :: DENSITY1 = 1000.0_DP !standard density parameter from LS12 [kg/m3]
       real(DP), parameter   :: MU_BAR = 0.37_DP !0.385#0.37#0.3333# 3.978 # 1/3 material parameter for hydrodynamic planet-size bodies (LS12)
       real(DP), parameter   :: BETA = 2.85_DP !slope of sfd for remnants from LS12 2.85
+      real(DP), parameter   :: ETA = -1.50_DP !! LS12 eq. (44)
       real(DP), parameter   :: C1 = 2.43_DP  !! Kokubo & Genda (2010) eq. (3)
       real(DP), parameter   :: C2 = -0.0408_DP !! Kokubo & Genda (2010) eq. (3)
       real(DP), parameter   :: C3 = 1.86_DP !! Kokubo & Genda (2010) eq. (3)
@@ -108,7 +110,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb
       real(DP)           :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta
       real(DP)           :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat
       real(DP)           :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat
-      real(DP)           :: Mint, Mtot
+      real(DP)           :: Mint, Mtot, Mtmp
       real(DP)           :: Rp, rhill 
       real(DP)           :: Mresidual
       real(DP)           :: U_binding
@@ -157,7 +159,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb
       Qr = mu*(Vimp**2) / Mtot / 2.0_DP
 
       !calculate mass largest remnant Mlr 
-      Mlr = (1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot  ! [kg] # LS12 eq (5)
+      Mlr = max((1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot, min_mfrag) ! [kg] # LS12 eq (5)
 
       !calculate Vsupercat
       Qr_supercat = SUPERCAT_QRATIO * Qrd_pstar ! See LS12 Section 4.1 
@@ -196,7 +198,7 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb
                Qloss = 0.0_DP
             else
                Mlr = m1
-               Mslr = calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star)
+               Mslr = max(calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star), min_mfrag)
                regime = COLLRESOLVE_REGIME_HIT_AND_RUN !hit and run
                Qloss = (c_star + 1.0_DP) * U_binding ! Qr 
             end if 
@@ -207,20 +209,26 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb
                Mslr = 0.0_DP
                Qloss = 0.0_DP
             else 
-               Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA)  ! LS12 eq (37)
+               Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag)  ! LS12 eq (37)
                regime = COLLRESOLVE_REGIME_DISRUPTION !disruption
                Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_erosion
             end if 
          else if (Vimp > Vsupercat) then 
-            Mlr = Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(-1.5_DP)   !LS12 eq (44)
-            Mslr = Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA)  !LS12 eq (37)
+            Mlr = max(Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(ETA), min_mfrag)   !LS12 eq (44)
+            Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag)  !LS12 eq (37)
             regime = COLLRESOLVE_REGIME_SUPERCATASTROPHIC ! supercatastrophic
             Qloss = (c_star + 1.0_DP) * U_binding ! Qr - Qr_supercat
          else 
             write(*,*) "Error no regime found in symba_regime"
          end if 
       end if 
 
+      if (Mslr > Mlr) then ! The second-largest fragment is actually larger than the largest, so we will swap them
+         Mtmp = Mlr
+         Mlr = Mslr
+         Mslr = Mtmp
+      end if
+
       Mresidual = Mtot - Mlr - Mslr
       if (Mresidual < 0.0_DP) then ! prevents final masses from going negative
          Mlr = Mlr + Mresidual

src/fraggle/fraggle_set.f90

@@ -27,60 +27,106 @@ module subroutine fraggle_set_budgets_fragments(self, colliders)
    end subroutine fraggle_set_budgets_fragments
 
 
-   module subroutine fraggle_set_mass_dist_fragments(self, colliders)
+   module subroutine fraggle_set_mass_dist_fragments(self, colliders, param)
       !! author: David A. Minton
       !!
       !! Sets the mass of fragments based on the mass distribution returned by the regime calculation.
       !! This subroutine must be run after the the setup rourtine has been run on the fragments
       implicit none
       ! Arguments
-      class(fraggle_fragments), intent(inout) :: self      !! Fraggle fragment system object
-      class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object
+      class(fraggle_fragments),     intent(inout) :: self      !! Fraggle fragment system object
+      class(fraggle_colliders),     intent(inout) :: colliders !! Fraggle collider system object
+      class(swiftest_parameters),   intent(in)    :: param     !! Current Swiftest run configuration parameters
       ! Internals
-      integer(I4B)              :: i, jproj, jtarg
+      integer(I4B)              :: i, jproj, jtarg, nfrag, istart
       real(DP), dimension(2)    :: volume
       real(DP), dimension(NDIM) :: Ip_avg
+      real(DP) :: mfrag, mremaining, min_mfrag
+      real(DP), parameter :: BETA = 2.85_DP
+      integer(I4B), parameter :: NFRAGMAX = 100  !! Maximum number of fragments that can be generated
+      integer(I4B), parameter :: NFRAGMIN = 7 !! 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
 
-      associate(frag => self, nfrag => self%nbody)
+      associate(frag => self)
          ! Get mass weighted mean of Ip and density
          volume(1:2) = 4._DP / 3._DP * PI * colliders%radius(1:2)**3
-         frag%density(1:nfrag) = frag%mtot / sum(volume(:))
          Ip_avg(:) = (colliders%mass(1) * colliders%Ip(:,1) + colliders%mass(2) * colliders%Ip(:,2)) / frag%mtot
-         do i = 1, nfrag
-            frag%Ip(:, i) = Ip_avg(:)
-         end do
+         if (colliders%mass(1) > colliders%mass(2)) then
+            jtarg = 1
+            jproj = 2
+         else
+            jtarg = 2
+            jproj = 1
+         end if
+
+         select type(param)
+         class is (symba_parameters)
+            min_mfrag = (param%min_GMfrag / param%GU) 
+            nfrag = ceiling(NFRAG_SIZE_MULTIPLIER  * log(frag%mtot / min_mfrag))
+            nfrag = max(min(nfrag, NFRAGMAX), NFRAGMIN)
+         class default
+            min_mfrag = 0.0_DP
+            nfrag = NFRAGMAX
+         end select
 
          select case(frag%regime)
          case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
-            ! Assign the first and second-largest masses that came out of the regime determination subroutine to the first two fragments
-            frag%mass(1:2) = frag%mass_dist(1:2)
+            istart = 2
+         case(COLLRESOLVE_REGIME_HIT_AND_RUN)
+            istart = 1
+         case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) 
+            call frag%setup(1, param)
+            frag%mass(1) = frag%mass_dist(1)
+            frag%radius(1) = colliders%radius(jtarg)
+            frag%density(1) = frag%mass_dist(1) / volume(jtarg)
+            frag%Ip(:, 1) = colliders%Ip(:,1)
+            return
+         case default
+            write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",frag%regime
+         end select
 
-            ! Distribute remaining mass among the remaining fragments
-            do i = 3, nfrag
-               frag%mass(i) = (frag%mtot - sum(frag%mass(1:2))) / (nfrag - 2) 
-            end do
-            frag%radius(1:nfrag) = (3 * frag%mass(1:nfrag) / (4 * PI * frag%density(1:nfrag)))**(1.0_DP / 3.0_DP)
+         i = istart 
+         mremaining = max(frag%mtot - sum(frag%mass_dist(1:istart)), 0.0_DP)
+         do while (i <= NFRAGMAX)
+            mfrag = (i - istart + 1)**(-3._DP / BETA) * frag%mass_dist(istart)
+            if (mremaining - mfrag < 0.0_DP) exit
+            mremaining = mremaining - mfrag
+            i = i + 1
+         end do
+
+         call frag%setup(nfrag, param)
+         frag%mass(1:istart) = frag%mass_dist(1:istart)
+         mremaining = max(frag%mtot - sum(frag%mass_dist(1:istart)), 0.0_DP)
+         do i = istart + 1, nfrag
+            mfrag = (i - istart + 1)**(-3._DP / BETA) * frag%mass_dist(istart)
+            mfrag = min(mfrag, mremaining)
+            frag%mass(i) = mfrag
+            mremaining = mremaining - mfrag
+         end do
+         select case(frag%regime)
          case(COLLRESOLVE_REGIME_HIT_AND_RUN) 
-            if (colliders%mass(1) > colliders%mass(2)) then
-               jtarg = 1
-               jproj = 2
-            else
-               jtarg = 2
-               jproj = 1
-            end if
             frag%mass(1) = frag%mass_dist(1)
             frag%radius(1) = colliders%radius(jtarg)
             frag%density(1) = frag%mass_dist(1) / volume(jtarg)
-
-            frag%mass(2:nfrag) = (frag%mtot - frag%mass(1)) / (nfrag - 1) 
-            frag%mass(nfrag) = frag%mass(nfrag) + (frag%mtot - sum(frag%mass(:)))
-            frag%radius(2:nfrag) = (3 * frag%mass(2:nfrag) / (4 * PI * frag%density(2:nfrag)))**(1.0_DP / 3.0_DP)
-         case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) 
-            frag%mass(1) = frag%mtot
-            frag%radius(1) = (3 * sum(volume(:)) / (4 * PI))**(1._DP / 3._DP)
+            frag%Ip(:, 1) = colliders%Ip(:,1)
+            istart = 2
          case default
-            write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",frag%regime
+            istart = 1
          end select
+         if (mremaining > 0.0_DP) then
+            ! Distribute remaining mass among the fragments
+            mfrag = 1._DP + mremaining / sum(frag%mass(istart:nfrag))
+            do i = istart, nfrag
+               frag%mass(i) = frag%mass(i) * mfrag
+            end do
+            mremaining = frag%mtot - sum(frag%mass(1:nfrag))
+            frag%mass(nfrag) = frag%mass(nfrag) + mremaining
+         end if
+         frag%density(istart:nfrag) = frag%mtot / sum(volume(:))
+         frag%radius(istart:nfrag) = (3 * frag%mass(istart:nfrag) / (4 * PI * frag%density(istart:nfrag)))**(1.0_DP / 3.0_DP)
+         do i = istart, nfrag
+            frag%Ip(:, i) = Ip_avg(:)
+         end do
 
       end associate
 
@@ -181,7 +227,6 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders)
 
       associate(frag => self)
          ! Find the center of mass of the collisional system	
-         frag%mtot = sum(colliders%mass(:)) 
          frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot
          frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot
 

src/modules/fraggle_classes.f90

@@ -173,7 +173,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param)
          class(fraggle_colliders),     intent(inout) :: self   !! Fraggle colliders object
          class(fraggle_fragments),     intent(inout) :: frag   !! Fraggle fragment system object
          class(swiftest_nbody_system), intent(in)    :: system !! Swiftest nbody system object
-         class(swiftest_parameters),   intent(in)    :: param  !! Current Swiftest run configuration parameter
+         class(swiftest_parameters),   intent(in)    :: param  !! Current Swiftest run configuration parameters
       end subroutine fraggle_regime_colliders
 
       module subroutine fraggle_set_budgets_fragments(self, colliders)
@@ -188,10 +188,11 @@ module subroutine fraggle_set_coordinate_system(self, colliders)
          class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object
       end subroutine fraggle_set_coordinate_system
 
-      module subroutine fraggle_set_mass_dist_fragments(self, colliders)
+      module subroutine fraggle_set_mass_dist_fragments(self, colliders, param)
          implicit none
-         class(fraggle_fragments), intent(inout) :: self      !! Fraggle fragment system object
-         class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object
+         class(fraggle_fragments),     intent(inout) :: self      !! Fraggle fragment system object
+         class(fraggle_colliders),     intent(inout) :: colliders !! Fraggle collider system object
+         class(swiftest_parameters),   intent(in)    :: param     !! Current Swiftest run configuration parameters
       end subroutine fraggle_set_mass_dist_fragments
 
       module subroutine fraggle_set_natural_scale_factors(self, colliders)

src/symba/symba_collision.f90

@@ -1,8 +1,6 @@
 submodule (symba_classes) s_symba_collision
    use swiftest
 
-   integer(I4B), parameter :: NFRAG_DISRUPT = 12
-   integer(I4B), parameter :: NFRAG_SUPERCAT = 20
 contains
 
    module function symba_collision_casedisruption(system, param, colliders, frag)  result(status)
@@ -26,17 +24,14 @@ module function symba_collision_casedisruption(system, param, colliders, frag)
       select case(frag%regime)
       case(COLLRESOLVE_REGIME_DISRUPTION)
          message = "Disruption between"
-         nfrag = NFRAG_DISRUPT 
       case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
          message = "Supercatastrophic disruption between"
-         nfrag = NFRAG_SUPERCAT
       end select
       call symba_collision_collider_message(system%pl, colliders%idx, message)
       call fraggle_io_log_one_message(message)
 
       ! Collisional fragments will be uniformly distributed around the pre-impact barycenter
-      call frag%setup(nfrag, param)
-      call frag%set_mass_dist(colliders)
+      call frag%set_mass_dist(colliders, param)
 
       ! Generate the position and velocity distributions of the fragments
       call frag%generate_fragments(colliders, system, param, lfailure)
@@ -53,6 +48,7 @@ module function symba_collision_casedisruption(system, param, colliders, frag)
          end select
       else
          ! Populate the list of new bodies
+         nfrag = frag%nbody
          write(message, *) nfrag
          call fraggle_io_log_one_message("Generating " // trim(adjustl(message)) // " fragments")
          select case(frag%regime)
@@ -106,10 +102,8 @@ module function symba_collision_casehitandrun(system, param, colliders, frag)  r
          nfrag = 0
          lpure = .true.
       else ! Imperfect hit and run, so we'll keep the largest body and destroy the other
-         nfrag = NFRAG_DISRUPT - 1
          lpure = .false.
-         call frag%setup(nfrag, param)
-         call frag%set_mass_dist(colliders)
+         call frag%set_mass_dist(colliders, param)
 
          ! Generate the position and velocity distributions of the fragments
          call frag%generate_fragments(colliders, system, param, lpure)
@@ -118,6 +112,7 @@ module function symba_collision_casehitandrun(system, param, colliders, frag)  r
             call fraggle_io_log_one_message("Should have been a pure hit and run instead")
             nfrag = 0
          else
+            nfrag = frag%nbody
             write(message, *) nfrag
             call fraggle_io_log_one_message("Generating " // trim(adjustl(message)) // " fragments")
          end if
@@ -174,8 +169,7 @@ module function symba_collision_casemerge(system, param, colliders, frag)  resul
       select type(pl => system%pl)
       class is (symba_pl)
 
-         call frag%setup(1, param)
-         call frag%set_mass_dist(colliders)
+         call frag%set_mass_dist(colliders, param)
          ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1))
          frag%id(1) = pl%id(ibiggest)
          frag%xb(:,1) = frag%xbcom(:)
@@ -723,7 +717,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status)
             associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody)
                ! Add the colliders%idx bodies to the subtraction list
                ncolliders = colliders%ncoll
-               nfrag   = size(frag%mass(:))
+               nfrag = frag%nbody
 
                ! Setup new bodies
                allocate(plnew, mold=pl)