Rearranged some of the methods in the collision system to keep consis…

…tency when computing spins

src/collision/collision_generate.f90

@@ -382,6 +382,7 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail
       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)
@@ -487,9 +488,7 @@ module subroutine collision_generate_simple_rot_vec(collider)
       class(collision_basic), intent(inout) :: collider !! Fraggle collision system object
       ! Internals
       real(DP), dimension(NDIM) :: Lresidual, Lbefore, Lafter, Lspin, rot
-      real(DP) :: ke_residual, rotmag_old, rotmag_new, vmag_old, vmag_new
       integer(I4B) :: i, loop
-      integer(I4B) :: MAXLOOP = 10
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
 
@@ -505,19 +504,10 @@ module subroutine collision_generate_simple_rot_vec(collider)
 
          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(:)
-
-         if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then 
-            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
+         fragments%rot(:,2:nfrag) = .unit. fragments%rot(:,2:nfrag) * .mag. fragments%rot(:,1)
+         call fragments%set_spins()
 
       end associate
 
@@ -576,12 +566,9 @@ module subroutine collision_generate_simple_vel_vec(collider)
          mass_vscale(:) = sqrt(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(:))
-
             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(:)
             else

src/collision/collision_module.f90

@@ -114,10 +114,13 @@ module collision
       real(DP), dimension(NDIM)                              :: Lspin     !! Spin angular momentum vector of all fragments
       real(DP)                                               :: ke_orbit  !! Orbital kinetic energy of all fragments
       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 :: reset                => collision_util_reset_fragments      !! Deallocates all allocatable arrays and sets everything else to 0
       procedure :: get_angular_momentum => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments
       procedure :: get_kinetic_energy   => collision_util_get_kinetic_energy   !! Calcualtes the current kinetic energy of the fragments
+      procedure :: set_spins            => collision_util_set_spins            !! Calcualtes the spins of all fragments from the angular momentum budget and residual
       final     ::                         collision_final_fragments           !! Finalizer deallocates all allocatables
    end type collision_fragments
 
@@ -149,6 +152,7 @@ module collision
       procedure :: construct_temporary_system => collision_util_construct_temporary_system !! Constructs temporary n-body nbody_system in order to compute pre- and post-impact energy and momentum
       procedure :: get_energy_and_momentum    => collision_util_get_energy_momentum        !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.)
       procedure :: reset                      => collision_util_reset_system               !! Deallocates all allocatables
+      procedure :: set_budgets                => collision_util_set_budgets                !! Sets the energy and momentum budgets of the fragments based on the collider value
       procedure :: set_coordinate_system      => collision_util_set_coordinate_collider    !! Sets the coordinate system of the collisional system
       procedure :: generate                   => collision_generate_basic                  !! Merges the impactors to make a single final body
       procedure :: hitandrun                  => collision_generate_hitandrun              !! Merges the impactors to make a single final body
@@ -393,7 +397,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_add_fragments_to_collider(self, nbody_system, param)
          implicit none
          class(collision_basic),  intent(in)    :: self         !! Collision system object
@@ -410,27 +413,25 @@ module subroutine collision_util_construct_temporary_system(self, nbody_system,
          class(base_parameters),   allocatable, intent(out)   :: tmpparam     !! Output temporary configuration run parameters
       end subroutine collision_util_construct_temporary_system 
 
-
       module subroutine collision_util_get_angular_momentum(self) 
          implicit none
-         class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object
+         class(collision_fragments(*)), intent(inout) :: self !! Fragment system object
       end subroutine collision_util_get_angular_momentum
 
       module subroutine collision_util_get_kinetic_energy(self) 
          implicit none
-         class(collision_fragments(*)), intent(inout) :: self !! Fraggle fragment system object
+         class(collision_fragments(*)), intent(inout) :: self !! Fragment system object
       end subroutine collision_util_get_kinetic_energy
 
       module subroutine collision_util_reset_fragments(self)
          implicit none
          class(collision_fragments(*)), intent(inout) :: self
       end subroutine collision_util_reset_fragments
 
-      module subroutine collision_util_set_mass_dist(self, param)
+      module subroutine collision_util_set_budgets(self)
          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
+         class(collision_basic), intent(inout) :: self !! Collision system object
+      end subroutine  collision_util_set_budgets
 
       module subroutine collision_util_set_coordinate_collider(self)
          implicit none
@@ -442,6 +443,17 @@ module subroutine collision_util_set_coordinate_impactors(self)
          class(collision_impactors), intent(inout) :: self      !! collisional system
       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(base_parameters),             intent(in)    :: param !! Current Swiftest run configuration parameters
+      end subroutine collision_util_set_mass_dist
+
+      module subroutine collision_util_set_spins(self)
+         implicit none
+         class(collision_fragments(*)), intent(inout)  :: self !! Collision fragment system object
+      end subroutine collision_util_set_spins
+
       module subroutine collision_util_setup_collider(self, nbody_system)
          implicit none
          class(collision_basic),   intent(inout) :: self         !! Encounter collision system object
@@ -554,7 +566,7 @@ subroutine collision_final_plpl(self)
          !! Finalizer will deallocate all allocatables
          implicit none
          ! Arguments
-         type(collision_list_plpl),  intent(inout) :: self !! Fraggle encountar storage object
+         type(collision_list_plpl),  intent(inout) :: self !! PL-PL collision list object
 
          call self%dealloc()
 
@@ -567,7 +579,7 @@ subroutine collision_final_pltp(self)
          !! Finalizer will deallocate all allocatables
          implicit none
          ! Arguments
-         type(collision_list_pltp),  intent(inout) :: self !! Fraggle encountar storage object
+         type(collision_list_pltp),  intent(inout) :: self !! PL-TP collision list object
 
          call self%dealloc()
 
@@ -581,7 +593,7 @@ subroutine collision_final_snapshot(self)
          !! Finalizer will deallocate all allocatables
          implicit none
          ! Arguments
-         type(collision_snapshot),  intent(inout) :: self !! Fraggle encountar storage object
+         type(collision_snapshot),  intent(inout) :: self !! Collsion snapshot object
 
          call encounter_final_snapshot(self%encounter_snapshot)
 

src/collision/collision_util.f90

@@ -175,7 +175,6 @@ module subroutine collision_util_get_idvalues_snapshot(self, idvals)
    end subroutine collision_util_get_idvalues_snapshot
 
 
-
    module subroutine collision_util_get_angular_momentum(self) 
       !! Author: David A. Minton
       !!
@@ -417,6 +416,31 @@ module subroutine collision_util_reset_system(self)
    end subroutine collision_util_reset_system
 
 
+   module subroutine collision_util_set_budgets(self)
+      !! author: David A. Minton
+      !!
+      !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum
+      implicit none
+      ! Arguments
+      class(collision_basic), intent(inout) :: self !! Fraggle collision system object
+      ! Internals
+      real(DP) :: dEtot
+      real(DP), dimension(NDIM) :: dL
+
+      associate(impactors => self%impactors, fragments => self%fragments)
+
+         dEtot = self%Etot(1)
+         dL(:) = self%Ltot(:,1)
+
+         fragments%L_budget(:) = -dL(:)
+         fragments%ke_budget = -(dEtot - impactors%Qloss)
+
+      end associate
+
+      return
+   end subroutine collision_util_set_budgets
+
+
    module subroutine collision_util_set_coordinate_collider(self)
       !! author: David A. Minton
       !!
@@ -675,6 +699,31 @@ module subroutine collision_util_set_mass_dist(self, param)
    end subroutine collision_util_set_mass_dist
 
 
+   module subroutine collision_util_set_spins(self)
+      !! Author: David A. Minton
+      !!
+      !! Distributes any residual angular momentum into the spins of the n>1 fragments
+      implicit none
+      ! Arguments
+      class(collision_fragments(*)), intent(inout)  :: self !! Collision fragment system object
+      ! Internals
+      integer(I4B) :: i
+      real(DP), dimension(NDIM) :: Lresidual
+
+      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
+
+      return
+   end subroutine collision_util_set_spins
+
+
    module subroutine collision_util_setup_collider(self, nbody_system)
       !! author: David A. Minton
       !!

src/fraggle/fraggle_generate.f90

@@ -197,7 +197,7 @@ subroutine fraggle_generate_minimize(collider, lfailure)
                end do
 
                ! Set spins in order to conserve angular momentum
-               call fraggle_generate_set_spins(fragments)
+               call collision_util_set_spins(fragments)
 
                if (.not.lfailure) exit
                tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution
@@ -238,7 +238,7 @@ function E_objective_function(val_input) result(fval)
 
                   call collision_util_shift_vector_to_origin(tmp_frag%mass, tmp_frag%vc)
                   ! Set spins in order to conserve angular momentum
-                  call fraggle_generate_set_spins(tmp_frag)
+                  call collision_util_set_spins(tmp_frag)
 
                   ! Get the current kinetic energy of the system
                   call tmp_frag%get_kinetic_energy()
@@ -261,34 +261,4 @@ end function E_objective_function
 
    end subroutine fraggle_generate_minimize
 
-
-   subroutine fraggle_generate_set_spins(fragments)
-      !! Author: David A. Minton
-      !!
-      !! Distributes any residual angular momentum into the spins of the n>1 fragments
-      implicit none
-      ! Arguments
-      class(fraggle_fragments(*)), intent(inout)  :: fragments !! Fraggle fragment system object
-      ! Internals
-      integer(I4B) :: i
-      real(DP), dimension(NDIM) :: Lresidual
-
-      call fragments%get_angular_momentum()
-      Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:))
-
-      ! Distribute residual angular momentum amongst the fragments
-      if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then 
-         do i = 2,fragments%nbody
-            fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / ((fragments%nbody - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) 
-         end do
-      end if
-
-      return
-   end subroutine fraggle_generate_set_spins
-
-
-
-
-
-
 end submodule s_fraggle_generate

src/fraggle/fraggle_module.f90

@@ -20,8 +20,6 @@ module fraggle
       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
-      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 :: 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)
@@ -39,7 +37,6 @@ module fraggle
       real(DP) :: Lscale = 1.0_DP  !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale)
    contains
       procedure :: disrupt                    => fraggle_generate_disrupt                !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum.
-      procedure :: set_budgets                => fraggle_util_set_budgets                !! Sets the energy and momentum budgets of the fragments based on the collider value
       procedure :: set_natural_scale          => fraggle_util_set_natural_scale_factors  !! Scales dimenional quantities to ~O(1) with respect to the collisional system.  
       procedure :: set_original_scale         => fraggle_util_set_original_scale_factors !! Restores dimenional quantities back to the original system units
       procedure :: setup_fragments            => fraggle_util_setup_fragments_system     !! Initializer for the fragments of the collision system. 
@@ -84,11 +81,6 @@ module subroutine fraggle_util_reset_system(self)
          class(collision_fraggle), intent(inout) :: self  !! Collision system object
       end subroutine fraggle_util_reset_system
 
-      module subroutine fraggle_util_set_budgets(self)
-         implicit none
-         class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object
-      end subroutine  fraggle_util_set_budgets
-
       module subroutine fraggle_util_set_natural_scale_factors(self)
          implicit none
          class(collision_fraggle), intent(inout) :: self  !! Fraggle collision system object

src/fraggle/fraggle_util.f90

@@ -87,34 +87,6 @@ module subroutine fraggle_util_reset_system(self)
    end subroutine fraggle_util_reset_system
 
 
-   module subroutine fraggle_util_set_budgets(self)
-      !! author: David A. Minton
-      !!
-      !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum
-      implicit none
-      ! Arguments
-      class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object
-      ! Internals
-      real(DP) :: dEtot
-      real(DP), dimension(NDIM) :: dL
-
-      associate(impactors => self%impactors)
-         select type(fragments => self%fragments)
-         class is (fraggle_fragments(*))
-
-            dEtot = self%Etot(1)
-            dL(:) = self%Ltot(:,1)
-
-            fragments%L_budget(:) = -dL(:)
-            fragments%ke_budget = -(dEtot - impactors%Qloss)
-
-         end select
-      end associate
-
-      return
-   end subroutine fraggle_util_set_budgets
-
-
    module subroutine fraggle_util_set_natural_scale_factors(self)
       !! author: David A. Minton
       !!