diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index b4141185a..13d97b4e1 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -705,10 +705,6 @@ module subroutine collision_util_set_spins(self) end do end if - ! Test to see if we were successful - call self%get_angular_momentum() - Lresidual(:) = self%L_budget(:) - (self%Lorbit(:) + self%Lspin(:)) - return end subroutine collision_util_set_spins diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 5205335b5..4ed5efac1 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -160,14 +160,14 @@ 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-5_DP - real(DP), parameter :: TOL_INIT = 1e-8_DP - integer(I4B), parameter :: MAXLOOP = 50 + real(DP), parameter :: TOL_INIT = 1e-5_DP + integer(I4B), parameter :: MAXLOOP = 100 real(DP), dimension(collider%fragments%nbody) :: input_v real(DP), dimension(:), allocatable :: output_v type(lambda_obj) :: Efunc real(DP) :: tol, fval integer(I4B) :: loop,i, nelem + logical :: lfirst_Efunc associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) select type(fragments => collider%fragments) @@ -175,35 +175,24 @@ subroutine fraggle_generate_minimize(collider, lfailure) nelem = nfrag lfailure = .false. - ! Find the local kinetic energy minimum for the nbody_system that conserves linear and angular momentum + ! Find the local kinetic energy minimum for the nbody_system that conserves linear and angular momentum Efunc = lambda_obj(E_objective_function) tol = TOL_INIT fragments%v_r_unit(:,:) = .unit. fragments%vc(:,:) fragments%vmag(:) = .mag. fragments%vc(:,1:nfrag) - 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) - - do concurrent(i=1:nfrag) - fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i) - end do - - ! Set spins in order to conserve angular momentum - call fragments%set_spins() - - if (.not.lfailure) exit - tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution + input_v(:) = fragments%vmag(1:nfrag) + lfirst_Efunc = .true. + fval = E_objective_function(input_v) + + call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v) + fragments%vmag(1:nfrag) = output_v(1:nfrag) + do concurrent(i=1:nfrag) + fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i) end do - + ! Set spins in order to conserve angular momentum + call fragments%set_spins() call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) end select end associate @@ -226,7 +215,8 @@ function E_objective_function(val_input) result(fval) integer(I4B) :: i type(fraggle_fragments(:)), allocatable :: tmp_frag real(DP) :: deltaE - + real(DP), save :: deltaE_scale = 1.0_DP + associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) select type(fragments => collider%fragments) class is (fraggle_fragments(*)) @@ -243,13 +233,17 @@ 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)) / (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 - fval = deltaE**8 + + ! The result will be scaled so such that the initial deltaE is 1.0. This helps improve the success of the + ! minimizer, by keeping values from getting to large + if (lfirst_Efunc) then + deltaE_scale = deltaE + deltaE = 1.0_DP + lfirst_Efunc = .false. else - fval = deltaE**2 + deltaE = deltaE/deltaE_scale end if + fval = deltaE**2 deallocate(tmp_frag) end select diff --git a/src/misc/minimizer_module.f90 b/src/misc/minimizer_module.f90 index 8beab3879..a0d56ae3a 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-8_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