Refactored name of disruption model

src/collision/collision_generate.f90

@@ -191,7 +191,7 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t)
    end subroutine collision_generate_hitandrun
 
 
-   module subroutine collision_generate_simple(self, nbody_system, param, t)
+   module subroutine collision_generate_disruption(self, nbody_system, param, t)
       !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!
       !! Create the fragments resulting from a non-catastrophic disruption collision
@@ -263,7 +263,7 @@ module subroutine collision_generate_simple(self, nbody_system, param, t)
       end select
       end select
       return
-   end subroutine collision_generate_simple
+   end subroutine collision_generate_disruption
 
 
    module subroutine collision_generate_merge(self, nbody_system, param, t)
@@ -384,15 +384,15 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail
       logical, optional,                   intent(out)   :: lfailure
       ! Internals
 
-      call collision_generate_simple_pos_vec(self)
-      call collision_generate_simple_rot_vec(self)
-      call collision_generate_simple_vel_vec(self)
+      call collision_generate_disruption_pos_vec(self)
+      call collision_generate_disruption_rot_vec(self)
+      call collision_generate_disruption_vel_vec(self)
 
       return
    end subroutine collision_generate_disrupt
 
 
-   module subroutine collision_generate_simple_pos_vec(collider)
+   module subroutine collision_generate_disruption_pos_vec(collider)
       !! 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.
@@ -477,10 +477,10 @@ module subroutine collision_generate_simple_pos_vec(collider)
       end associate
 
       return
-   end subroutine collision_generate_simple_pos_vec
+   end subroutine collision_generate_disruption_pos_vec
 
 
-   module subroutine collision_generate_simple_rot_vec(collider)
+   module subroutine collision_generate_disruption_rot_vec(collider)
       !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!
       !! Calculates the spins of a collection of fragments such that they conserve angular momentum 
@@ -519,10 +519,10 @@ module subroutine collision_generate_simple_rot_vec(collider)
       end associate
 
       return
-   end subroutine collision_generate_simple_rot_vec
+   end subroutine collision_generate_disruption_rot_vec
 
 
-   module subroutine collision_generate_simple_vel_vec(collider)
+   module subroutine collision_generate_disruption_vel_vec(collider)
       !! Author:  David A. Minton
       !!
       !! Generates an initial velocity distribution. For disruptions, the velocity magnitude is set to be
@@ -644,7 +644,7 @@ module subroutine collision_generate_simple_vel_vec(collider)
 
       end associate
       return
-   end subroutine collision_generate_simple_vel_vec
+   end subroutine collision_generate_disruption_vel_vec
 
 
 end submodule s_collision_generate

src/collision/collision_module.f90

@@ -167,7 +167,7 @@ module collision
 
    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 :: generate      => collision_generate_disruption    !! 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
@@ -258,28 +258,28 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
          real(DP),                 intent(in)    :: t            !! The time of the collision
       end subroutine collision_generate_merge
 
-      module subroutine collision_generate_simple(self, nbody_system, param, t)
+      module subroutine collision_generate_disruption(self, nbody_system, param, t)
          implicit none
          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
+      end subroutine collision_generate_disruption
 
-      module subroutine collision_generate_simple_pos_vec(collider)
+      module subroutine collision_generate_disruption_pos_vec(collider)
          implicit none
          class(collision_disruption), intent(inout) :: collider !! Collision system object
-      end subroutine collision_generate_simple_pos_vec 
+      end subroutine collision_generate_disruption_pos_vec 
 
-      module subroutine collision_generate_simple_rot_vec(collider)
+      module subroutine collision_generate_disruption_rot_vec(collider)
          implicit none
          class(collision_basic), intent(inout) :: collider !! Collision system object
-      end subroutine collision_generate_simple_rot_vec 
+      end subroutine collision_generate_disruption_rot_vec 
 
-      module subroutine collision_generate_simple_vel_vec(collider)
+      module subroutine collision_generate_disruption_vel_vec(collider)
          implicit none
          class(collision_disruption), intent(inout) :: collider !! Collision system object
-      end subroutine collision_generate_simple_vel_vec
+      end subroutine collision_generate_disruption_vel_vec
 
       module subroutine collision_io_collider_message(pl, collidx, collider_message)
          implicit none

src/collision/collision_util.f90

@@ -700,11 +700,14 @@ module subroutine collision_util_set_spins(self)
       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)) 
+         do i = 1,self%nbody
+            self%rot(:,i) = self%rot(:,i) + Lresidual(:) / (self%nbody * self%mass(i) * self%radius(i)**2 * self%Ip(3,i)) 
          end do
       end if
 
+      call self%get_angular_momentum()
+      Lresidual(:) = self%L_budget(:) - (self%Lorbit(:) + self%Lspin(:))
+
       return
    end subroutine collision_util_set_spins
 

src/fraggle/fraggle_generate.f90

@@ -93,9 +93,9 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
             call self%set_budgets()
 
             ! Minimize difference between energy/momentum and budgets
-            call fraggle_generate_minimize(self, lfailure_local)
-            call fraggle_generate_tan_vel(self, lfailure_local)
-            ! call fraggle_generate_rad_vel(self, lfailure_local)
+            ! call fraggle_generate_minimize(self, lfailure_local)
+            ! call fraggle_generate_tan_vel(self, lfailure_local)
+            call fraggle_generate_rad_vel(self, lfailure_local)
 
             call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
             dEtot = self%Etot(2) - self%Etot(1)
@@ -188,6 +188,7 @@ subroutine fraggle_generate_minimize(collider, lfailure)
             fval = E_objective_function(input_v)
 
             call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v)
+            fval = E_objective_function(output_v)
             fragments%vmag(1:nfrag) = output_v(1:nfrag)
             do concurrent(i=1:nfrag)
                fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i)
@@ -287,7 +288,7 @@ subroutine fraggle_generate_tan_vel(collider, lfailure)
       type(lambda_obj_err)                 :: objective_function
       real(DP), dimension(NDIM)            :: Li, L_remainder, L_frag_tot
       character(len=STRMAX)                :: message
-      real(DP), dimension(:), allocatable :: output_v
+      real(DP), dimension(:), allocatable :: v_t_output
       logical                             :: lfirst_func
 
       associate(impactors => collider%impactors,  nfrag => collider%fragments%nbody)
@@ -310,10 +311,11 @@ subroutine fraggle_generate_tan_vel(collider, lfailure)
 
                ! ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum
                objective_function = lambda_obj(tangential_objective_function, lfailure)
+               fval = tangential_objective_function(v_t_output(:), lfailure)
 
-               call minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure, output_v)
-               fval = tangential_objective_function(output_v(:), lfailure)
-               fragments%v_t_mag(7:nfrag) = output_v(:)
+               call minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure, v_t_output)
+               fval = tangential_objective_function(v_t_initial(:), lfailure)
+               fragments%v_t_mag(7:nfrag) = v_t_output(:)
                ! Now that the KE-minimized values of the i>6 fragments are found, calculate the momentum-conserving solution for tangential velociteis
                v_t_initial(7:nfrag) = fragments%v_t_mag(7:nfrag)
 
@@ -474,33 +476,35 @@ subroutine fraggle_generate_rad_vel(collider, lfailure)
       real(DP), parameter                   :: TOL_INIT = 1e-14_DP
       real(DP), parameter                   :: VNOISE_MAG = 1e-10_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure
       integer(I4B), parameter               :: MAXLOOP = 100
-      real(DP)                              :: ke_radial, tol 
+      real(DP)                              :: ke_radial, tol, fval
       integer(I4B)                          :: i
       real(DP), dimension(:), allocatable   :: v_r_initial
       real(DP), dimension(collider%fragments%nbody)       :: vnoise
       type(lambda_obj)                      :: objective_function
       character(len=STRMAX)                 :: message
-      real(DP), dimension(:), allocatable :: output_v
+      real(DP), dimension(:), allocatable :: v_r_output
 
       associate(impactors => collider%impactors,  nfrag => collider%fragments%nbody)
          select type(fragments => collider%fragments)
          class is (fraggle_fragments(*))
             ! Set the "target" ke for the radial component
             ke_radial = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit
+
+            do i = 1, nfrag
+               fragments%v_t_mag(i) = dot_product(fragments%vc(:,i), fragments%v_t_unit(:,i))
+               fragments%v_r_mag(i) = dot_product(fragments%vc(:,i), fragments%v_r_unit(:,i))
+            end do
 
             allocate(v_r_initial, source=fragments%v_r_mag)
-            ! Initialize radial velocity magnitudes with a random value that related to equipartition of kinetic energy with some noise
-            call random_number(vnoise(1:nfrag))
-            vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1.0_DP) 
-            v_r_initial(1:nfrag) = sqrt(abs(2 * ke_radial) / (fragments%mass(1:nfrag) * nfrag)) * vnoise(1:nfrag)
 
             ! Initialize the lambda function using a structure constructor that calls the init method
             ! Minimize the ke objective function using the BFGS optimizer
             objective_function = lambda_obj(radial_objective_function)
             tol = TOL_INIT
             do while(tol < TOL_MIN)
-               call minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lfailure, output_v)
-               fragments%v_r_mag = output_v
+               fval = radial_objective_function(v_r_initial)
+               call minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lfailure, v_r_output)
+               fragments%v_r_mag = v_r_output
                if (.not.lfailure) exit
                tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution
                v_r_initial(:) = fragments%v_r_mag(:)
@@ -515,7 +519,7 @@ subroutine fraggle_generate_rad_vel(collider, lfailure)
                                  fragments%v_t_mag(1:nfrag), fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) 
             do i = 1, nfrag
                fragments%vc(:, i) = fragments%vb(:, i) - impactors%vbcom(:)
-               fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vb(:, i), fragments%vb(:, i))
+               fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:, i), fragments%vc(:, i))
             end do
             fragments%ke_orbit = 0.5_DP * fragments%ke_orbit
 
@@ -559,6 +563,7 @@ function radial_objective_function(v_r_mag_input) result(fval)
                   allocate(v_shift, mold=fragments%vb)
                   v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%v_r_unit, fragments%v_t_mag, fragments%v_t_unit, fragments%mass, impactors%vbcom) 
                   do i = 1,fragments%nbody
+                     v_shift(:,i) = v_shift(:,i) - impactors%vbcom(:)
                      rotmag2 = fragments%rot(1,i)**2 + fragments%rot(2,i)**2 + fragments%rot(3,i)**2
                      vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2
                      kearr(i) = fragments%mass(i) * (fragments%Ip(3, i) * fragments%radius(i)**2 * rotmag2 + vmag2)