Moved the fragment position and initial velocity distribution code in…

…to the "simple disruption" model. This will allow these simpler distributions to be tested before adding in the energy and momentum constraints that Fraggle uses.

examples/Fragmentation/Fragmentation_Movie.py

@@ -206,8 +206,8 @@ def data_stream(self, frame=0):
         # Set fragmentation parameters
         minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body
         gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades
-        sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
+        sim.set_parameter(collision_model="simple", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
         sim.run(dt=1e-3, tstop=1.0e-3, istep_out=1, dump_cadence=0)
 
         print("Generating animation")
-        anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1)
+        anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1)

src/collision/collision_generate.f90

@@ -123,7 +123,7 @@ end subroutine collision_generate_hitandrun
    module subroutine collision_generate_merge(self, nbody_system, param, t)
       !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!
-      !! Merge massive bodies in any collisionals ystem.
+      !! Merge massive bodies in any collisional system.
       !! 
       !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90
       !!
@@ -221,11 +221,146 @@ end subroutine collision_generate_merge
 
 
    module subroutine collision_generate_simple_disruption(self, nbody_system, param, t)
+      !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
+      !!
+      !! Generates a simple fragment position and velocity distribution based on the collision 
+      !! regime. It makes no attempt to constrain the energy of the collision
       implicit none
+      ! Arguments
       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
+      ! Internals
+      real(DP) :: r_max_start
+
+      associate(impactors => self%impactors, fragments => self%fragments,  nfrag => self%fragments%nbody)
+         r_max_start = 1.1_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1))
+         call collision_generate_simple_pos_vec(self, r_max_start)
+         call self%set_coordinate_system()
+         call collision_generate_simple_vel_vec(self)
+      end associate
+
+      return
    end subroutine collision_generate_simple_disruption
 
+
+   module subroutine collision_generate_simple_pos_vec(collider, r_max_start)
+      !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
+      !!
+      !! Initializes the position vectors of the fragments around the center of mass based on the collision style.
+      !! For hit and run with disruption, the fragments are generated in a random cloud around the smallest of the two colliders (body 2)
+      !! For disruptive collisions, the fragments are generated in a random cloud around the impact point. Bodies are checked for overlap and
+      !! regenerated if they overlap.
+      implicit none
+      ! Arguments
+      class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object
+      real(DP),                           intent(in)    :: r_max_start    !! The maximum radial distance of fragments for disruptive collisions
+      ! Internals
+      real(DP)  :: dis, rad, r_max, fdistort
+      logical, dimension(:), allocatable :: loverlap
+      integer(I4B) :: i, j
+      logical :: lnoncat, lhitandrun
+
+      associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
+         allocate(loverlap(nfrag))
+
+         lnoncat = (impactors%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point
+         lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) ! Disruptive hit and runs have their own fragment distribution
+
+         ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies
+         ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down)
+         r_max = r_max_start
+         rad = sum(impactors%radius(:))
+
+         ! This is a factor that will "distort" the shape of the fragment cloud in the direction of the impact velocity 
+         fdistort = .mag. (impactors%y_unit(:) .cross. impactors%v_unit(:)) 
+
+         ! We will treat the first two fragments of the list as special cases. They get initialized at the original positions of the impactor bodies
+         fragments%rc(:, 1) = impactors%rb(:, 1) - impactors%rbcom(:) 
+         fragments%rc(:, 2) = impactors%rb(:, 2) - impactors%rbcom(:)
+         call random_number(fragments%rc(:,3:nfrag))
+         loverlap(:) = .true.
+         do while (any(loverlap(3:nfrag)))
+            if (lhitandrun) then ! For a hit-and-run with disruption, the fragment cloud size is based on the radius of the disrupted body
+               r_max = 2 * impactors%radius(2) 
+            else ! For disruptions, the the fragment cloud size is based on the mutual collision system
+               r_max = r_max + 0.1_DP * rad
+            end if
+            do i = 3, nfrag
+               if (loverlap(i)) then 
+                  call random_number(fragments%rc(:,i))
+                  fragments%rc(:,i) = 2 * (fragments%rc(:, i) - 0.5_DP)  
+                  fragments%rc(:, i) = fragments%rc(:,i) + fdistort * impactors%v_unit(:) 
+                  fragments%rc(:, i) = r_max * fragments%rc(:, i)  
+                  fragments%rc(:, i) = fragments%rc(:, i) + (impactors%rbimp(:) - impactors%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM
+                  if (lnoncat .and. dot_product(fragments%rc(:,i), impactors%y_unit(:)) < 0.0_DP) fragments%rc(:, i) = -fragments%rc(:, i) ! Make sure the fragment cloud points away from the impact point
+               end if
+            end do
+
+            ! Check for any overlapping bodies.
+            loverlap(:) = .false.
+            do j = 1, nfrag
+               do i = j + 1, nfrag
+                  dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i))
+                  loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) 
+               end do
+            end do
+         end do
+         call collision_util_shift_vector_to_origin(fragments%mass, fragments%rc)
+         call collider%set_coordinate_system()
+
+         do concurrent(i = 1:nfrag)
+            fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:)
+         end do
+
+         impactors%rbcom(:) = 0.0_DP
+         do concurrent(i = 1:nfrag)
+            impactors%rbcom(:) = impactors%rbcom(:) + fragments%mass(i) * fragments%rb(:,i) 
+         end do
+         impactors%rbcom(:) = impactors%rbcom(:) / fragments%mtot
+      end associate
+
+      return
+   end subroutine collision_generate_simple_pos_vec
+
+
+   module subroutine collision_generate_simple_vel_vec(collider)
+      !! Author:  David A. Minton
+      !!
+      !! Generates an initial "guess" for the velocitity vectors 
+      implicit none
+      ! Arguments
+      class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object
+      ! Internals
+      integer(I4B) :: i, j
+      logical :: lcat
+      real(DP), dimension(NDIM) :: vimp_unit
+      real(DP), dimension(NDIM,collider%fragments%nbody) :: vnoise
+      real(DP), parameter :: VNOISE_MAG = 0.10_DP
+      real(DP) :: vmag
+
+      associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
+         lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) 
+
+         ! "Bounce" the first two bodies
+         do i = 1,2
+            fragments%vc(:,i) = impactors%vb(:,i) - impactors%vbcom(:)  - 2 * impactors%vbimp(:)
+         end do
+
+         vmag = .mag.impactors%vbcom(:) 
+         vimp_unit(:) = .unit. impactors%vbimp(:)
+         call random_number(vnoise)
+         vnoise = (2 * vnoise - 1.0_DP) * vmag
+         do i = 3, nfrag
+            vimp_unit(:) = .unit. (fragments%rc(:,i) + impactors%rbcom(:) - impactors%rbimp(:))
+            fragments%vc(:,i) = vmag * vimp_unit(:) + vnoise(:,i)
+         end do
+      end associate
+      return
+   end subroutine collision_generate_simple_vel_vec
+
+
+
+
 end submodule s_collision_generate

src/collision/collision_module.f90

@@ -233,6 +233,17 @@ module subroutine collision_generate_simple_disruption(self, nbody_system, param
          real(DP),                 intent(in)    :: t            !! The time of the collision
       end subroutine collision_generate_simple_disruption
 
+      module subroutine collision_generate_simple_pos_vec(collider, r_max_start)
+         implicit none
+         class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object
+         real(DP),                           intent(in)    :: r_max_start    !! The maximum radial distance of fragments for disruptive collisions
+      end subroutine collision_generate_simple_pos_vec 
+
+      module subroutine collision_generate_simple_vel_vec(collider)
+         implicit none
+         class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object
+      end subroutine collision_generate_simple_vel_vec
+
       module subroutine collision_io_collider_message(pl, collidx, collider_message)
          implicit none
          class(base_object),            intent(in)    :: pl               !! Swiftest massive body object
@@ -349,32 +360,6 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec)
          integer(I4B),               intent(in)    :: irec   !! Current recursion level
       end subroutine collision_resolve_pltp
 
-      module subroutine collision_util_set_coordinate_collider(self)
-         implicit none
-         class(collision_merge), intent(inout) :: self      !! collisional system
-      end subroutine collision_util_set_coordinate_collider
-
-      module subroutine collision_util_set_coordinate_impactors(self)
-         implicit none
-         class(collision_impactors), intent(inout) :: self      !! collisional system
-      end subroutine collision_util_set_coordinate_impactors
-
-      module subroutine collision_util_setup_collider(self, nbody_system)
-         implicit none
-         class(collision_merge),   intent(inout) :: self         !! Encounter collision system object
-         class(base_nbody_system), intent(in)    :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots
-      end subroutine collision_util_setup_collider
-
-      module subroutine collision_util_setup_impactors_collider(self)
-         implicit none
-         class(collision_merge), intent(inout) :: self   !! Encounter collision system object
-      end subroutine collision_util_setup_impactors_collider
-
-      module subroutine collision_util_setup_fragments_collider(self, nfrag)
-         implicit none
-         class(collision_merge), intent(inout) :: self  !! Encounter collision system object
-         integer(I4B),            intent(in)    :: nfrag !! Number of fragments to create
-      end subroutine collision_util_setup_fragments_collider
 
       module subroutine collision_util_add_fragments_to_collider(self, nbody_system, param)
          implicit none
@@ -403,6 +388,39 @@ module subroutine collision_util_set_mass_dist(self, param)
          class(base_parameters),             intent(in)    :: param !! Current Swiftest run configuration parameters
       end subroutine collision_util_set_mass_dist
 
+      module subroutine collision_util_set_coordinate_collider(self)
+         implicit none
+         class(collision_merge), intent(inout) :: self      !! collisional system
+      end subroutine collision_util_set_coordinate_collider
+
+      module subroutine collision_util_set_coordinate_impactors(self)
+         implicit none
+         class(collision_impactors), intent(inout) :: self      !! collisional system
+      end subroutine collision_util_set_coordinate_impactors
+
+      module subroutine collision_util_setup_collider(self, nbody_system)
+         implicit none
+         class(collision_merge),   intent(inout) :: self         !! Encounter collision system object
+         class(base_nbody_system), intent(in)    :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots
+      end subroutine collision_util_setup_collider
+
+      module subroutine collision_util_setup_impactors_collider(self)
+         implicit none
+         class(collision_merge), intent(inout) :: self   !! Encounter collision system object
+      end subroutine collision_util_setup_impactors_collider
+
+      module subroutine collision_util_setup_fragments_collider(self, nfrag)
+         implicit none
+         class(collision_merge), intent(inout) :: self  !! Encounter collision system object
+         integer(I4B),           intent(in)    :: nfrag !! Number of fragments to create
+      end subroutine collision_util_setup_fragments_collider
+
+      module subroutine collision_util_shift_vector_to_origin(m_frag, vec_frag)
+         implicit none
+         real(DP), dimension(:),   intent(in)    :: m_frag    !! Fragment masses
+         real(DP), dimension(:,:), intent(inout) :: vec_frag  !! Fragment positions or velocities in the center of mass frame
+      end subroutine
+
       module subroutine collision_util_get_idvalues_snapshot(self, idvals)
          implicit none
          class(collision_snapshot),               intent(in)  :: self   !! Collision snapshot object

src/collision/collision_util.f90

@@ -459,8 +459,8 @@ module subroutine collision_util_set_coordinate_impactors(self)
          ! Find the point of impact between the two bodies
          impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * impactors%y_unit(:)
 
-         ! The "bounce" unit vector is the projection of 
-         impactors%vbimp(:) = dot_product(delta_v(:),impactors%x_unit(:)) * impactors%x_unit(:)
+         ! The "bounce" unit vector is the projection of the velocity vector into the distance vector
+         impactors%vbimp(:) = dot_product(delta_v(:),impactors%y_unit(:)) * impactors%y_unit(:)
       end associate
 
       return
@@ -656,6 +656,36 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag)
    end subroutine collision_util_setup_fragments_collider
 
 
+   module subroutine collision_util_shift_vector_to_origin(m_frag, vec_frag)
+      !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
+      !!
+      !! Adjusts the position or velocity of the fragments as needed to align them with the origin
+      implicit none
+      ! Arguments
+      real(DP), dimension(:),   intent(in)    :: m_frag    !! Fragment masses
+      real(DP), dimension(:,:), intent(inout) :: vec_frag  !! Fragment positions or velocities in the center of mass frame
+
+      ! Internals
+      real(DP), dimension(NDIM) :: mvec_frag, COM_offset
+      integer(I4B) :: i, nfrag
+      real(DP) :: mtot
+
+      mvec_frag(:) = 0.0_DP
+      mtot = sum(m_frag)
+      nfrag = size(m_frag)
+
+      do i = 1, nfrag
+         mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i)
+      end do
+      COM_offset(:) = -mvec_frag(:) / mtot
+      do i = 1, nfrag 
+         vec_frag(:, i) = vec_frag(:, i) + COM_offset(:)
+      end do
+
+      return
+   end subroutine collision_util_shift_vector_to_origin
+
+
    subroutine collision_util_save_snapshot(collision_history, snapshot)
       !! author: David A. Minton
       !!