From a91a7cd0463f3635f92e1abbb608f47458f19156 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 11 Feb 2023 14:44:08 -0500 Subject: [PATCH] Improved what happens when Fraggle fails to find a solution. It now restructures by merging the two largest bodies and reducing the number of fragments by 1. It keeps doing this until it runs out of fragments, then converts to merge if it never found a solution. --- examples/Fragmentation/Fragmentation_Movie.py | 12 +- src/fraggle/fraggle_generate.f90 | 190 ++++++++---------- src/fraggle/fraggle_module.f90 | 9 + src/fraggle/fraggle_util.f90 | 64 ++++++ 4 files changed, 162 insertions(+), 113 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 0917abef4..8850f0621 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -59,9 +59,9 @@ } vel_vectors = {"disruption_headon" : [np.array([ 0.00, 6.280005, 0.0]), - np.array([ 0.00, 3.28, 0.0])], + np.array([ 0.00, 3.90, 0.0])], "disruption_off_axis" : [np.array([ 0.00, 6.280005, 0.0]), - np.array([ 0.05, 3.28, 0.0])], + np.array([ 0.05, 3.90, 0.0])], "supercatastrophic_headon": [np.array([ 0.00, 6.28, 0.0]), np.array([ 0.00, 4.28, 0.0])], "supercatastrophic_off_axis": [np.array([ 0.00, 6.28, 0.0]), @@ -140,13 +140,13 @@ def encounter_combiner(sim): enc = sim.encounters[keep_vars].load() # Remove any encounter data at the same time steps that appear in the data to prevent duplicates - t_not_duplicate = ~enc['time'].isin(data['time']) - enc = enc.where(t_not_duplicate,drop=True) + t_not_duplicate = ~data['time'].isin(enc['time']) + data = data.sel(time=t_not_duplicate) tgood=enc.time.where(~np.isnan(enc.time),drop=True) enc = enc.sel(time=tgood) # The following will combine the two datasets along the time dimension, sort the time dimension, and then fill in any time gaps with interpolation - ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time") + ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time", method="akima") # Rename the merged Target body so that their data can be combined tname=[n for n in ds['name'].data if names[0] in n] @@ -336,7 +336,7 @@ def vec_props(self, c): sim.add_body(name=names, Gmass=body_Gmass[style], radius=body_radius[style], rh=pos_vectors[style], vh=vel_vectors[style], rot=rot_vectors[style]) # Set fragmentation parameters - minimum_fragment_gmass = 0.05 * body_Gmass[style][1] + minimum_fragment_gmass = 0.01 * body_Gmass[style][1] gmtiny = 0.10 * body_Gmass[style][1] sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 0e5bfd874..00a7c5f3f 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -57,18 +57,9 @@ module subroutine fraggle_generate(self, nbody_system, param, t) 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. 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 merge.") - call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge - return - end if + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Treating this as a merge.") + call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge + return end if associate (fragments => self%fragments) @@ -517,23 +508,21 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals - real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-3_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount) + real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-5_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount) real(DP) :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy) integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best logical :: lhitandrun, lsupercat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new, dL_metric real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, rn, dL1_mag, dE_conv - integer(I4B), dimension(collider%fragments%nbody) :: vsign - real(DP), dimension(collider%fragments%nbody) :: vscale - real(DP), dimension(collider%fragments%nbody) :: dLi_mag + integer(I4B), dimension(:), allocatable :: vsign + real(DP), dimension(:), allocatable :: vscale + real(DP), dimension(:), allocatable :: dLi_mag real(DP), parameter :: L_ROT_VEL_RATIO = 0.5_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have real(DP), parameter :: hitandrun_vscale = 0.25_DP real(DP) :: vmin_guess real(DP) :: vmax_guess - real(DP) :: delta_v, GC - integer(I4B), parameter :: MAXINNER = 100 - integer(I4B), parameter :: MAXOUTER = 10 + integer(I4B), parameter :: MAXINNER = 20 integer(I4B), parameter :: MAXANGMTM = 1000 class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -543,64 +532,61 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) lsupercat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - GC = impactors%Gmass(1) / impactors%mass(1) - allocate(collider_local, source=collider) - associate(fragments => collider_local%fragments) - - ! The fragments will be divided into two "clouds" based on identified origin body. - ! These clouds will collectively travel like two impactors bouncing off of each other. - where(fragments%origin_body(:) == 1) - vsign(:) = -1 - elsewhere - vsign(:) = 1 - end where - - ! Hit and run collisions should only affect the runners, not the target - if (lhitandrun) then - istart = 2 - else - istart = 1 - end if + ! Hit and run collisions should only affect the runners, not the target + if (lhitandrun) then + istart = 2 + else + istart = 1 + end if - ! The minimum fragment velocity will be set by the escape velocity - vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) - vesc = sqrt(2 * sum(fragments%Gmass(istart:nfrag)) / sum(fragments%radius(istart:nfrag))) - if (lhitandrun) then - vmin_guess = .mag.impactors%vc(:,2) - vimp * hitandrun_vscale - vmax_guess = .mag.impactors%vc(:,2) + vimp * hitandrun_vscale - else - vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) - vmin_guess = 1.001_DP * vesc - vmax_guess = vimp - end if + ! The minimum fragment velocity will be set by the escape velocity + vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) + + E_residual_best = huge(1.0_DP) + lfailure = .false. + dE_metric = huge(1.0_DP) + dE_best = huge(1.0_DP) + nsteps_best = 0 + nsteps = 0 + outer: do try = 1, nfrag - istart - 1 + associate(fragments => collider_local%fragments) + if (allocated(vsign)) deallocate(vsign); allocate(vsign(fragments%nbody)) + if (allocated(vscale)) deallocate(vscale); allocate(vscale(fragments%nbody)) + if (allocated(dLi_mag)) deallocate(dLi_mag); allocate(dLi_mag(fragments%nbody)) + + if (lhitandrun) then + vesc = sqrt(2 * sum(fragments%Gmass(istart:fragments%nbody)) / sum(fragments%radius(istart:fragments%nbody))) + vmin_guess = .mag.impactors%vc(:,2) - vimp * hitandrun_vscale + vmax_guess = .mag.impactors%vc(:,2) + vimp * hitandrun_vscale + else + vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) + vmin_guess = 1.001_DP * vesc + vmax_guess = vimp + end if + + ! The fragments will be divided into two "clouds" based on identified origin body. + ! These clouds will collectively travel like two impactors bouncing off of each other. + where(fragments%origin_body(:) == 1) + vsign(:) = -1 + elsewhere + vsign(:) = 1 + end where - E_residual_best = huge(1.0_DP) - lfailure = .false. - dE_metric = huge(1.0_DP) - dE_best = huge(1.0_DP) - nsteps_best = 0 - nsteps = 0 - outer: do try = 1, MAXOUTER ! 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) + do concurrent(i = 1:fragments%nbody) rimp(:) = fragments%rc(:,i) - impactors%rcimp(:) vscale(i) = .mag. rimp(:) / sum(impactors%radius(1:2)) end do - ! Set the velocity scale factor to span from vmin/vesc to vmax/vesc - if (nfrag == 2) then - vscale(:) = 1.0_DP - else - vscale(:) = vscale(:)/minval(vscale(:)) - fscale = log(vmax_guess - vmin_guess + 1.0_DP) / log(maxval(vscale(:))) - vscale(:) = vscale(:)**fscale + vmin_guess - 1.0_DP - end if + vscale(:) = vscale(:) / minval(vscale(:)) + fscale = log(vmax_guess - vmin_guess + 1.0_DP) / log(maxval(vscale(:))) + vscale(:) = vscale(:)**fscale + vmin_guess - 1.0_DP ! Set the velocities of all fragments using all of the scale factors determined above if (istart > 1) fragments%vc(:,1) = impactors%vc(:,1) * impactors%mass(1) / fragments%mass(1) - do concurrent(i = istart:nfrag) + do concurrent(i = istart:fragments%nbody) j = fragments%origin_body(i) vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) if (lhitandrun) then @@ -621,7 +607,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu E_residual = huge(1.0_DP) inner: do loop = 1, MAXINNER nsteps = nsteps + 1 - mfrag = sum(fragments%mass(istart:nfrag)) + mfrag = sum(fragments%mass(istart:fragments%nbody)) ! Try to put residual angular momentum into the spin, but if this would go past the spin barrier, then put it into velocity shear instead call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") @@ -629,14 +615,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu L_residual_unit(:) = .unit. L_residual(:) ! Use equipartition of spin kinetic energy to distribution spin angular momentum - do concurrent(i = istart:nfrag) + do concurrent(i = istart:fragments%nbody) dLi_mag(i) = ((fragments%mass(i) / fragments%mass(istart)) * & (fragments%radius(i) / fragments%radius(istart))**2 * & (fragments%Ip(3,i) / fragments%Ip(3,istart)))**(1.5_DP) end do - dL1_mag = .mag.L_residual(:) / sum(dLi_mag(istart:nfrag)) + dL1_mag = .mag.L_residual(:) / sum(dLi_mag(istart:fragments%nbody)) - do i = istart,nfrag + do i = istart,fragments%nbody dL(:) = -dL1_mag * dLi_mag(i) * L_residual_unit(:) drot(:) = L_ROT_VEL_RATIO * dL(:) / (fragments%mass(i) * fragments%Ip(3,i) * fragments%radius(i)**2) rot_new(:) = fragments%rot(:,i) + drot(:) @@ -665,8 +651,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu if (all(dL_metric(:) <= 1.0_DP)) exit angmtm - do i = istart, nfrag - dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:nfrag)) + do i = istart, fragments%nbody + dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody)) call collision_util_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) @@ -697,14 +683,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu allocate(collider%fragments, source=fragments) dE_metric = abs(E_residual) / impactors%Qloss end if - if ((dE_best < 0.0_DP) .and. (dE_metric <= ENERGY_SUCCESS_METRIC * try)) exit outer ! As the tries increase, we relax the success metric. What was once a failure might become a success - if (dE_conv < ENERGY_SUCCESS_METRIC * try) exit inner + if ((dE_best < 0.0_DP) .and. (dE_metric <= ENERGY_SUCCESS_METRIC)) exit outer + if (dE_conv < ENERGY_SUCCESS_METRIC) exit inner 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/try,0.0_DP)**2 + 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) @@ -718,45 +704,35 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) call fragments%set_coordinate_system() end do inner + end associate - call fragments%reset() - call fraggle_generate_pos_vec(collider_local, nbody_system, param, lfailure) - if (lfailure) exit - call fraggle_generate_rot_vec(collider_local, nbody_system, param) - - ! Increase the spatial size factor to get a less dense cloud - collider_local%fail_scale = collider_local%fail_scale * 1.001_DP + call collider_local%restructure(nbody_system, param, lfailure) + if (lfailure) exit outer + end do outer + lfailure = lfailure .or. (dE_best > 0.0_DP) - ! Bring the minimum and maximum velocities closer together - delta_v = (vmax_guess - vmin_guess) / 16.0_DP - vmin_guess = vmin_guess + delta_v - vmax_guess = vmax_guess - delta_v - end do outer - lfailure = lfailure .or. (dE_best > 0.0_DP) + 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.") - 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) - 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 + 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) - write(message,*) "dE/Qloss = ",-dE_best / impactors%Qloss - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - write(message,*) nsteps_best - call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Best solution came after " // trim(adjustl(message)) // " steps.") + 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) + write(message,*) "dE/Qloss = ",-dE_best / impactors%Qloss + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + write(message,*) nsteps_best + call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Best solution came after " // trim(adjustl(message)) // " steps.") - end associate end associate return end subroutine fraggle_generate_vel_vec diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index fda3ebe1b..37905a188 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -22,6 +22,7 @@ module fraggle procedure :: generate => fraggle_generate !! A simple disruption models that does not constrain energy loss in collisions procedure :: hitandrun => fraggle_generate_hitandrun !! Generates either a pure hit and run, or one in which the runner is disrupted procedure :: set_mass_dist => fraggle_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type + procedure :: restructure => fraggle_util_restructure !! Restructures the fragment distribution after a failure to converge on a solution end type collision_fraggle interface @@ -73,6 +74,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu logical, intent(out) :: lfailure !! Did the velocity computation fail? end subroutine fraggle_generate_vel_vec + module subroutine fraggle_util_restructure(self, nbody_system, param, lfailure) + implicit none + class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + logical, intent(out) :: lfailure !! Did the computation fail? + end subroutine fraggle_util_restructure + module subroutine fraggle_util_set_mass_dist(self, param) implicit none class(collision_fraggle), intent(inout) :: self !! Fraggle collision object diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index e8647d589..18fc34bf2 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -11,6 +11,70 @@ use swiftest contains + module subroutine fraggle_util_restructure(self, nbody_system, param, lfailure) + !! author: David A. Minton + !! + !! Restructures the fragment distribution after a failure to converge on a solution. + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object + class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + logical, intent(out) :: lfailure !! Did the computation fail? + ! Internals + class(collision_fragments), allocatable :: new_fragments + integer(I4B) :: i,nnew, nold + real(DP) :: volume + + associate(old_fragments => self%fragments) + lfailure = .false. + ! Merge the second-largest fragment into the first and shuffle the masses up + nold = old_fragments%nbody + if (nold <= 3) then + lfailure = .true. + return + end if + nnew = nold - 1 + allocate(new_fragments, mold=old_fragments) + call new_fragments%setup(nnew) + + ! Merge the first and second bodies + volume = 4._DP / 3._DP * PI * sum(old_fragments%radius(1:2)**3) + new_fragments%mass(1) = sum(old_fragments%mass(1:2)) + new_fragments%Gmass(1) =sum(old_fragments%Gmass(1:2)) + new_fragments%density(1) = new_fragments%mass(1) / volume + new_fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD) + do concurrent(i = 1:NDIM) + new_fragments%Ip(i,1) = sum(old_fragments%mass(1:2) * old_fragments%Ip(i,1:2)) + end do + new_fragments%Ip(:,1) = new_fragments%Ip(:,1) / new_fragments%mass(1) + new_fragments%origin_body(1) = old_fragments%origin_body(1) + + ! Copy the n>2 old fragments to the n>1 new fragments + new_fragments%mass(2:nnew) = old_fragments%mass(3:nold) + new_fragments%Gmass(2:nnew) = old_fragments%Gmass(3:nold) + new_fragments%density(2:nnew) = old_fragments%density(3:nold) + new_fragments%radius(2:nnew) = old_fragments%radius(3:nold) + do concurrent(i = 1:NDIM) + new_fragments%Ip(i,2:nnew) = old_fragments%Ip(i,3:nold) + end do + new_fragments%origin_body(2:nnew) = old_fragments%origin_body(3:nold) + + new_fragments%mtot = old_fragments%mtot + end associate + + deallocate(self%fragments) + call move_alloc(new_fragments, self%fragments) + + call fraggle_generate_pos_vec(self, nbody_system, param, lfailure) + if (lfailure) return + call fraggle_generate_rot_vec(self, nbody_system, param) + + ! Increase the spatial size factor to get a less dense cloud + self%fail_scale = self%fail_scale * 1.001_DP + + end subroutine fraggle_util_restructure + module subroutine fraggle_util_set_mass_dist(self, param) !! author: David A. Minton !!