From 05bb7284e4ad4aa631d3f29569000e0cb04f324e Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 4 Jun 2021 13:32:18 -0400 Subject: [PATCH] In the middle of troubleshooting failure to find tangential velocity. Added write statments to objective function --- src/symba/symba_frag_pos.f90 | 47 ++++++++++++++++++++++-------------- 1 file changed, 29 insertions(+), 18 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index ac93679f5..73bb44035 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -510,14 +510,12 @@ subroutine set_fragment_tangential_velocities(lerr) logical, intent(out) :: lerr ! Internals integer(I4B) :: i - real(DP) :: L_orb_mag, f_spin - real(DP), parameter :: TOL = 2e-8_DP + real(DP) :: L_orb_mag, f_spin, rbar + real(DP), parameter :: TOL = 1e-8_DP real(DP), dimension(:), allocatable :: v_t_initial type(lambda_obj) :: spinfunc type(lambda_obj_err) :: objective_function - real(DP), dimension(NDIM) :: L_frag_spin - real(DP), dimension(1) :: alpha - real(DP), dimension(1), parameter :: alpha_init = [1.0] + real(DP), dimension(NDIM) :: L_frag_spin, r_cross_L ! Initialize the fragments with 0 velocity and spin so we can divide up the balance between the tangential, radial, and spin components while conserving momentum lerr = .false. @@ -545,12 +543,15 @@ subroutine set_fragment_tangential_velocities(lerr) ! Divide up the pre-impact spin angular momentum equally amongst the fragments and calculate the spin kinetic energy do i = 1, nfrag rot_frag(:,i) = L_frag_spin(:) / (nfrag * m_frag(i) * Ip_frag(3, i) * rad_frag(i)**2) - v_t_initial(i) = L_orb_mag / (m_frag(i) * rmag(i) * nfrag) + call util_crossproduct(v_r_unit(:, i), z_col_unit, r_cross_L(:)) + rbar = rmag(i) * norm2(r_cross_L(:)) + v_t_mag(i) = L_orb_mag / (m_frag(i) * rbar * nfrag) end do - objective_function = lambda_obj(tangential_objective_function, lerr) - v_t_mag(1:nfrag) = util_minimize_bfgs(objective_function, nfrag, v_t_initial(1:nfrag), TOL, lerr) v_t_initial(:) = v_t_mag(:) v_t_mag(1:nfrag) = solve_fragment_tangential_velocities(v_t_mag_input=v_t_initial(7:nfrag), lerr=lerr) + v_t_initial(:) = v_t_mag(:) + objective_function = lambda_obj(tangential_objective_function, lerr) + v_t_mag(1:nfrag) = util_minimize_bfgs(objective_function, nfrag, v_t_initial(1:nfrag), TOL, lerr) ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) 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(:)) ke_frag = 0.0_DP @@ -566,7 +567,7 @@ subroutine set_fragment_tangential_velocities(lerr) lerr = (ke_radial < 0.0_DP) write(*,*) 'Tangential' write(*,*) 'ke_frag_budget: ',ke_frag_budget - write(*,*) 'ke_frag : ',ke_frag + write(*,*) 'ke_frag : ',ke_frag write(*,*) 'ke_frag_spin : ',ke_frag_spin write(*,*) 'ke_radial : ',ke_radial @@ -587,27 +588,37 @@ function tangential_objective_function(v_t_mag_input, lerr) result(fval) ! Internals integer(I4B) :: i real(DP), dimension(:,:), allocatable :: v_shift + real(DP), dimension(:), allocatable :: v_t_new real(DP), dimension(NDIM) :: L_frag, L - real(DP) :: dL + real(DP) :: dL, keo, kes lerr = .false. allocate(v_shift, mold=vb_frag) - v_shift(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_mag_input, v_t_unit, m_frag, vcom) + allocate(v_t_new, mold=v_t_mag_input) - ke_frag = 0.0_DP - ke_frag_spin = 0.0_DP + v_t_new(1:nfrag) = solve_fragment_tangential_velocities(v_t_mag_input=v_t_mag_input(7:nfrag), lerr=lerr) + v_shift(:,:) = vmag_to_vb(v_r_mag, v_r_unit, v_t_new, v_t_unit, m_frag, vcom) + + keo = 0.0_DP + kes = 0.0_DP L_frag(:) = 0.0_DP + write(*,*) 'tan obj' do i = 1, nfrag - ke_frag = ke_frag + m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) - ke_frag_spin = ke_frag_spin + m_frag(i) * Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) + keo = keo + m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) + kes = kes + m_frag(i) * Ip_frag(3, i) * rad_frag(i)**2 * dot_product(rot_frag(:, i), rot_frag(:, i)) call util_crossproduct(x_frag(:, i), v_shift(:, i), L(:)) L_frag(:) = L_frag(:) + L(:) * m_frag(i) + write(*,*) i,v_t_new(i) end do dL = norm2(L_frag(:) - L_frag_orb(:)) / norm2(L_frag_orb(:)) - ke_frag = 0.5_DP * ke_frag - ke_frag_spin = 0.5_DP * ke_frag_spin - fval = (ke_frag + ke_frag_spin) / ke_frag_budget * dL + keo = 0.5_DP * keo + kes = 0.5_DP * kes + fval = (keo + kes) / ke_frag_budget + write(*,*) 'ke_frag : ',keo + write(*,*) 'ke_frag_spin : ',kes + write(*,*) 'dL : ',dL + write(*,*) 'fval: ',fval lerr = .false. return end function tangential_objective_function