Fixed problems with angular momentum convergence. Residual rotation w…

…as not being recomputed correctly

examples/Fragmentation/Fragmentation_Movie.py

@@ -114,6 +114,16 @@
          "merge_spinner"             : 5.0e-3,
          }
 
+nfrag_reduction = {"disruption_headon" : 10.0,
+         "disruption_off_axis"         : 10.0,
+         "supercatastrophic_headon"    : 10.0,
+         "supercatastrophic_off_axis"  : 10.0,
+         "hitandrun_disrupt"           : 1.0,
+         "hitandrun_pure"              : 1.0,
+         "merge"                       : 1.0,
+         "merge_spinner"               : 1.0,
+         }
+
 density = 3000 * swiftest.AU2M**3 / swiftest.MSun
 GU = swiftest.GMSun * swiftest.YR2S**2 / swiftest.AU2M**3
 body_radius = body_Gmass.copy()
@@ -338,7 +348,7 @@ def vec_props(self, c):
         # Set fragmentation parameters
         minimum_fragment_gmass = 0.01 * body_Gmass[style][1] 
         gmtiny = 0.10 * body_Gmass[style][1] 
-        sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=10.0)
+        sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style])
         sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0)
 
         print("Generating animation")

src/fraggle/fraggle_generate.f90

@@ -533,7 +533,7 @@ 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
-      real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-2_DP    !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount)
+      real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-3_DP    !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount)
       real(DP)  :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy)
       integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best, posloop
       logical :: lhitandrun, lsupercat
@@ -542,13 +542,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       integer(I4B), dimension(:), allocatable :: vsign
       real(DP), dimension(:), allocatable :: vscale
       real(DP), dimension(:), allocatable :: dLi_mag
-      real(DP), parameter :: L_ROT_VEL_RATIO = 0.5_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments
       ! 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), parameter     :: hitandrun_vscale = 0.25_DP 
       real(DP)                :: vmin_guess 
       real(DP)                :: vmax_guess 
-      integer(I4B), parameter :: MAXLOOP = 20
-      integer(I4B), parameter :: MAXTRY  = 10
+      integer(I4B), parameter :: MAXLOOP = 100
+      integer(I4B), parameter :: MAXTRY  = 100
       integer(I4B), parameter :: MAXANGMTM = 1000
       class(collision_fraggle), allocatable :: collider_local
       character(len=STRMAX) :: message
@@ -653,12 +652,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
 
                   do i = istart,fragments%nbody
                      dL(:) = -dL1_mag * dLi_mag(i) * L_residual_unit(:)
-                     drot(:) = L_ROT_VEL_RATIO * dL(:) / (fragments%mass(i) * fragments%Ip(3,i) * fragments%radius(i)**2)
+                     drot(:) = dL(:) / (fragments%mass(i) * fragments%Ip(3,i) * fragments%radius(i)**2)
                      rot_new(:) = fragments%rot(:,i) + drot(:)
                      if (.mag.rot_new(:) < collider_local%max_rot) then
                         fragments%rot(:,i) = rot_new(:)
                         fragments%rotmag(i) = .mag.fragments%rot(:,i)
-                     else ! We would break the spin barrier here. Put less into spin and more into velocity shear. 
+                     else ! We would break the spin barrier here. Add a random component of rotation that is less than what would break the limit. The rest will go in velocity shear
                         call random_number(drot)
                         call random_number(rn)
                         drot(:) = (rn * collider_local%max_rot - fragments%rotmag(i)) * 2 * (drot(:) - 0.5_DP)
@@ -669,7 +668,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                            fragments%rot(:,i) = fragments%rotmag(i) * .unit. fragments%rot(:,i)
                         end if
                      end if
-                     L_residual(:) = L_residual(:) - drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 
+                     L_residual(:) = L_residual(:) + drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 
                   end do
 
                   ! Put any remaining residual into velocity shear
@@ -682,7 +681,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                      if (all(dL_metric(:) <= 1.0_DP)) exit angmtm
 
                      do i = istart, fragments%nbody
-                        dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody))
+                        dL(:) = -L_residual(:) / (fragments%nbody - istart + 1) 
                         call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i))
                         call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)  
                         fragments%vmag(i) = .mag.fragments%vc(:,i)