More refinements to Fraggle attempting to get better success rate

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="disruption", 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

@@ -231,7 +231,7 @@ module subroutine collision_generate_simple(self, nbody_system, param, t)
             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 impactors%set_coordinate_system()
             call self%disrupt(nbody_system, param, t)
             call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
 

src/collision/collision_util.f90

@@ -423,17 +423,11 @@ module subroutine collision_util_set_budgets(self)
       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)
+         fragments%L_budget(:) = self%Ltot(:,1)
+         fragments%ke_budget = self%Etot(1) - impactors%Qloss
 
       end associate
 

src/fraggle/fraggle_generate.f90

@@ -101,7 +101,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
 
             lfailure_local = (dEtot > 0.0_DP) 
             if (lfailure_local) then
-               write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL
+               write(message, *) "dEtot: ",dEtot, "dEtot/Qloss", dEtot / impactors%Qloss
                call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to energy gain: " // &
                                                         trim(adjustl(message)))
                !cycle
@@ -161,7 +161,7 @@ subroutine fraggle_generate_minimize(collider, lfailure)
       logical,                  intent(out)   :: lfailure  !! Logical flag indicating whether this step fails or succeeds
       ! Internals
       real(DP), parameter :: TOL_MIN = 1.0e-5_DP
-      real(DP), parameter :: TOL_INIT = 1e-6_DP
+      real(DP), parameter :: TOL_INIT = 1e-8_DP
       integer(I4B), parameter :: MAXLOOP = 50
       real(DP), dimension(collider%fragments%nbody) :: input_v
       real(DP), dimension(:), allocatable :: output_v
@@ -181,13 +181,13 @@ subroutine fraggle_generate_minimize(collider, lfailure)
 
             fragments%v_r_unit(:,:) = .unit. fragments%vc(:,:)
             fragments%vmag(:) = .mag. fragments%vc(:,1:nfrag) 
-            fragments%rot(:,1:nfrag) = fragments%rot(:,1:nfrag) * 1e-12_DP
             do loop = 1, 3 !while(tol < TOL_MIN)
 
                input_v(:) = fragments%vmag(1:nfrag)
                fval = E_objective_function(input_v)
                call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v)
                fval = E_objective_function(output_v)
+               lfailure = lfailure .and. (fval > tol)
                input_v(:) = output_v(:)
 
                fragments%vmag(1:nfrag) = output_v(1:nfrag)
@@ -242,7 +242,7 @@ function E_objective_function(val_input) result(fval)
 
                   ! Get the current kinetic energy of the system
                   call tmp_frag%get_kinetic_energy()
-                  deltaE = tmp_frag%ke_budget - (tmp_frag%ke_orbit + tmp_frag%ke_spin)
+                  deltaE = (tmp_frag%ke_budget - (tmp_frag%ke_orbit + tmp_frag%ke_spin)) / (tmp_frag%ke_budget)
 
                   ! Use the deltaE as the basis of our objective function, with a higher penalty for having excess kinetic energy compared with having a deficit
                   if (deltaE < 0.0_DP) then

src/misc/minimizer_module.f90

@@ -128,7 +128,7 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1)
          real(DP), dimension(:), intent(out), allocatable :: x1
          ! Internals
          integer(I4B) ::  i, j, k, l, conv
-         real(DP), parameter     :: graddelta = 1e-2_DP !! Delta x for gradient calculations
+         real(DP), parameter     :: graddelta = 1e-4_DP !! Delta x for gradient calculations
          real(DP), dimension(N) :: S               !! Direction vectors 
          real(DP), dimension(N,N) :: H             !! Approximated inverse Hessian matrix 
          real(DP), dimension(N) :: grad1           !! gradient of f