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 
Improved objective function and allow now for some of the velocities to be unknown but others fixed based on the initial guess.
  • Loading branch information
@daminton
daminton committed May 13, 2021
1 parent 1f87eaf commit 2d03140
Show file tree
Hide file tree
Showing 3 changed files with 47 additions and 54 deletions.
44 changes: 17 additions & 27 deletions src/modules/module_symba.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -130,18 +131,20 @@ 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
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_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
Expand All @@ -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

Expand All @@ -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
33 changes: 18 additions & 15 deletions src/symba/symba_frag_pos.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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

Expand Down
24 changes: 12 additions & 12 deletions src/util/util_minimize_bfgs.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -80,16 +80,16 @@ 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
return
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

Expand Down

0 comments on commit 2d03140

Please sign in to comment.