Added gravitational binding energy to the energy bookeeping.

src/collision/collision_generate.f90

@@ -155,7 +155,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
       ! Internals
       integer(I4B)                              :: i, j, k, ibiggest
       real(DP), dimension(NDIM)                 :: Lspin_new
-      real(DP)                                  :: volume
+      real(DP)                                  :: volume, G
       character(len=STRMAX) :: message
 
       select type(nbody_system)
@@ -185,7 +185,9 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
 
                   ! Compute the physical properties of the new body after the merge.
                   volume = 4._DP / 3._DP * PI * sum(impactors%radius(:)**3)
+                  G = nbody_system%collider%impactors%Gmass(1) / nbody_system%collider%impactors%mass(1)
                   fragments%mass(1) = impactors%mass_dist(1)
+                  fragments%Gmass(1) = G * fragments%mass(1)
                   fragments%density(1) = fragments%mass(1) / volume
                   fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD)
                   if (param%lrotation) then

src/collision/collision_io.f90

@@ -231,6 +231,9 @@ module subroutine collision_io_netcdf_initialize_output(self, param)
                call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname,      nc%out_type,&
                   [                                                nc%stage_dimid,  nc%event_dimid], nc%PE_varid),      "collision_io_netcdf_initialize_output nf90_def_var PE_varid"  )
 
+               call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname,      nc%out_type,&
+                  [                                                nc%stage_dimid,  nc%event_dimid], nc%BE_varid),      "collision_io_netcdf_initialize_output nf90_def_var BE_varid"  )
+
                call netcdf_io_check( nf90_def_var(nc%id, nc%L_orb_varname,   nc%out_type, &
                   [              nc%space_dimid,                   nc%stage_dimid,  nc%event_dimid], nc%L_orb_varid),   "collision_io_netcdf_initialize_output nf90_def_var L_orb_varid"  )
 
@@ -332,6 +335,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param)
                call netcdf_io_check( nf90_put_var(nc%id, nc%ke_orb_varid,  collider%ke_orbit(:), start=[   1, eslot], count=[      2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_orb_varid before" )
                call netcdf_io_check( nf90_put_var(nc%id, nc%ke_spin_varid, collider%ke_spin(:),  start=[   1, eslot], count=[      2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_spin_varid before" )
                call netcdf_io_check( nf90_put_var(nc%id, nc%pe_varid,      collider%pe(:),       start=[   1, eslot], count=[      2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" )
+               call netcdf_io_check( nf90_put_var(nc%id, nc%be_varid,      collider%be(:),       start=[   1, eslot], count=[      2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" )
                call netcdf_io_check( nf90_put_var(nc%id, nc%L_orb_varid,   collider%Lorbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_orb_varid before" )
                call netcdf_io_check( nf90_put_var(nc%id, nc%Lspin_varid,   collider%Lspin(:,:),  start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Lspin_varid before" )
             end if

src/collision/collision_module.f90

@@ -98,6 +98,7 @@ module collision
       real(DP),                  dimension(NDIM,nbody)       :: vb           !! Barycentric velocity
       real(DP),                  dimension(NDIM,nbody)       :: rot          !! rotation vectors of fragments
       real(DP),                  dimension(NDIM,nbody)       :: Ip           !! Principal axes moment of inertia for fragments
+      real(DP),                  dimension(nbody)            :: Gmass        !! G*mass of fragments
       real(DP),                  dimension(nbody)            :: mass         !! masses of fragments
       real(DP),                  dimension(nbody)            :: radius       !! Radii  of fragments
       real(DP),                  dimension(nbody)            :: density      !! Radii  of fragments
@@ -117,14 +118,16 @@ module collision
       real(DP), dimension(NDIM,nbody)                        :: Lspin        !! Spin angular momentum vector of each individual fragment
       real(DP)                                               :: ke_orbit_tot !! Orbital kinetic energy of all fragments
       real(DP)                                               :: ke_spin_tot  !! Spin kinetic energy of all fragments
-      real(DP)                                               :: ke_budget    !! Kinetic energy budget for computing fragment trajectories
+      real(DP)                                               :: pe           !! Potential energy of all fragments
+      real(DP)                                               :: be           !! Binding energy of all fragments
+      real(DP)                                               :: E_budget    !! Kinetic energy budget for computing fragment trajectories
       real(DP), dimension(NDIM)                              :: L_budget     !! Angular momentum budget for computing fragment trajectories
       real(DP),                  dimension(nbody)            :: ke_orbit     !! Orbital kinetic energy of each individual fragment
       real(DP),                  dimension(nbody)            :: ke_spin      !! Spin kinetic energy of each individual fragment
    contains
       procedure :: reset                 => collision_util_reset_fragments      !! Deallocates all allocatable arrays and sets everything else to 0
       procedure :: get_angular_momentum  => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments
-      procedure :: get_kinetic_energy    => collision_util_get_kinetic_energy   !! Calcualtes the current kinetic energy of the fragments
+      procedure :: get_energy    => collision_util_get_energy   !! Calcualtes the current kinetic energy of the fragments
       procedure :: set_coordinate_system => collision_util_set_coordinate_fragments !! Sets the coordinate system of the fragments
       final     ::                          collision_final_fragments           !! Finalizer deallocates all allocatables
    end type collision_fragments
@@ -148,6 +151,7 @@ module collision
       real(DP), dimension(2)      :: ke_orbit !! Before/after orbital kinetic energy
       real(DP), dimension(2)      :: ke_spin  !! Before/after spin kinetic energy
       real(DP), dimension(2)      :: pe       !! Before/after potential energy
+      real(DP), dimension(2)      :: be       !! Before/after binding energy
       real(DP), dimension(2)      :: Etot     !! Before/after total nbody_system energy
    contains
       procedure :: setup                      => collision_util_setup_collider             !! Initializer for the encounter collision system and the before/after snapshots
@@ -385,10 +389,10 @@ module subroutine collision_util_get_angular_momentum(self)
          class(collision_fragments(*)), intent(inout) :: self !! Fragment system object
       end subroutine collision_util_get_angular_momentum
 
-      module subroutine collision_util_get_kinetic_energy(self) 
+      module subroutine collision_util_get_energy(self) 
          implicit none
          class(collision_fragments(*)), intent(inout) :: self !! Fragment system object
-      end subroutine collision_util_get_kinetic_energy
+      end subroutine collision_util_get_energy
 
       module subroutine collision_util_reset_fragments(self)
          implicit none

src/collision/collision_resolve.f90

@@ -548,10 +548,6 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
                      call collider%generate(nbody_system, param, t)
 
                      call nbody_system%get_energy_and_momentum(param)
-                     collider%pe(2) = nbody_system%pe
-                     dpe = collider%pe(2) - collider%pe(1) 
-                     nbody_system%Ecollisions = nbody_system%Ecollisions - dpe 
-                     nbody_system%Euntracked  = nbody_system%Euntracked + dpe
 
                      call collision_history%take_snapshot(param,nbody_system, t, "after") 
 

src/collision/collision_util.f90

@@ -200,15 +200,16 @@ module subroutine collision_util_get_angular_momentum(self)
    end subroutine collision_util_get_angular_momentum
 
 
-   module subroutine collision_util_get_kinetic_energy(self) 
+   module subroutine collision_util_get_energy(self) 
       !! Author: David A. Minton
       !!
-      !! Calculates the current kinetic energy of the fragments
+      !! Calculates the current energy of the fragments
       implicit none
       ! Argument
       class(collision_fragments(*)), intent(inout)  :: self !! Fragment system object
       ! Internals
-      integer(I4B) :: i
+      integer(I4B) :: i,j
+      real(DP), dimension(self%nbody) :: pepl
 
       associate(fragments => self, nfrag => self%nbody)
 
@@ -222,10 +223,19 @@ module subroutine collision_util_get_kinetic_energy(self)
          fragments%ke_orbit_tot = sum(fragments%ke_orbit(:))
          fragments%ke_spin_tot = sum(fragments%ke_spin(:))
 
+         fragments%pe = 0.0_DP
+         do i = 1, nfrag
+            do concurrent(j = i+1:nfrag)
+               pepl(j) = - (fragments%Gmass(i) * fragments%mass(j)) / .mag.(fragments%rc(:,i) - fragments%rc(:,j))
+            end do
+            fragments%pe = fragments%pe + sum(pepl(i+1:nfrag))
+         end do
+         fragments%be = -sum(3*fragments%Gmass(:)*fragments%mass(:)/(5*fragments%radius(:)))
+
       end associate
 
       return
-   end subroutine collision_util_get_kinetic_energy
+   end subroutine collision_util_get_energy
 
 
    module subroutine collision_util_get_energy_momentum(self,  nbody_system, param, lbefore)
@@ -242,9 +252,10 @@ module subroutine collision_util_get_energy_momentum(self,  nbody_system, param,
       class(base_parameters),   intent(inout) :: param   !! Current swiftest run configuration parameters
       logical,                  intent(in)    :: lbefore !! Flag indicating that this the "before" state of the nbody_system, with impactors included and fragments excluded or vice versa
       ! Internals
-      integer(I4B)  :: stage,i
+      integer(I4B)  :: stage,i,j
       real(DP), dimension(NDIM) :: Lorbit, Lspin
-      real(DP) :: ke_orbit, ke_spin
+      real(DP) :: ke_orbit, ke_spin, pe, be
+      real(DP), dimension(self%fragments%nbody) :: pepl
 
       select type(nbody_system)
       class is (swiftest_nbody_system)
@@ -254,7 +265,6 @@ module subroutine collision_util_get_energy_momentum(self,  nbody_system, param,
 
 
             if (lbefore) then
-
                Lorbit(:) = sum(impactors%Lorbit(:,:), dim=2)
                Lspin(:) = sum(impactors%Lspin(:,:), dim=2)
                ke_orbit = 0.0_DP
@@ -265,15 +275,20 @@ module subroutine collision_util_get_energy_momentum(self,  nbody_system, param,
                end do
                ke_orbit = ke_orbit / 2
                ke_spin = ke_spin / 2
+
+               pe = -(impactors%Gmass(1) * impactors%mass(2)) / .mag.(impactors%rc(:,2) - impactors%rc(:,1))
+               be = -sum(3*impactors%Gmass(:)*impactors%mass(:)/(5*impactors%radius(:)))
 
             else
                call fragments%get_angular_momentum()
                Lorbit(:) = fragments%Lorbit_tot(:) 
                Lspin(:) = fragments%Lspin_tot(:) 
 
-               call fragments%get_kinetic_energy()
+               call fragments%get_energy()
                ke_orbit = fragments%ke_orbit_tot
                ke_spin = fragments%ke_spin_tot
+               pe = fragments%pe
+               be = fragments%be
 
             end if 
             ! Calculate the current fragment energy and momentum balances
@@ -287,7 +302,9 @@ module subroutine collision_util_get_energy_momentum(self,  nbody_system, param,
             self%Ltot(:,stage) = Lorbit(:) + Lspin(:)
             self%ke_orbit(stage) = ke_orbit
             self%ke_spin(stage) = ke_spin
-            self%Etot(stage) = ke_orbit + ke_spin
+            self%pe(stage) = pe
+            self%be(stage) = be
+            self%Etot(stage) = ke_orbit + ke_spin + pe + be
          end associate
       end select
       end select
@@ -429,7 +446,7 @@ module subroutine collision_util_set_budgets(self)
       associate(impactors => self%impactors, fragments => self%fragments)
 
          fragments%L_budget(:) = self%Ltot(:,1)
-         fragments%ke_budget = self%Etot(1) - impactors%Qloss
+         fragments%E_budget = self%Etot(1) - impactors%Qloss
 
       end associate
 
@@ -616,7 +633,7 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag)
       self%fragments%L_budget(:) = 0.0_DP
       self%fragments%ke_orbit_tot = 0.0_DP
       self%fragments%ke_spin_tot = 0.0_DP
-      self%fragments%ke_budget = 0.0_DP
+      self%fragments%E_budget = 0.0_DP
 
       return
    end subroutine collision_util_setup_fragments_collider

src/fraggle/fraggle_generate.f90

@@ -373,7 +373,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
       integer(I4B) :: i, j, loop, try, istart, n, ndof
       logical :: lhitandrun, lsupercat
       real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual, vshear, vunit
-      real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof, ke_tot, ke_residual_min
+      real(DP) :: vmag, vesc, rotmag, E_residual, ke_per_dof, ke_tot, E_residual_min
       integer(I4B), dimension(collider%fragments%nbody) :: vsign
       real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail
       integer(I4B), parameter :: MAXLOOP = 100
@@ -442,30 +442,30 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
             istart = 1
          end if
          call fragments%set_coordinate_system()
-         ke_residual_min = -huge(1.0_DP)
+         E_residual_min = -huge(1.0_DP)
          outer: do try = 1, MAXTRY
             do loop = 1, MAXLOOP
-               call fragments%get_kinetic_energy()
-               ke_residual = fragments%ke_budget - (fragments%ke_orbit_tot + fragments%ke_spin_tot)
-               if ((abs(ke_residual) < abs(ke_residual_min)) .or. ((ke_residual >= 0.0_DP) .and. (ke_residual_min < 0.0_DP))) then ! This is our best case so far. Save it for posterity
+               call fragments%get_energy()
+               E_residual = fragments%E_budget - (fragments%ke_orbit_tot + fragments%ke_spin_tot + fragments%pe + fragments%be)
+               if ((abs(E_residual) < abs(E_residual_min)) .or. ((E_residual >= 0.0_DP) .and. (E_residual_min < 0.0_DP))) then ! This is our best case so far. Save it for posterity
                   if (allocated(collider%fragments)) deallocate(collider%fragments)
                   allocate(collider%fragments, source=fragments)
-                  ke_residual_min = ke_residual
-                  if ((ke_residual > 0.0_DP) .and. (ke_residual < TOL * fragments%ke_budget)) exit outer
+                  E_residual_min = E_residual
+                  if ((E_residual > 0.0_DP) .and. (E_residual < TOL * fragments%E_budget)) exit outer
                end if
                ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape
                ke_avail(:) = fragments%ke_orbit(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%rmag(:)
                ke_tot = 0.0_DP
-               ke_per_dof = -ke_residual
+               ke_per_dof = -E_residual
                do i = 1, 2*(nfrag - istart + 1)
-                  n = count(ke_avail(istart:nfrag) > -ke_residual/i) 
-                  if (ke_residual < 0.0_DP) n = n + count(fragments%ke_spin(istart:nfrag) > -ke_residual/i)
+                  n = count(ke_avail(istart:nfrag) > -E_residual/i) 
+                  if (E_residual < 0.0_DP) n = n + count(fragments%ke_spin(istart:nfrag) > -E_residual/i)
                   if (abs(n * ke_per_dof) > ke_tot) then
-                     ke_per_dof = -ke_residual/i
+                     ke_per_dof = -E_residual/i
                      ke_tot = n * ke_per_dof
                      ndof = i
-                     if (abs(ke_tot) > abs(ke_residual)) then
-                        ke_tot = -ke_residual
+                     if (abs(ke_tot) > abs(E_residual)) then
+                        ke_tot = -E_residual
                         ke_per_dof = ke_tot/n
                         exit
                      end if
@@ -511,7 +511,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure)
             call fraggle_generate_rot_vec(collider)
             collider%fail_scale = collider%fail_scale*1.01_DP
          end do outer
-         lfailure = ke_residual < 0.0_DP
+         lfailure = E_residual < 0.0_DP
 
          do concurrent(i = 1:nfrag)
             fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:)