From 8fd92935efac32b19f8b6b5bc7ffd015e09dfbb0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 8 Feb 2023 09:44:37 -0500 Subject: [PATCH] More adjustments to angular momentum conservation --- src/collision/collision_resolve.f90 | 49 +++++++++++++++++------------ src/collision/collision_util.f90 | 15 ++++----- src/fraggle/fraggle_generate.f90 | 15 +++++---- src/swiftest/swiftest_util.f90 | 1 + 4 files changed, 47 insertions(+), 33 deletions(-) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index dc3d58a8d..3b5402162 100644 --- a/src/collision/collision_resolve.f90 +++ b/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 diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 7adbe2ec5..2b468c6f7 100644 --- a/src/collision/collision_util.f90 +++ b/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 diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 0d6fa5078..e4763519d 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/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) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 74a795856..5e271605c 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/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)