From 7805c9365dac6b6ed892a62368e0a91e7e0f453b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 27 Dec 2022 18:02:22 -0500 Subject: [PATCH] Fraggle sort of works again. The constraints are not met as well, but it does *something* at least --- src/collision/collision_generate.f90 | 16 ++++++++++++---- src/collision/collision_util.f90 | 4 +++- src/fraggle/fraggle_generate.f90 | 23 +++++++++++++---------- src/misc/minimizer_module.f90 | 2 +- 4 files changed, 29 insertions(+), 16 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 978dfdcff..165ee42a8 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -328,6 +328,8 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) ! The fragment trajectory will be the barycentric trajectory fragments%rb(:,1) = impactors%rbcom(:) fragments%vb(:,1) = impactors%vbcom(:) + fragments%rc(:,1) = 0.0_DP + fragments%vc(:,1) = 0.0_DP ! Get the energy of the system after the collision call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) @@ -485,7 +487,9 @@ module subroutine collision_generate_simple_rot_vec(collider) class(collision_basic), intent(inout) :: collider !! Fraggle collision system object ! Internals real(DP), dimension(NDIM) :: Lresidual, Lbefore, Lafter, Lspin, rot - integer(I4B) :: i + real(DP) :: ke_residual, rotmag_old, rotmag_new, vmag_old, vmag_new + integer(I4B) :: i, loop + integer(I4B) :: MAXLOOP = 10 associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) @@ -537,7 +541,7 @@ module subroutine collision_generate_simple_vel_vec(collider) real(DP), dimension(2) :: vimp real(DP) :: vmag, vdisp integer(I4B), dimension(collider%fragments%nbody) :: vsign - real(DP), dimension(collider%fragments%nbody) :: vscale + real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) @@ -554,7 +558,6 @@ module subroutine collision_generate_simple_vel_vec(collider) vdisp = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) else vdisp = 2 * sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) - !vmag = abs(dot_product(impactors%vb(:,2) - impactors%vb(:,1), impactors%y_unit(:))) end if vimp(:) = .mag.impactors%vc(:,:) @@ -566,13 +569,18 @@ module subroutine collision_generate_simple_vel_vec(collider) end do vscale(:) = vscale(:)/maxval(vscale(:)) + ! Give the fragment velocities a random value that is somewhat scaled with fragment mass + call random_number(mass_vscale) + mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2 + mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) + mass_vscale(:) = sqrt(2*mass_vscale(:) / maxval(mass_vscale(:))) fragments%vc(:,1) = .mag.impactors%vc(:,1) * impactors%bounce_unit(:) do concurrent(i = 2:nfrag) j = fragments%origin_body(i) vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rb(:,j) + impactors%rbcom(:)) - vmag = .mag.impactors%vc(:,j) * vscale(i) + vmag = .mag.impactors%vc(:,j) * vscale(i) * mass_vscale(i) if (lhitandrun) then fragments%vc(:,i) = vmag * 0.5_DP * impactors%bounce_unit(:) * vsign(i) + vrot(:) diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index f182b6dfb..94dce4f6d 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -216,7 +216,7 @@ module subroutine collision_util_get_kinetic_energy(self) do i = 1, nfrag fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) - fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) + fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) ) end do fragments%ke_orbit = fragments%ke_orbit / 2 @@ -263,6 +263,8 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, ke_orbit = ke_orbit + impactors%mass(i) * dot_product(impactors%vc(:,i), impactors%vc(:,i)) ke_spin = ke_spin + impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * dot_product(impactors%rot(:,i), impactors%rot(:,i)) end do + ke_orbit = ke_orbit / 2 + ke_spin = ke_spin / 2 else call fragments%get_angular_momentum() Lorbit(:) = fragments%Lorbit(:) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 2738039e3..8a22b41e4 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -99,12 +99,12 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur dEtot = self%Etot(2) - self%Etot(1) dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1)) - lfailure_local = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) + lfailure_local = (dEtot > 0.0_DP) if (lfailure_local) then write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // & + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to energy gain: " // & trim(adjustl(message))) - cycle + !cycle lfailure_local = .false. exit end if @@ -160,8 +160,8 @@ subroutine fraggle_generate_minimize(collider, lfailure) class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds ! Internals - real(DP), parameter :: TOL_MIN = 1.0e-6_DP - real(DP), parameter :: TOL_INIT = 1e-8_DP + real(DP), parameter :: TOL_MIN = 1.0e-5_DP + real(DP), parameter :: TOL_INIT = 1e-6_DP integer(I4B), parameter :: MAXLOOP = 50 real(DP), dimension(collider%fragments%nbody) :: input_v real(DP), dimension(:), allocatable :: output_v @@ -179,11 +179,11 @@ subroutine fraggle_generate_minimize(collider, lfailure) Efunc = lambda_obj(E_objective_function) tol = TOL_INIT + 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 while(tol < TOL_MIN) - fragments%v_r_unit(:,:) = .unit. fragments%vc(:,:) - fragments%vmag(:) = .mag. fragments%vc(:,:) - input_v(:) = fragments%vmag(1:nfrag) fval = E_objective_function(input_v) call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v) @@ -196,13 +196,15 @@ subroutine fraggle_generate_minimize(collider, lfailure) fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i) end do - call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) ! Set spins in order to conserve angular momentum call fraggle_generate_set_spins(fragments) if (.not.lfailure) exit tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution end do + + + call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) end select end associate @@ -244,10 +246,11 @@ function E_objective_function(val_input) result(fval) ! 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 - fval = deltaE**4 + fval = deltaE**8 else fval = deltaE**2 end if + deallocate(tmp_frag) end select end associate diff --git a/src/misc/minimizer_module.f90 b/src/misc/minimizer_module.f90 index 8beab3879..e17cee3e9 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-4_DP !! Delta x for gradient calculations + real(DP), parameter :: graddelta = 1e-2_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