Restructured a bit so that the fragment SFD generator is in the main …

…collision module, not Fraggle

src/collision/collision_generate.f90

@@ -219,7 +219,7 @@ end subroutine collision_generate_bounce
 
    module subroutine collision_generate_simple(self, nbody_system, param, t)
       implicit none
-      class(collision_simple),  intent(inout) :: self         !! Simple fragment system object 
+      class(collision_simple_disruption),  intent(inout) :: self         !! Simple fragment system object 
       class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
       class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters 
       real(DP),                 intent(in)    :: t            !! The time of the collision

src/collision/collision_module.f90

@@ -147,11 +147,12 @@ module collision
       final     ::             collision_final_bounce    !! Finalizer will deallocate all allocatables
    end type collision_bounce
 
-   type, extends(collision_merge) :: collision_simple
+   type, extends(collision_merge) :: collision_simple_disruption
    contains 
-      procedure :: generate => collision_generate_simple !! If a collision would result in a disruption [TODO: SOMETHING LIKE CHAMBERS 2012]
-      final     ::             collision_final_simple !! Finalizer will deallocate all allocatables
-   end type collision_simple
+      procedure :: generate      => collision_generate_simple    !! If a collision would result in a disruption [TODO: SOMETHING LIKE CHAMBERS 2012]
+      procedure :: set_mass_dist => collision_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type
+      final     ::                  collision_final_simple       !! Finalizer will deallocate all allocatables
+   end type collision_simple_disruption
 
 
    !! NetCDF dimension and variable names for the enounter save object
@@ -224,7 +225,7 @@ end subroutine collision_generate_bounce
 
       module subroutine collision_generate_simple(self, nbody_system, param, t)
          implicit none
-         class(collision_simple),  intent(inout) :: self         !! Simple fragment nbody_system object 
+         class(collision_simple_disruption),  intent(inout) :: self         !! Simple fragment nbody_system object 
          class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
          class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters 
          real(DP),                 intent(in)    :: t            !! The time of the collision
@@ -389,6 +390,12 @@ module subroutine collision_util_reset_fragments(self)
          class(collision_fragments(*)), intent(inout) :: self
       end subroutine collision_util_reset_fragments
 
+      module subroutine collision_util_set_mass_dist(self, param)
+         implicit none
+         class(collision_simple_disruption), intent(inout) :: self  !! Simple disruption collision object
+         class(base_parameters),             intent(in)    :: param !! Current Swiftest run configuration parameters
+      end subroutine collision_util_set_mass_dist
+
       module subroutine collision_util_get_idvalues_snapshot(self, idvals)
          implicit none
          class(collision_snapshot),               intent(in)  :: self   !! Collision snapshot object
@@ -554,7 +561,7 @@ subroutine collision_final_simple(self)
          !! Finalizer will deallocate all allocatables
          implicit none
          ! Arguments
-         type(collision_simple),  intent(inout) :: self !!  Collision system object
+         type(collision_simple_disruption),  intent(inout) :: self !!  Collision system object
 
          call self%reset()
          if (allocated(self%impactors)) deallocate(self%impactors)

src/collision/collision_setup.f90

@@ -58,11 +58,11 @@ module subroutine collision_setup_fragments_system(self, nfrag)
 
       if (allocated(self%fragments)) deallocate(self%fragments)
       allocate(collision_fragments(nfrag) :: self%fragments)
+      self%fragments%nbody = nfrag
 
       return
    end subroutine collision_setup_fragments_system
 
-
 end submodule s_collision_setup
 
 

src/collision/collision_util.f90

@@ -435,6 +435,142 @@ module subroutine collision_util_set_coordinate_system(self)
    end subroutine collision_util_set_coordinate_system
 
 
+   module subroutine collision_util_set_mass_dist(self, 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 routine has been run on the fragments
+      !!
+      implicit none
+      ! Arguments
+      class(collision_simple_disruption), intent(inout) :: self  !! Fraggle collision system object
+      class(base_parameters),             intent(in)    :: param !! Current Swiftest run configuration parameters
+      ! Internals
+      integer(I4B)              :: i, jproj, jtarg, nfrag, istart
+      real(DP), dimension(2)    :: volume
+      real(DP), dimension(NDIM) :: Ip_avg
+      real(DP) :: mfrag, mremaining, min_mfrag, mtot
+      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
+      integer(I4B), parameter :: iMlr = 1
+      integer(I4B), parameter :: iMslr = 2
+      integer(I4B), parameter :: iMrem = 3
+
+      associate(impactors => self%impactors)
+         ! Get mass weighted mean of Ip and density
+         volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3
+         mtot = sum(impactors%mass(:))
+         Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / mtot
+
+         if (impactors%mass(1) > impactors%mass(2)) then
+            jtarg = 1
+            jproj = 2
+         else
+            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.
+            ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on
+            ! 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)
+               min_mfrag = (param%min_GMfrag / param%GU) 
+               ! 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
+
+            i = iMrem
+            mremaining = impactors%mass_dist(iMrem)
+            do while (i <= nfrag)
+               mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr)
+               if (mremaining - mfrag < 0.0_DP) exit
+               mremaining = mremaining - mfrag
+               i = i + 1
+            end do
+            if (i < nfrag) nfrag = max(i, NFRAGMIN)  ! The sfd would actually give us fewer fragments than our maximum
+            call self%setup_fragments(nfrag)
+
+         case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) 
+
+            call self%setup_fragments(1)
+            select type(fragments => self%fragments)
+            class is (collision_fragments(*))
+               fragments%mass(1) = impactors%mass_dist(1)
+               fragments%radius(1) = impactors%radius(jtarg)
+               fragments%density(1) = impactors%mass_dist(1) / volume(jtarg)
+               if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1)
+            end select
+            return
+         case default
+            write(*,*) "collision_util_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime
+         end select
+
+         select type(fragments => self%fragments)
+         class is (collision_fragments(*))
+            fragments%mtot = mtot
+
+            ! Make the first two bins the same as the Mlr and Mslr values that came from collision_regime
+            fragments%mass(1) = impactors%mass_dist(iMlr) 
+            fragments%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
+               mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr)
+               fragments%mass(i) = mfrag
+               mremaining = mremaining - 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 remains the second largest
+            end if
+            mfrag = 1._DP + mremaining / sum(fragments%mass(istart:nfrag))
+            fragments%mass(istart:nfrag) = fragments%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(fragments%mass(1:nfrag))
+            fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining
+
+            ! Compute physical properties of the new fragments
+            select case(impactors%regime)
+            case(COLLRESOLVE_REGIME_HIT_AND_RUN)  ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment
+               fragments%radius(1) = impactors%radius(jtarg)
+               fragments%density(1) = impactors%mass_dist(iMlr) / volume(jtarg)
+               fragments%Ip(:, 1) = impactors%Ip(:,1)
+               istart = 2
+            case default
+               istart = 1
+            end select
+
+            fragments%density(istart:nfrag) = fragments%mtot / sum(volume(:))
+            fragments%radius(istart:nfrag) = (3 * fragments%mass(istart:nfrag) / (4 * PI * fragments%density(istart:nfrag)))**(1.0_DP / 3.0_DP)
+            do i = istart, nfrag
+               fragments%Ip(:, i) = Ip_avg(:)
+            end do
+
+         end select
+      end associate
+
+      return
+   end subroutine collision_util_set_mass_dist
+
+
    subroutine collision_util_save_snapshot(collision_history, snapshot)
       !! author: David A. Minton
       !!

src/fraggle/fraggle_module.f90

@@ -15,10 +15,8 @@ module fraggle
    implicit none
    public
 
-
    !> Class definition for the variables that describe a collection of fragments by Fraggle barycentric coordinates
    type, extends(collision_fragments) :: fraggle_fragments
-
       real(DP), dimension(nbody) :: v_r_mag   !! Array of radial direction velocity magnitudes of individual fragments 
       real(DP), dimension(nbody) :: v_t_mag   !! Array of tangential direction velocity magnitudes of individual fragments
       real(DP), dimension(nbody) :: v_n_mag   !! Array of normal direction velocity magnitudes of individual fragments
@@ -28,16 +26,15 @@ module fraggle
       real(DP)                   :: ke_spin   !! Spin kinetic energy of all fragments
       real(DP)                   :: ke_budget !! Kinetic energy budget for computing fragment trajectories
       real(DP), dimension(NDIM)  :: L_budget  !! Angular momentum budget for computing fragment trajectories
-
    contains
-      procedure :: get_angular_momentum  => fraggle_util_get_angular_momentum  !! Calcualtes the current angular momentum of the fragments
-      procedure :: reset                 => fraggle_util_reset_fragments       !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate)
-      procedure :: restructure           => fraggle_util_restructure           !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints
-      final     ::                          fraggle_final_fragments       !! Finalizer will deallocate all allocatables
+      procedure :: get_angular_momentum  => fraggle_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments
+      procedure :: reset                 => fraggle_util_reset_fragments      !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate)
+      procedure :: restructure           => fraggle_util_restructure          !! Restructure the inputs after a failed attempt failed to find a set of positions and velocities that satisfy the energy and momentum constraints
+      final     ::                          fraggle_final_fragments           !! Finalizer will deallocate all allocatables
    end type fraggle_fragments
 
 
-   type, extends(collision_merge) :: collision_fraggle
+   type, extends(collision_simple_disruption) :: collision_fraggle
       ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1
       real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor
       real(DP) :: mscale = 1.0_DP !! Mass scale factor
@@ -48,13 +45,12 @@ module fraggle
    contains
       procedure :: generate                   => fraggle_generate_system                 !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum.
       procedure :: set_budgets                => fraggle_set_budgets                     !! Sets the energy and momentum budgets of the fragments based on the collider value
-      procedure :: set_mass_dist              => fraggle_set_mass_dist                   !! Sets the distribution of mass among the fragments depending on the regime type
       procedure :: set_natural_scale          => fraggle_set_natural_scale_factors       !! Scales dimenional quantities to ~O(1) with respect to the collisional system.  
       procedure :: set_original_scale         => fraggle_set_original_scale_factors      !! Restores dimenional quantities back to the original system units
       procedure :: setup_fragments            => fraggle_setup_fragments_system          !! Initializer for the fragments of the collision system. 
       procedure :: construct_temporary_system => fraggle_util_construct_temporary_system !! Constructs temporary n-body system in order to compute pre- and post-impact energy and momentum
       procedure :: reset                      => fraggle_util_reset_system               !! Deallocates all allocatables
-      final     ::                               fraggle_final_system               !! Finalizer will deallocate all allocatables
+      final     ::                               fraggle_final_system                    !! Finalizer will deallocate all allocatables
    end type collision_fraggle  
 
 
@@ -98,12 +94,6 @@ module subroutine fraggle_set_budgets(self)
          class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object
       end subroutine  fraggle_set_budgets
 
-      module subroutine fraggle_set_mass_dist(self, param)
-         implicit none
-         class(collision_fraggle),        intent(inout) :: self  !! Fraggle collision system object
-         class(swiftest_parameters),   intent(in)    :: param !! Current Swiftest run configuration parameters
-      end subroutine fraggle_set_mass_dist
-
       module subroutine fraggle_set_natural_scale_factors(self)
          implicit none
          class(collision_fraggle), intent(inout) :: self  !! Fraggle collision system object
@@ -176,7 +166,6 @@ end function fraggle_util_vmag_to_vb
    end interface
 
    contains
-
 
       subroutine fraggle_final_fragments(self)
          !! author: David A. Minton