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 b276d7813..003a7918e 100644 Binary files a/python/.DS_Store and b/python/.DS_Store differ 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