From c0afe56121f590a2f0fb10300a80a9e2149534b1 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 13 Feb 2023 10:53:55 -0500 Subject: [PATCH 1/3] Put in minimum number of fragment limit --- src/fraggle/fraggle_util.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 98ebe0654..2f498a0cf 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -94,7 +94,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) real(DP), parameter :: BETA = 2.85_DP integer(I4B), parameter :: MASS_NOISE_FACTOR = 5 !! The number of digits of random noise that get added to the minimum mass value to prevent identical masses from being generated in a single run integer(I4B), parameter :: NFRAGMAX_UNSCALED = 3000 !! Maximum number of fragments that can be generated - integer(I4B), parameter :: NFRAGMIN = 1 !! Minimum number of fragments that can be generated (set by the fraggle_generate algorithm for constraining momentum and energy) + integer(I4B), parameter :: NFRAGMIN = 3 !! Minimum number of fragments that can be generated integer(I4B), parameter :: iMlr = 1 integer(I4B), parameter :: iMslr = 2 integer(I4B), parameter :: iMrem = 3 @@ -136,7 +136,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) mfrag = max((nfrag - 1)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) mremaining = mremaining - mfrag end do - nfrag = min(ceiling(nfrag / param%nfrag_reduction), nfragmax) + nfrag = max(min(ceiling(nfrag / param%nfrag_reduction), nfragmax), NFRAGMIN) call self%setup_fragments(nfrag) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) From 506bbb85657438e841eecda5260a92e19d6ec0b2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 13 Feb 2023 11:08:24 -0500 Subject: [PATCH 2/3] Improvements to hit and run case --- examples/Fragmentation/Fragmentation_Movie.py | 2 +- src/fraggle/fraggle_generate.f90 | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 5afb9484d..4d8c45d70 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -337,7 +337,7 @@ def vec_props(self, c): # Set fragmentation parameters minimum_fragment_gmass = 0.01 * body_Gmass[style][1] - gmtiny = 0.10 * body_Gmass[style][1] + gmtiny = 0.010 * body_Gmass[style][1] sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=10.0) 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 75ab1521b..c696fd4c8 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -536,7 +536,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu 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), dimension(NDIM) :: vimp_unit, rimp, vrot, vdisp, 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(:), allocatable :: vsign real(DP), dimension(:), allocatable :: vscale @@ -581,7 +581,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu 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))) + vesc = sqrt(2 * sum(fragments%Gmass(istart:fragments%nbody)) / impactors%radius(2)) vmin_guess = .mag.impactors%vc(:,2) * (1.0_DP - hitandrun_vscale) vmax_guess = .mag.impactors%vc(:,2) * (1.0_DP + hitandrun_vscale) else @@ -615,7 +615,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu j = fragments%origin_body(i) vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) if (lhitandrun) then - fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) + vrot(:) + vdisp(:) = .unit.(fragments%rc(:,i) - impactors%rc(:,2)) * vesc + fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) + vrot(:) + vdisp(:) else vmag = vscale(i) rimp(:) = fragments%rc(:,i) - impactors%rcimp(:) From 704eb59d20574e0dfb9af28b5a72cdced64f6619 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 13 Feb 2023 11:17:11 -0500 Subject: [PATCH 3/3] Fixed bug that occurs when L_residual_best is not set --- examples/Fragmentation/Fragmentation_Movie.py | 2 +- src/fraggle/fraggle_generate.f90 | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 4d8c45d70..5afb9484d 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -337,7 +337,7 @@ def vec_props(self, c): # Set fragmentation parameters minimum_fragment_gmass = 0.01 * body_Gmass[style][1] - gmtiny = 0.010 * 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, nfrag_reduction=10.0) 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 c696fd4c8..62e182fbb 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -569,6 +569,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) E_residual_best = huge(1.0_DP) + L_residual_best(:) = 0.0_DP lfailure = .false. dE_metric = huge(1.0_DP) dE_best = huge(1.0_DP) @@ -639,6 +640,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) L_residual_unit(:) = .unit. L_residual(:) + if (nsteps == 1) L_residual_best(:) = L_residual(:) ! Use equipartition of spin kinetic energy to distribution spin angular momentum do concurrent(i = istart:fragments%nbody)