Improved energy balancing due to collisions with central body by allowing radius to change and also to keep track of energy terms that can no longer be accounted for (such as potential energy of a body that has collided with the central body)

src/io/io_conservation_report.f90

@@ -64,18 +64,12 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
             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(*,*) 'Comparisons with last time: '
-               write(*,*) 'ke_orb : ',ke_orbit - ke_orb_last, (ke_orbit - ke_orb_last) / abs(Eorbit_last)
-               write(*,*) 'ke_spin: ',ke_spin - ke_spin_last, (ke_spin - ke_spin_last) / abs(Eorbit_last)
-               write(*,*) 'pe     : ',pe - pe_last, (pe - pe_last) / abs(Eorbit_last)
-               write(*,*) 'Etot   : ',Eorbit - Eorbit_last, (Eorbit - Eorbit_last) / abs(Eorbit_last)
-               write(*,*)
             end if
+
          end if
          ke_orb_last = ke_orbit
          ke_spin_last = ke_spin
          pe_last = pe
-         Ltot_last(:) = Ltot_now(:)
          Eorbit_last = Eorbit
       end associate
       return

src/main/swiftest_symba.f90

@@ -146,6 +146,7 @@ program swiftest_symba
 
       ! Save initial mass and angular momentum of the central body
       symba_plA%helio%swiftest%Mcb_initial = symba_plA%helio%swiftest%mass(1)
+      symba_plA%helio%swiftest%Rcb_initial = symba_plA%helio%swiftest%radius(1)
       symba_plA%helio%swiftest%Lcb_initial(:) = symba_plA%helio%swiftest%Ip(3,1) * symba_plA%helio%swiftest%mass(1) * &
                                                 symba_plA%helio%swiftest%radius(1)**2 * symba_plA%helio%swiftest%rot(:,1)
 
@@ -189,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_before - Eorbit_after   ! 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,12 +877,13 @@ END SUBROUTINE symba_discard_sun_pl
      END INTERFACE
 
      INTERFACE
-            subroutine symba_discard_conserve_mtm(swiftest_plA, ipl, lescape)
+            subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape)
             use swiftest_globals
             use swiftest_data_structures
             implicit none
-            integer(I4B), intent(in)    :: ipl
+            type(user_input_parameters), intent(inout) :: param
             type(swiftest_pl), intent(inout) :: swiftest_plA
+            integer(I4B), intent(in)    :: ipl
             logical, intent(in)        :: lescape
             end subroutine
       END INTERFACE
@@ -1093,7 +1094,7 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA,
          type(symba_merger),          intent(inout) :: mergeadd_list 
          logical(LGT),                intent(in)    :: ldiscard_pl, ldiscard_tp 
          real(DP),                    intent(in)    :: mtiny
-         type(user_input_parameters), intent(in)    :: param
+         type(user_input_parameters), intent(inout) :: param
          logical, dimension(:), allocatable, intent(inout) :: discard_l_pl
          integer(I4B), dimension(:), allocatable, intent(inout) :: discard_stat_list
       end subroutine symba_rearray

src/modules/swiftest_data_structures.f90

@@ -44,11 +44,13 @@ module swiftest_data_structures
       real(DP),     dimension(:,:),   allocatable :: Ip     !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. 
       real(DP),     dimension(:,:),   allocatable :: rot    !! Body rotation vector in inertial coordinate frame 
       real(DP),     dimension(NDIM)               :: Lcb_initial  !! Initial angular momentum of the central body
-      real(DP),     dimension(NDIM)               :: dLcb = (/0.0_DP, 0.0_DP, 0.0_DP/) !! Change in angular momentum of the central body
-      real(DP)                                    :: Mcb_initial !! Initial mass and change in mass of the central body
+      real(DP),     dimension(NDIM)               :: dLcb = [0.0_DP, 0.0_DP, 0.0_DP] !! Change in angular momentum of the central body
+      real(DP),     dimension(NDIM)               :: Lescape = [0.0_DP, 0.0_DP, 0.0_DP] !! Angular momentum of bodies that escaped the system (used for bookeeping)
+      real(DP)                                    :: Mcb_initial !! Initial mass of the central body
       real(DP)                                    :: dMcb = 0.0_DP !! Change in mass of the central body
-      real(DP),     dimension(NDIM)               :: Lescape = (/0.0_DP, 0.0_DP, 0.0_DP/)!! Angular momentum of bodies that escaped the system (used for bookeeping)
       real(DP)                                    :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping)
+      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
       integer(I4B), dimension(:,:), allocatable   :: k_plpl
       integer(I8B)                                :: num_plpl_comparisons

src/symba/symba_discard_conserve_mtm.f90

@@ -1,45 +1,59 @@
-subroutine symba_discard_conserve_mtm(swiftest_plA, ipl, lescape)
+subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape)
    !! author: David A. Minton
    !! 
    !! Conserves system momentum when a body is lost from the system or collides with central body
 
    use swiftest
-   use module_interfaces, EXCEPT_THIS_ONE => symba_discard_conserve_mtm
+   use module_interfaces, except_this_one => symba_discard_conserve_mtm
    implicit none
-
-   integer(I4B), intent(in)    :: ipl
+   ! Arguments
    type(swiftest_pl), intent(inout) :: swiftest_plA
+   type(user_input_parameters), intent(inout) :: param
+   integer(I4B), intent(in)    :: ipl
    logical, intent(in)         :: lescape
-
+   ! Internals
    real(DP), dimension(NDIM) :: Lpl, Lcb, xcom, vcom
+   reaL(DP)                  :: pe, ke
 
    associate(npl => swiftest_plA%nbody, xb => swiftest_plA%xb, vb => swiftest_plA%vb, &
-      rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, rad => swiftest_plA%radius, mass => swiftest_plA%mass, &
+      rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, radius => swiftest_plA%radius, mass => swiftest_plA%mass, &
       statipl => swiftest_plA%status(ipl), xbpl => swiftest_plA%xb(:,ipl), xbcb => swiftest_plA%xb(:,1))
 
       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(:) + rad(ipl)**2 * Ip(3,ipl) * rot(:, 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(:) 
 
       ! Add planet mass to central body accumulator
       if (lescape) then
+         ! Add the potential and kinetic energy of the lost body to the records
+         pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1))
+         ke = 0.5_DP * mass(ipl) * dot_product(vb(:, ipl), vb(:, ipl))
+         swiftest_plA%Elost_bodies = swiftest_plA%Elost_bodies + pe + ke
          swiftest_plA%Mescape = swiftest_plA%Mescape + mass(ipl)
       else
+         ! Add the potential energy of the lost body to the records
+         pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1))
+         ke = 0.0_DP
+         swiftest_plA%Elost_bodies = swiftest_plA%Elost_bodies + pe
          swiftest_plA%dMcb = swiftest_plA%dMcb + mass(ipl)
+         swiftest_plA%dRcb = swiftest_plA%dRcb + 1.0_DP / 3.0_DP * (radius(ipl) / radius(1))**3 - 2.0_DP / 9.0_DP * (radius(ipl) / radius(1))**6
          ! Update mass of central body to be consistent with its total mass
          mass(1) = swiftest_plA%Mcb_initial + swiftest_plA%dMcb
+         radius(1) = swiftest_plA%Rcb_initial + swiftest_plA%dRcb
+         param%rmin = radius(1)
       end if
+      swiftest_plA%Ecollisions  = swiftest_plA%Ecollisions - (ke + pe)
 
       ! Add planet angular momentum to central body accumulator
       swiftest_plA%dLcb(:) = Lpl(:) + Lcb(:) + swiftest_plA%dLcb(:)
 
       ! Update rotation of central body to by consistent with its angular momentum 
-      swiftest_plA%rot(:,1) = (swiftest_plA%Lcb_initial(:) + swiftest_plA%dLcb(:)) / (Ip(3, 1) * mass(1) * rad(1)**2)        
+      swiftest_plA%rot(:,1) = (swiftest_plA%Lcb_initial(:) + swiftest_plA%dLcb(:)) / (Ip(3, 1) * mass(1) * radius(1)**2)        
 
       ! Update position and velocity of central body
       xb(:, 1) = xcom(:)

src/symba/symba_rearray.f90

@@ -22,7 +22,7 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA,
    type(symba_merger),          intent(inout) :: mergeadd_list 
    logical,                     intent(in)    :: ldiscard_pl, ldiscard_tp 
    real(DP),                    intent(in)    :: mtiny
-   type(user_input_parameters), intent(in)    :: param
+   type(user_input_parameters), intent(inout) :: param
    logical, dimension(:), allocatable, intent(inout) :: discard_l_pl
    integer(I4B), dimension(:), allocatable, intent(inout) :: discard_stat_list
 
@@ -50,7 +50,7 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA,
             cycle
          end if
          ! Resolve the discard
-         call symba_discard_conserve_mtm(symba_plA%helio%swiftest, discard_index_list(i), lescape)
+         call symba_discard_conserve_mtm(param, symba_plA%helio%swiftest, discard_index_list(i), lescape)
          ! Flip the main status flag to the discard state
       end do
       symba_plA%helio%swiftest%status(discard_index_list(:)) = discard_stat_list(:)