Fixed more bugs and trying new schemes for convergence

src/collision/collision_util.f90

@@ -248,11 +248,11 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par
                      fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * fragments%rot(:,i)
                   end do
                   call swiftest_util_get_potential_energy(nfrag, [(.true., i = 1, nfrag)], constraint_system%cb%Gmass, fragments%Gmass, fragments%mass, fragments%rb, fragments%pe)
-                  fragments%be = sum(-3*fragments%Gmass(:)*fragments%mass(:)/(5*fragments%radius(:)))
-                  fragments%L_orbit_tot(:) = sum(fragments%L_orbit(:,:),dim=2)
-                  fragments%L_spin_tot(:) = sum(fragments%L_spin(:,:),dim=2)
-                  fragments%ke_orbit_tot = sum(fragments%ke_orbit(:))
-                  fragments%ke_spin_tot = sum(fragments%ke_spin(:))
+                  fragments%be = sum(-3*fragments%Gmass(1:nfrag)*fragments%mass(1:nfrag)/(5*fragments%radius(1:nfrag)))
+                  fragments%L_orbit_tot(:) = sum(fragments%L_orbit(:,1:nfrag),dim=2)
+                  fragments%L_spin_tot(:) = sum(fragments%L_spin(:,1:nfrag),dim=2)
+                  fragments%ke_orbit_tot = sum(fragments%ke_orbit(1:nfrag))
+                  fragments%ke_spin_tot = sum(fragments%ke_spin(1:nfrag))
                end if
             end select
 
@@ -755,7 +755,7 @@ module subroutine collision_util_shift_vector_to_origin(m_frag, vec_frag)
 
       mvec_frag(:) = 0.0_DP
       mtot = sum(m_frag)
-      nfrag = size(m_frag)
+      nfrag = count(m_frag > tiny(0.0_DP))
 
       do i = 1, nfrag
          mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i)

src/encounter/encounter_io.f90

@@ -40,10 +40,10 @@ module subroutine encounter_io_netcdf_dump(self, param)
                end if
             end do
 
+            nc%max_tslot = nc%max_tslot + maxval(self%tmap(1:self%iframe))
             call nc%close()
-            call self%reset()
             ! Update the time slot tracker
-            nc%max_tslot = nc%max_tslot + maxval(self%tmap(1:self%iframe))
+            call self%reset()
          end if
       end select
 

src/fraggle/fraggle_generate.f90

@@ -457,9 +457,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       real(DP)                :: vmin_guess = 1.5_DP 
       real(DP)                :: vmax_guess = 10.0_DP
       real(DP)                :: delta_v, volume
-      integer(I4B), parameter :: MAXLOOP = 20
-      integer(I4B), parameter :: MAXTRY  = 100
-      real(DP), parameter :: SUCCESS_METRIC = 1.0e-3_DP
+      integer(I4B), parameter :: MAXLOOP = 100
+      real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP
       class(collision_fraggle), allocatable :: collider_local
       character(len=STRMAX) :: message
 
@@ -497,7 +496,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
             lfailure = .false.
             dE_metric = huge(1.0_DP)
 
-            outer: do try = 1, nfrag - 1
+            outer: do try = 1, nfrag - 2
                ! 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)
@@ -510,7 +509,6 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                fscale = log(vmax_guess - vmin_guess + 1.0_DP) / log(maxval(vscale(:)))
                vscale(:) = vscale(:)**fscale + vmin_guess - 1.0_DP
 
-
                ! Set the velocities of all fragments using all of the scale factors determined above
                do concurrent(i = 1:nfrag)
                   j = fragments%origin_body(i)
@@ -526,7 +524,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                      ! Add more velocity dispersion to disruptions vs hit and runs.
                      vmag = vesc * vscale(i) 
                      rimp(:) = fragments%rc(:,i) - impactors%rbimp(:)
-                     vimp_unit(:) = .unit. rimp(:)
+                     vimp_unit(:) = .unit. (rimp(:) + 0.5_DP * impactors%bounce_unit(:))
                      fragments%vc(:,i) = vmag * vscale(i) * vimp_unit(:) + vrot(:) 
                   end if
                end do
@@ -541,6 +539,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1)
                   fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) 
                   fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) 
+                  fragments%rotmag(:) = .mag.fragments%rot(:,:)
 
                   call collider_local%get_energy_and_momentum(nbody_system, param, phase="after")
                   L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) 
@@ -573,7 +572,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   fragments%vc(:,:) = fscale * fragments%vc(:,:)
 
                   f_spin = 1.0_DP - f_orbit
-                  ke_remove = min(f_spin * E_residual, fragments%ke_spin_tot)
+                  ke_remove = min(f_spin * E_residual, 0.01_DP*fragments%ke_spin_tot)
                   ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot)
                   where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) 
                   do concurrent(i = 1:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP)))
@@ -600,6 +599,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                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
 
                call fragments%reset()