From b7e645cf993056669298894d0ed5e18b055fbaf0 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 13 May 2021 17:34:10 -0400 Subject: [PATCH] Hooked up the BFGS optimizer to the fragment velocity code --- src/modules/module_symba.f90 | 37 +++++----- src/symba/symba_frag_pos.f90 | 138 ++++++++++++----------------------- 2 files changed, 64 insertions(+), 111 deletions(-) diff --git a/src/modules/module_symba.f90 b/src/modules/module_symba.f90 index 74ced777d..26284b3a1 100644 --- a/src/modules/module_symba.f90 +++ b/src/modules/module_symba.f90 @@ -105,10 +105,10 @@ 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 real(DP), dimension(:,:), allocatable :: v_r_unit - real(DP), dimension(NDIM) :: tau, Gam - real(DP) :: Beta, Lam + real(DP), dimension(NDIM) :: L_lin_tan + real(DP) :: T_rad contains generic :: init => ke_objective_func_init procedure :: eval => ke_objective_func_eval @@ -117,33 +117,34 @@ MODULE module_symba end type symba_vel_lambda_obj abstract interface - function abstract_objective_func(v_r_mag, m_frag, v_r_unit, tau, Gam, Beta, Lam) result(fnorm) + function abstract_objective_func(v_r_mag, 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, m_frag, tau, Gam - real(DP), dimension(:,:), intent(in) :: v_r_unit - real(DP), intent(in) :: Beta, Lam - real(DP) :: fnorm + 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 + real(DP), intent(in) :: T_rad !! Target radial kinetic energ + real(DP) :: fnorm !! The objective function result: norm of the vector composed of the tangential momentum and energy end function end interface contains - subroutine ke_objective_func_init(self, lambda, m_frag, v_r_unit, tau, Gam, Beta, Lam) + subroutine ke_objective_func_init(self, lambda, 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 - real(DP), dimension(:), intent(in) :: m_frag, tau, Gam - real(DP), dimension(:,:), intent(in) :: v_r_unit - real(DP), intent(in) :: Beta, Lam + procedure(abstract_objective_func) :: lambda + 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 + real(DP), intent(in) :: T_rad !! Target radial kinetic ener self%ke_objective_func_ptr => lambda allocate(self%m_frag, source=m_frag) allocate(self%v_r_unit, source=v_r_unit) - self%tau(:) = tau(:) - self%Gam(:) = Gam(:) - self%Beta = Beta - self%Lam = Lam + self%L_lin_tan(:) = L_lin_tan(:) + self%T_rad = T_rad end subroutine ke_objective_func_init subroutine ke_objective_func_destroy(self) @@ -163,7 +164,7 @@ function ke_objective_func_eval(self, x) result(fnorm) real(DP) :: fnorm if (associated(self%ke_objective_func_ptr)) then - fnorm = self%ke_objective_func_ptr(x, self%m_frag, self%v_r_unit, self%tau, self%Gam, self%Beta, self%Lam) + fnorm = self%ke_objective_func_ptr(x, 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 diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 396a0bf5d..cde64d629 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -324,7 +324,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! - !! Adjust the fragment velocities to set the fragment orbital kinent + !! Adjust the fragment velocities to set the fragment orbital kinetic energy. !! It will check that we don't end up with negative energy (bound system). If so, we'll set the fragment velocities to !! zero in the center of mass frame and indicate the the fragmentation should instead by a merger. !! It takes in the initial "guess" of velocities and solve for the a scaling factor applied to the radial component wrt the @@ -340,13 +340,17 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr logical, intent(out) :: lmerge ! Internals - real(DP) :: mtot !! Total mass of fragments - real(DP) :: Lambda !! Sum of the radial kinetic energy of all fragments - real(DP), dimension(NDIM) :: tau !! Sum of the tangential momentum vector of all fragments - integer(I4B) :: i, nfrag - 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) :: mtot !! Total mass of fragments + real(DP) :: T_rad !! Sum of the radial kinetic energy of all fragments + real(DP), dimension(NDIM) :: L_lin_tan !! Sum of the tangential momentum vector of all fragments + 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 + integer(I4B), parameter :: MAXITER = 500 + real(DP), parameter :: TOL = 1e-5_DP + real(DP) :: vmult + type(symba_vel_lambda_obj) :: ke_objective_func nfrag = size(m_frag) mtot = sum(m_frag) @@ -362,19 +366,27 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr v_t(:,i) = v_frag(:, i) - dot_product(v_frag(:,i), v_r_unit(:, i)) * v_r_unit(:, i) end do - tau = 0.0_DP - Lambda = ke_target - 0.5_DP * mtot * dot_product(vcom(:), vcom(:)) + L_lin_tan = 0.0_DP + T_rad = ke_target - 0.5_DP * mtot * dot_product(vcom(:), vcom(:)) do i = 1, nfrag - Lambda = Lambda - 0.5_DP * m_frag(i) * dot_product(v_t(:, i), v_t(:, i)) - tau(:) = tau(:) + m_frag(i) * v_t(:, i) + T_rad = T_rad - 0.5_DP * m_frag(i) * dot_product(v_t(:, i), v_t(:, i)) + L_lin_tan(:) = L_lin_tan(:) + m_frag(i) * v_t(:, i) end do - if (Lambda > 0.0_DP) then + if (T_rad > 0.0_DP) then lmerge = .false. - v_r_mag(:) = symba_frag_pos_fragment_velocity(nfrag, m_frag, v_r_unit, Lambda, tau) + + 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) + + ! 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) + ! Minimize error using the BFGS optimizer + neval = util_minimize_bfgs(ke_objective_func, nfrag, v_r_mag, TOL) + do i = 1, nfrag v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t(:, i) end do - call symba_frag_pos_com_adjust(vcom, m_frag, v_frag) ! Temporary until we can get the fragment velocity constraints fixed else ! No solution exists for this case, so we need to indicate that this should be a merge ! This may happen due to setting the tangential velocities too high when setting the angular momentum constraint @@ -385,93 +397,33 @@ 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_fragment_velocity(nfrag, m_frag, v_r_unit, Lambda, tau) result(v_r_mag) + function symba_frag_pos_ke_objective_function(v_r_mag, 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 - integer(I4B), intent(in) :: nfrag - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial velocity unit vector for each fragment - real(DP), intent(in) :: Lambda !! Sum of the radial kinetic energy of all fragments - real(DP), dimension(:), intent(in) :: tau !! Sum of the tangential momentum vector of all fragments + 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) :: L_lin_tan !! Tangential component of linear momentum + real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial unit vectors + real(DP), intent(in) :: T_rad !! Target radial kinetic energy ! Result - real(DP), dimension(nfrag) :: v_r_mag !! Radial velocity magnitudes that satisfy the kinetic energy and momntum constraint + 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 - integer(I4B) :: i, j - real(DP), dimension(NDIM) :: Gam !! Sum of the radial momentum vector of i>4 fragments - real(DP) :: Beta !! Sum of the radial kinetic energy of i>4 fragments - real(DP), dimension(4) :: v_r_mag_01, v_r_mag_02 !! Two initial value guesses for the radial velocity magnitude of the first four fragments - real(DP), dimension(4) :: f_vec_01, f_vec_02, f_vec_const !! Equation vectors for initial value guesses 1 and 2 - integer(I4B), parameter :: MAXITER = 500 - real(DP), parameter :: TOL = 1e-9_DP - real(DP) :: err_rel, vmult - - ! Initialize the fragment radial velocities with random values. The first 4 serve as guesses that get updated with the secant method solver. - ! We shift the random variate to the range 0.5, 1.5 to prevent any zero values for radial velocities - call random_number(v_r_mag(:)) - - vmult = 1.0_DP - do j = 1, MAXITER - v_r_mag(:) = vmult * sqrt(8 * Lambda / nfrag / m_frag(:)) * (v_r_mag(:) + 0.1_DP) - ! Set up the constant values (those invovling i>4 fragments) - Gam = 0.0_DP - Beta = 0.0_DP - do i = 5, nfrag - Gam = Gam + m_frag(i) * v_r_mag(i) * v_r_unit(:,i) - Beta = Beta + 0.5_DP * m_frag(i) * v_r_mag(i)**2 - end do - if (Beta - Lambda > 0.0_DP) then - vmult = 0.9_DP * vmult - else - exit - end if - end do - - f_vec_const(1:3) = Gam(:) + tau(:) - f_vec_const(4) = Beta - Lambda - - ! The secant method requires two guesses, so we will use small values to start it off - v_r_mag_01(:) = 0.0_DP - v_r_mag_02(:) = sum(v_r_mag(:)) / (nfrag - 4) - - f_vec_01(:) = f_vec_const(:) - - do j = 1, MAXITER - f_vec_02(:) = f_vec_const(:) - do i = 1, 4 - if (j == 1) then - f_vec_01(1:3) = f_vec_01(1:3) + m_frag(i) * v_r_mag_01(i) * v_r_unit(:, i) - f_vec_01(4) = f_vec_01(4) + 0.5_DP * m_frag(i) * v_r_mag_01(i)**2 - end if - f_vec_02(1:3) = f_vec_02(1:3) + m_frag(i) * v_r_mag_02(i) * v_r_unit(:,i) - f_vec_02(4) = f_vec_02(4) + 0.5_DP * m_frag(i) * v_r_mag_02(i)**2 - end do - err_rel = norm2((f_vec_02(:) - f_vec_01(:)) / f_vec_01(:)) - - write(*,*) 'Iteration: ',j - write(*,*) 'v_r_mag_01: ',v_r_mag_01(:) - write(*,*) 'v_r_mag_02: ',v_r_mag_02(:) - write(*,*) 'f_vec_01: ',f_vec_01(:) - write(*,*) 'f_vec_02: ',f_vec_02(:) - write(*,*) 'err_rel: ',err_rel - - do i = 1, 4 - v_r_mag(i) = v_r_mag_02(i) - f_vec_02(i) * (v_r_mag_02(i) - v_r_mag_01(i)) / (f_vec_02(i) - f_vec_01(i)) - end do - - if (err_rel < TOL) exit - - write(*,*) 'v_r_mag: ',v_r_mag(1:4) + integer(I4B), parameter :: Neqs = 4 + real(DP), dimension(:), allocatable :: f_vec - v_r_mag_01(:) = v_r_mag_02(:) - v_r_mag_02(:) = v_r_mag(1:4) - f_vec_01(:) = f_vec_02(:) + 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 - end do + fnorm = norm2(f_vec(:)) return - end function symba_frag_pos_fragment_velocity + end function symba_frag_pos_ke_objective_function subroutine symba_frag_pos_com_adjust(vec_com, m_frag, vec_frag) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton