Adjusted parameters after failure to reduce the mass in fragments whe…

…n there is too much residual energy.

src/collision/collision_module.f90

@@ -151,6 +151,7 @@ module collision
       integer(I4B)                            :: status              !! Status flag to pass to the collision list once the collision has been resolved
       integer(I4B)                            :: collision_id        !! ID number of this collision event
       integer(I4B)                            :: maxid_collision = 0 !! The current maximum collision id number
+      real(DP)                                :: min_mfrag           !! Minimum fragment mass
 
       ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1
       real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor

src/collision/collision_util.f90

@@ -921,6 +921,8 @@ module subroutine collision_util_set_natural_scale_factors(self)
          fragments%Gmass(:)  = fragments%Gmass(:)  / (collider%dscale**3/collider%tscale**2)
          fragments%radius(:) = fragments%radius(:) / collider%dscale
          impactors%Qloss     = impactors%Qloss     / collider%Escale
+
+         collider%min_mfrag = collider%min_mfrag / collider%mscale
       end associate
 
       return
@@ -986,6 +988,7 @@ module subroutine collision_util_set_original_scale_factors(self)
          collider%pe(:)        = collider%pe(:)        * collider%Escale
          collider%be(:)        = collider%be(:)        * collider%Escale
          collider%te(:)        = collider%te(:)        * collider%Escale
+         collider%min_mfrag    = collider%min_mfrag    * collider%mscale
 
          collider%mscale = 1.0_DP
          collider%dscale = 1.0_DP

src/fraggle/fraggle_generate.f90

@@ -468,18 +468,19 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       class(swiftest_parameters),   intent(inout) :: param        !! Current run configuration parameters 
       logical,                      intent(out)   :: lfailure     !! Did the velocity computation fail?
       ! Internals
-      integer(I4B) :: i, j, loop, try, istart, nfrag, nlast, nsteps
+      integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps
       logical :: lhitandrun, lsupercat
       real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual
-      real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric
+      real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric, dM
       integer(I4B), dimension(collider%fragments%nbody) :: vsign
-      real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove
+      real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove, volume
       ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have
-      real(DP)                :: vmin_guess = 1.1_DP 
-      real(DP)                :: vmax_guess = 10.0_DP
-      real(DP)                :: delta_v, volume
-      integer(I4B), parameter :: MAXLOOP = 200
-      integer(I4B), parameter :: MAXTRY = 1000
+      real(DP)                :: vmin_guess = 1.01_DP 
+      real(DP)                :: vmax_guess = 8.0_DP
+      real(DP)                :: delta_v, GC
+      integer(I4B), parameter :: MAXLOOP = 100
+      integer(I4B), parameter :: MAXTRY = 20
+      real(DP), parameter     :: mass_reduction_ratio = 0.1_DP ! Ratio of difference between first and second fragment mass to remove from the largest fragment in case of a failure
       real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP
       class(collision_fraggle), allocatable :: collider_local
       character(len=STRMAX) :: message
@@ -489,9 +490,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
          lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) 
          lsupercat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) 
 
+         GC = impactors%Gmass(1) / impactors%mass(1)
+
          allocate(collider_local, source=collider)
          associate(fragments => collider_local%fragments)
 
+            volume(:) = 4.0_DP / 3.0_DP * PI * (fragments%radius(:))**3
+            fragments%density(:) = fragments%mass(:) / volume(:)
+
             ! The fragments will be divided into two "clouds" based on identified origin body. 
             ! These clouds will collectively travel like two impactors bouncing off of each other. 
             where(fragments%origin_body(:) == 1)
@@ -517,8 +523,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
             E_residual_best = huge(1.0_DP)
             lfailure = .false.
             dE_metric = huge(1.0_DP)
+            dE_best = huge(1.0_DP)
 
-            outer: do try = 1, nfrag
+            outer: do try = 1, MAXTRY
                ! Scale the magnitude of the velocity by the distance from the impact point
                ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct  
                do concurrent(i = 1:nfrag)
@@ -582,7 +589,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   dE = collider_local%te(2) - collider_local%te(1) 
                   E_residual = dE + impactors%Qloss
 
-                  if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (E_residual_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity
+                  if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (dE_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity
                      E_residual_best = E_residual
                      dE_best = dE
 
@@ -620,24 +627,23 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   call fragments%set_coordinate_system()
 
                end do
-               if (dE_best < 0.0_DP) exit outer
                ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters
-               if (fragments%nbody == 2) exit outer
-               ! Reduce the number of fragments by one
-               nlast = fragments%nbody
-               fragments%Ip(:,1) = fragments%mass(1) * fragments%Ip(:,1) + fragments%mass(nlast) * fragments%Ip(:,nlast)
-               fragments%mass(1) = fragments%mass(1) + fragments%mass(nlast)
-               fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1)
-               fragments%Gmass(1) = fragments%Gmass(1) + fragments%mass(nlast)
-               volume = 4.0_DP / 3.0_DP * PI * ((fragments%radius(1))**3 + (fragments%radius(nlast))**3)
-               fragments%density(1) = fragments%mass(1) / volume
-               fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD)
-               fragments%Ip(:,nlast) = 0.0_DP
-               fragments%mass(nlast) = 0.0_DP
-               fragments%Gmass(nlast) = 0.0_DP
-               fragments%radius(nlast) = 0.0_DP
-               fragments%status(nlast) = INACTIVE
-               fragments%nbody = nlast - 1
+               ! Reduce the fragment masses and add it to the largest remenant and try again
+               if (any(fragments%mass(2:nfrag) > collider%min_mfrag)) then
+                  do i = 2, nfrag
+                     if (fragments%mass(i) > collider%min_mfrag) then
+                        dM = min(mass_reduction_ratio * fragments%mass(i), fragments%mass(i) - collider%min_mfrag)
+                        fragments%mass(i) = fragments%mass(i) - dM
+                        fragments%mass(1) = fragments%mass(1) + dM
+                     end if
+                  end do
+               else
+                  exit outer
+               end if
+               fragments%Gmass(:) = GC * fragments%mass(:)
+
+               volume(:) = fragments%mass(:) / fragments%density(:)
+               fragments%radius(:) = (3._DP * volume(:) / (4._DP * PI))**(THIRD)
 
                call fragments%reset()
                call fraggle_generate_pos_vec(collider_local)

src/fraggle/fraggle_util.f90

@@ -25,7 +25,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
       integer(I4B)              :: i, j, jproj, jtarg, nfrag, istart
       real(DP), dimension(2)    :: volume
       real(DP), dimension(NDIM) :: Ip_avg
-      real(DP) :: mfrag, mremaining, min_mfrag, mtot, mcumul, G
+      real(DP) :: mfrag, mremaining, mtot, mcumul, G
       real(DP), dimension(:), allocatable :: mass
       real(DP), parameter :: BETA = 2.85_DP
       integer(I4B), parameter :: NFRAGMAX = 100  !! Maximum number of fragments that can be generated
@@ -37,7 +37,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
       integer(I4B), dimension(:), allocatable :: ind
       logical :: flipper
 
-      associate(impactors => self%impactors)
+      associate(impactors => self%impactors, min_mfrag => self%min_mfrag)
          ! Get mass weighted mean of Ip and density
          volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3
          mtot = sum(impactors%mass(:))

src/swiftest/swiftest_util.f90

@@ -1869,21 +1869,21 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
                call plplenc_old%copy(nbody_system%plpl_encounter)
             end if
 
-            ! ! Re-build the encounter list
-            ! ! Be sure to get the level info if this is a SyMBA nbody_system
-            ! select type(nbody_system)
-            ! class is (symba_nbody_system)
-            ! select type(pl)
-            ! class is (symba_pl)
-            ! select type(tp)
-            ! class is (symba_tp)
-            !    lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) 
-            !    if (tp%nbody > 0) then
-            !       lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec)
-            !    end if
-            ! end select
-            ! end select
-            ! end select
+            ! Re-build the encounter list
+            ! Be sure to get the level info if this is a SyMBA nbody_system
+            select type(nbody_system)
+            class is (symba_nbody_system)
+            select type(pl)
+            class is (symba_pl)
+            select type(tp)
+            class is (symba_tp)
+               lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) 
+               if (tp%nbody > 0) then
+                  lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec)
+               end if
+            end select
+            end select
+            end select
 
             ! Re-index the encounter list as the index values may have changed
             if (nenc_old > 0) then