Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Browse files Browse the repository at this point in the history 
Hooked up the BFGS optimizer to the fragment velocity code
  • Loading branch information
@daminton
daminton committed May 13, 2021
1 parent ffdf0ad commit b7e645c
Show file tree
Hide file tree
Showing 2 changed files with 64 additions and 111 deletions.
37 changes: 19 additions & 18 deletions src/modules/module_symba.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down
138 changes: 45 additions & 93 deletions src/symba/symba_frag_pos.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit b7e645c

Please sign in to comment.