From 2d031401b079e59982a6c612666f5b455220bce0 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 13 May 2021 18:26:46 -0400 Subject: [PATCH] Improved objective function and allow now for some of the velocities to be unknown but others fixed based on the initial guess. --- src/modules/module_symba.f90 | 44 +++++++++++++-------------------- src/symba/symba_frag_pos.f90 | 33 ++++++++++++++----------- src/util/util_minimize_bfgs.f90 | 24 +++++++++--------- 3 files changed, 47 insertions(+), 54 deletions(-) diff --git a/src/modules/module_symba.f90 b/src/modules/module_symba.f90 index 26284b3a1..536ef4ae1 100644 --- a/src/modules/module_symba.f90 +++ b/src/modules/module_symba.f90 @@ -95,8 +95,8 @@ MODULE module_symba ! number of fragments real(DP), dimension(:,:), allocatable :: xb ! barycentric position real(DP), dimension(:,:), allocatable :: vb ! barycentric velocity - real(DP), dimension(:), allocatable :: mass ! mass - real(DP), dimension(:), allocatable :: radius ! radius + real(DP), dimension(:), allocatable :: mass ! mass + real(DP), dimension(:), allocatable :: radius ! radius real(DP), dimension(:,:), allocatable :: IP ! moment of intertia real(DP), dimension(:,:), allocatable :: rot ! rotation type(swiftest_particle_info), dimension(:), allocatable :: info @@ -105,7 +105,7 @@ MODULE module_symba type, public, extends(lambda_obj) :: symba_vel_lambda_obj procedure(abstract_objective_func), pointer, nopass :: ke_objective_func_ptr => null() - real(DP), dimension(:), allocatable :: m_frag + real(DP), dimension(:), allocatable :: m_frag, v_r_mag real(DP), dimension(:,:), allocatable :: v_r_unit real(DP), dimension(NDIM) :: L_lin_tan real(DP) :: T_rad @@ -117,10 +117,11 @@ MODULE module_symba end type symba_vel_lambda_obj abstract interface - function abstract_objective_func(v_r_mag, m_frag, v_r_unit, L_lin_tan, T_rad) result(fnorm) + function abstract_objective_func(v_r_mag_unknowns, v_r_mag_knowns, m_frag, v_r_unit, L_lin_tan, T_rad) result(fnorm) ! Template for the kinetic energy constraint function used for minimizing import DP - real(DP), dimension(:), intent(in) :: v_r_mag !! Radial velocity magnitude + real(DP), dimension(:), intent(in) :: v_r_mag_unknowns !! Unknown radial velocity magnitudes + real(DP), dimension(:), intent(in) :: v_r_mag_knowns !! Known Radial velocity magnitude real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial unit vectors real(DP), dimension(:), intent(in) :: L_lin_tan !! Tangential component of linear momentum @@ -130,11 +131,12 @@ function abstract_objective_func(v_r_mag, m_frag, v_r_unit, L_lin_tan, T_rad) re end interface contains - subroutine ke_objective_func_init(self, lambda, m_frag, v_r_unit, L_lin_tan, T_rad) + subroutine ke_objective_func_init(self, lambda, v_r_mag, m_frag, v_r_unit, L_lin_tan, T_rad) implicit none ! Arguments class(symba_vel_lambda_obj), intent(out) :: self - procedure(abstract_objective_func) :: lambda + procedure(abstract_objective_func) :: lambda !! The lambda function + real(DP), dimension(:), intent(in) :: v_r_mag !! Radial velocity magnitude real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial unit vectors real(DP), dimension(:), intent(in) :: L_lin_tan !! Tangential component of linear momentum @@ -142,6 +144,7 @@ subroutine ke_objective_func_init(self, lambda, m_frag, v_r_unit, L_lin_tan, T_r self%ke_objective_func_ptr => lambda allocate(self%m_frag, source=m_frag) + allocate(self%v_r_mag, source=v_r_mag) allocate(self%v_r_unit, source=v_r_unit) self%L_lin_tan(:) = L_lin_tan(:) self%T_rad = T_rad @@ -152,6 +155,7 @@ subroutine ke_objective_func_destroy(self) type(symba_vel_lambda_obj) :: self if (allocated(self%m_frag)) deallocate(self%m_frag) if (allocated(self%v_r_unit)) deallocate(self%v_r_unit) + if (allocated(self%v_r_mag)) deallocate(self%v_r_mag) if (associated(self%ke_objective_func_ptr)) nullify(self%ke_objective_func_ptr) end subroutine ke_objective_func_destroy @@ -162,29 +166,15 @@ function ke_objective_func_eval(self, x) result(fnorm) real(DP), dimension(:), intent(in) :: x ! Result real(DP) :: fnorm + ! Internals + integer(I4B) :: nfrag, nunknown if (associated(self%ke_objective_func_ptr)) then - fnorm = self%ke_objective_func_ptr(x, self%m_frag, self%v_r_unit, self%L_lin_tan, self%T_rad) + nfrag = size(self%v_r_mag) + nunknown = size(x) + fnorm = self%ke_objective_func_ptr(x, self%v_r_mag(nunknown + 1:nfrag), self%m_frag, self%v_r_unit, self%L_lin_tan, self%T_rad) else error stop "KE Objective function was not initialized." end if end function ke_objective_func_eval -END MODULE module_symba -!********************************************************************************************************************************** -! -! Author(s) : David E. Kaufmann -! -! Revision Control System (RCS) Information -! -! Source File : $RCSfile$ -! Full Path : $Source$ -! Revision : $Revision$ -! Date : $Date$ -! Programmer : $Author$ -! Locked By : $Locker$ -! State : $State$ -! -! Modification History: -! -! $Log$ -!********************************************************************************************************************************** +END MODULE module_symba \ No newline at end of file diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index ac3f5bd7c..6f44f21c4 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -346,9 +346,9 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr integer(I4B) :: i, nfrag, neval real(DP), dimension(:,:), allocatable :: v_r_unit, v_t real(DP), dimension(NDIM) :: v_t_unit, h_unit, L_orb_frag - real(DP), dimension(:), allocatable :: v_r_mag + real(DP), dimension(:), allocatable :: v_r_mag, v_r_mag_sol integer(I4B), parameter :: MAXITER = 500 - real(DP), parameter :: TOL = 1e-5_DP + real(DP), parameter :: TOL = 1e-12_DP real(DP) :: vmult type(symba_vel_lambda_obj) :: ke_objective_func @@ -375,14 +375,13 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr if (T_rad > 0.0_DP) then lmerge = .false. - vmult = 1.0_DP call random_number(v_r_mag(:)) - v_r_mag(:) = vmult * sqrt(2 * T_rad / nfrag / m_frag(:)) * (v_r_mag(:) + 0.1_DP) + v_r_mag(:) = sqrt(2 * T_rad / nfrag / m_frag(:)) * (v_r_mag(:) + 0.1_DP) ! Initialize the lambda function with all the parameters that stay constant during the minimization - call ke_objective_func%init(symba_frag_pos_ke_objective_function, m_frag, v_r_unit, L_lin_tan, T_rad) + call ke_objective_func%init(symba_frag_pos_ke_objective_function, v_r_mag, m_frag, v_r_unit, L_lin_tan, T_rad) ! Minimize error using the BFGS optimizer - neval = util_minimize_bfgs(ke_objective_func, nfrag, v_r_mag, TOL) + neval = util_minimize_bfgs(ke_objective_func, 4, v_r_mag(1:4), TOL) do i = 1, nfrag v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t(:, i) @@ -397,11 +396,12 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr return end subroutine symba_frag_pos_kinetic_energy - function symba_frag_pos_ke_objective_function(v_r_mag, m_frag, v_r_unit, L_lin_tan, T_rad) result(fnorm) + function symba_frag_pos_ke_objective_function(v_r_mag_unknowns, v_r_mag_knowns, m_frag, v_r_unit, L_lin_tan, T_rad) result(fnorm) ! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value implicit none ! Arguments - real(DP), dimension(:), intent(in) :: v_r_mag !! Radial velocity magnitude + real(DP), dimension(:), intent(in) :: v_r_mag_unknowns !! Radial velocity magnitude + real(DP), dimension(:), intent(in) :: v_r_mag_knowns !! Radial velocity magnitude real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses real(DP), dimension(:), intent(in) :: L_lin_tan !! Tangential component of linear momentum real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial unit vectors @@ -410,15 +410,18 @@ function symba_frag_pos_ke_objective_function(v_r_mag, m_frag, v_r_unit, L_lin_t real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy !! Minimizing this brings us closer to our objective ! Internals - real(DP), dimension(4) :: f_vec + real(DP), dimension(4) :: f_vec, f_vec_0 - f_vec(1) = sum(m_frag(:) * v_r_mag(:) * v_r_unit(1,:)) - f_vec(2) = sum(m_frag(:) * v_r_mag(:) * v_r_unit(2,:)) - f_vec(3) = sum(m_frag(:) * v_r_mag(:) * v_r_unit(3,:)) - f_vec(1:3) = f_vec(1:3) + L_lin_tan(1:3) - f_vec(4) = sum(m_frag(:) * v_r_mag(:)**2) - T_rad + f_vec_0(1:3) = L_lin_tan(1:3) + f_vec_0(4) = -T_rad + f_vec(1) = sum(m_frag(:) * [v_r_mag_unknowns(:), v_r_mag_knowns(:)] * v_r_unit(1,:)) + f_vec(2) = sum(m_frag(:) * [v_r_mag_unknowns(:), v_r_mag_knowns(:)] * v_r_unit(2,:)) + f_vec(3) = sum(m_frag(:) * [v_r_mag_unknowns(:), v_r_mag_knowns(:)] * v_r_unit(3,:)) + f_vec(4) = sum(m_frag(:) * [v_r_mag_unknowns(:), v_r_mag_knowns(:)] **2) - fnorm = norm2(f_vec(:)) + f_vec(:) = f_vec(:) + f_vec_0(:) + + fnorm = norm2(f_vec(:) / f_vec_0(:)) return diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index a98bffc1d..2e1d01304 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -26,14 +26,14 @@ function util_minimize_bfgs(f, N, x1, eps) result(fnum) integer(I4B) :: fnum ! Internals integer(I4B) :: i, j, k, l, conv, num - integer(I4B), parameter :: MAXLOOP = 2000 !! Maximum number of loops before method is determined to have failed - real(DP), dimension(N) :: S ! Direction vectors + integer(I4B), parameter :: MAXLOOP = 2000 !! Maximum number of loops before method is determined to have failed + real(DP), dimension(N) :: S !! Direction vectors real(DP), dimension(N) :: x0 - real(DP), dimension(N) :: Snorm ! normalized direction - real(DP), dimension(N,N) :: H ! Approximated inverse Hessian matrix - real(DP), dimension(N) :: grad ! gradient of f - real(DP), dimension(N) :: grad0 ! old value of gradient - real(DP) :: astar ! 1 - D minimized value + real(DP), dimension(N) :: Snorm !! normalized direction + real(DP), dimension(N,N) :: H !! Approximated inverse Hessian matrix + real(DP), dimension(N) :: grad !! gradient of f + real(DP), dimension(N) :: grad0 !! old value of gradient + real(DP) :: astar !! 1D minimized value real(DP), dimension(N) :: y, P real(DP), dimension(N,N) :: PP, PyH, HyP real(DP) :: yHy, Py @@ -45,11 +45,11 @@ function util_minimize_bfgs(f, N, x1, eps) result(fnum) do i = 1, MAXLOOP xmag = norm2(x1(:)) grad0(:) = grad(:) - fnum = fnum + gradf(f, N, x1, grad, eps) + fnum = fnum + gradf(f, N, x1(1:N), grad, eps) if (i > 1) then ! set up factors for H matrix update y(:) = grad(:) - grad0(:) - P(:) = x1(:) - x0(:) + P(:) = x1(1:N) - x0(:) Py = sum(P(:) * y(:)) yHy = 0._DP do k = 1, N @@ -80,7 +80,7 @@ function util_minimize_bfgs(f, N, x1, eps) result(fnum) if (conv == 0) return ! normalize gradient Snorm(:) = S(:) / norm2(S) - num = fnum + minimize1D(f, x1, Snorm, N, eps, astar) + num = fnum + minimize1D(f, x1(1:N), Snorm, N, eps, astar) if (num == fnum) then write(*,*) "Exiting BFGS" fnum = 0 @@ -88,8 +88,8 @@ function util_minimize_bfgs(f, N, x1, eps) result(fnum) end if fnum = num ! Get new x values - x0(:) = x1(:) - x1(:) = x1(:) + astar * Snorm(:) + x0(:) = x1(1:N) + x1(1:N) = x1(1:N) + astar * Snorm(:) end do return