Merge remote-tracking branch 'origin/debug' into debug

src/collision/collision_resolve.f90

@@ -142,7 +142,6 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa
             lflag = .true.
 
             ! Shift the impactors so that they are not overlapping
-
             rlim = sum(impactors%radius(1:2))
             vrel = impactors%vb(:,2) - impactors%vb(:,1)
             rrel = impactors%rb(:,2) - impactors%rb(:,1)
@@ -366,6 +365,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
             plnew%mass(1:nfrag) = fragments%mass(1:nfrag)
             plnew%Gmass(1:nfrag) = param%GU * fragments%mass(1:nfrag)
             plnew%radius(1:nfrag) = fragments%radius(1:nfrag)
+            if (allocated(pl%a)) call plnew%xv2el(cb)
 
             if (param%lrotation) then
                plnew%Ip(:, 1:nfrag) = fragments%Ip(:, 1:nfrag)
@@ -520,13 +520,14 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
       real(DP),                   intent(in)    :: dt     !! Current simulation step size
       integer(I4B),               intent(in)    :: irec   !! Current recursion level
       ! Internals
-      real(DP) :: E_before, E_after, dLmag
-      real(DP), dimension(NDIM) :: L_orbit_before, L_orbit_after, L_spin_before, L_spin_after, L_before, L_after, dL_orbit, dL_spin, dL
+      real(DP) :: E_before, E_after, mnew
+      real(DP), dimension(NDIM) ::L_before, L_after, dL
       logical :: lplpl_collision
       character(len=STRMAX) :: timestr, idstr
       integer(I4B), dimension(2) :: idx_parent       !! Index of the two bodies considered the "parents" of the collision
       logical  :: lgoodcollision
-      integer(I4B) :: i, loop, ncollisions
+      integer(I4B) :: i, j, nnew, loop, ncollisions
+      integer(I4B), dimension(:), allocatable :: idnew
       integer(I4B), parameter :: MAXCASCADE = 1000
 
       select type (nbody_system)
@@ -549,8 +550,6 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
             if (param%lenergy) then
                call nbody_system%get_energy_and_momentum(param)
                E_before = nbody_system%te
-               L_orbit_before(:) = nbody_system%L_orbit(:)
-               L_spin_before(:) = nbody_system%L_spin(:)
                L_before(:) = nbody_system%L_total(:)
             end if
 
@@ -597,6 +596,10 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
                   call plpl_collision%setup(0_I8B)
 
                   if ((nbody_system%pl_adds%nbody == 0) .and. (nbody_system%pl_discards%nbody == 0)) exit
+                  if (allocated(idnew)) deallocate(idnew)
+                  nnew = nbody_system%pl_adds%nbody
+                  allocate(idnew, source=nbody_system%pl_adds%id)
+                  mnew = sum(nbody_system%pl_adds%mass(:))
 
                   ! Rearrange the arrays: Remove discarded bodies, add any new bodies, re-sort, and recompute all indices and encounter lists
                   call pl%rearray(nbody_system, param)
@@ -605,6 +608,26 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
                   call nbody_system%pl_discards%setup(0, param)
                   call nbody_system%pl_adds%setup(0, param)
 
+                  if (param%lenergy) then
+                     call nbody_system%get_energy_and_momentum(param)
+                     L_after(:) = nbody_system%L_total(:)
+                     dL = L_after(:) - L_before(:)
+
+                     ! Add some velocity torque to the new bodies to remove residual angular momentum difference
+                     do j = 1, nnew
+                        i = findloc(pl%id,idnew(j),dim=1)
+                        if (i == 0) cycle
+                        call collision_util_velocity_torque(-dL * pl%mass(i)/mnew, pl%mass(i), pl%rb(:,i), pl%vb(:,i)) 
+                     end do
+
+                     call nbody_system%get_energy_and_momentum(param)
+                     E_after = nbody_system%te
+                     nbody_system%E_collisions = nbody_system%E_collisions + (E_after - E_before)
+                     L_after(:) = nbody_system%L_total(:)
+                     dL = L_after(:) - L_before(:)
+                  end if
+
+
                   ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plpl_collision 
                   call self%collision_check(nbody_system, param, t, dt, irec, lplpl_collision)
 
@@ -617,21 +640,6 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
                end associate
             end do
 
-            if (param%lenergy) then
-               call nbody_system%get_energy_and_momentum(param)
-               E_after = nbody_system%te
-               nbody_system%E_collisions = nbody_system%E_collisions + (E_after - E_before)
-
-               L_orbit_after(:) = nbody_system%L_orbit(:)
-               L_spin_after(:) = nbody_system%L_spin(:)
-               L_after(:) = nbody_system%L_total(:)
-
-               dL_orbit = L_orbit_after(:) - L_orbit_before(:)
-               dL_spin = L_spin_after(:) - L_spin_before(:)
-               dL = L_after(:) - L_before(:)
-               dLmag = .mag.dL
-            end if
-
          end associate
       end select
       end select

src/collision/collision_util.f90

@@ -116,6 +116,8 @@ module subroutine collision_util_construct_constraint_system(collider, nbody_sys
 
             ! Remove spins and velocities from all bodies other than the new fragments so that we can isolate the kinetic energy and momentum of the collision system, but still be able to compute
             ! the potential energy correctly
+            tmpsys%cb%Gmass = 0.0_DP
+            tmpsys%cb%mass = 0.0_DP
             tmpsys%cb%rot(:) = 0.0_DP
             tmpsys%pl%rot(:,:) = 0.0_DP
             tmpsys%pl%vb(:,:) = 0.0_DP
@@ -1039,17 +1041,16 @@ module subroutine collision_util_velocity_torque(dL, mass, r, v)
       real(DP), dimension(:), intent(in)    :: r    !! Position of body wrt system center of mass
       real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass
       ! Internals
-      real(DP), dimension(NDIM) :: dL_unit, r_unit, r_lever, vapply
-      real(DP) :: rmag, r_lever_mag
+      real(DP), dimension(NDIM) :: dL_unit, r_unit, vapply
+      real(DP) :: rmag, vmag, vapply_mag
 
       dL_unit(:) = .unit. dL
       r_unit(:) = .unit.r(:)
       rmag = .mag.r(:)
-      ! Project the position vector onto the plane defined by the angular momentum vector and the origin to get the "lever arm" distance
-      r_lever(:) = dL_unit(:) .cross. (r(:) .cross. dL_unit(:))
-      r_lever_mag = .mag.r_lever(:)
-      if ((r_lever_mag > epsilon(1.0_DP)) .and. (rmag > epsilon(1.0_DP))) then
-         vapply(:) = (dL(:) .cross. r(:)) / (mass * rmag**2)
+      vmag = .mag.v(:)
+      vapply_mag = .mag.(dL(:)/(rmag * mass))
+      if ((vapply_mag / vmag) > epsilon(1.0_DP)) then
+         vapply(:) = vapply_mag * (dL_unit(:) .cross. r_unit(:))
          v(:) = v(:) + vapply(:)
       end if
 

src/fraggle/fraggle_generate.f90

@@ -26,7 +26,7 @@ module subroutine fraggle_generate(self, nbody_system, param, t)
       real(DP),                 intent(in)    :: t            !! Time of collision
       ! Internals
       integer(I4B)          :: i, ibiggest, nfrag
-      real(DP), dimension(NDIM) :: L_residual, vbcom_orig, dvb
+      real(DP), dimension(NDIM) :: L_residual, vbcom_orig
       character(len=STRMAX) :: message 
       logical               :: lfailure
 
@@ -37,10 +37,6 @@ module subroutine fraggle_generate(self, nbody_system, param, t)
          associate(impactors => self%impactors, status => self%status, maxid => nbody_system%maxid)
             ! Set the coordinate system of the impactors
             call impactors%set_coordinate_system()
-
-
-            vbcom_orig(:) = impactors%vbcom(:)
-
             select case (impactors%regime) 
             case (COLLRESOLVE_REGIME_HIT_AND_RUN)
                call self%hitandrun(nbody_system, param, t)
@@ -78,30 +74,30 @@ module subroutine fraggle_generate(self, nbody_system, param, t)
                      fragments%radius(1:2) = impactors%radius(1:2)
                      fragments%rb(:,1:2) = impactors%rb(:,1:2)
                      fragments%vb(:,1:2) = impactors%vb(:,1:2)
-                     do concurrent(i = 1:2)
-                        fragments%rc(:,i) = fragments%rb(:,i) - impactors%rbcom(:)
-                        fragments%vc(:,i) = fragments%vb(:,i) - impactors%vbcom(:)
-                     end do
                      fragments%Ip(:,1:2) = impactors%Ip(:,1:2)
                      fragments%rot(:,1:2) = impactors%rot(:,1:2)
                      fragments%mtot = sum(fragments%mass(1:2))
                   end associate
                end if
             end if
 
-            ! Get the energy and momentum of the system before and after the collision
-            call self%get_energy_and_momentum(nbody_system, param, phase="before")
-            call self%get_energy_and_momentum(nbody_system, param, phase="after")
-            L_residual(:) = (self%L_total(:,2) - self%L_total(:,1))
 
             associate (fragments => self%fragments)
+               ! Get the energy and momentum of the system before and after the collision
+               call self%get_energy_and_momentum(nbody_system, param, phase="before")
                nfrag = fragments%nbody
+               do concurrent(i = 1:2)
+                  fragments%rc(:,i) = fragments%rb(:,i) - impactors%rbcom(:)
+                  fragments%vc(:,i) = fragments%vb(:,i) - impactors%vbcom(:)
+               end do
+               call self%get_energy_and_momentum(nbody_system, param, phase="after")
+               L_residual(:) = (self%L_total(:,2) - self%L_total(:,1))
 
                ! Put any residual angular momentum into orbital velocity
+               vbcom_orig = impactors%vbcom(:)
                call collision_util_velocity_torque(-L_residual(:), fragments%mtot, impactors%rbcom(:), impactors%vbcom(:))
-               dvb(:) = impactors%vbcom(:) - vbcom_orig(:)
-               do concurrent(i = 1:nfrag)
-                  fragments%vb(:,i) = fragments%vb(:,i) + dvb(:)
+               do i=1,nfrag
+                  fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:)
                end do
 
                select case(impactors%regime)
@@ -320,13 +316,15 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
       real(DP), dimension(collider%fragments%nbody) :: mass_rscale, phi, theta, u
       integer(I4B) :: i, j, loop, istart
       logical :: lsupercat, lhitandrun
-      integer(I4B), parameter :: MAXLOOP = 100
+      integer(I4B), parameter :: MAXLOOP = 200
       real(DP), parameter :: rdistance_scale_factor = 1.0_DP ! Scale factor to apply to distance scaling of cloud centers in the event of overlap
                                                              ! The distance is chosen to be close to the original locations of the impactors
                                                              ! but far enough apart to prevent a collisional cascade between fragments 
       real(DP), parameter :: cloud_size_scale_factor = 3.0_DP ! Scale factor to apply to the size of the cloud relative to the distance from the impact point. 
-                                                               ! A larger value puts more space between fragments initially
-
+                                                              ! A larger value puts more space between fragments initially
+      real(DP) :: rbuffer  ! Body radii are inflated by this scale factor to prevent secondary collisions 
+      real(DP), parameter :: rbuffer_max = 1.2
+      rbuffer = 1.05_DP 
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody, &
          pl => nbody_system%pl, tp => nbody_system%tp, npl => nbody_system%pl%nbody, ntp => nbody_system%tp%nbody)
@@ -358,8 +356,8 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
                fragment_cloud_center(:,1) = impactors%rc(:,1) 
                fragment_cloud_center(:,2) = impactors%rc(:,2) + rdistance * impactors%bounce_unit(:)
             else ! Keep the first and second bodies at approximately their original location, but so as not to be overlapping
-               fragment_cloud_center(:,1) = impactors%rc(:,1)
-               fragment_cloud_center(:,2) = impactors%rcimp(:) + 1.001_DP * max(fragments%radius(2),impactors%radius(2)) * impactors%y_unit(:)
+               fragment_cloud_center(:,1) = impactors%rcimp(:) - rbuffer * max(fragments%radius(1),impactors%radius(1)) * impactors%y_unit(:)
+               fragment_cloud_center(:,2) = impactors%rcimp(:) + rbuffer * max(fragments%radius(2),impactors%radius(2)) * impactors%y_unit(:)
                fragment_cloud_radius(:) = cloud_size_scale_factor * rdistance
             end if
             fragments%rc(:,1) = fragment_cloud_center(:,1)
@@ -379,7 +377,7 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
 
                ! Make a random cloud
                phi(i) = TWOPI * phi(i)
-               theta(i) = acos( 2 * theta(i) - 1.0_DP)
+               theta(i) = acos(2 * theta(i) - 1.0_DP)
                ! Scale the cloud size
                fragments%rmag(i) = fragment_cloud_radius(j) * mass_rscale(i) * u(i)**(THIRD)
 
@@ -418,17 +416,18 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
                ! Check for overlaps between fragments
                do i = j + 1, nfrag
                   dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i))
-                  loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) 
-                  loverlap(j) = loverlap(j) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) 
+                  loverlap(i) = loverlap(i) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) 
+                  loverlap(j) = loverlap(j) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) 
                end do
                ! Check for overlaps with existing bodies that are not involved in the collision 
                do i = 1, npl
                   if (any(impactors%id(:) == i)) cycle
                   dis = .mag. (fragments%rc(:,j) - (pl%rb(:,i) / collider%dscale - impactors%rbcom(:)))
-                  loverlap(j) = loverlap(j) .or. (dis <= (pl%radius(i) / collider%dscale + fragments%radius(j))) 
+                  loverlap(j) = loverlap(j) .or. (dis <= rbuffer * (pl%radius(i) / collider%dscale + fragments%radius(j))) 
                end do
             end do
             rdistance = rdistance * collider%fail_scale
+            rbuffer = min(rbuffer * collider%fail_scale, rbuffer_max)
          end do
 
          lfailure = any(loverlap(:))
@@ -631,11 +630,6 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   call fragments%set_coordinate_system()
                   call collider_local%get_energy_and_momentum(nbody_system, param, phase="after")
 
-                  ke_avail = 0.0_DP
-                  do i = fragments%nbody, 1, -1
-                     ke_avail = ke_avail + 0.5_DP * fragments%mass(i) * max(fragments%vmag(i) - vesc,0.0_DP)**2
-                  end do
-
                   dE = collider_local%te(2) - collider_local%te(1) 
                   E_residual_last = E_residual
                   E_residual = dE + impactors%Qloss
@@ -661,6 +655,11 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   end if
 
                   ! Remove a constant amount of velocity from the bodies so we don't shift the center of mass and screw up the momentum 
+                  ke_avail = 0.0_DP
+                  do i = fragments%nbody, 1, -1
+                     ke_avail = ke_avail + 0.5_DP * fragments%mass(i) * max(fragments%vmag(i) - vesc,0.0_DP)**2
+                  end do
+
                   ke_remove = min(E_residual, ke_avail)
                   fscale = sqrt((max(fragments%ke_orbit_tot - ke_remove, 0.0_DP))/fragments%ke_orbit_tot)
                   fragments%vc(:,:) = fscale * fragments%vc(:,:)
@@ -688,17 +687,20 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
             end do outer
             lfailure = (dE_best > 0.0_DP) 
 
-            call collision_util_velocity_torque(-L_residual_best(:), fragments%mtot, impactors%rbcom, impactors%vbcom)
-
-            do concurrent(i = 1:collider%fragments%nbody)
-               collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:)
-            end do
-
             write(message, *) nsteps
             if (lfailure) then
                call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle velocity calculation failed to converge after " // trim(adjustl(message)) // " steps. The best solution found had:")
             else 
                call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Fraggle velocity calculation converged after " // trim(adjustl(message)) // " steps.")
+
+               call collider%get_energy_and_momentum(nbody_system, param, phase="after")
+               L_residual(:) = (collider%L_total(:,2) - collider%L_total(:,1))
+               call collision_util_velocity_torque(-L_residual(:), collider%fragments%mtot, impactors%rbcom, impactors%vbcom)
+
+               do concurrent(i = 1:collider%fragments%nbody)
+                  collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:)
+               end do
+
             end if
             write(message,*) "dL/|L0|  = ",(L_residual_best(:))/.mag.collider_local%L_total(:,1) 
             call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)

src/swiftest/swiftest_util.f90

@@ -2661,6 +2661,7 @@ module subroutine swiftest_util_set_rhill(self,cb)
       class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
 
       if (self%nbody == 0) return
+      if (cb%Gmass <= tiny(1.0_DP)) return
 
       call self%xv2el(cb) 
       where(self%a(1:self%nbody) > 0.0_DP)