Improved energy and momentum bookkeeping

examples/symba_energy_momentum/escape.in

@@ -0,0 +1,18 @@
+3
+1 39.47841760435743
+0.0 0.0 0.0
+0.0 0.0 0.0
+0.4 0.4 0.4 !Ip 
+0.0 0.0 0.0 !rot
+2 1e-07 0.0009
+7e-06
+99.9 0.0 0.0
+100.00 10.00 0.0
+0.4 0.4 0.4 !Ip 
+0.0 0.0 0.0 !rot
+3 1e-08 0.0004
+3.25e-06
+1.0 4.20E-05 0.0
+0.00 -6.28 0.0
+0.4 0.4 0.4 !Ip 
+0.0 0.0 0.0 !rot

examples/symba_energy_momentum/param.escape.in

@@ -0,0 +1,33 @@
+T0             0.0e0
+TSTOP          0.00100           ! simulation length in seconds = 100 years
+DT             0.00010              ! stepsize in seconds
+PL_IN          escape.in
+TP_IN          tp.in
+IN_TYPE        ASCII
+ISTEP_OUT      1                 ! output cadence every year
+BIN_OUT        bin.escape.dat
+PARTICLE_FILE  particle.escape.dat
+OUT_TYPE       REAL8             ! double precision real output
+OUT_FORM       XV                ! osculating element output
+OUT_STAT       REPLACE
+ISTEP_DUMP     1                 ! system dump cadence
+J2             0.0               ! no J2 term
+J4             0.0               ! no J4 term
+CHK_CLOSE      yes               ! check for planetary close encounters
+CHK_RMIN       0.005
+CHK_RMAX       1e2
+CHK_EJECT      -1.0              ! ignore this check
+CHK_QMIN       -1.0              ! ignore this check
+!CHK_QMIN_COORD HELIO            ! commented out here
+!CHK_QMIN_RANGE 1.0 1000.0       ! commented out here
+ENC_OUT        enc.escape.dat
+EXTRA_FORCE    no                ! no extra user-defined forces
+BIG_DISCARD    no                ! output all planets if anything discarded
+RHILL_PRESENT  yes               ! Hill's sphere radii in input file
+MTINY          1.0e-16
+FRAGMENTATION  yes
+MU2KG          1.98908e30
+TU2S           3.1556925e7
+DU2M           1.49598e11
+ENERGY         yes
+ROTATION       yes

examples/symba_energy_momentum/param.merger.in

@@ -1,8 +1,6 @@
 !
 ! Parameter file for Mars in SI units
 !
-NPLMAX         3                 ! not used
-NTPMAX         0                 ! not used
 T0             0.0e0
 TSTOP          0.10              ! simulation length in seconds = 100 years
 DT             0.01              ! stepsize in seconds
@@ -11,6 +9,7 @@ TP_IN          tp.in
 IN_TYPE        ASCII
 ISTEP_OUT      1                 ! output cadence every year
 BIN_OUT        bin.merger.dat
+PARTICLE_FILE  particle.merger.dat
 OUT_TYPE       REAL8             ! double precision real output
 OUT_FORM       XV                ! osculating element output
 OUT_STAT       REPLACE

examples/symba_energy_momentum/sun.in

@@ -1,4 +1,4 @@
-2
+3
 1 39.47841760435743
 0.0 0.0 0.0
 0.0 0.0 0.0
@@ -10,7 +10,7 @@
 -10.00 2.00 0.0
 0.4 0.4 0.4 !Ip 
 0.0 0.0 0.0 !rot
-3 1e-08 0.0004
+3 1e-05 0.0004
 3.25e-06
 1.0 4.20E-05 0.0
 0.00 -6.28 0.0

src/io/io_conservation_report.f90

@@ -32,7 +32,8 @@ 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, &
-                  mass => symba_plA%helio%swiftest%mass, dMcb => symba_plA%helio%swiftest%dMcb, Mcb_initial => symba_plA%helio%swiftest%Mcb_initial)
+                Eescape => symba_plA%helio%swiftest%Eescape, 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
                open(unit = egyiu, file = energy_file, form = "formatted", status = "old", action = "write", position = "append")
@@ -58,8 +59,8 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
             Lmag_now = norm2(Ltot_now)
             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 - (Eorbit_orig + Ecollisions)) / abs(Eorbit_orig)
+            Ecoll_error = Ecollisions / abs(Eorbit_orig)
+            Etotal_error = (Eorbit - Ecollisions - Eorbit_orig - Eescape) / 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

src/main/swiftest_symba.f90

@@ -21,7 +21,7 @@ program swiftest_symba
    logical                       :: lfrag_add, ldiscard_pl, ldiscard_tp
    integer(I4B)                  :: nplm, ntp, ntp0, nsppl, nsptp, iout, idump, iloop, i
    integer(I4B)                  :: nplplenc, npltpenc, nmergeadd, nmergesub
-   real(DP)                      :: t, tfrac, tbase,  msys
+   real(DP)                      :: t, tfrac, tbase,  msys, Eorbit_orig
    real(DP)                      :: Ecollision, Eorbit_before, Eorbit_after, ke_orbit_before, ke_spin_before, ke_orbit_after, ke_spin_after, pe_before, pe_after
    real(DP), dimension(NDIM)     :: Ltot
    character(STRMAX)             :: inparfile
@@ -153,6 +153,7 @@ program swiftest_symba
       if (param%lenergy) then
          call coord_h2b(npl, symba_plA%helio%swiftest, msys)
          call io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, lterminal=.false.) 
+         call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_orig, Ltot)
       end if
       write(*, *) " *************** Main Loop *************** "
 
@@ -175,7 +176,6 @@ program swiftest_symba
          call symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, mergeadd_list, nmergeadd, param)
          ldiscard_pl = ldiscard_pl .or. lfrag_add
 
-         ! These next two blocks should be streamlined/improved but right now we treat discards separately from collisions so it has to be this way
          if (ldiscard_pl .or. ldiscard_tp) then
             if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot)
             call symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, &
@@ -190,7 +190,7 @@ program swiftest_symba
 
             if (param%lenergy)  then
                call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot)
-               Ecollision = Eorbit_before - Eorbit_after   ! Energy change resulting in this collisional event Total running energy offset from collision in this step
+               Ecollision = Eorbit_after - Eorbit_before   ! Energy change resulting in this collisional event Total running energy offset from collision in this step
                symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision
             end if
             !if (ntp > 0) call util_dist_index_pltp(nplm, ntp, symba_plA, symba_tpA)

src/modules/module_interfaces.f90

@@ -877,14 +877,14 @@ END SUBROUTINE symba_discard_sun_pl
      END INTERFACE
 
      INTERFACE
-            subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape)
+            subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body)
             use swiftest_globals
             use swiftest_data_structures
             implicit none
             type(user_input_parameters), intent(inout) :: param
             type(swiftest_pl), intent(inout) :: swiftest_plA
             integer(I4B), intent(in)    :: ipl
-            logical, intent(in)        :: lescape
+            logical, intent(in)        :: lescape_body
             end subroutine
       END INTERFACE
 

src/modules/swiftest_data_structures.f90

@@ -52,6 +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
       integer(I4B), dimension(:,:), allocatable   :: k_plpl
       integer(I8B)                                :: num_plpl_comparisons
    contains

src/symba/symba_discard_conserve_mtm.f90

@@ -1,4 +1,4 @@
-subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape)
+subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body)
    !! author: David A. Minton
    !! 
    !! Conserves system momentum when a body is lost from the system or collides with central body
@@ -10,63 +10,72 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape)
    type(swiftest_pl), intent(inout) :: swiftest_plA
    type(user_input_parameters), intent(inout) :: param
    integer(I4B), intent(in)    :: ipl
-   logical, intent(in)         :: lescape
+   logical, intent(in)         :: lescape_body
    ! Internals
    real(DP), dimension(NDIM) :: Lpl, Lcb, xcom, vcom
    real(DP)                  :: pe, ke_orbit, ke_spin
+   integer(I4B)              :: i
 
 
    associate(npl => swiftest_plA%nbody, xb => swiftest_plA%xb, vb => swiftest_plA%vb, &
       rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, radius => swiftest_plA%radius, mass => swiftest_plA%mass, &
       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)
-
-      xcom(:) = (mass(ipl) * xb(:, ipl) + mass(1) * xb(:,1)) / (mass(1) + mass(ipl))
-      vcom(:) = (mass(ipl) * vb(:, ipl) + mass(1) * vb(:,1)) / (mass(1) + mass(ipl))
-
-      call util_crossproduct(xb(:,ipl) - xcom(:), vb(:,ipl) - vcom(:), Lpl)
-      Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl))
-
-      call util_crossproduct(xb(:, 1) - xcom(:), vb(:, 1) - vcom(:), Lcb)
-      Lcb(:) = mass(1) * Lcb(:) 
+      Mcb_initial => swiftest_plA%Mcb_initial, dMcb => swiftest_plA%dMcb, &
+      Mescape => swiftest_plA%Mescape, Lescape => swiftest_plA%Lescape, Eescape => swiftest_plA%Eescape)
 
       ! Add the potential and kinetic energy of the lost body to the records
       pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1))
       ke_orbit = 0.5_DP * mass(ipl) * dot_product(vb(:, ipl), vb(:, ipl)) 
       ke_spin  = 0.5_DP * mass(ipl) * radius(ipl)**2 * Ip(3, ipl) * dot_product(rot(:, ipl), rot(:, ipl))
 
       ! Add the pre-collision ke of the central body to the records
-      ke_spin  = 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1)) + ke_spin
       ! Add planet mass to central body accumulator
-      if (lescape) then
+      if (lescape_body) then
          Mescape = Mescape + mass(ipl)
-         ke_orbit = 0.0_DP
+         call util_crossproduct(xb(:,ipl), vb(:,ipl), Lpl)
+         Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl))
+         Lescape(:) = Lescape(:)  + Lpl(:) 
+         do i = 2, npl
+            if (i == ipl) cycle
+            pe = pe - mass(i) * mass(ipl) / norm2(xb(:, ipl) - xb(:, i))
+         end do
       else
+         xcom(:) = (mass(ipl) * xb(:, ipl) + mass(1) * xb(:,1)) / (mass(1) + mass(ipl))
+         vcom(:) = (mass(ipl) * vb(:, ipl) + mass(1) * vb(:,1)) / (mass(1) + mass(ipl))
+
+         call util_crossproduct(xb(:,ipl) - xcom(:), vb(:,ipl) - vcom(:), Lpl)
+         Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl))
+
+         call util_crossproduct(xb(:, 1) - xcom(:), vb(:, 1) - vcom(:), Lcb)
+         Lcb(:) = mass(1) * Lcb(:) 
          ke_orbit = 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1)) + ke_orbit
+         ke_spin = ke_spin + 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1))
          ! Update mass of central body to be consistent with its total mass
          dMcb = dMcb + mass(ipl)
          dRcb = dRcb + 1.0_DP / 3.0_DP * (radius(ipl) / radius(1))**3 - 2.0_DP / 9.0_DP * (radius(ipl) / radius(1))**6
          mass(1) = Mcb_initial + dMcb
          radius(1) = Rcb_initial + dRcb
          param%rmin = radius(1)
-         ! Update position and velocity of central body
+         ! Add planet angular momentum to central body accumulator
+         dLcb(:) = Lpl(:) + Lcb(:) + dLcb(:)
+         ! Update rotation of central body to by consistent with its angular momentum 
+         rot(:,1) = (Lcb_initial(:) + dLcb(:)) / (Ip(3, 1) * mass(1) * radius(1)**2)        
+         ke_spin  = ke_spin - 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1)) 
          xb(:, 1) = xcom(:)
          vb(:, 1) = vcom(:)
          ke_orbit = ke_orbit - 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1)  ) 
       end if
-
-      ! Add planet angular momentum to central body accumulator
-      dLcb(:) = Lpl(:) + Lcb(:) + dLcb(:)
-
-      ! Update rotation of central body to by consistent with its angular momentum 
-      rot(:,1) = (Lcb_initial(:) + dLcb(:)) / (Ip(3, 1) * mass(1) * radius(1)**2)        
-
-      ! Remove the new kinetic energy of the central body from the records
-      ke_spin  = ke_spin - 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1)) 
+      ! Update position and velocity of central body
 
       ! We must do this for proper book-keeping, since we can no longer track this body's contribution to energy directly
-      Ecollisions  = Ecollisions - 2 * (ke_orbit + ke_spin + pe)
+      if (lescape_body) then
+         Ecollisions  = Ecollisions + ke_orbit + ke_spin + pe
+         Eescape  = Eescape - (ke_orbit + ke_spin + pe)
+      else
+         Ecollisions  = Ecollisions + pe 
+         Eescape = Eescape - pe
+      end if
 
       ! Update the heliocentric coordinates of everything else
       call coord_b2h(npl, swiftest_plA)