From d24de6d0e470ce1abb91439837f286342e92ea3c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 10 Aug 2021 15:07:31 -0400 Subject: [PATCH] Added total energy tracker and update the collisional energy bookkeeping term after discards --- src/symba/symba_discard.f90 | 12 +++++++++++- src/util/util_get_energy_momentum.f90 | 2 ++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 5f6d3926a..253fb2700 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -271,7 +271,6 @@ subroutine symba_discard_peri_pl(pl, system, param) pl%lfirst = lfirst_orig return - end subroutine symba_discard_peri_pl @@ -285,6 +284,8 @@ module subroutine symba_discard_pl(self, system, param) class(symba_pl), intent(inout) :: self !! SyMBA test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + real(DP) :: Eorbit_before, Eorbit_after select type(system) class is (symba_nbody_system) @@ -309,8 +310,17 @@ module subroutine symba_discard_pl(self, system, param) end if if (any(pl%ldiscard(:))) then + if (param%lenergy) then + call system%get_energy_and_momentum(param) + Eorbit_before = system%te + end if call symba_discard_nonplpl_conservation(self, system, param) call pl%rearray(system, param) + if (param%lenergy) then + call system%get_energy_and_momentum(param) + Eorbit_after = system%te + system%Ecollisions = Eorbit_after - Eorbit_before + end if end if end associate diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 700ecbe40..fa7cda43d 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -129,6 +129,8 @@ module subroutine util_get_energy_momentum_system(self, param) system%Lspin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), lstatus(1:npl)) system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), lstatus(1:npl)) end if + + system%te = system%ke_orbit + system%ke_spin + system%pe end associate return