From e9461d4f9d11b46b4402c83be9313ff1a93ecddd Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 31 Jan 2023 12:10:07 -0500 Subject: [PATCH] Fixed up a bunch of issues with Fraggle that was preventing it from finding solutions when rotations were high. --- src/collision/collision_generate.f90 | 9 +----- src/collision/collision_module.f90 | 7 +++++ src/collision/collision_util.f90 | 23 ++++++++++++++ src/fraggle/fraggle_generate.f90 | 46 +++++++++++++++++++++++----- 4 files changed, 69 insertions(+), 16 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 7bcf7bdbe..4ef6ec215 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -70,14 +70,7 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) nimp = size(impactors%id(:)) do i = 1, nimp j = impactors%id(i) - rcom(:) = pl%rb(:,j) - impactors%rbcom(:) - vcom(:) = pl%vb(:,j) - impactors%vbcom(:) - rnorm(:) = .unit. rcom(:) - ! Do the reflection - vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:) - ! Shift the positions so that the collision doesn't immediately occur again - pl%rb(:,j) = pl%rb(:,j) + 0.5_DP * pl%radius(j) * rnorm(:) - pl%vb(:,j) = impactors%vbcom(:) + vcom(:) + call collision_util_bounce_one(pl%rb(:,j),pl%vb(:,j),impactors%rbcom(:),impactors%vbcom(:),pl%radius(j)) self%status = DISRUPTED pl%status(j) = ACTIVE pl%ldiscard(j) = .false. diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 21b606202..b15aade01 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -397,6 +397,13 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine collision_util_add_fragments_to_collider + module subroutine collision_util_bounce_one(r,v,rcom,vcom,radius) + implicit none + real(DP), dimension(:), intent(inout) :: r,v + real(DP), dimension(:), intent(in) :: rcom,vcom + real(DP), intent(in) :: radius + end subroutine collision_util_bounce_one + module subroutine collision_util_construct_constraint_system(collider, nbody_system, param, constraint_system, phase) implicit none class(collision_basic), intent(inout) :: collider !! Collision system object diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index a73238ef7..d1e2c99e5 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -60,6 +60,29 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p return end subroutine collision_util_add_fragments_to_collider + module subroutine collision_util_bounce_one(r,v,rcom,vcom,radius) + !! Author: David A. Minton + !! + !! Performs a "bounce" operation on a single body by reversing its velocity in a center of mass frame. + implicit none + ! Arguments + real(DP), dimension(:), intent(inout) :: r,v + real(DP), dimension(:), intent(in) :: rcom,vcom + real(DP), intent(in) :: radius + ! Internals + real(DP), dimension(NDIM) :: rrel, vrel, rnorm + + rrel(:) = r(:) - rcom(:) + vrel(:) = v(:) - vcom(:) + rnorm(:) = .unit. rrel(:) + ! Do the reflection + vrel(:) = vrel(:) - 2 * dot_product(vrel(:),rnorm(:)) * rnorm(:) + ! Shift the positions so that the collision doesn't immediately occur again + r(:) = r(:) + 0.5_DP * radius * rnorm(:) + v(:) = vcom(:) + vrel(:) + + return + end subroutine collision_util_bounce_one module subroutine collision_util_construct_constraint_system(collider, nbody_system, param, constraint_system, phase) !! Author: David A. Minton diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 1f90fbe15..ba3d6dd01 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -54,10 +54,34 @@ module subroutine fraggle_generate(self, nbody_system, param, t) end select call self%set_mass_dist(param) call self%disrupt(nbody_system, param, t, lfailure) + if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Treating this as a hit and run.") - call self%hitandrun(nbody_system, param, t) - return + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Simplifying the collisional model.") + + impactors%mass_dist(1) = impactors%mass(1) + impactors%mass_dist(2) = max(0.5_DP * impactors%mass(2), self%min_mfrag) + impactors%mass_dist(3) = impactors%mass(2) - impactors%mass_dist(2) + impactors%regime = COLLRESOLVE_REGIME_DISRUPTION + call self%set_mass_dist(param) + call self%disrupt(nbody_system, param, t, lfailure) + + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Treating this as a bounce.") + call collision_util_bounce_one(impactors%rb(:,1),impactors%vb(:,1),impactors%rbcom(:),impactors%vbcom(:),impactors%radius(1)) + call collision_util_bounce_one(impactors%rb(:,2),impactors%vb(:,2),impactors%rbcom(:),impactors%vbcom(:),0.0_DP) + call impactors%set_coordinate_system() + call self%setup_fragments(2) + associate (fragments => self%fragments) + fragments%mass(1:2) = impactors%mass(1:2) + fragments%Gmass(1:2) = impactors%Gmass(1:2) + fragments%radius(1:2) = impactors%radius(1:2) + fragments%rb(:,1:2) = impactors%rb(:,1:2) + fragments%vb(:,1:2) = impactors%vb(:,1:2) + fragments%Ip(:,1:2) = impactors%Ip(:,1:2) + fragments%rot(:,1:2) = impactors%rot(:,1:2) + end associate + + end if end if associate (fragments => self%fragments) @@ -572,11 +596,6 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu 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. - dL(:) = -L_residual(:) * fragments%mass(i) / fragments%mtot + drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 - call fraggle_generate_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) - if (i >= istart) then call random_number(drot) drot(:) = (0.5_DP * collider_local%max_rot - fragments%rotmag(i)) * 2 * (drot(:) - 0.5_DP) @@ -586,12 +605,21 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fragments%rotmag(i) = 0.5_DP * collider%max_rot fragments%rot(:,i) = fragments%rotmag(i) * .unit. fragments%rot(:,i) end if + else + drot(:) = 0.0_DP end if + dL(:) = -L_residual(:) * fragments%mass(i) / fragments%mtot + drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 + call fraggle_generate_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) + end if end do end do angmtm + call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) + call fragments%set_coordinate_system() call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") ke_avail = 0.0_DP do i = 1, fragments%nbody @@ -627,6 +655,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu fscale = sqrt((max(fragments%ke_orbit_tot - ke_remove, 0.0_DP))/fragments%ke_orbit_tot) fragments%vc(:,:) = fscale * fragments%vc(:,:) fragments%vmag(:) = .mag.fragments%vc(:,:) + fragments%rc(:,:) = 1.0_DP / fscale * fragments%rc(:,:) + fragments%rmag(:) = .mag.fragments%rc(:,:) ! Update the unit vectors and magnitudes for the fragments based on their new orbits and rotations call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)