From 71fcd33c787ad412cf8a629a25f9ef19bb1e5866 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 13 Jan 2023 14:29:33 -0500 Subject: [PATCH] Found a set of criteria that can handle the weird edge case I've been throwing at it. It reduces the number of fragments if there is a failure. --- src/fraggle/fraggle_generate.f90 | 47 +++++++++++++++++++++----------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index c2de18f73..846a63b57 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -103,6 +103,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur real(DP), dimension(NDIM) :: dL character(len=STRMAX) :: message real(DP), parameter :: fail_scale_initial = 1.001_DP + integer(I4B) :: nfrag_start ! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we ! can simply fail the attempt and try again. So we need to turn off any floating point exception halting modes temporarily @@ -114,9 +115,10 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - associate(impactors => self%impactors, fragments => self%fragments, nfrag => self%fragments%nbody, pl => nbody_system%pl) + associate(impactors => self%impactors, pl => nbody_system%pl) - write(message,*) nfrag + nfrag_start = self%fragments%nbody + write(message,*) nfrag_start call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") if (param%lflatten_interactions) then @@ -135,6 +137,10 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call fraggle_generate_rot_vec(self, nbody_system, param) call fraggle_generate_vel_vec(self, nbody_system, param, lfailure_local) + if (self%fragments%nbody /= nfrag_start) then + write(message,*) self%fragments%nbody + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle found a solution with " // trim(adjustl(message)) // " fragments" ) + end if call self%get_energy_and_momentum(nbody_system, param, phase="after") dL = self%L_total(:,2)- self%L_total(:,1) @@ -407,6 +413,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) fragments%rot(:,i) = impactors%rot(:,2) fragments%vc(:,i) = impactors%vc(:,2) end do + fragments%rotmag(:) = .mag.fragments%rot(:,:) call collider%get_energy_and_momentum(nbody_system, param, phase="after") dKE = 0.5_DP * (collider%pe(2) - collider%pe(1) + collider%be(2) - collider%be(1) - impactors%Qloss) KE_final = max(KE_init + dKE,0.0_DP) @@ -416,9 +423,9 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) Lbefore(:) = mass_init * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1)) Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:)) - L_spin(:) = impactors%L_spin(:,1) + (Lbefore(:) - Lafter(:)) - fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) - fragments%rotmag(1) = .mag.fragments%rot(:,1) + L_spin(:) = impactors%L_spin(:,1) + random_scale_factor * (Lbefore(:) - Lafter(:)) + !fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) + !fragments%rotmag(1) = .mag.fragments%rot(:,1) ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random do i = 2,nfrag call random_number(rotdir) @@ -450,7 +457,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu integer(I4B) :: i, j, loop, try, istart, nfrag, nlast logical :: lhitandrun, lsupercat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual - real(DP) :: vmag, vesc, dE, E_residual, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric + real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have @@ -458,6 +465,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP) :: vmax_guess = 10.0_DP real(DP) :: delta_v, volume integer(I4B), parameter :: MAXLOOP = 100 + integer(I4B), parameter :: MAXTRY = 1000 real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -496,7 +504,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 - 2 + outer: do try = 1, maxtry ! 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) @@ -524,7 +532,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(:) + 0.5_DP * impactors%bounce_unit(:)) + vimp_unit(:) = .unit. (rimp(:) + vsign(i) * impactors%bounce_unit(:)) fragments%vc(:,i) = vmag * vscale(i) * vimp_unit(:) + vrot(:) end if end do @@ -532,15 +540,22 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Every time the collision-frame velocities are altered, we need to be sure to shift everything back to the center-of-mass frame call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) call fragments%set_coordinate_system() + ke_min = 0.5_DP * fragments%mtot * vesc**2 do loop = 1, MAXLOOP call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + ke_avail = max(fragments%ke_orbit_tot - ke_min, 0.0_DP) ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1) - do i = 1, fragments%nbody - fragments%L_spin(:,i) = fragments%L_spin(:,i) - L_residual(:) * fragments%mass(i) / fragments%mtot - fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) - end do + if (ke_avail < epsilon(1.0_DP)) then + do i = 1, fragments%nbody + fragments%L_spin(:,i) = fragments%L_spin(:,i) - L_residual(:) * fragments%mass(i) / fragments%mtot + fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) + end do + else + 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)) + end if fragments%rotmag(:) = .mag.fragments%rot(:,:) call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") @@ -562,8 +577,6 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if if ((dE_best < 0.0_DP) .and. (dE_metric <= SUCCESS_METRIC * try)) exit outer ! As the tries increase, we relax the success metric. What was once a failure might become a success - ke_avail = fragments%ke_orbit_tot - 0.5_DP * fragments%mtot * vesc**2 - ! Remove a constant amount of velocity from the bodies so we don't shift the center of mass and screw up the momentum f_spin = (fragments%ke_spin_tot )/ (fragments%ke_spin_tot + fragments%ke_orbit_tot) f_orbit = 1.0_DP - f_spin @@ -574,11 +587,12 @@ 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, 0.01_DP*fragments%ke_spin_tot) + ke_remove = min(f_spin * E_residual, 0.9_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))) fscale = sqrt((fragments%ke_spin(i) - ke_rot_remove(i))/fragments%ke_spin(i)) + fragments%rotmag(i) = fscale * fragments%rotmag(i) fragments%rot(:,i) = fscale * fragments%rot(:,i) end do @@ -588,13 +602,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end do ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters + if (fragments%nbody == 2) exit ! Reduce the number of fragments by one nlast = fragments%nbody fragments%Ip(:,1) = fragments%mass(1) * fragments%Ip(:,1) + fragments%mass(nlast) * fragments%Ip(:,nlast) fragments%mass(1) = fragments%mass(1) + fragments%mass(nlast) fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1) fragments%Gmass(1) = fragments%Gmass(1) + fragments%mass(nlast) - volume = 4.0_DP / 3.0_DP * PI * (fragments%radius(1) + fragments%radius(nlast)) + volume = 4.0_DP / 3.0_DP * PI * ((fragments%radius(1))**3 + (fragments%radius(nlast))**3) fragments%density(1) = fragments%mass(1) / volume fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD) fragments%Ip(:,nlast) = 0.0_DP