Made more improvements to the "simple" model (now called "disruption")

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

@@ -131,7 +131,7 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t)
             ! The simple disruption model (and its extended types allow for non-pure hit and run. 
             !For the basic merge model, all hit and runs are pure
             select type(self) 
-            class is (collision_simple_disruption)
+            class is (collision_disruption)
                if (impactors%mass_dist(2) > 0.9_DP * impactors%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Pure hit and run. No new fragments generated.")
                   nfrag = 0
@@ -198,7 +198,7 @@ module subroutine collision_generate_simple(self, nbody_system, param, t)
       !! 
       implicit none
       ! Arguments
-      class(collision_simple_disruption), intent(inout) :: self
+      class(collision_disruption), intent(inout) :: self
       class(base_nbody_system),           intent(inout) :: nbody_system !! Swiftest nbody system object
       class(base_parameters),             intent(inout) :: param        !! Current run configuration parameters with SyMBA additions
       real(DP),                           intent(in)    :: t            !! Time of collision
@@ -230,6 +230,8 @@ module subroutine collision_generate_simple(self, nbody_system, param, t)
             end select
             call self%set_mass_dist(param) 
             call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)
+            call self%set_budgets()
+            call self%set_coordinate_system()
             call self%disrupt(nbody_system, param, t)
             call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
 
@@ -372,21 +374,19 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail
       !! 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
+      !! regime. 
       implicit none
       ! Arguments
-      class(collision_simple_disruption),  intent(inout) :: self         !! Simple fragment system object 
+      class(collision_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
       logical, optional,                   intent(out)   :: lfailure
       ! Internals
 
-      call self%set_budgets()
       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 collision_generate_simple_vel_vec(self)
 
       return
    end subroutine collision_generate_disrupt
@@ -401,7 +401,7 @@ module subroutine collision_generate_simple_pos_vec(collider)
       !! regenerated if they overlap.
       implicit none
       ! Arguments
-      class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object
+      class(collision_disruption), intent(inout) :: collider !! Fraggle collision system object
       ! Internals
       real(DP)  :: dis
       real(DP), dimension(NDIM,2) :: fragment_cloud_center
@@ -461,6 +461,7 @@ module subroutine collision_generate_simple_pos_vec(collider)
             rdistance = rdistance * fail_scale
             fragment_cloud_radius(:) = fragment_cloud_radius(:) * fail_scale
          end do
+
          call collision_util_shift_vector_to_origin(fragments%mass, fragments%rc)
          call collider%set_coordinate_system()
 
@@ -487,27 +488,33 @@ 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, Lbefore, Lafter, Lspin, rot
-      integer(I4B) :: i, loop
+      real(DP), dimension(NDIM) :: Lbefore, Lafter, Lspin, rotdir
+      real(DP) :: v_init, v_final, mass_init, mass_final, rotmag
+      integer(I4B) :: i
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
 
-         fragments%rot(:,:) = 0.0_DP
-         ! 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))
+         ! Torque the first body based on the change in angular momentum betwen the pre- and post-impact system assuming a simple bounce
+         mass_init = impactors%mass(2)
+         mass_final = sum(fragments%mass(2:nfrag))
+         v_init = .mag.(impactors%vb(:,2) - impactors%vb(:,1))
+         v_final = sqrt(mass_init / mass_final * v_init**2 - 2 * impactors%Qloss / mass_final)
 
-         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
+         Lbefore(:) = mass_init * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1))
+
+         Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:))
          Lspin(:) = impactors%Lspin(:,1) + (Lbefore(:) - Lafter(:))
-
          fragments%rot(:,1) = Lspin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1))
 
-         ! Randomize the rotational vector direction of the n>1 fragments and distribute the residual momentum amongst them
-         call random_number(fragments%rot(:,2:nfrag))
-         fragments%rot(:,2:nfrag) = .unit. fragments%rot(:,2:nfrag) * .mag. fragments%rot(:,1)
-         call fragments%set_spins()
+         ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random
+         do concurrent(i = 2:nfrag)
+            call random_number(rotdir)
+            rotdir = rotdir - 0.5_DP
+            rotdir = .unit. rotdir
+            call random_number(rotmag)
+            rotmag = rotmag * .mag. (Lbefore(:) - Lafter(:)) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i))
+            fragments%rot(:,i) = rotmag * rotdir
+         end do
 
       end associate
 
@@ -523,13 +530,13 @@ module subroutine collision_generate_simple_vel_vec(collider)
       !! 2x the escape velocity of the smallest of the two bodies.
       implicit none
       ! Arguments
-      class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object
+      class(collision_disruption), intent(inout) :: collider !! Fraggle collision system object
       ! Internals
       integer(I4B) :: i,j
       logical :: lhitandrun, lnoncat
-      real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot
+      real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear
       real(DP), dimension(2) :: vimp
-      real(DP) :: vmag, vdisp
+      real(DP) :: vmag, vdisp, Lmag
       integer(I4B), dimension(collider%fragments%nbody) :: vsign
       real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale
 
@@ -563,11 +570,11 @@ module subroutine collision_generate_simple_vel_vec(collider)
          call random_number(mass_vscale)
          mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2
          mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP)
-         mass_vscale(:) = sqrt(2*mass_vscale(:) / maxval(mass_vscale(:)))
+         mass_vscale(:) = 2*mass_vscale(:) / maxval(mass_vscale(:))
          fragments%vc(:,1) = .mag.impactors%vc(:,1) * impactors%bounce_unit(:) 
          do concurrent(i = 2:nfrag)
             j = fragments%origin_body(i)
-            vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rb(:,j) + impactors%rbcom(:))
+            vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j))
             vmag = .mag.impactors%vc(:,j) * vscale(i) * mass_vscale(i)
             if (lhitandrun) then
                fragments%vc(:,i) = vmag * 0.5_DP * impactors%bounce_unit(:) * vsign(i) + vrot(:)
@@ -578,6 +585,18 @@ module subroutine collision_generate_simple_vel_vec(collider)
                fragments%vc(:,i) = vmag * 0.5_DP * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:)
             end if
          end do
+         call collider%set_coordinate_system()
+         ! Check for any residual angular momentum, and if there is any, put it into velocity shear of the fragments
+         call fragments%get_angular_momentum()
+         if (all(fragments%L_budget(:) / (fragments%Lorbit(:) + fragments%Lspin(:)) > epsilon(1.0_DP))) then
+            Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:))
+            do concurrent(i = 3:nfrag)
+               vshear(:) = Lresidual(:) / (nfrag - 2) / (fragments%mass(i) * fragments%rc(:,i) .cross. fragments%v_t_unit(:,i))
+               vshear(:) = .mag.vshear(:) * fragments%v_t_unit(:,i)
+               fragments%vc(:,i) = fragments%vc(:,i) + vshear(:)
+            end do
+         end if
+
          do concurrent(i = 1:nfrag)
             fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:)
          end do
@@ -588,6 +607,9 @@ module subroutine collision_generate_simple_vel_vec(collider)
          end do
          impactors%vbcom(:) = impactors%vbcom(:) / fragments%mtot
 
+         ! Distribute any remaining angular momentum into fragments pin
+         call fragments%set_spins()
+
       end associate
       return
    end subroutine collision_generate_simple_vel_vec

src/collision/collision_module.f90

@@ -165,13 +165,13 @@ module collision
       final     ::             collision_final_bounce    !! Finalizer will deallocate all allocatables
    end type collision_bounce
 
-   type, extends(collision_basic) :: collision_simple_disruption
+   type, extends(collision_basic) :: collision_disruption
    contains 
       procedure :: generate      => collision_generate_simple    !! A simple disruption models that does not constrain energy loss in collisions
       procedure :: disrupt       => collision_generate_disrupt   !! Disrupt the colliders into the fragments
       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
+   end type collision_disruption
 
 
    !! NetCDF dimension and variable names for the enounter save object
@@ -235,7 +235,7 @@ end subroutine collision_generate_bounce
 
       module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfailure)
          implicit none
-         class(collision_simple_disruption), intent(inout) :: self
+         class(collision_disruption), intent(inout) :: self
          class(base_nbody_system),           intent(inout) :: nbody_system !! Swiftest nbody system object
          class(base_parameters),             intent(inout) :: param        !! Current run configuration parameters with SyMBA additions
          real(DP),                           intent(in)    :: t            !! Time of collision
@@ -260,15 +260,15 @@ end subroutine collision_generate_merge
 
       module subroutine collision_generate_simple(self, nbody_system, param, t)
          implicit none
-         class(collision_simple_disruption),  intent(inout) :: self         !! Simple fragment nbody_system object 
+         class(collision_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
       end subroutine collision_generate_simple
 
       module subroutine collision_generate_simple_pos_vec(collider)
          implicit none
-         class(collision_simple_disruption), intent(inout) :: collider !! Collision system object
+         class(collision_disruption), intent(inout) :: collider !! Collision system object
       end subroutine collision_generate_simple_pos_vec 
 
       module subroutine collision_generate_simple_rot_vec(collider)
@@ -278,7 +278,7 @@ end subroutine collision_generate_simple_rot_vec
 
       module subroutine collision_generate_simple_vel_vec(collider)
          implicit none
-         class(collision_simple_disruption), intent(inout) :: collider !! Collision system object
+         class(collision_disruption), intent(inout) :: collider !! Collision system object
       end subroutine collision_generate_simple_vel_vec
 
       module subroutine collision_io_collider_message(pl, collidx, collider_message)
@@ -445,7 +445,7 @@ end subroutine collision_util_set_coordinate_impactors
 
       module subroutine collision_util_set_mass_dist(self, param)
          implicit none
-         class(collision_simple_disruption), intent(inout) :: self  !! Simple disruption collision object
+         class(collision_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
 
@@ -641,7 +641,7 @@ subroutine collision_final_simple(self)
          !! Finalizer will deallocate all allocatables
          implicit none
          ! Arguments
-         type(collision_simple_disruption),  intent(inout) :: self !!  Collision system object
+         type(collision_disruption),  intent(inout) :: self !!  Collision system object
 
          call self%reset()
          if (allocated(self%impactors)) deallocate(self%impactors)

src/collision/collision_util.f90

@@ -432,8 +432,8 @@ module subroutine collision_util_set_budgets(self)
          dEtot = self%Etot(1)
          dL(:) = self%Ltot(:,1)
 
-         fragments%L_budget(:) = -dL(:)
-         fragments%ke_budget = -(dEtot - impactors%Qloss)
+         fragments%L_budget(:) = dL(:)
+         fragments%ke_budget = (dEtot - impactors%Qloss)
 
       end associate
 
@@ -448,28 +448,20 @@ module subroutine collision_util_set_coordinate_collider(self)
       implicit none
       ! Arguments
       class(collision_basic),    intent(inout) :: self      !! Collisional nbody_system
-      ! Internals
-      integer(I4B) :: i
-      real(DP), dimension(NDIM, self%fragments%nbody) :: L_sigma
 
       associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody)
          call impactors%set_coordinate_system() 
 
-         if ((.not.allocated(self%fragments)) .or. (nfrag == 0) .or. (.not.any(fragments%rc(:,:) > 0.0_DP))) return
+         if (.not.allocated(self%fragments)) return
+         if ((nfrag == 0) .or. (.not.any(fragments%rc(:,:) > 0.0_DP))) return
 
          fragments%rmag(:) = .mag. fragments%rc(:,:)
+         fragments%vmag(:) = .mag. fragments%vc(:,:)
 
-         ! Randomize the tangential velocity direction. 
-         ! This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, otherwise we can get an ill-conditioned nbody_system
-         call random_number(L_sigma(:,:)) 
-         do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP)
-            fragments%v_n_unit(:, i) = impactors%z_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP)
-         end do
-
          ! Define the radial, normal, and tangential unit vectors for each individual fragment
          fragments%v_r_unit(:,:) = .unit. fragments%rc(:,:) 
-         fragments%v_n_unit(:,:) = .unit. fragments%v_n_unit(:,:) 
-         fragments%v_t_unit(:,:) = .unit. (fragments%v_n_unit(:,:) .cross. fragments%v_r_unit(:,:))
+         fragments%v_t_unit(:,:) = .unit. fragments%vc(:,:) 
+         fragments%v_n_unit(:,:) = .unit. (fragments%v_r_unit(:,:) .cross. fragments%v_t_unit(:,:))
 
       end associate
 
@@ -546,7 +538,7 @@ module subroutine collision_util_set_mass_dist(self, param)
       !!
       implicit none
       ! Arguments
-      class(collision_simple_disruption), intent(inout) :: self  !! Fraggle collision system object
+      class(collision_disruption), intent(inout) :: self  !! Fraggle collision system object
       class(base_parameters),             intent(in)    :: param !! Current Swiftest run configuration parameters
       ! Internals
       integer(I4B)              :: i, j, jproj, jtarg, nfrag, istart
@@ -712,14 +704,17 @@ module subroutine collision_util_set_spins(self)
 
       call self%get_angular_momentum()
       Lresidual(:) = self%L_budget(:) - (self%Lorbit(:) + self%Lspin(:))
-
       ! Distribute residual angular momentum amongst the fragments
       if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then 
          do i = 2,self%nbody
             self%rot(:,i) = self%rot(:,i) + Lresidual(:) / ((self%nbody - 1) * self%mass(i) * self%radius(i)**2 * self%Ip(3,i)) 
          end do
       end if
 
+      ! Test to see if we were successful
+      call self%get_angular_momentum()
+      Lresidual(:) = self%L_budget(:) - (self%Lorbit(:) + self%Lspin(:))
+
       return
    end subroutine collision_util_set_spins
 
@@ -772,6 +767,30 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag)
       if (allocated(self%fragments)) deallocate(self%fragments)
       allocate(collision_fragments(nfrag) :: self%fragments)
       self%fragments%nbody = nfrag
+      self%fragments%nbody = nfrag
+      self%fragments%status(:) = ACTIVE
+      self%fragments%rh(:,:) = 0.0_DP
+      self%fragments%vh(:,:) = 0.0_DP
+      self%fragments%rb(:,:) = 0.0_DP
+      self%fragments%vb(:,:) = 0.0_DP
+      self%fragments%rc(:,:) = 0.0_DP
+      self%fragments%vc(:,:) = 0.0_DP
+      self%fragments%rot(:,:) = 0.0_DP
+      self%fragments%Ip(:,:) = 0.0_DP
+      self%fragments%v_r_unit(:,:) = 0.0_DP
+      self%fragments%v_t_unit(:,:) = 0.0_DP
+      self%fragments%v_n_unit(:,:) = 0.0_DP
+      self%fragments%mass(:) = 0.0_DP
+      self%fragments%radius(:) = 0.0_DP
+      self%fragments%density(:) = 0.0_DP
+      self%fragments%rmag(:) = 0.0_DP
+      self%fragments%vmag(:) = 0.0_DP
+      self%fragments%Lorbit(:) = 0.0_DP
+      self%fragments%Lspin(:) = 0.0_DP
+      self%fragments%L_budget(:) = 0.0_DP
+      self%fragments%ke_orbit = 0.0_DP
+      self%fragments%ke_spin = 0.0_DP
+      self%fragments%ke_budget = 0.0_DP
 
       return
    end subroutine collision_util_setup_fragments_collider