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

Commit

Permalink
symba_frag_pos_fragment_velocity using secant method
Browse files Browse the repository at this point in the history
  • Loading branch information
@cwishard
cwishard committed May 12, 2021
1 parent b67327f commit 8336612
Show file tree
Hide file tree
Showing 3 changed files with 85 additions and 17 deletions.
11 changes: 5 additions & 6 deletions Makefile.Defines
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,13 +69,12 @@ GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak
GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries

#FFLAGS = -init=snan,arrays -traceback -no-wrap-margin -O3 -g -CB -nogen-interfaces -no-pie -fp-speculation=safe $(SIMDVEC) $(PAR) #$(HEAPARR)
#FFLAGS = $(IDEBUG) $(HEAPARR)
FFLAGS = -init=snan,arrays -no-wrap-margin -fp-model strict -fp-model no-except -traceback -g $(OPTIMIZE) -O3 $(SIMDVEC) $(PAR) $(HEAPARR)
FORTRAN = ifort
#AR = xiar
#FFLAGS = $(IDEBUG)
#FFLAGS = -init=snan,arrays -no-wrap-margin -fp-model strict -fp-model no-except $(OPTIMIZE) -O3 $(SIMDVEC) $(PAR) #-g -traceback -CB
#FORTRAN = ifort

#FORTRAN = gfortran
#FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM)
FORTRAN = gfortran
FFLAGS = -ffree-line-length-none -O3 #$(GDEBUG) $(GMEM)
AR = ar

# DO NOT include in CFLAGS the "-c" option to compile object only
Expand Down
Binary file modified python/.DS_Store
View file
Binary file not shown.
91 changes: 80 additions & 11 deletions src/symba/symba_frag_pos.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -338,6 +338,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr
! 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
Expand All @@ -357,13 +358,15 @@ 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(:))
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)
end do
if (Lambda > 0.0_DP) then
lmerge = .false.
call symba_frag_pos_fragment_velocity(m_frag, v_r_unit, Lambda, v_r_mag)
call symba_frag_pos_fragment_velocity(m_frag, x_frag, v_r_unit, Lambda, v_r_mag, Tau)
do i = 1, nfrag
v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t(:, i)
end do
Expand All @@ -378,23 +381,89 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr
return
end subroutine symba_frag_pos_kinetic_energy

subroutine symba_frag_pos_fragment_velocity(m_frag, v_r_unit, Lambda, v_r_mag)
subroutine symba_frag_pos_fragment_velocity(m_frag, x_frag, v_r_unit, Lambda, v_r_mag, Tau)
implicit none
! Arguments
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(out) :: v_r_mag !! Radial velocity magnitude (the intial guess values for i=1:4, and the final values for i=5:nfrag)
real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses
real(DP), dimension(:,:), intent(in) :: x_frag !! Fragment position in the center of mass frame
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(out) :: v_r_mag !! Radial velocity magnitude (the intial guess values for i=1:4, and the final values for i=5:nfrag)
real(DP), dimension(NDIM) :: Tau !! Sum of the tangential momentum vector of all fragments
! Internals
integer(I4B) :: i, j, k, nfrag
real(DP) :: Gam !! Sum of the radial momentum vector of i>4 fragments
real(DP), dimension(NDIM) :: tau !! Sum of the tangential momentum vector of all fragments
integer(I4B) :: i, j, k, nfrag
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, v_vec_02 !! Equation vectors for initial value guesses 1 and 2
real(DP) :: f_1_01, f_1_02 !! Equation 1 for initial value guesses 1 and 2 : angular momentum in the x direction
real(DP) :: f_2_01, f_2_02 !! Equation 2 for initial value guesses 1 and 2 : angular momentum in the y direction
real(DP) :: f_3_01, f_3_02 !! Equation 3 for initial value guesses 1 and 2 : angular momentum in the z direction
real(DP) :: f_4_01, f_4_02 !! Equation 4 for initial value guesses 1 and 2 : kinetic energy
real(DP) :: KE_after !! Radial kinetic energy of the four fragments given their new velocities
integer(I4B), parameter :: MAXITER = 50

nfrag = size(m_frag(:))

Gam = 0.0_DP
Beta = 0.0_DP

do i = 4, nfrag
Gam = Gam + m_frag(i) * v_r_mag(i) * v_r_unit(:,i)
Beta = Beta + m_frag(i) * v_r_mag(i)**2
end do

f_1_01 = Gam(1) + Tau(1)
f_1_02 = Gam(1) + Tau(1)

f_2_01 = Gam(2) + Tau(2)
f_2_02 = Gam(2) + Tau(2)

f_3_01 = Gam(3) + Tau(3)
f_3_02 = Gam(3) + Tau(3)

f_4_01 = Beta - Lambda
f_4_02 = Beta - Lambda

! Our initial guess for the first 4 fragments and the values of will be based on an equipartition of KE with some random variation
! 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(:))
v_r_mag(:) = sqrt(2 * Lambda / nfrag / m_frag(:)) * (1.0_DP + 0.02_DP * (v_r_mag(:) - 0.5_DP))
call random_number(v_r_mag_01(:))
v_r_mag_01(:) = sqrt(2 * Lambda / nfrag / m_frag(:)) * (v_r_mag_01(:) + 0.5_DP)
call random_number(v_r_mag_02(:))
v_r_mag_02(:) = sqrt(2 * Lambda / nfrag / m_frag(:)) * (v_r_mag_02(:) + 0.5_DP)

do i = 1, MAXITER
do i = 1, 4
f_1_01 = f_1_01 + m_frag(i) * v_r_mag_01(i) * x_frag(1,i)
f_1_02 = f_1_02 + m_frag(i) * v_r_mag_02(i) * x_frag(1,i)

f_2_01 = f_2_01 + m_frag(i) * v_r_mag_01(i) * x_frag(2,i)
f_2_02 = f_2_02 + m_frag(i) * v_r_mag_02(i) * x_frag(2,i)

f_3_01 = f_3_01 + m_frag(i) * v_r_mag_01(i) * x_frag(3,i)
f_3_02 = f_3_02 + m_frag(i) * v_r_mag_02(i) * x_frag(3,i)

f_4_01 = f_4_01 + m_frag(i) * v_r_mag_01(i)**2
f_4_02 = f_4_02 + m_frag(i) * v_r_mag_02(i)**2
end do

f_vec_01 = (/f_1_01, f_2_01, f_3_01, f_4_01/)
f_vec_02 = (/f_1_02, f_2_02, f_3_02, f_4_02/)

KE_after = 0.0_DP

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))
KE_after = KE_after + 0.5_DP * m_frag(i) * v_r_mag(i)**2
end do

if (KE_after <= Lambda) then
exit
else
v_r_mag_01(:) = v_r_mag_02(:)
v_r_mag_02(:) = v_r_mag(:)
end if
end do

return

Expand Down

0 comments on commit 8336612

Please sign in to comment.