From ec271484d624b3f3e92847505579030267988428 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 16 Aug 2021 13:04:55 -0400 Subject: [PATCH] Getting the new fragmentation code more similar to the old one for comparison --- src/fragmentation/fragmentation.f90 | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 0dd7ddd64..31080d44b 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -132,7 +132,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, call calculate_system_energy(linclude_fragments=.true.) ke_frag_budget = -dEtot - Qloss - call define_coordinate_system() call set_fragment_tan_vel(lfailure) ke_avg_deficit = ke_avg_deficit - ke_radial subtry = subtry + 1 @@ -324,8 +323,10 @@ subroutine define_coordinate_system() x_col_unit(:) = y_col_unit(:) .cross. z_col_unit(:) rmag(:) = .mag. x_frag(:,:) + if (.not.any(rmag(:) > 0.0_DP)) return + call random_number(L_sigma(:,:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector, - ! otherwise we can get an ill-conditioned system + ! otherwise we can get an ill-conditioned system do concurrent(i = 1:nfrag, rmag(i) > 0.0_DP) v_r_unit(:, i) = x_frag(:, i) / rmag(i) v_h_unit(:, i) = z_col_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) @@ -690,15 +691,14 @@ subroutine set_fragment_tan_vel(lerr) ! This will be done using a linear solver that solves for the tangential velocities of the first 6 fragments, constrained by the linear and angular momentum vectors, ! which is embedded in a non-linear minimizer that will adjust the tangential velocities of the remaining i>6 fragments to minimize kinetic energy for a given momentum solution ! The initial conditions fed to the minimizer for the fragments will be the remaining angular momentum distributed between the fragments. - Li(:) = L_remainder(:) / nfrag - do concurrent (i = 1:nfrag) - v_t_initial(i) = norm2(Li(:)) / (m_frag(i) * norm2(x_frag(:,i))) - end do - - v_t_mag(:) = v_t_initial(:) - vb_frag(:,1:nfrag) = vmag_to_vb(v_r_mag(1:nfrag), v_r_unit(:,1:nfrag), v_t_mag(1:nfrag), v_t_unit(:,1:nfrag), m_frag(1:nfrag), vcom(:)) - do concurrent(i = 1:nfrag) - v_frag(:, i) = vb_frag(:, i) - vcom(:) + ! Li(:) = L_remainder(:) / nfrag + ! do concurrent (i = 1:nfrag) + ! v_t_initial(i) = norm2(Li(:)) / (m_frag(i) * norm2(x_frag(:,i))) + ! end do + do i = 1, nfrag + v_t_initial(i) = norm2(L_remainder(:)) / ((nfrag - i + 1) * m_frag(i) * norm2(x_frag(:,i))) + Li(:) = m_frag(i) * (x_frag(:,i) .cross. (v_t_initial(i) * v_t_unit(:, i))) + L_remainder(:) = L_remainder(:) - Li(:) end do ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum