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

Commit

Permalink
Added new 'COMPACT' output format that will eventually be parsed by t…
Browse files Browse the repository at this point in the history 
…he Python side
  • Loading branch information
@daminton
daminton committed Nov 21, 2022
1 parent 6f90ed3 commit adb4ea3
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 17 deletions.
69 changes: 59 additions & 10 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,44 @@

contains

module subroutine io_compact_output(self, param, timer)
!! author: David Minton
!!
!! Generates the terminal output displayed when display_style is set to COMPACT. This is used by the Python driver to
!! make nice-looking progress reports.
implicit none
! Arguments
class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Input colleciton of user-defined parameters
class(*), intent(in) :: timer !! Object used for computing elapsed wall time (must be unlimited polymorphic because the walltimer module requires swiftest_classes)
! Internals
character(*), parameter :: COMPACTFMT = '("ILOOP",I,";T",ES23.16,";WT",ES23.16,";IWT",ES23.16,";WTPS",ES23.16,";NPL",I,";NTP",I,";"$)'
character(*), parameter :: SYMBAFMT = '(";NPLM",I,$)'
character(*), parameter :: EGYFMT = '("LTOTERR",ES24.16,";ETOTERR",ES24.16,";MTOTERR",ES24.16,";KEOERR",ES24.16,";KESERR",ES24.16,";PEERR",ES24.16' &
// '";EORBERR",ES24.16,";ECOLERR",ES24.16,";EUNTRERR",ES24.16,";LSPINERR",ES24.16,";LESCERR",ES24.16' &
// '";MESCERR",ES24.16,$)'
select type(timer)
class is (walltimer)
if (.not. timer%main_is_started) then ! This is the start of a new run
write(*,*) "START"
write(*,COMPACTFMT) 0,param%t,0.0,0.0,0.0,self%pl%nbody, self%tp%nbody
else
write(*,COMPACTFMT) param%iloop,param%t, timer%wall_main, timer%wall_step, timer%wall_per_substep,self%pl%nbody, self%tp%nbody
end if
select type(pl => self%pl)
class is (symba_pl)
write(*,SYMBAFMT) pl%nplm
end select
if (param%lenergy) then
write(*,EGYFMT) self%Ltot_error, self%Etot_error, self%Mtot_error, self%ke_orbit_error, self%ke_spin_error, self%pe_error, &
self%Eorbit_error, self%Ecoll_error, self%Euntracked_error, self%Lspin_error, self%Lescape_error, self%Mescape_error
end if
write(*,*)
end select
return
end subroutine io_compact_output


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
!!
Expand All @@ -24,19 +62,18 @@ module subroutine io_conservation_report(self, param, lterminal)
! 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) :: Eorbit_error, Etotal_error, Ecoll_error
real(DP) :: Eorbit_error, Etot_error, Ecoll_error
real(DP) :: GMtot_now
real(DP) :: Lerror, Merror
character(len=STRMAX) :: errmsg
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 &
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, display_unit => param%display_unit)

call pl%vb2vh(cb)
call pl%xh2xb(cb)

Expand All @@ -51,6 +88,9 @@ module subroutine io_conservation_report(self, param, lterminal)
GMtot_now = system%GMtot + system%GMescape

if (param%lfirstenergy) then
system%ke_orbit_orig = ke_orbit_now
system%ke_spin_orig = ke_spin_now
system%pe_orig = pe_now
system%Eorbit_orig = Eorbit_now
system%GMtot_orig = GMtot_now
system%Lorbit_orig(:) = Lorbit_now(:)
Expand All @@ -60,12 +100,21 @@ module subroutine io_conservation_report(self, param, lterminal)
end if

if (.not.param%lfirstenergy) then
Lerror = norm2(Ltot_now(:) - system%Ltot_orig(:)) / norm2(system%Ltot_orig(:))
Eorbit_error = (Eorbit_now - system%Eorbit_orig) / abs(system%Eorbit_orig)
Ecoll_error = system%Ecollisions / abs(system%Eorbit_orig)
Etotal_error = (Eorbit_now - system%Ecollisions - system%Eorbit_orig - system%Euntracked) / abs(system%Eorbit_orig)
Merror = (GMtot_now - system%GMtot_orig) / system%GMtot_orig
if (lterminal) write(display_unit, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror
system%ke_orbit_error = (ke_orbit_now - system%ke_orbit_orig) / abs(system%Eorbit_orig)
system%ke_spin_error = (ke_spin_now - system%ke_spin_orig) / abs(system%Eorbit_orig)
system%pe_error = (pe_now - system%pe_orig) / abs(system%Eorbit_orig)
system%Eorbit_error = (Eorbit_now - system%Eorbit_orig) / abs(system%Eorbit_orig)
system%Ecoll_error = system%Ecollisions / abs(system%Eorbit_orig)
system%Euntracked_error = system%Euntracked / abs(system%Eorbit_orig)
system%Etot_error = (Eorbit_now - system%Ecollisions - system%Eorbit_orig - system%Euntracked) / abs(system%Eorbit_orig)

system%Lorbit_error = norm2(Lorbit_now(:) - system%Lorbit_orig(:)) / norm2(system%Ltot_orig(:))
system%Lspin_error = norm2(Lspin_now(:) - system%Lspin_orig(:)) / norm2(system%Ltot_orig(:))
system%Lescape_error = norm2(system%Lescape(:)) / norm2(system%Ltot_orig(:))
system%Ltot_error = norm2(Ltot_now(:) - system%Ltot_orig(:)) / norm2(system%Ltot_orig(:))
system%Mescape_error = system%GMescape / system%GMtot_orig
system%Mtot_error = (GMtot_now - system%GMtot_orig) / system%GMtot_orig
if (lterminal) write(display_unit, EGYTERMFMT) system%Ltot_error, system%Ecoll_error, system%Etot_error,system%Mtot_error
if (abs(Merror) > 100 * epsilon(Merror)) then
write(*,*) "Severe error! Mass not conserved! Halting!"
! Save the frame of data to the bin file in the slot just after the present one for diagnostics
Expand Down
17 changes: 11 additions & 6 deletions src/main/swiftest_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ program swiftest_driver
character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters
character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD"
integer(I4B) :: ierr !! I/O error code
integer(I8B) :: iloop !! Loop counter
integer(I8B) :: idump !! Dump cadence counter
integer(I8B) :: iout !! Output cadence counter
integer(I8B) :: ioutput_t0 !! The output frame counter at time 0
Expand All @@ -33,13 +32,15 @@ program swiftest_driver
type(walltimer) :: integration_timer !! Object used for computing elapsed wall time
real(DP) :: tfrac
type(progress_bar) :: pbar !! Object used to print out a progress bar
character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
'"; Number of active pl, tp = ", I5, ", ", I5)'
character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
'"; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)'
character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
'"; Number of active pl, tp = ", I6, ", ", I6)'
character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
'"; Number of active plm, pl, tp = ", I6, ", ", I6, ", ", I6)'
character(*), parameter :: pbarfmt = '("Time = ", ES12.5," of ",ES12.5)'
character(len=64) :: pbarmessage

character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)'


call io_get_args(integrator, param_file_name, display_style)

Expand Down Expand Up @@ -68,6 +69,7 @@ program swiftest_driver
t0 => param%t0, &
dt => param%dt, &
tstop => param%tstop, &
iloop => param%iloop, &
istep_out => param%istep_out, &
istep_dump => param%istep_dump, &
ioutput => param%ioutput, &
Expand All @@ -93,13 +95,14 @@ program swiftest_driver
if (istep_out > 0) call nbody_system%write_frame(param)
end if


write(display_unit, *) " *************** Main Loop *************** "
if (param%lrestart .and. param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.)
if (display_style == "PROGRESS") then
call pbar%reset(nloops)
write(pbarmessage,fmt=pbarfmt) t0, tstop
call pbar%update(1,message=pbarmessage)
else if (display_style == "COMPACT") then
call nbody_system%compact_output(param,integration_timer)
end if
do iloop = 1, nloops
!> Step the system forward in time
Expand Down Expand Up @@ -136,6 +139,8 @@ program swiftest_driver
if (display_style == "PROGRESS") then
write(pbarmessage,fmt=pbarfmt) t, tstop
call pbar%update(1,message=pbarmessage)
else if (display_style == "COMPACT") then
call nbody_system%compact_output(param,integration_timer)
end if
end if
end if
Expand Down
30 changes: 29 additions & 1 deletion src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,11 +35,12 @@ module swiftest_classes
integer(I4B) :: integrator = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator used
character(STRMAX) :: param_file_name = "param.in" !! The default name of the parameter input file
integer(I4B) :: maxid = -1 !! The current maximum particle id number
integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number
integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number
real(DP) :: t0 = -1.0_DP !! Integration start time
real(DP) :: t = -1.0_DP !! Integration current time
real(DP) :: tstop = -1.0_DP !! Integration stop time
real(DP) :: dt = -1.0_DP !! Time step
integer(I8B) :: iloop = 0_I8B !! Main loop counter
integer(I8B) :: ioutput = 0_I8B !! Output counter
character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body
character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies
Expand Down Expand Up @@ -370,6 +371,9 @@ module swiftest_classes
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) :: Ltot = 0.0_DP !! System angular momentum vector
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) :: Eorbit_orig = 0.0_DP !! Initial orbital energy
real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass
real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector
Expand All @@ -379,6 +383,22 @@ module swiftest_classes
real(DP) :: GMescape = 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
real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies

! Energy, momentum, and mass errors (used in error reporting)
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) :: Eorbit_error = 0.0_DP
real(DP) :: Ecoll_error = 0.0_DP
real(DP) :: Euntracked_error = 0.0_DP
real(DP) :: Etot_error = 0.0_DP
real(DP) :: Lorbit_error = 0.0_DP
real(DP) :: Lspin_error = 0.0_DP
real(DP) :: Lescape_error = 0.0_DP
real(DP) :: Ltot_error = 0.0_DP
real(DP) :: Mtot_error = 0.0_DP
real(DP) :: Mescape_error = 0.0_DP

logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated
!! separately from massive bodies. Massive body variables are saved at half steps, and passed to
!! the test particles
Expand All @@ -388,6 +408,7 @@ module swiftest_classes

! 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 :: compact_output => io_compact_output !! Prints out out terminal output when display_style is set to COMPACT
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 :: get_old_t_final_bin => io_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output.
Expand Down Expand Up @@ -589,6 +610,13 @@ pure module subroutine gr_vh2pv_body(self, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine gr_vh2pv_body

module subroutine io_compact_output(self, param, timer)
implicit none
class(swiftest_nbody_system), intent(in) :: self !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Input colleciton of user-defined parameters
class(*), intent(in) :: timer !! Object used for computing elapsed wall time
end subroutine io_compact_output

module subroutine io_conservation_report(self, param, lterminal)
implicit none
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
Expand Down

0 comments on commit adb4ea3

Please sign in to comment.