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

@@ -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

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

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

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

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
 

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