From 8336612ba9449e91bfdd5d2075b2499167fc6c14 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Wed, 12 May 2021 10:03:31 -0400 Subject: [PATCH] symba_frag_pos_fragment_velocity using secant method --- Makefile.Defines | 11 ++--- python/.DS_Store | Bin 8196 -> 8196 bytes src/symba/symba_frag_pos.f90 | 91 ++++++++++++++++++++++++++++++----- 3 files changed, 85 insertions(+), 17 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index bf1270235..6366989d7 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -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 diff --git a/python/.DS_Store b/python/.DS_Store index b276d7813dc28ba7949c14764160013cc55a275a..003a7918e1f539f333acc9f83d8cb191664f5032 100644 GIT binary patch literal 8196 zcmeHMU2GLa6h3GB!|nnyaDnpI_CjgITA+Wnm=?YL1FJv_E%XO!*?V`PtJ~drckjKF zQqu=RjH1y8G%6;(5CbGUm|%P`#s_13ASM$0eK0ZL&y(@ZXf%3eW_tz7eN~h0Br|8u zoO5RPocYdd&n*C8N5-rNSOWkmbxzr08t&4#ocDQ630F#pBH053V1ZAg4cQ3Rq0EXWA(`yo!9(?m#TB^6f(4Wa}fT1?BL z(VXf4VG|E&BBZmD3RBusRu34OVpw9JG$(nSs}oIxbXHPn4k*n5!zW`{p`ds=*~NuA zU_#Ppj3N+4V15L6?XFYJo}-%m*sb+@kK?71j+Z1kPpMYwTmBM_W(91LIe#*=xtg4G zY(J>0`x2F+;*!!uszz<8HrhY#<@y67$oL&b`Y`AEj9@$&&Idft$k z8yL20rG>F=TBLGh&~{BP*WqPc(--^nOM!UM7+XI*z4d|SmfCoIV?*6cZG5`EVMkqU zd`Erb%#5a%RoCzAIXW?Qa{AQF)8~XofZ^2w%Osy4-z?;><7YdizeS~U{(hDBS6>`E zG^h^`*(uNMw@+GQY?)8)_7K~bwzA~fV9Ikc6RtnTmZuCSm2r%~YIhv2><ykyn7%Bn3}8(P}h zyZ7$DzVx1DOkWXOIheJ5JLy;lbG8|b4;y~Ua!uPkK0+~c1N)e5r8V}pI^Y@;mN1Hy z#n!5S4a#nYL(@%X)R}UX2;* zfKQ&uEo$^yjrrQZlXRZabh}=!GE+O~7)i^KOHF#S${y1Ov(iDi_@LgZF2`DJDCJE| z&};j0`Jq@sWuK^z^Jnid{GfNza{qa&E>#`ol6KRiMYe1ABR!ACVt2-2EvUfEV!+ zUcnFWBm5MF%($vSKxST+e*wqpl&;vRx#pMd6J9K|s_N&tNV9fGGz@SGxe zp28VCjnCq9_&lD$v-m2$hOgrrcn;6wMG3e(0DEo+AYLj2z-)iobv*Yt$>+d(yXgYt zilonLc>Uj9`1k(>7@5fBC;|&70$AMD-_=3>7P=XEtsSR+lsa!*Z&p%aLW3yBiHdTZ mX#O9DG>?-iw@HL_R#IA__OJgCP~b1+==&di{{h)N6z;c4$jkty!$1-?5m#d*gakGaO#%To*@O@g6W#m)fz9mBK&G6X2|Kg9 zBr#MhJ>p;agH@JQ7UjucdBTGpylCOUDz&6a^x(x(4_>@_u<&*FYzh+g;6Y;*(^dVx z?tbr0cfa{&dfzM|1p0D%J0W$15D6F3PeDELOfy+7o|y@HP*pE z5P={9K?H&b1Q7@#@PCK^p4qHOCC+`J4P+33AOiPf1lao_f|X&+$0dRM)q#Si0Ep!{ z%nPi5oV4U^KnU_a3>(#3HUREpP>Ms9sMHQ zoIp%qAcF`55m*}mmb)vMvllRDKYM%rK4dxRl;xyQ&mO6RSetM0f~^uhse(J}`&>`W zTBhqoqu)cQth{1Fr6gn7ET2xyIE93#c{#UF%bsGr5zU)P`R!51u~ynzYTPsmNxC)e z*q&zEMwUC9x`8HB6Q->@g+3=|>n^`eLMZq|Nm{e8(B2y9>|AV(EVOrYL|Y?W(dgo$ zEY&u(9~zp@%*`)6w|Ma~-wXl&twO3qe?`7sYO9u;kBhB_=89Xrxw2Jdqk3fUsFI{x z^0fevQg=tnn~g1o$GTlj&-=C*yEJ-s(J-&aMex2VcV z436!YXH6qJRxovM=7i>^4O=(ua}#;fHB*)`YHJyzy6m(hbFZ3ORh`VG%=0j6REZkGis$ zYR_Ip9kaBQVNqQc^X%_X)QP-^7(M0>c6BT2WZKDOa4T*#wEDw`4y)=JHuIC3>kZEu zHm>SRfBHUL)d{BS)pfQesWH|p?k?YMV=O9qGM+!%$sRmyj*=<-iZx5-$R+YJd4;?~ z-X|ZCYveO>oqR>UCf|^6$&ch`@(cNu{6_vDf00{&1Of$Wpcd-jL3jk3VHdQ(e&~V& za1eT60FJ?NcnXGL5~d*qY0x19d6u!Pj58-%tu9GBVCMae(s z;1#k=K0*F{M{bgzkbA!)@BT#2RYEoATs<_v4%msjYk}R+29H4}bR+lT(2v|3gb{ce z5-Wl2X(>3U=e+0_J66#p~KT!Jn|N1X9SV0hhAOind1faY>(cg!cFZHTsxpoBW zG*(ua-I74=f`XZ_ah{GNuKZz$?TD|7`M4wyHCX=h4*|jX&whUf=RarPy(|6(QseFW diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 0774d3f74..6c2a0e2b3 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -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 @@ -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 @@ -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