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

Commit

Permalink
Added energy and momentum report back in
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 5, 2021
1 parent 5dee8d3 commit fb21ff7
Show file tree
Hide file tree
Showing 4 changed files with 299 additions and 10 deletions.
86 changes: 86 additions & 0 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,92 @@
use swiftest
contains

module subroutine io_conservation_report(self, param, lterminal)
!! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
!!
!! Reports the current state of energy, mass, and angular momentum conservation in a run
implicit none
! Arguments
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
class(swiftest_parameters), intent(inout) :: param !! Input colleciton of user-defined parameters
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), dimension(NDIM), save :: Ltot_last, Lorbit_last, Lspin_last
real(DP), save :: ke_orbit_last, ke_spin_last, pe_last, Eorbit_last
real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now
real(DP) :: Eorbit_error, Etotal_error, Ecoll_error
real(DP) :: Mtot_now, Merror
real(DP) :: Lmag_now, Lerror
character(len=*), parameter :: EGYFMT = '(ES23.16,10(",",ES23.16,:))' ! Format code for all simulation output
character(len=*), parameter :: EGYHEADER = '("t,Eorbit,Ecollisions,Lx,Ly,Lz,Mtot")'
integer(I4B), parameter :: EGYIU = 72
character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5 &
"; DEcollisions/|E0| = ", ES12.5, &
"; D(Eorbit+Ecollisions)/|E0| = ", ES12.5, &
"; DM/M0 = ", ES12.5)'

associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, Ecollisions => self%Ecollisions, Lescape => self%Lescape, Mescape => self%Mescape, &
Euntracked => self%Euntracked, Eorbit_orig => param%Eorbit_orig, Mtot_orig => param%Mtot_orig, &
Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), &
lfirst => param%lfirstenergy)
if (lfirst) then
if (param%out_stat == "OLD") then
open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "old", action = "write", position = "append")
else
open(unit = EGYIU, file = ENERGY_FILE, form = "formatted", status = "replace", action = "write")
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)
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
if (lfirst) then
Eorbit_orig = Eorbit_now
Mtot_orig = Mtot_now
Lorbit_orig(:) = Lorbit_now(:)
Lspin_orig(:) = Lspin_now(:)
Ltot_orig(:) = Ltot_now(:)
Lmag_orig = norm2(Ltot_orig(:))
lfirst = .false.
end if

write(EGYIU,EGYFMT) param%t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now
flush(EGYIU)
if (.not.lfirst .and. lterminal) then
Lmag_now = norm2(Ltot_now)
Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig
Eorbit_error = (Eorbit_now - Eorbit_orig) / abs(Eorbit_orig)
Ecoll_error = Ecollisions / abs(Eorbit_orig)
Etotal_error = (Eorbit_now - Ecollisions - Eorbit_orig - Euntracked) / abs(Eorbit_orig)
Merror = (Mtot_now - Mtot_orig) / Mtot_orig
write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror
if (Ecoll_error > 0.0_DP) then
write(*,*) 'Something has gone wrong! Collisional energy should not be positive!'
write(*,*) 'dke_orbit: ',(ke_orbit_now - ke_orbit_last) / abs(Eorbit_orig)
write(*,*) 'dke_spin : ',(ke_spin_now - ke_spin_last) / abs(Eorbit_orig)
write(*,*) 'dpe : ',(pe_now - pe_last) / abs(Eorbit_orig)
write(*,*)
end if
if (Lerror > 1e-6) then
write(*,*) 'Something has gone wrong! Angular momentum is too high!'
write(*,*) 'Lerror: ', Lerror
end if
end if
ke_orbit_last = ke_orbit_now
ke_spin_last = ke_spin_now
pe_last = pe_now
Eorbit_last = Eorbit_now
Lorbit_last(:) = Lorbit_now(:)
Lspin_last(:) = Lspin_now(:)
Ltot_last(:) = Ltot_now(:)
end associate
return

end subroutine io_conservation_report


module subroutine io_dump_param(self, param_file_name)
!! author: David A. Minton
!!
Expand Down
62 changes: 52 additions & 10 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,9 @@ module swiftest_classes
real(QP) :: DU2M = -1.0_QP !! Converts distance unit to centimeters
real(DP) :: GU = -1.0_DP !! Universal gravitational constant in the system units
real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units
character(STRMAX) :: ennergy_out = "" !! Name of output energy and momentum report file

!Logical flags to turn on or off various features of the code
! Logical flags to turn on or off various features of the code
logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually)
logical :: lextra_force = .false. !! User defined force function turned on
logical :: lbig_discard = .false. !! Save big bodies on every discard
Expand All @@ -56,6 +57,16 @@ module swiftest_classes
logical :: lrotation = .false. !! Include rotation states of big bodies
logical :: ltides = .false. !! Include tidal dissipation

! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values)
real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy
real(DP) :: Mtot_orig = 0.0_DP !! Initial system mass
real(DP) :: Lmag_orig = 0.0_DP !! Initial total angular momentum magnitude
real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector
real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum
real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector
logical :: lfirstenergy = .true. !! This is the first time computing energe
logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step

! Future features not implemented or in development
logical :: lgr = .false. !! Turn on GR
logical :: lyarkovsky = .false. !! Turn on Yarkovsky effect
Expand Down Expand Up @@ -259,7 +270,7 @@ module swiftest_classes
class(swiftest_pl), allocatable :: pl !! Massive body data structure
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) :: Gmtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion
real(DP) :: ke = 0.0_DP !! System kinetic energy
real(DP) :: pe = 0.0_DP !! System potential energy
real(DP) :: te = 0.0_DP !! System total energy
Expand All @@ -276,14 +287,16 @@ module swiftest_classes
procedure(abstract_step_system), deferred :: step

! Concrete classes that are common to the basic integrator (only test particles considered for discard)
procedure :: discard => discard_system !! Perform a discard step on the system
procedure :: dump => io_dump_system !! Dump the state of the system to a file
procedure :: read_frame => io_read_frame_system !! Read in a frame of input data from file
procedure :: write_discard => io_write_discard !! Write out information about discarded test particles
procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file
procedure :: initialize => setup_initialize_system !! Initialize the system from input files
procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides.
procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies.
procedure :: discard => discard_system !! Perform a discard step on the system
procedure :: conservation_report => io_conservation_report !! Compute energy and momentum and print out the change with time
procedure :: dump => io_dump_system !! Dump the state of the system to a file
procedure :: read_frame => io_read_frame_system !! Read in a frame of input data from file
procedure :: write_discard => io_write_discard !! Write out information about discarded test particles
procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file
procedure :: initialize => setup_initialize_system !! Initialize the system from input files
procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides.
procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies.
procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum
end type swiftest_nbody_system

type :: swiftest_encounter
Expand Down Expand Up @@ -487,6 +500,13 @@ module pure subroutine gr_vh2pv_body(self, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine gr_vh2pv_body

module subroutine io_conservation_report(self, param, lterminal)
implicit none
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
class(swiftest_parameters), intent(inout) :: param !! Input colleciton of user-defined parameters
logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen
end subroutine io_conservation_report

module subroutine io_dump_param(self, param_file_name)
implicit none
class(swiftest_parameters),intent(in) :: self !! Output collection of parameters
Expand Down Expand Up @@ -662,6 +682,17 @@ module subroutine obl_acc_tp(self, system)
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
end subroutine obl_acc_tp

module subroutine obl_pot(npl, Mcb, Mpl, j2rp2, j4rp4, xh, irh, oblpot)
implicit none
integer(I4B), intent(in) :: npl
real(DP), intent(in) :: Mcb
real(DP), dimension(:), intent(in) :: Mpl
real(DP), intent(in) :: j2rp2, j4rp4
real(DP), dimension(:), intent(in) :: irh
real(DP), dimension(:, :), intent(in) :: xh
real(DP), intent(out) :: oblpot
end subroutine obl_pot

module subroutine orbel_el2xv_vec(self, cb)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest body object
Expand Down Expand Up @@ -983,6 +1014,17 @@ 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)
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)
implicit none
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
Expand Down
40 changes: 40 additions & 0 deletions src/obl/obl.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,5 +106,45 @@ module subroutine obl_acc_tp(self, system)

end subroutine obl_acc_tp

module subroutine obl_pot(npl, Mcb, Mpl, j2rp2, j4rp4, xh, irh, oblpot)
!! author: David A. Minton
!!
!! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body
!! Returned value does not include monopole term or terms higher than J4
!!
!! Reference: MacMillan, W. D. 1958. The Theory of the Potential, (Dover Publications), 363.
!!
!! Adapted from David E. Kaufmann's Swifter routine: obl_pot.f90
!! Adapted from Hal Levison's Swift routine obl_pot.f
implicit none
! Arguments
integer(I4B), intent(in) :: npl
real(DP), intent(in) :: Mcb
real(DP), dimension(:), intent(in) :: Mpl
real(DP), intent(in) :: j2rp2, j4rp4
real(DP), dimension(:), intent(in) :: irh
real(DP), dimension(:, :), intent(in) :: xh
real(DP), intent(out) :: oblpot

! Internals
integer(I4B) :: i
real(DP) :: rinv2, t0, t1, t2, t3, p2, p4, mu

oblpot = 0.0_DP
mu = Mcb
do i = 1, npl
rinv2 = irh(i)**2
t0 = mu * Mpl(i) * rinv2 * irh(i)
t1 = j2rp2
t2 = xh(3, i) * xh(3, i) * rinv2
t3 = j4rp4 * rinv2
p2 = 0.5_DP * (3 * t2 - 1.0_DP)
p4 = 0.125_DP * ((35 * t2 - 30.0_DP) * t2 + 3.0_DP)
oblpot = oblpot + t0 * (t1 * p2 + t3 * p4)
end do

return
end subroutine obl_pot


end submodule s_obl
Loading

0 comments on commit fb21ff7

Please sign in to comment.