From 6f79fa573fe24618f00115ca24fb9f43bc38e022 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 12 May 2021 14:32:36 -0400 Subject: [PATCH] Applied some pre-conditioning checks on the initial guess values of everything, but it always leads to divergence. --- src/symba/symba_frag_pos.f90 | 37 ++++++++++++++++++++++-------------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 1dffa270c..fc4d50a08 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -362,7 +362,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr 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) + tau(:) = tau(:) + m_frag(i) * v_t(:, i) end do if (Lambda > 0.0_DP) then lmerge = .false. @@ -400,29 +400,38 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd 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 = 50 real(DP), parameter :: TOL = 1e-8_DP + real(DP) :: 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(:)) - v_r_mag(:) = sqrt(2 * Lambda / nfrag / m_frag(:)) * (v_r_mag(:) + 0.5_DP) - ! Our initial guess for the first 4 fragments and the values of will be based on an equipartition of KE with some random variation - v_r_mag_01(:) = v_r_mag(1:4) - ! The secant method requires two guesses, so we will use a small random variate to update the initial guesses - call random_number(v_r_mag_02(:)) - v_r_mag_02(:) = v_r_mag_01 * (1._DP + 2e-1_DP * (v_r_mag_02 - 0.5_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 + m_frag(i) * v_r_mag(i)**2 + vmult = 1.0_DP + do j = 1, MAXITER + v_r_mag(:) = vmult * sqrt(2 * 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 + m_frag(i) * v_r_mag(i)**2 + end do + if (Beta - Lambda > 0.0_DP) then + vmult = 0.8_DP * vmult + else + exit + end if end do f_vec_const(1:3) = Gam(:) + tau(:) f_vec_const(4) = 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 + v_r_mag_01(:) = 0.0_DP + ! The secant method requires two guesses, so we will use a small random variate to update the initial guesses + !call random_number(v_r_mag_02(:)) + v_r_mag_02(:) = 1e-2_DP * sum(v_r_mag(:)) / (nfrag - 4) + f_vec_01(:) = f_vec_const(:) do j = 1, MAXITER