diff --git a/src/io/io.f90 b/src/io/io.f90 index 10bb911f8..d7b899475 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -39,7 +39,12 @@ module subroutine io_conservation_report(self, param, lterminal) write(EGYIU,EGYHEADER) end if end if - call system%get_energy_and_momentum(param, ke_orbit_now, ke_spin_now, pe_now, Lorbit_now, Lspin_now) + call system%get_energy_and_momentum(param) + ke_orbit_now = system%ke_orbit + ke_spin_now = system%ke_spin + pe_now = system%pe + Lorbit_now = system%Lorbit + Lspin_now = system%Lspin Eorbit_now = ke_orbit_now + ke_spin_now + pe_now Ltot_now(:) = Lorbit_now(:) + Lspin_now(:) + Lescape(:) Mtot_now = cb%mass + sum(pl%mass(1:npl)) + system%Mescape diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a4042eb27..10e60a491 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -276,10 +276,12 @@ module swiftest_classes class(swiftest_tp), allocatable :: tp !! Test particle data structure class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure real(DP) :: Gmtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion - real(DP) :: ke = 0.0_DP !! System kinetic energy + real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy + real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy real(DP) :: pe = 0.0_DP !! System potential energy real(DP) :: te = 0.0_DP !! System total energy - real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector + real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! System orbital angular momentum vector + real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! System spin angular momentum vector real(DP), dimension(NDIM) :: Lescape = 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) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions @@ -1019,15 +1021,10 @@ module subroutine util_resize_tp(self, nnew) integer(I4B), intent(in) :: nnew !! New size neded end subroutine util_resize_tp - module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin, pe, Lorbit, Lspin) + module subroutine util_get_energy_momentum_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(out) :: ke_orbit !! Orbital kinetic energy - real(DP), intent(out) :: ke_spin !! Spin kinetic energy - real(DP), intent(out) :: pe !! Potential energy - real(DP), dimension(:), intent(out) :: Lorbit !! Orbital angular momentum - real(DP), dimension(:), intent(out) :: Lspin !! Spin angular momentum end subroutine util_get_energy_momentum_system module subroutine util_set_beg_end_pl(self, xbeg, xend, vbeg) diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index edb641907..33b5c07c4 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -133,6 +133,7 @@ module subroutine setup_initialize_system(self, param) call self%tp%set_mu(self%cb) call self%pl%eucl_index() if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) + !if (param%lfirstenergy) then return end subroutine setup_initialize_system diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 6ab835e36..5f6d3926a 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -112,7 +112,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) Ltot(:) = Ltot(:) - Lpl(:) end do Ltot(:) = Ltot(:) - cb%mass * cb%xb(:) .cross. cb%vb(:) - system%Lescape(:) = system%Lescape(:) + system%Ltot(:) + system%Lescape(:) = system%Lescape(:) + Ltot(:) if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) else diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 69936e1b2..38701229d 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -1,7 +1,7 @@ submodule (swiftest_classes) s_util_get_energy_momentum use swiftest contains - module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin, pe, Lorbit, Lspin) + module subroutine util_get_energy_momentum_system(self, param) !! author: David A. Minton !! !! Compute total system angular momentum vector and kinetic, potential and total system energy @@ -12,11 +12,6 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(out) :: ke_orbit !! Orbital kinetic energy - real(DP), intent(out) :: ke_spin !! Spin kinetic energy - real(DP), intent(out) :: pe !! Potential energy - real(DP), dimension(:), intent(out) :: Lorbit !! Orbital angular momentum - real(DP), dimension(:), intent(out) :: Lspin !! Spin angular momentum ! Internals integer(I4B) :: i, j integer(I8B) :: k @@ -28,11 +23,12 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin logical, dimension(self%pl%nplpl) :: lstatpl logical, dimension(self%pl%nbody) :: lstatus - Lorbit(:) = 0.0_DP - Lspin(:) = 0.0_DP - ke_orbit = 0.0_DP - ke_spin = 0.0_DP associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) + system%Lorbit(:) = 0.0_DP + system%Lspin(:) = 0.0_DP + system%ke_orbit = 0.0_DP + system%ke_spin = 0.0_DP + kepl(:) = 0.0_DP Lplorbitx(:) = 0.0_DP Lplorbity(:) = 0.0_DP @@ -56,6 +52,7 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin ! Kinetic energy from orbit and spin kepl(i) = pl%mass(i) * v2 end do + if (param%lrotation) then do i = 1, npl rot2 = dot_product(pl%rot(:,i), pl%rot(:,i)) @@ -90,10 +87,10 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin end associate end do - ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) - if (param%lrotation) ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) + system%ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) + if (param%lrotation) system%ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) - pe = sum(pepl(:), lstatpl(:)) + sum(pecb(2:npl), lstatus(2:npl)) + system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(2:npl), lstatus(2:npl)) ! Potential energy from the oblateness term if (param%loblatecb) then @@ -102,17 +99,16 @@ module subroutine util_get_energy_momentum_system(self, param, ke_orbit, ke_spin irh(i) = 1.0_DP / norm2(pl%xh(:,i)) end do call obl_pot(npl, cb%mass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) - pe = pe + oblpot + system%pe = system%pe + oblpot end if - Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) - Lorbit(2) = sum(Lplorbity(1:npl), lstatus(1:npl)) - Lorbit(3) = sum(Lplorbitz(1:npl), lstatus(1:npl)) - - Lspin(1) = sum(Lplspinx(1:npl), lstatus(1:npl)) - Lspin(2) = sum(Lplspiny(1:npl), lstatus(1:npl)) - Lspin(3) = sum(Lplspinz(1:npl), lstatus(1:npl)) + system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) + system%Lorbit(2) = sum(Lplorbity(1:npl), lstatus(1:npl)) + system%Lorbit(3) = sum(Lplorbitz(1:npl), lstatus(1:npl)) + system%Lspin(1) = sum(Lplspinx(1:npl), lstatus(1:npl)) + system%Lspin(2) = sum(Lplspiny(1:npl), lstatus(1:npl)) + system%Lspin(3) = sum(Lplspinz(1:npl), lstatus(1:npl)) end associate return