More adjustments to angular momentum conservation

src/collision/collision_resolve.f90

@@ -365,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)
@@ -519,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)
@@ -548,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
 
@@ -596,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)
@@ -604,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)
 
@@ -616,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

@@ -686,17 +686,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)