From 613e5d4996dcbdbe6ae9f2d93f6db28e465be73d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 31 Dec 2022 21:08:06 -0500 Subject: [PATCH] Added gravitational binding energy to the energy bookeeping. --- src/collision/collision_generate.f90 | 4 +- src/collision/collision_io.f90 | 4 ++ src/collision/collision_module.f90 | 12 +++-- src/collision/collision_resolve.f90 | 4 -- src/collision/collision_util.f90 | 39 +++++++++++----- src/fraggle/fraggle_generate.f90 | 28 ++++++------ src/fraggle/fraggle_util.f90 | 67 +++++++++++++++------------- src/netcdf_io/netcdf_io_module.f90 | 22 ++++----- src/swiftest/swiftest_io.f90 | 19 ++++++-- src/swiftest/swiftest_module.f90 | 3 ++ src/swiftest/swiftest_util.f90 | 4 +- src/symba/symba_discard.f90 | 12 +++-- 12 files changed, 132 insertions(+), 86 deletions(-) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 202cc5297..b946a387e 100644 --- a/src/collision/collision_generate.f90 +++ b/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 diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 3cf031b22..a798004d1 100644 --- a/src/collision/collision_io.f90 +++ b/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 diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index fed0191b4..3fa91a220 100644 --- a/src/collision/collision_module.f90 +++ b/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 diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index b0fff79dd..0036a8a00 100644 --- a/src/collision/collision_resolve.f90 +++ b/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") diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 3659bbc63..7f073eecf 100644 --- a/src/collision/collision_util.f90 +++ b/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 diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index a68cf7e4f..cce8f2af1 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/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(:) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index d405fc2b5..057dfe2b7 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -72,7 +72,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) integer(I4B) :: i, j, jproj, jtarg, nfrag, istart real(DP), dimension(2) :: volume real(DP), dimension(NDIM) :: Ip_avg - real(DP) :: mfrag, mremaining, min_mfrag, mtot, mcumul + real(DP) :: mfrag, mremaining, min_mfrag, mtot, mcumul, G real(DP), parameter :: BETA = 2.85_DP integer(I4B), parameter :: NFRAGMAX = 100 !! Maximum number of fragments that can be generated integer(I4B), parameter :: NFRAGMIN = 7 !! Minimum number of fragments that can be generated (set by the fraggle_generate algorithm for constraining momentum and energy) @@ -86,6 +86,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! Get mass weighted mean of Ip and density volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3 mtot = sum(impactors%mass(:)) + G = impactors%Gmass(1) / impactors%mass(1) Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / mtot if (impactors%mass(1) > impactors%mass(2)) then @@ -133,6 +134,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) select type(fragments => self%fragments) class is (collision_fragments(*)) fragments%mass(1) = impactors%mass_dist(1) + fragments%Gmass(1) = G * impactors%mass_dist(1) fragments%radius(1) = impactors%radius(jtarg) fragments%density(1) = impactors%mass_dist(1) / volume(jtarg) if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1) @@ -171,6 +173,8 @@ module subroutine fraggle_util_set_mass_dist(self, param) mremaining = fragments%mtot - sum(fragments%mass(1:nfrag)) fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining + fragments%Gmass(:) = G * fragments%mass(:) + ! Compute physical properties of the new fragments select case(impactors%regime) case(COLLRESOLVE_REGIME_HIT_AND_RUN) ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment @@ -190,25 +194,25 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! For catastrophic impacts, we will assign each of the n>2 fragments to one of the two original bodies so that the fragment cloud occupies ! roughly the same space as both original bodies. For all other disruption cases, we use body 2 as the center of the cloud. - fragments%origin_body(1) = 1 - fragments%origin_body(2) = 2 - if (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) then - mcumul = fragments%mass(1) - flipper = .true. - j = 2 - do i = 1, nfrag - if (flipper .and. (mcumul < impactors%mass(1))) then - flipper = .false. - j = 1 - else - j = 2 - flipper = .true. - end if - fragments%origin_body(i) = j - end do - else - fragments%origin_body(3:nfrag) = 2 - end if + fragments%origin_body(1) = 1 + fragments%origin_body(2) = 2 + if (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) then + mcumul = fragments%mass(1) + flipper = .true. + j = 2 + do i = 1, nfrag + if (flipper .and. (mcumul < impactors%mass(1))) then + flipper = .false. + j = 1 + else + j = 2 + flipper = .true. + end if + fragments%origin_body(i) = j + end do + else + fragments%origin_body(3:nfrag) = 2 + end if end select @@ -263,10 +267,11 @@ module subroutine fraggle_util_set_natural_scale_factors(self) impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do - fragments%mtot = fragments%mtot / collider%mscale - fragments%mass = fragments%mass / collider%mscale - fragments%radius = fragments%radius / collider%dscale - impactors%Qloss = impactors%Qloss / collider%Escale + fragments%mtot = fragments%mtot / collider%mscale + fragments%mass(:) = fragments%mass(:) / collider%mscale + fragments%Gmass(:) = fragments%Gmass(:) / (collider%dscale**3/collider%tscale**2) + fragments%radius(:) = fragments%radius(:) / collider%dscale + impactors%Qloss = impactors%Qloss / collider%Escale end associate return @@ -310,12 +315,13 @@ module subroutine fraggle_util_set_original_scale_factors(self) impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do - fragments%mtot = fragments%mtot * collider%mscale - fragments%mass = fragments%mass * collider%mscale - fragments%radius = fragments%radius * collider%dscale - fragments%rot = fragments%rot / collider%tscale - fragments%rc = fragments%rc * collider%dscale - fragments%vc = fragments%vc * collider%vscale + fragments%mtot = fragments%mtot * collider%mscale + fragments%mass(:) = fragments%mass(:) * collider%mscale + fragments%Gmass(:) = fragments%Gmass(:) * (collider%dscale**3/collider%tscale**2) + fragments%radius(:) = fragments%radius(:) * collider%dscale + fragments%rot(:,:) = fragments%rot(:,:) / collider%tscale + fragments%rc(:,:) = fragments%rc(:,:) * collider%dscale + fragments%vc(:,:) = fragments%vc(:,:) * collider%vscale do i = 1, fragments%nbody fragments%rb(:, i) = fragments%rc(:, i) + impactors%rbcom(:) @@ -330,6 +336,7 @@ module subroutine fraggle_util_set_original_scale_factors(self) collider%ke_orbit(:) = collider%ke_orbit(:) * collider%Escale collider%ke_spin(:) = collider%ke_spin(:) * collider%Escale collider%pe(:) = collider%pe(:) * collider%Escale + collider%be(:) = collider%be(:) * collider%Escale collider%Etot(:) = collider%Etot(:) * collider%Escale collider%mscale = 1.0_DP diff --git a/src/netcdf_io/netcdf_io_module.f90 b/src/netcdf_io/netcdf_io_module.f90 index 7199451b0..3926d2eea 100644 --- a/src/netcdf_io/netcdf_io_module.f90 +++ b/src/netcdf_io/netcdf_io_module.f90 @@ -96,24 +96,26 @@ module netcdf_io integer(I4B) :: k2_varid !! ID for the Love number variable character(NAMELEN) :: q_varname = "Q" !! name of the energy dissipation variable integer(I4B) :: Q_varid !! ID for the energy dissipation variable - character(NAMELEN) :: ke_orb_varname = "KE_orb" !! name of the nbody_system orbital kinetic energy variable - integer(I4B) :: KE_orb_varid !! ID for the nbody_system orbital kinetic energy variable - character(NAMELEN) :: ke_spin_varname = "KE_spin" !! name of the nbody_system spin kinetic energy variable - integer(I4B) :: KE_spin_varid !! ID for the nbody_system spin kinetic energy variable - character(NAMELEN) :: pe_varname = "PE" !! name of the nbody_system potential energy variable - integer(I4B) :: PE_varid !! ID for the nbody_system potential energy variable + character(NAMELEN) :: ke_orb_varname = "KE_orb" !! name of the system orbital kinetic energy variable + integer(I4B) :: KE_orb_varid !! ID for the system orbital kinetic energy variable + character(NAMELEN) :: ke_spin_varname = "KE_spin" !! name of the system spin kinetic energy variable + integer(I4B) :: KE_spin_varid !! ID for the system spin kinetic energy variable + character(NAMELEN) :: pe_varname = "PE" !! name of the system potential energy variable + integer(I4B) :: PE_varid !! ID for the system potential energy variable + character(NAMELEN) :: be_varname = "BE" !! name of the system binding energy variable + integer(I4B) :: BE_varid !! ID for the system binding energy variable character(NAMELEN) :: L_orb_varname = "L_orb" !! name of the orbital angular momentum vector variable - integer(I4B) :: L_orb_varid !! ID for the nbody_system orbital angular momentum vector variable + integer(I4B) :: L_orb_varid !! ID for the system orbital angular momentum vector variable character(NAMELEN) :: Lspin_varname = "Lspin" !! name of the spin angular momentum vector variable - integer(I4B) :: Lspin_varid !! ID for the nbody_system spin angular momentum vector variable + integer(I4B) :: Lspin_varid !! ID for the system spin angular momentum vector variable character(NAMELEN) :: L_escape_varname = "L_escape" !! name of the escaped angular momentum vector variable integer(I4B) :: L_escape_varid !! ID for the escaped angular momentum vector variable character(NAMELEN) :: Ecollisions_varname = "Ecollisions" !! name of the escaped angular momentum y variable integer(I4B) :: Ecollisions_varid !! ID for the energy lost in collisions variable character(NAMELEN) :: Euntracked_varname = "Euntracked" !! name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) integer(I4B) :: Euntracked_varid !! ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) - character(NAMELEN) :: GMescape_varname = "GMescape" !! name of the G*Mass of bodies that escape the nbody_system - integer(I4B) :: GMescape_varid !! ID for the G*Mass of bodies that escape the nbody_system + character(NAMELEN) :: GMescape_varname = "GMescape" !! name of the G*Mass of bodies that escape the system + integer(I4B) :: GMescape_varid !! ID for the G*Mass of bodies that escape the system character(NAMELEN) :: origin_type_varname = "origin_type" !! name of the origin type variable (Initial Conditions, Disruption, etc.) integer(I4B) :: origin_type_varid !! ID for the origin type character(NAMELEN) :: origin_time_varname = "origin_time" !! name of the time of origin variable diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 5f1a90276..fe4ecba18 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -121,7 +121,7 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen ! Internals real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now - real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now + real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now, be_now real(DP) :: GMtot_now character(len=STRMAX) :: errmsg integer(I4B), parameter :: EGYIU = 72 @@ -139,9 +139,10 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) ke_orbit_now = nbody_system%ke_orbit ke_spin_now = nbody_system%ke_spin pe_now = nbody_system%pe + be_now = nbody_system%be Lorbit_now(:) = nbody_system%Lorbit(:) Lspin_now(:) = nbody_system%Lspin(:) - Eorbit_now = ke_orbit_now + ke_spin_now + pe_now + Eorbit_now = ke_orbit_now + ke_spin_now + pe_now + be_now Ltot_now(:) = nbody_system%Ltot(:) + nbody_system%Lescape(:) GMtot_now = nbody_system%GMtot + nbody_system%GMescape @@ -149,6 +150,7 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) nbody_system%ke_orbit_orig = ke_orbit_now nbody_system%ke_spin_orig = ke_spin_now nbody_system%pe_orig = pe_now + nbody_system%be_orig = be_now nbody_system%Eorbit_orig = Eorbit_now nbody_system%GMtot_orig = GMtot_now nbody_system%Lorbit_orig(:) = Lorbit_now(:) @@ -161,6 +163,7 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) nbody_system%ke_orbit_error = (ke_orbit_now - nbody_system%ke_orbit_orig) / abs(nbody_system%Eorbit_orig) nbody_system%ke_spin_error = (ke_spin_now - nbody_system%ke_spin_orig) / abs(nbody_system%Eorbit_orig) nbody_system%pe_error = (pe_now - nbody_system%pe_orig) / abs(nbody_system%Eorbit_orig) + nbody_system%be_error = (be_now - nbody_system%be_orig) / abs(nbody_system%Eorbit_orig) nbody_system%Eorbit_error = (Eorbit_now - nbody_system%Eorbit_orig) / abs(nbody_system%Eorbit_orig) nbody_system%Ecoll_error = nbody_system%Ecollisions / abs(nbody_system%Eorbit_orig) nbody_system%Euntracked_error = nbody_system%Euntracked / abs(nbody_system%Eorbit_orig) @@ -512,7 +515,7 @@ module function swiftest_io_netcdf_get_old_t_final_system(self, param) result(ol real(DP), dimension(:), allocatable :: vals real(DP), dimension(1) :: rtemp real(DP), dimension(NDIM) :: rot0, Ip0, Lnow - real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig + real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, BE_orig associate (nc => param%system_history%nc, cb => self%cb) call nc%open(param) @@ -534,10 +537,13 @@ module function swiftest_io_netcdf_get_old_t_final_system(self, param) result(ol call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_old_t_final_system PE_varid" ) PE_orig = rtemp(1) + call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_old_t_final_system BE_varid" ) + BE_orig = rtemp(1) + call netcdf_io_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[1]), "netcdf_io_get_old_t_final_system Ecollisions_varid" ) call netcdf_io_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[1]), "netcdf_io_get_old_t_final_system Euntracked_varid" ) - self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + self%Ecollisions + self%Euntracked + self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + BE_orig + self%Ecollisions + self%Euntracked call netcdf_io_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_old_t_final_system L_orb_varid" ) call netcdf_io_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_old_t_final_system Lspin_varid" ) @@ -698,6 +704,7 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type, nc%time_dimid, nc%KE_orb_varid), "netcdf_io_initialize_output nf90_def_var KE_orb_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type, nc%time_dimid, nc%KE_spin_varid), "netcdf_io_initialize_output nf90_def_var KE_spin_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, nc%time_dimid, nc%PE_varid), "netcdf_io_initialize_output nf90_def_var PE_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type, nc%time_dimid, nc%BE_varid), "netcdf_io_initialize_output nf90_def_var PE_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_orb_varid), "netcdf_io_initialize_output nf90_def_var L_orb_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%Lspin_varid), "netcdf_io_initialize_output nf90_def_var Lspin_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%L_escape_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_escape_varid), "netcdf_io_initialize_output nf90_def_var L_escape_varid" ) @@ -861,6 +868,7 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) status = nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%KE_orb_varid) status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) + status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid) status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) @@ -1182,6 +1190,8 @@ module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar KE_spin_varid" ) status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar PE_varid" ) + status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, self%be, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar BE_varid" ) status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_orb_varid" ) status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) @@ -1606,6 +1616,7 @@ module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) call netcdf_io_check( nf90_put_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_orb_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_spin_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var PE_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%BE_varid, self%be, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var BE_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_orb_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var Lspin_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_escape_varid" ) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index a492bdda6..39a6b65e5 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -331,6 +331,7 @@ module swiftest real(DP) :: ke_orbit = 0.0_DP !! nbody_system orbital kinetic energy real(DP) :: ke_spin = 0.0_DP !! nbody_system spin kinetic energy real(DP) :: pe = 0.0_DP !! nbody_system potential energy + real(DP) :: be = 0.0_DP !! nbody_system binding energy of all bodies real(DP) :: te = 0.0_DP !! nbody_system total energy real(DP) :: oblpot = 0.0_DP !! nbody_system potential energy due to oblateness of the central body real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! nbody_system orbital angular momentum vector @@ -339,6 +340,7 @@ module swiftest real(DP) :: ke_orbit_orig = 0.0_DP !! Initial orbital kinetic energy real(DP) :: ke_spin_orig = 0.0_DP !! Initial spin kinetic energy real(DP) :: pe_orig = 0.0_DP !! Initial potential energy + real(DP) :: be_orig = 0.0_DP !! Initial binding energy real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy real(DP) :: GMtot_orig = 0.0_DP !! Initial nbody_system mass real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector @@ -353,6 +355,7 @@ module swiftest real(DP) :: ke_orbit_error = 0.0_DP real(DP) :: ke_spin_error = 0.0_DP real(DP) :: pe_error = 0.0_DP + real(DP) :: be_error = 0.0_DP real(DP) :: Eorbit_error = 0.0_DP real(DP) :: Ecoll_error = 0.0_DP real(DP) :: Euntracked_error = 0.0_DP diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index e3e0fccee..45b91aa73 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1322,7 +1322,9 @@ module subroutine swiftest_util_get_energy_momentum_system(self, param) nbody_system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), pl%lmask(1:npl)) end if - nbody_system%te = nbody_system%ke_orbit + nbody_system%ke_spin + nbody_system%pe + nbody_system%be = sum(-3*pl%Gmass(1:npl)*pl%mass(1:npl)/(5*pl%radius(1:npl)), pl%lmask(1:npl)) + + nbody_system%te = nbody_system%ke_orbit + nbody_system%ke_spin + nbody_system%pe + nbody_system%be nbody_system%Ltot(:) = nbody_system%Lorbit(:) + nbody_system%Lspin(:) end associate diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index cdeeeb0fe..216ab5f28 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -117,14 +117,15 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body logical, intent(in) :: lescape_body ! Internals real(DP), dimension(NDIM) :: Lpl, Ltot, Lcb, xcom, vcom - real(DP) :: pe, ke_orbit, ke_spin + real(DP) :: pe, be, ke_orbit, ke_spin integer(I4B) :: i, oldstat select type(cb => nbody_system%cb) class is (symba_cb) - ! Add the potential and kinetic energy of the lost body to the records + ! Add the potential, binding, and kinetic energy of the lost body to the records pe = -cb%Gmass * pl%mass(ipl) / norm2(pl%rb(:, ipl) - cb%rb(:)) + be = -3*pl%Gmass(ipl) * pl%mass(ipl) / (5 * pl%radius(ipl)) ke_orbit = 0.5_DP * pl%mass(ipl) * dot_product(pl%vb(:, ipl), pl%vb(:, ipl)) if (param%lrotation) then ke_spin = 0.5_DP * pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * dot_product(pl%rot(:, ipl), pl%rot(:, ipl)) @@ -194,11 +195,8 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, 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 - nbody_system%Ecollisions = nbody_system%Ecollisions + ke_orbit + ke_spin + pe - nbody_system%Euntracked = nbody_system%Euntracked - (ke_orbit + ke_spin + pe) - else - nbody_system%Ecollisions = nbody_system%Ecollisions + pe - nbody_system%Euntracked = nbody_system%Euntracked - pe + nbody_system%Ecollisions = nbody_system%Ecollisions + ke_orbit + ke_spin + pe + be + nbody_system%Euntracked = nbody_system%Euntracked - (ke_orbit + ke_spin + pe + be) end if end select