Renabled most of the Fraggle success criteria. I think there is an is…

…sue with the units that still needs to be resolved, because the results look weird
MintonGroup · Dec 30, 2022 · 44f066d · 44f066d
1 parent 4eecf34
commit 44f066d
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 13 deletions.
diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90
@@ -12,7 +12,7 @@
    use symba
 
    real(DP), parameter     :: FRAGGLE_LTOL = 1e-4_DP !10 * epsilon(1.0_DP)
-   real(DP), parameter     :: FRAGGLE_ETOL = 1e-1_DP
+   real(DP), parameter     :: FRAGGLE_ETOL = 1e-12_DP
 
 contains
 
@@ -54,11 +54,7 @@ module subroutine fraggle_generate(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%set_budgets()
-            call impactors%set_coordinate_system()
             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 
@@ -106,7 +102,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
        ! Internals
       integer(I4B)                         :: try
       real(DP)                             :: dEtot, dLmag
-      integer(I4B), parameter              :: MAXTRY = 10
+      integer(I4B), parameter              :: MAXTRY = 100
       logical                              :: lk_plpl, lfailure_local
       logical, dimension(size(IEEE_ALL))   :: fpe_halting_modes, fpe_quiet_modes
       logical, dimension(size(IEEE_USUAL)) :: fpe_flag 
@@ -156,23 +152,22 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
             ! Use the disruption collision model to generate initial conditions
             ! Compute the "before" energy/momentum and the budgets
             call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)
+            call self%set_budgets()
             call fraggle_generate_pos_vec(self)
             call fraggle_generate_rot_vec(self)
             call fraggle_generate_vel_vec(self)
             call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
 
-            call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
             dEtot = self%Etot(2) - self%Etot(1)
             dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1))
 
-            lfailure_local = (dEtot > 0.0_DP) 
+            lfailure_local = (dEtot > FRAGGLE_ETOL) 
             if (lfailure_local) then
                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
+               cycle
                lfailure_local = .false.
-               exit
             end if
 
             lfailure_local = ((abs(dLmag) / (.mag.self%Ltot(:,1))) > FRAGGLE_LTOL) 
@@ -296,8 +291,6 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t)
    end subroutine fraggle_generate_hitandrun
 
 
-
-
    module subroutine fraggle_generate_pos_vec(collider)
       !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!
@@ -510,11 +503,15 @@ module subroutine fraggle_generate_vel_vec(collider)
             ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)
             ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape
             ke_avail(:) = fragments%ke_orbit_frag(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%vmag(:)
-            ke_per_dof = -ke_residual/(nfrag - istart + 1)
+            ke_per_dof = -ke_residual/(2 * (nfrag - istart + 1))
             do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof)
                fragments%vmag(i) = sqrt(2 * (fragments%ke_orbit_frag(i) - ke_per_dof)/fragments%mass(i))
                fragments%vc(:,i) = fragments%vmag(i) * fragments%v_unit(:,i)
             end do
+            do concurrent(i = istart:nfrag, fragments%ke_spin_frag(i) > ke_per_dof)
+               fragments%rotmag(i) = sqrt(2 * (fragments%ke_spin_frag(i) - ke_per_dof)/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)))
+               fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i)
+            end do
             call fragments%get_kinetic_energy()
             ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)
 

diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90
@@ -250,6 +250,7 @@ module subroutine fraggle_util_set_natural_scale_factors(self)
          impactors%rc(:,:)     = impactors%rc(:,:)     / collider%dscale
          impactors%vc(:,:)     = impactors%vc(:,:)     / collider%vscale
          impactors%mass(:)     = impactors%mass(:)     / collider%mscale
+         impactors%Gmass(:)    = impactors%Gmass(:)   / 0.5_DP * collider%vscale**2 * collider%dscale 
          impactors%Mcb         = impactors%Mcb         / collider%mscale
          impactors%radius(:)   = impactors%radius(:)   / collider%dscale
          impactors%Lspin(:,:)  = impactors%Lspin(:,:)  / collider%Lscale
@@ -293,6 +294,7 @@ module subroutine fraggle_util_set_original_scale_factors(self)
          impactors%bounce_unit(:) = impactors%bounce_unit(:) * collider%vscale
 
          impactors%mass      = impactors%mass      * collider%mscale
+         impactors%Gmass(:)  = impactors%Gmass(:)  * 0.5_DP * collider%vscale**2 * collider%dscale 
          impactors%Mcb       = impactors%Mcb       * collider%mscale
          impactors%mass_dist = impactors%mass_dist * collider%mscale
          impactors%radius    = impactors%radius    * collider%dscale