From 698ce177a0500322399221d9a69ddb60481a54c7 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 19 May 2021 12:26:44 -0400 Subject: [PATCH] Rearranged the fragment generation code for clarity and fixed numerous issues arising due to the internal unit conversion. Code is now passing some of the tests for the first time (head on disruption and supercatastrophic), but not others (off-axis cases still fail) --- src/symba/symba_frag_pos.f90 | 124 ++++++++++++++++++-------------- src/util/util_minimize_bfgs.f90 | 4 +- 2 files changed, 72 insertions(+), 56 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 3e7096bae..1218f42de 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -116,9 +116,9 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad real(DP), dimension(:,:), intent(inout) :: xb_frag, vb_frag, rot_frag logical, intent(out) :: lmerge ! Answers the question: Should this have been a merger instead? ! Internals - real(DP), parameter :: f_spin = 0.20_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin - real(DP) :: mscale, rscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system - real(DP) :: mtot + real(DP), parameter :: f_spin = 0.00_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin + real(DP) :: mscale = 1.0_DP, rscale = 1.0_DP, vscale = 1.0_DP, tscale = 1.0_DP, Lscale = 1.0_DP, Escale = 1.0_DP! Scale factors that reduce quantities to O(~1) in the collisional system + real(DP) :: mtot real(DP), dimension(NDIM) :: xcom, vcom integer(I4B) :: i, j, nfrag, fam_size logical, dimension(:), allocatable :: lexclude @@ -129,14 +129,23 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad real(DP), dimension(NDIM) :: Ltot_after real(DP) :: Etot_before, ke_orb_before, ke_spin_before, pe_before, Lmag_before real(DP) :: Etot_after, ke_orb_after, ke_spin_after, pe_after, Lmag_after - real(DP), dimension(NDIM) :: L_frag_spin, L_frag_tot, L_frag_orb - real(DP) :: ke_family, ke_target, ke_frag + real(DP), dimension(NDIM) :: L_frag_spin, L_frag_tot, L_frag_orb, L_offset + real(DP) :: ke_family, ke_target, ke_frag, ke_offset real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit character(len=*), parameter :: fmtlabel = "(A14,10(ES9.2,1X,:))" + integer(I4B), dimension(:), allocatable :: seed + integer(I4B) :: nseed + + ! Temporary until we get random number stuff settled for the whole project + call random_seed(size=nseed) + allocate(seed(nseed)) + seed = [(12*i, i=1, nseed)] + call random_seed(put=seed) allocate(x_frag, source=xb_frag) allocate(v_frag, source=vb_frag) nfrag = size(m_frag) + mtot = sum(m_frag) associate(npl => symba_plA%helio%swiftest%nbody, status => symba_plA%helio%swiftest%status) allocate(lexclude(npl)) @@ -147,32 +156,31 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad end where end associate - call set_scale_factors() - call define_coordinate_system() - call calculate_system_energy(ke_orb_before, ke_spin_before, pe_before, Ltot_before, linclude_fragments=.false.) Lmag_before = norm2(Ltot_before(:)) Etot_before = ke_orb_before + ke_spin_before + pe_before - call define_pre_collisional_family() - call set_fragment_position_vectors() - write(*, "(' ---------------------------------------------------------------------------')") write(*, "(' Energy normalized by |Etot_before|')") - write(*, "(' | T_orb T_spin T pe Etot Ltot')") write(*, "(' ---------------------------------------------------------------------------')") write(*,fmtlabel) ' Qloss |',-Qloss / abs(Etot_before) + + call set_scale_factors() + call define_coordinate_system() + call define_pre_collisional_family() + write(*, "(' ---------------------------------------------------------------------------')") write(*, "(' First pass to get angular momentum ')") - write(*, "(' ---------------------------------------------------------------------------')") - + call set_fragment_position_vectors() call set_fragment_tangential_velocities() - call calculate_system_energy(ke_orb_after, ke_spin_after, pe_after, Ltot_after, linclude_fragments=.true.) Etot_after = ke_orb_after + ke_spin_after + pe_after Lmag_after = norm2(Ltot_after(:)) + write(*, "(' ---------------------------------------------------------------------------')") + write(*, "(' | T_orb T_spin T pe Etot Ltot')") + write(*, "(' ---------------------------------------------------------------------------')") write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), & (ke_spin_after - ke_spin_before)/ abs(Etot_before), & (ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), & @@ -180,23 +188,21 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad (Etot_after - Etot_before) / abs(Etot_before), & norm2(Ltot_after - Ltot_before) / Lmag_before - call set_fragment_radial_velocities(lmerge) write(*, "(' ---------------------------------------------------------------------------')") write(*, "(' Second pass to get energy ')") write(*, "(' ---------------------------------------------------------------------------')") - write(*, "(' ---------------------------------------------------------------------------')") - write(*,fmtlabel) ' T_family |',ke_family / abs(Etot_before) - write(*,fmtlabel) ' T_frag targ |',ke_target / abs(Etot_before) - write(*, "(' ---------------------------------------------------------------------------')") - write(*,fmtlabel) ' T_frag calc |',ke_frag / abs(Etot_before) - write(*,fmtlabel) ' residual |',(ke_frag - ke_target) / abs(Etot_before) + + call set_fragment_radial_velocities(lmerge) call restore_scale_factors() call calculate_system_energy(ke_orb_after, ke_spin_after, pe_after, Ltot_after, linclude_fragments=.true.) Etot_after = ke_orb_after + ke_spin_after + pe_after Lmag_after = norm2(Ltot_after(:)) - + write(*,fmtlabel) ' T_frag calc |',ke_frag / abs(Etot_before) + write(*,fmtlabel) ' residual |',(ke_frag - ke_target) / abs(Etot_before) + write(*, "(' ---------------------------------------------------------------------------')") + write(*, "(' | T_orb T_spin T pe Etot Ltot')") write(*, "(' ---------------------------------------------------------------------------')") write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), & (ke_spin_after - ke_spin_before)/ abs(Etot_before), & @@ -246,6 +252,13 @@ subroutine set_scale_factors() m_frag = m_frag / mscale rad_frag = rad_frag / rscale Qloss = Qloss / Escale + + ke_orb_before = ke_orb_before / Escale + ke_spin_before = ke_spin_before / Escale + pe_before = pe_before * rscale / mscale**2 + Etot_before = ke_orb_before + ke_spin_before + pe_before + Ltot_before(:) = Ltot_before(:) / Lscale + Lmag_before = norm2(Ltot_before(:)) return end subroutine set_scale_factors @@ -281,8 +294,9 @@ subroutine restore_scale_factors() Lmag_before = Lmag_before * Lscale ke_orb_before = ke_orb_before * Escale ke_spin_before = ke_spin_before * Escale - pe_before = pe_before * Escale - Etot_before = Etot_before * Escale + pe_before = pe_before * mscale**2 / rscale + Etot_before = ke_orb_before + ke_spin_before + pe_before + Lmag_before = norm2(Ltot_before(:)) Qloss = Qloss * Escale mscale = 1.0_DP rscale = 1.0_DP @@ -414,11 +428,6 @@ subroutine calculate_system_energy(ke_orb, ke_spin, pe, Ltot, linclude_fragments ltmp(1:npl) = lexclude(1:npl) ltmp(npl+1:npl_new) = .false. call move_alloc(ltmp, lexclude) - - ke_frag = 0._DP - do i = 1, nfrag - ke_frag = ke_frag + 0.5_DP * m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) - end do end if where (lexclude(1:npl)) @@ -433,6 +442,17 @@ subroutine calculate_system_energy(ke_orb, ke_spin, pe, Ltot, linclude_fragments deallocate(symba_plwksp%helio%swiftest%k_plpl) nplm = count(symba_plA%helio%swiftest%mass > param%mtiny) if (lk_plpl) call util_dist_index_plpl(npl, nplm, symba_plA) + + ! Calculate the current fragment energy and momentum balances + if (linclude_fragments) then + ke_frag = 0._DP + do i = 1, nfrag + ke_frag = ke_frag + 0.5_DP * m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) + end do + ke_target = ke_family + (ke_spin_before - ke_spin) + (pe_before - pe) - Qloss + ke_offset = ke_frag - ke_target + L_offset(:) = Ltot_before(:) - Ltot(:) + end if end associate return end subroutine calculate_system_energy @@ -483,19 +503,19 @@ subroutine set_fragment_position_vectors() allocate(v_r_unit(NDIM,nfrag)) allocate(v_t_unit(NDIM,nfrag)) - ! Re-normalize position and velocity vectors by the fragment number so that for our initial guess we weight each - ! fragment position by the mass and assume equipartition of energy for the velocity - r_max = 1.5_DP * sum(rad_frag(:)) / PI + ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies + ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) + r_max = 4 * sum(rad_frag(:)) / PI - ! We will treat the first fragment of the list as a special case. - x_frag(:, 1) = -z_col_unit(:) + ! We will treat the first fragment of the list as a special case. It gets initialized the maximum distance along the original impactor distance vector. + ! This is done because in a regular disruption, the first body is the largest fragment. + x_frag(:, 1) = -y_col_unit(:) * r_max call random_number(x_frag(:,2:nfrag)) x_frag(:, :) = x_frag(:, :) * r_max call shift_vector_to_origin(m_frag, x_frag) do i = 1, nfrag - xb_frag(:,i) = x_frag(:,i) + xcom(:) rmag(i) = norm2(x_frag(:, i)) v_r_unit(:, i) = x_frag(:, i) / rmag(i) 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, @@ -503,6 +523,7 @@ subroutine set_fragment_position_vectors() L(:) = z_col_unit(:) + 2e-3_DP * (L_sigma(:) - 0.5_DP) L(:) = L(:) / norm2(L(:)) call util_crossproduct(L(:), v_r_unit(:, i), v_t_unit(:, i)) + xb_frag(:,i) = x_frag(:,i) + xcom(:) end do return @@ -539,7 +560,7 @@ subroutine set_fragment_tangential_velocities() call util_crossproduct(v_r_unit(:, i), v_t_unit(:, i), L(:)) A(4:6, i) = m_frag(i) * rmag(i) * L(:) else !For the remining bodies, distribute the angular momentum equally amongs them - v_t_mag(i) = L_orb_mag / (m_frag(i) * rmag(i) * nfrag) + v_t_mag(i) = L_orb_mag / (m_frag(i) * rmag(i) * nfrag) v_frag(:, i) = v_t_mag(i) * v_t_unit(:, i) L_lin_others(:) = L_lin_others(:) + m_frag(i) * v_frag(:, i) call util_crossproduct(x_frag(:, i), v_frag(:, i), L(:)) @@ -574,19 +595,18 @@ subroutine set_fragment_radial_velocities(lmerge) use symba_frag_pos_lambda_implementation implicit none ! Arguments - logical, intent(out) :: lmerge + logical, intent(out) :: lmerge ! Internals - real(DP), parameter :: TOL = epsilon(1.0_DP) + real(DP), parameter :: TOL = 1e-10_DP real(DP), dimension(:), allocatable :: vflat logical :: lerr integer(I4B) :: i - real(DP) :: ke_tangential real(DP), dimension(:), allocatable :: v_r_initial real(DP), dimension(:,:), allocatable :: v_r ! Set the "target" ke_orb_after (the value of the orbital kinetic energy that the fragments ought to have) - ke_target = ke_family + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss - if (ke_target < 0.0_DP) then + + if ((ke_target < 0.0_DP) .or. (ke_offset > 0.0_DP)) then lmerge = .true. return end if @@ -594,7 +614,7 @@ subroutine set_fragment_radial_velocities(lmerge) ! Initialize radial velocity magnitudes with a random value: allocate(v_r_initial, source=v_r_mag) call random_number(v_r_initial(:)) - v_r_initial(:) = 3 * v_r_initial(:) * sqrt(2 * ke_target / nfrag / m_frag(:)) + v_r_initial(:) = v_r_initial(:) * sqrt(2 * ke_target / nfrag / m_frag(:)) ! Initialize the lambda function using a structure constructor that calls the init method ! Minimize error using the BFGS optimizer v_r_mag(:) = util_minimize_bfgs(ke_constraint(ke_objective_function, v_r_unit, v_t_mag, v_t_unit, x_frag, m_frag, L_frag_orb, ke_target), nfrag, v_r_initial, TOL, lerr) @@ -608,7 +628,7 @@ subroutine set_fragment_radial_velocities(lmerge) ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) allocate(v_r, mold=v_frag) do i = 1, nfrag - v_r(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_r(:, i) = abs(v_r_mag(i)) * v_r_unit(:, i) end do call shift_vector_to_origin(m_frag, v_r) @@ -637,8 +657,8 @@ function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_f real(DP), dimension(:), intent(in) :: L_target !! Target orbital momentum real(DP), intent(in) :: ke_target !! Target kinetic energ ! Result - real(DP) :: fval !! The objective function result: norm of the vector composed of the tangential momentum and energy - !! Minimizing this brings us closer to our objective + real(DP) :: fval !! The objective function result, which is the square of the difference between the calculated fragment kinetic energy and our target + !! Minimizing this brings us closer to our objective ! Internals integer(I4B) :: i, nfrag, nsol real(DP), dimension(NDIM) :: L @@ -646,10 +666,11 @@ function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_f nfrag = size(m_frag) allocate(v_shift, mold=v_r_unit) - ! In order to keep satisfying the kinetic energy constraint, we must shift the origin of the radial component of the velocities to the center of mass + ! Make sure the velocity magnitude stays positive do i = 1, nfrag - v_shift(:,i) = v_r_mag(i) * v_r_unit(:, i) + v_shift(:,i) = abs(v_r_mag(i)) * v_r_unit(:, i) end do + ! In order to keep satisfying the kinetic energy constraint, we must shift the origin of the radial component of the velocities to the center of mass call shift_vector_to_origin(m_frag, v_shift) fval = -ke_target @@ -657,12 +678,7 @@ function ke_objective_function(v_r_mag, v_r_unit, v_t_mag, v_t_unit, x_frag, m_f v_shift(:, i) = v_shift(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:) fval = fval + 0.5_DP * m_frag(i) * dot_product(v_shift(:, i), v_shift(:, i)) end do - if (fval < 0.0_DP) then - fval = fval**8 - else - fval = fval**2 - end if - + fval = fval**2 ! This ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for return diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index 74b42a50e..9c0e3d73d 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -28,8 +28,8 @@ function util_minimize_bfgs(f, N, x0_d, eps_d, lerr) result(x1_d) real(DP), dimension(:), allocatable :: x1_d ! Internals integer(I4B) :: i, j, k, l, conv, num - integer(I4B), parameter :: MAXLOOP = 500 !! Maximum number of loops before method is determined to have failed - real(QP), parameter :: gradeps = 1e-8_QP !! Tolerance for gradient calculations + integer(I4B), parameter :: MAXLOOP = 10000 !! Maximum number of loops before method is determined to have failed + real(QP), parameter :: gradeps = 1e-5_QP !! Tolerance for gradient calculations real(QP), dimension(N) :: S !! Direction vectors real(QP), dimension(N) :: Snorm !! normalized direction real(QP), dimension(N,N) :: H !! Approximated inverse Hessian matrix