diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index f216fc242..ebad58734 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -232,7 +232,7 @@ def data_stream(self, frame=0): # Set fragmentation parameters minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades - sim.set_parameter(collision_model="disruption", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) sim.run(dt=1e-3, tstop=1.0e-3, istep_out=1, dump_cadence=0) print("Generating animation") diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index f9ce06541..eca87ed83 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -231,7 +231,7 @@ module subroutine collision_generate_simple(self, nbody_system, param, t) call self%set_mass_dist(param) call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) call self%set_budgets() - call self%set_coordinate_system() + call impactors%set_coordinate_system() call self%disrupt(nbody_system, param, t) call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 8151bdfa8..b4141185a 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -423,17 +423,11 @@ module subroutine collision_util_set_budgets(self) implicit none ! Arguments class(collision_basic), intent(inout) :: self !! Fraggle collision system object - ! Internals - real(DP) :: dEtot - real(DP), dimension(NDIM) :: dL associate(impactors => self%impactors, fragments => self%fragments) - dEtot = self%Etot(1) - dL(:) = self%Ltot(:,1) - - fragments%L_budget(:) = dL(:) - fragments%ke_budget = (dEtot - impactors%Qloss) + fragments%L_budget(:) = self%Ltot(:,1) + fragments%ke_budget = self%Etot(1) - impactors%Qloss end associate diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 4ccc9d12c..5205335b5 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -101,7 +101,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur lfailure_local = (dEtot > 0.0_DP) if (lfailure_local) then - write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL + write(message, *) "dEtot: ",dEtot, "dEtot/Qloss", dEtot / impactors%Qloss call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to energy gain: " // & trim(adjustl(message))) !cycle @@ -161,7 +161,7 @@ subroutine fraggle_generate_minimize(collider, lfailure) logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds ! Internals real(DP), parameter :: TOL_MIN = 1.0e-5_DP - real(DP), parameter :: TOL_INIT = 1e-6_DP + real(DP), parameter :: TOL_INIT = 1e-8_DP integer(I4B), parameter :: MAXLOOP = 50 real(DP), dimension(collider%fragments%nbody) :: input_v real(DP), dimension(:), allocatable :: output_v @@ -181,13 +181,13 @@ subroutine fraggle_generate_minimize(collider, lfailure) fragments%v_r_unit(:,:) = .unit. fragments%vc(:,:) fragments%vmag(:) = .mag. fragments%vc(:,1:nfrag) - fragments%rot(:,1:nfrag) = fragments%rot(:,1:nfrag) * 1e-12_DP do loop = 1, 3 !while(tol < TOL_MIN) input_v(:) = fragments%vmag(1:nfrag) fval = E_objective_function(input_v) call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v) fval = E_objective_function(output_v) + lfailure = lfailure .and. (fval > tol) input_v(:) = output_v(:) fragments%vmag(1:nfrag) = output_v(1:nfrag) @@ -242,7 +242,7 @@ function E_objective_function(val_input) result(fval) ! Get the current kinetic energy of the system call tmp_frag%get_kinetic_energy() - deltaE = tmp_frag%ke_budget - (tmp_frag%ke_orbit + tmp_frag%ke_spin) + deltaE = (tmp_frag%ke_budget - (tmp_frag%ke_orbit + tmp_frag%ke_spin)) / (tmp_frag%ke_budget) ! Use the deltaE as the basis of our objective function, with a higher penalty for having excess kinetic energy compared with having a deficit if (deltaE < 0.0_DP) then diff --git a/src/misc/minimizer_module.f90 b/src/misc/minimizer_module.f90 index e17cee3e9..8beab3879 100644 --- a/src/misc/minimizer_module.f90 +++ b/src/misc/minimizer_module.f90 @@ -128,7 +128,7 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1) real(DP), dimension(:), intent(out), allocatable :: x1 ! Internals integer(I4B) :: i, j, k, l, conv - real(DP), parameter :: graddelta = 1e-2_DP !! Delta x for gradient calculations + real(DP), parameter :: graddelta = 1e-4_DP !! Delta x for gradient calculations real(DP), dimension(N) :: S !! Direction vectors real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix real(DP), dimension(N) :: grad1 !! gradient of f