Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Streamlined the energy calc
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 5, 2021
1 parent aee2941 commit b4f1cb7
Show file tree
Hide file tree
Showing 5 changed files with 30 additions and 31 deletions.
7 changes: 6 additions & 1 deletion src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
13 changes: 5 additions & 8 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
1 change: 1 addition & 0 deletions src/setup/setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion src/symba/symba_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
38 changes: 17 additions & 21 deletions src/util/util_get_energy_momentum.f90
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down

0 comments on commit b4f1cb7

Please sign in to comment.