Started a simplified Fraggle. There are some issues with it (probably…

… related to energy/momentum budgets) but the infrastructure is in place.

examples/Fragmentation/Fragmentation_Movie.py

@@ -232,7 +232,7 @@ 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="simple", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
+        sim.set_parameter(collision_model="fraggle", 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")

src/collision/collision_generate.f90

@@ -141,7 +141,9 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t)
                   call self%set_mass_dist(param)
 
                   ! Generate the position and velocity distributions of the fragments
+                  call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)
                   call self%disrupt(nbody_system, param, t, lpure)
+                  call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
 
                   dpe = self%pe(2) - self%pe(1) 
                   nbody_system%Ecollisions = nbody_system%Ecollisions - dpe 
@@ -227,7 +229,9 @@ module subroutine collision_generate_simple(self, nbody_system, param, t)
                call util_exit(FAILURE)
             end select
             call self%set_mass_dist(param) 
+            call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)
             call self%disrupt(nbody_system, param, t)
+            call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
 
             dpe = self%pe(2) - self%pe(1) 
             nbody_system%Ecollisions = nbody_system%Ecollisions - dpe 
@@ -376,12 +380,11 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail
       logical, optional,                   intent(out)   :: lfailure
       ! Internals
 
-      call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)
       call collision_generate_simple_pos_vec(self)
       call self%set_coordinate_system()
       call collision_generate_simple_vel_vec(self)
       call collision_generate_simple_rot_vec(self)
-      call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
+
       return
    end subroutine collision_generate_disrupt
 
@@ -481,28 +484,34 @@ module subroutine collision_generate_simple_rot_vec(collider)
       ! Arguments
       class(collision_basic), intent(inout) :: collider !! Fraggle collision system object
       ! Internals
-      real(DP), dimension(NDIM) :: Lresidual
+      real(DP), dimension(NDIM) :: Lresidual, Lbefore, Lafter, Lspin, rot
       integer(I4B) :: i
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
 
          fragments%rot(:,:) = 0.0_DP
-         ! Keep the first two bodies spinning as before to start with
-         Lresidual(:) = 0.0_DP
-         do i = 1,2
-            fragments%rot(:,i) = impactors%Lspin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i))
-            Lresidual(:) = Lresidual(:) + impactors%Lorbit(:,i)
-         end do
+         ! Torque the first body based on the change in angular momentum betwen the pre- and post-impact system
+         Lbefore(:) = impactors%mass(2) * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1))
 
-         ! Compute the current orbital angular momentum
-         do i = 1, nfrag
-            Lresidual(:) = Lresidual(:) - fragments%mass(i) * (fragments%rc(:,i) .cross. fragments%vc(:,i))
+         Lafter(:) = 0.0_DP
+         do i = 2, nfrag
+            Lafter(:) = Lafter(:) + fragments%mass(i) * (fragments%rb(:,i) -  fragments%rb(:,1)) .cross. (fragments%vb(:,i) - fragments%vb(:,1))
          end do
+         Lspin(:) = impactors%Lspin(:,1) + (Lbefore(:) - Lafter(:))
+
+         fragments%rot(:,1) = Lspin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1))
+
+         Lresidual(:) = sum(impactors%Lspin(:,:) + impactors%Lorbit(:,:), dim=2) - Lspin(:)
+
+         ! Randomize the rotational vector direction of the n>1 fragments and distribute the residual momentum amongst them
+         call random_number(fragments%rot(:,2:nfrag))
+         rot(:) = Lresidual(:) / sum(fragments%mass(2:nfrag) * fragments%radius(2:nfrag)**2 * fragments%Ip(3,2:nfrag))
+
+         fragments%rot(:,2:nfrag) = .unit. fragments%rot(:,2:nfrag) * .mag. rot(:)
 
-         ! Distributed most of the remaining angular momentum amongst all the particles
          if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then 
-            do i = 1,nfrag
-               fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / (nfrag * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) 
+            do i = 2,nfrag
+               fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) 
             end do
          end if
 

src/collision/collision_module.f90

@@ -498,8 +498,7 @@ subroutine collision_final_fragments(self)
          ! Arguments
          type(collision_fragments(*)), intent(inout) :: self
 
-         call self%reset()
-
+         if (allocated(self%info)) deallocate(self%info)
          return
       end subroutine collision_final_fragments
 

src/collision/collision_util.f90

@@ -191,8 +191,7 @@ module subroutine collision_util_get_energy_momentum(self,  nbody_system, param,
       ! Internals
       class(base_nbody_system), allocatable, save :: tmpsys
       class(base_parameters), allocatable, save   :: tmpparam
-      integer(I4B)  :: npl_before, npl_after, stage
-
+      integer(I4B)  :: npl_before, npl_after, stage,i
       select type(nbody_system)
       class is (swiftest_nbody_system)
       select type(param)
@@ -226,15 +225,20 @@ module subroutine collision_util_get_energy_momentum(self,  nbody_system, param,
                   ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum
                   call self%add_fragments(tmpsys, tmpparam)
                   tmpsys%pl%status(impactors%id(1:impactors%ncoll)) = INACTIVE
-                  tmpsys%pl%status(npl_before+1:npl_after) = ACTIVE
+                  do concurrent(i = npl_before+1:npl_after)
+                     tmpsys%pl%status(i) = ACTIVE
+                     tmpsys%pl%rot(:,i) = 0.0_DP
+                     tmpsys%pl%vb(:,i) = impactors%vbcom(:)
+                  end do
                end select
             end if 
             select type(tmpsys)
             class is (swiftest_nbody_system)
 
                if (param%lflatten_interactions) call tmpsys%pl%flatten(param)
 
-               call tmpsys%get_energy_and_momentum(param) 
+               call tmpsys%get_energy_and_momentum(param)
+
 
                ! Calculate the current fragment energy and momentum balances
                if (lbefore) then
@@ -248,8 +252,7 @@ module subroutine collision_util_get_energy_momentum(self,  nbody_system, param,
                self%ke_orbit(stage) = tmpsys%ke_orbit
                self%ke_spin(stage) = tmpsys%ke_spin
                self%pe(stage) = tmpsys%pe
-               self%Etot(stage) = tmpsys%te 
-               if (stage == 2) self%Etot(stage) = self%Etot(stage) - (self%pe(2) - self%pe(1)) ! Gotta be careful with PE when number of bodies changes.
+               self%Etot(stage) = tmpsys%ke_orbit + tmpsys%ke_spin
             end select
          end associate
       end select