From 243ffcd0d0b24780391b1ed2f793b8ad1b0e938b Mon Sep 17 00:00:00 2001 From: David Minton Date: Sun, 23 May 2021 13:30:09 -0400 Subject: [PATCH] Made a number of improvements to energy bookkeeping after identifying source of energy gain from mergers (spurious addition of potential energy) and improved robustness of symba_frag_pos --- src/io/io_conservation_report.f90 | 8 +- src/modules/swiftest_data_structures.f90 | 2 +- src/symba/symba_casemerge.f90 | 98 ++++++++++++++---------- src/symba/symba_discard_conserve_mtm.f90 | 6 +- src/symba/symba_frag_pos.f90 | 29 ++++--- src/util/util_minimize_bfgs.f90 | 17 +++- 6 files changed, 98 insertions(+), 62 deletions(-) diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90 index 40c8580bd..94020a34b 100644 --- a/src/io/io_conservation_report.f90 +++ b/src/io/io_conservation_report.f90 @@ -32,7 +32,7 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, "; DM/M0 = ", ES12.5)' associate(Ecollisions => symba_plA%helio%swiftest%Ecollisions, Lescape => symba_plA%helio%swiftest%Lescape, Mescape => symba_plA%helio%swiftest%Mescape, & - Eescape => symba_plA%helio%swiftest%Eescape, mass => symba_plA%helio%swiftest%mass, dMcb => symba_plA%helio%swiftest%dMcb, & + Euntracked => symba_plA%helio%swiftest%Euntracked, mass => symba_plA%helio%swiftest%mass, dMcb => symba_plA%helio%swiftest%dMcb, & Mcb_initial => symba_plA%helio%swiftest%Mcb_initial) if (lfirst) then if (param%out_stat == "OLD") then @@ -60,11 +60,15 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig Eorbit_error = (Eorbit - Eorbit_orig) / abs(Eorbit_orig) Ecoll_error = Ecollisions / abs(Eorbit_orig) - Etotal_error = (Eorbit - Ecollisions - Eorbit_orig - Eescape) / abs(Eorbit_orig) + Etotal_error = (Eorbit - Ecollisions - Eorbit_orig - Euntracked) / abs(Eorbit_orig) Merror = (Mtot_now - Mtot_orig) / Mtot_orig write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror if (Ecoll_error > 0.0_DP) then write(*,*) 'Something has gone wrong! Collisional energy should not be positive!' + write(*,*) 'dke_orbit: ',(ke_orbit - ke_orb_last) / abs(Eorbit_orig) + write(*,*) 'dke_spin : ',(ke_spin - ke_spin_last) / abs(Eorbit_orig) + write(*,*) 'dpe : ',(pe - pe_last) / abs(Eorbit_orig) + write(*,*) end if end if ke_orb_last = ke_orbit diff --git a/src/modules/swiftest_data_structures.f90 b/src/modules/swiftest_data_structures.f90 index 5aebc5a87..9d4cef58c 100644 --- a/src/modules/swiftest_data_structures.f90 +++ b/src/modules/swiftest_data_structures.f90 @@ -52,7 +52,7 @@ module swiftest_data_structures real(DP) :: Rcb_initial !! Initial radius of the central body real(DP) :: dRcb = 0.0_DP!! Change in the radius of the central body real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions - real(DP) :: Eescape = 0.0_DP !! Energy gained from system due to escaped bodies + real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies integer(I4B), dimension(:,:), allocatable :: k_plpl integer(I8B) :: num_plpl_comparisons contains diff --git a/src/symba/symba_casemerge.f90 b/src/symba/symba_casemerge.f90 index 054ad8d2e..4c19a9e5a 100644 --- a/src/symba/symba_casemerge.f90 +++ b/src/symba/symba_casemerge.f90 @@ -23,59 +23,75 @@ function symba_casemerge (symba_plA, family, nmergeadd, mergeadd_list, x, v, mas ! Result integer(I4B) :: status ! Internals - integer(I4B) :: j, mergeid, ibiggest - real(DP) :: mass_new, radius_new, volume_new + integer(I4B) :: i, j, mergeid, ibiggest, nfamily + real(DP) :: mass_new, radius_new, volume_new, pe real(DP), dimension(NDIM) :: xcom, vcom, xc, vc, xcrossv real(DP), dimension(2) :: vol real(DP), dimension(NDIM) :: L_orb_old, L_spin_old real(DP), dimension(NDIM) :: L_spin_new, rot_new, Ip_new - status = MERGED + associate(Mpl => symba_plA%helio%swiftest%mass, id => symba_plA%helio%swiftest%id, info => symba_plA%helio%swiftest%info, & + xb => symba_plA%helio%swiftest%xb, Euntracked => symba_plA%helio%swiftest%Euntracked, Ecollisions => symba_plA%helio%swiftest%Ecollisions) - mass_new = sum(mass(:)) + status = MERGED - ! The merged body's name will be that of the largest of the two parents - ibiggest = maxloc(symba_plA%helio%swiftest%mass(family(:)), dim=1) - mergeid = symba_plA%helio%swiftest%id(family(ibiggest)) + mass_new = sum(mass(:)) - ! Merged body is created at the barycenter of the original bodies - xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mass_new - vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mass_new + ! The merged body's name will be that of the largest of the two parents + ibiggest = maxloc(Mpl(family(:)), dim=1) + mergeid = id(family(ibiggest)) - ! Get mass weighted mean of Ip and - Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mass_new - vol(:) = 4._DP / 3._DP * PI * radius(:)**3 - volume_new = sum(vol(:)) - radius_new = (3 * volume_new / (4 * PI))**(1._DP / 3._DP) - - L_spin_old(:) = L_spin(:,1) + L_spin(:,2) - L_orb_old(:) = 0.0_DP - ! Compute orbital angular momentum of pre-impact system - do j = 1, 2 - xc(:) = x(:, j) - xcom(:) - vc(:) = v(:, j) - vcom(:) - call utiL_crossproduct(xc(:), vc(:), xcrossv(:)) - L_orb_old(:) = L_orb_old(:) + mass(j) * xcrossv(:) - end do + ! Merged body is created at the barycenter of the original bodies + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mass_new + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mass_new - ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body - L_spin_new(:) = L_orb_old(:) + L_spin_old(:) + ! Get mass weighted mean of Ip and + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mass_new + vol(:) = 4._DP / 3._DP * PI * radius(:)**3 + volume_new = sum(vol(:)) + radius_new = (3 * volume_new / (4 * PI))**(1._DP / 3._DP) + + L_spin_old(:) = L_spin(:,1) + L_spin(:,2) + L_orb_old(:) = 0.0_DP + ! Compute orbital angular momentum of pre-impact system + do i = 1, 2 + xc(:) = x(:, i) - xcom(:) + vc(:) = v(:, i) - vcom(:) + call utiL_crossproduct(xc(:), vc(:), xcrossv(:)) + L_orb_old(:) = L_orb_old(:) + mass(i) * xcrossv(:) + end do - ! Assume prinicpal axis rotation on 3rd Ip axis - rot_new(:) = L_spin_new(:) / (Ip_new(3) * mass_new * radius_new**2) + ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body + L_spin_new(:) = L_orb_old(:) + L_spin_old(:) - ! Populate the list of new bodies - call symba_merger_size_check(mergeadd_list, nmergeadd + 1) - nmergeadd = nmergeadd + 1 - mergeadd_list%id(nmergeadd) = mergeid - mergeadd_list%status(nmergeadd) = status - mergeadd_list%xb(:,nmergeadd) = xcom(:) - mergeadd_list%vb(:,nmergeadd) = vcom(:) - mergeadd_list%mass(nmergeadd) = mass_new - mergeadd_list%radius(nmergeadd) = radius_new - mergeadd_list%Ip(:,nmergeadd) = Ip_new(:) - mergeadd_list%rot(:,nmergeadd) = rot_new(:) - mergeadd_list%info(nmergeadd) = symba_plA%helio%swiftest%info(family(ibiggest)) + ! Assume prinicpal axis rotation on 3rd Ip axis + rot_new(:) = L_spin_new(:) / (Ip_new(3) * mass_new * radius_new**2) + + ! Keep track of the component of potential energy due to the pre-impact family for book-keeping + nfamily = size(family(:)) + pe = 0.0_DP + do j = 1, nfamily + do i = j + 1, nfamily + pe = pe - Mpl(i) * Mpl(j) / norm2(xb(:, i) - xb(:, j)) + end do + end do + Ecollisions = Ecollisions + pe + Euntracked = Euntracked - pe + + ! Populate the list of new bodies + call symba_merger_size_check(mergeadd_list, nmergeadd + 1) + nmergeadd = nmergeadd + 1 + mergeadd_list%id(nmergeadd) = mergeid + mergeadd_list%status(nmergeadd) = status + mergeadd_list%xb(:,nmergeadd) = xcom(:) + mergeadd_list%vb(:,nmergeadd) = vcom(:) + mergeadd_list%mass(nmergeadd) = mass_new + mergeadd_list%radius(nmergeadd) = radius_new + mergeadd_list%Ip(:,nmergeadd) = Ip_new(:) + mergeadd_list%rot(:,nmergeadd) = rot_new(:) + mergeadd_list%info(nmergeadd) = info(family(ibiggest)) + + end associate return end function symba_casemerge diff --git a/src/symba/symba_discard_conserve_mtm.f90 b/src/symba/symba_discard_conserve_mtm.f90 index df478d598..d2846a68a 100644 --- a/src/symba/symba_discard_conserve_mtm.f90 +++ b/src/symba/symba_discard_conserve_mtm.f90 @@ -22,7 +22,7 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body) Lcb_initial => swiftest_plA%Lcb_initial, dLcb => swiftest_plA%dLcb, Ecollisions => swiftest_plA%Ecollisions, & Rcb_initial => swiftest_plA%Rcb_initial, dRcb => swiftest_plA%dRcb, & Mcb_initial => swiftest_plA%Mcb_initial, dMcb => swiftest_plA%dMcb, & - Mescape => swiftest_plA%Mescape, Lescape => swiftest_plA%Lescape, Eescape => swiftest_plA%Eescape) + Mescape => swiftest_plA%Mescape, Lescape => swiftest_plA%Lescape, Euntracked => swiftest_plA%Euntracked) ! Add the potential and kinetic energy of the lost body to the records pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1)) @@ -71,10 +71,10 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body) ! We must do this for proper book-keeping, since we can no longer track this body's contribution to energy directly if (lescape_body) then Ecollisions = Ecollisions + ke_orbit + ke_spin + pe - Eescape = Eescape - (ke_orbit + ke_spin + pe) + Euntracked = Euntracked - (ke_orbit + ke_spin + pe) else Ecollisions = Ecollisions + pe - Eescape = Eescape - pe + Euntracked = Euntracked - pe end if ! Update the heliocentric coordinates of everything else diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index b7bd75334..7079f50ae 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -113,7 +113,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi logical, intent(out) :: lmerge ! Answers the question: Should this have been a merger instead? real(DP), intent(inout) :: Qloss ! Internals - real(DP), parameter :: f_spin = 0.03_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin + real(DP), parameter :: f_spin = 0.20_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 @@ -231,19 +231,21 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi ! (ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), & ! (pe_after - pe_before) / abs(Etot_before), & ! dEtot, dLmag - else - write(*,*) 'Failed try ',try,': Could not find radial velocities.' + !else + !write(*,*) 'Failed try ',try,': Could not find radial velocities.' end if if (.not.lmerge) exit call restructure_failed_fragments() + try = try + 1 end do !write(*, "(' -------------------------------------------------------------------------------------')") !write(*, "(' Final diagnostic')") !write(*, "(' -------------------------------------------------------------------------------------')") if (lmerge) then - write(*,*) "symba_frag_pos can't find a solution, so this collision is flagged as a merge" + write(*,*) "symba_frag_pos failed after: ",try," tries" else + write(*,*) "symba_frag_pos succeeded after: ",try," tries" ! write(*, "(' dL_tot should be very small' )") ! write(*,fmtlabel) ' dL_tot |', dLmag ! write(*, "(' dE_tot should be negative and equal to Qloss' )") @@ -660,7 +662,7 @@ subroutine set_fragment_radial_velocities(lmerge) ! Arguments logical, intent(out) :: lmerge ! Internals - real(DP), parameter :: TOL = 1e-10_DP + real(DP), parameter :: TOL = 1e-9_DP real(DP), dimension(:), allocatable :: vflat logical :: lerr integer(I4B) :: i @@ -670,6 +672,8 @@ subroutine set_fragment_radial_velocities(lmerge) ! Set the "target" ke_orb_after (the value of the orbital kinetic energy that the fragments ought to have) if ((ke_target < 0.0_DP) .or. (ke_offset > 0.0_DP)) then + !if (ke_target < 0.0_DP) write(*,*) 'Negative ke_target: ',ke_target + !if (ke_offset > 0.0_DP) write(*,*) 'Positive ke_offset: ',ke_offset lmerge = .true. return end if @@ -678,8 +682,8 @@ subroutine set_fragment_radial_velocities(lmerge) ! Initialize radial velocity magnitudes with a random value that is approximately 10% of that found by distributing the kinetic energy equally allocate(v_r_sigma, source=v_r_mag) call random_number(v_r_sigma(1:nfrag)) - v_r_sigma(1:nfrag) = 1.0_DP + 2 * (v_r_sigma(1:nfrag) - 0.5_DP) * 1e-1_DP - v_r_initial(1:nfrag) = v_r_sigma(1:nfrag) * abs(ke_offset) / (m_frag(1:nfrag) * nfrag) + v_r_sigma(1:nfrag) = sqrt((1.0_DP + 2 * (v_r_sigma(1:nfrag) - 0.5_DP) * 1e-3_DP)) + v_r_initial(1:nfrag) = v_r_sigma(1:nfrag) * sqrt(abs(ke_offset) / (2 * m_frag(1:nfrag) * nfrag)) ! Initialize the lambda function using a structure constructor that calls the init method ! Minimize the ke objective function using the BFGS optimizer @@ -782,17 +786,18 @@ subroutine restructure_failed_fragments() integer(I4B) :: i integer(I4B), save :: iflip = 1 - if (iflip == 1) then - iflip = 2 - else - iflip = 1 + if (ke_offset < 0.0_DP) then + if (iflip == 1) then + iflip = 2 + else + iflip = 1 + end if m_frag(iflip) = m_frag(iflip) + m_frag(nfrag) rad_frag(iflip) = (rad_frag(iflip)**3 + rad_frag(nfrag)**3)**(1._DP / 3._DP) m_frag(nfrag) = 0.0_DP rad_frag(nfrag) = 0.0_DP nfrag = nfrag - 1 end if - try = try + 1 r_max_start = r_max_start + 1.0_DP ! The larger lever arm can help if the problem is in the angular momentum step end subroutine restructure_failed_fragments diff --git a/src/util/util_minimize_bfgs.f90 b/src/util/util_minimize_bfgs.f90 index 775d7b4ab..45000b5f3 100644 --- a/src/util/util_minimize_bfgs.f90 +++ b/src/util/util_minimize_bfgs.f90 @@ -111,6 +111,7 @@ function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) end do call ieee_get_flag(ieee_usual, fpe_flag) lerr = lerr .or. any(fpe_flag) + !if (any(fpe_flag)) write(*,*) 'BFGS did not converge due to fpe' !if (lerr) write(*,*) "BFGS did not converge!" call ieee_set_status(original_fpe_status) @@ -283,7 +284,7 @@ subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) integer(I4B) :: i integer(I4B), parameter :: MAXLOOP = 1000 ! maximum number of loops before method is determined to have failed real(DP), parameter :: eps = epsilon(lo) ! small number precision to test floating point equality - real(DP), parameter :: dela = 12.4423_DP ! arbitrary number to test if function is constant + real(DP), parameter :: dela = 2.0324_DP ! arbitrary number to test if function is constant ! set up initial bracket points lerr = .false. @@ -338,9 +339,10 @@ subroutine bracket(f, x0, S, N, gam, step, lo, hi, lerr) f0 = n2one(f, x0, S, N, a0) else ! all values equal stops if there is no minimum or takes RHS min if it exists ! test if function itself is constant - fcon = n2one(f, x0, S, N, a2 + dela) !add by an arbitrary number to see if constant + fcon = n2one(f, x0, S, N, a2 * dela) !change a2 by an arbitrary number to see if constant if (abs(f2 - fcon) < eps) then lerr = .true. + !write(*,*) 'bracket determined function is constant' return ! function is constant end if a3 = a0 + 0.5_DP * (a1 - a0) @@ -471,6 +473,9 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) errval = abs((astar - aold) / astar) call ieee_get_flag(ieee_usual, fpe_flag) if (any(fpe_flag)) then + !write(*,*) 'quadfit fpe' + !write(*,*) 'aold : ',aold + !write(*,*) 'astar: ',astar lerr = .true. exit end if @@ -491,7 +496,11 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) soln(:) = util_solve_linear_system(lhs, rhs, 3, lerr) call ieee_set_flag(ieee_all, .false.) ! Set all flags back to quiet call ieee_set_halting_mode(ieee_divide_by_zero, .false.) - if (lerr) exit + if (lerr) then + !write(*,*) 'quadfit fpe:' + !write(*,*) 'uttil_solve_linear_system failed' + exit + end if aold = astar if (soln(2) == soln(3)) then ! Handles the case where they are both 0. 0/0 is an unhandled exception astar = 0.5_DP @@ -500,6 +509,8 @@ subroutine quadfit(f, x0, S, N, eps, lo, hi, lerr) end if call ieee_get_flag(ieee_usual, fpe_flag) if (any(fpe_flag)) then + !write(*,*) 'quadfit fpe' + !write(*,*) 'soln(2:3): ',soln(2:3) lerr = .true. exit end if