diff --git a/src/io/io.f90 b/src/io/io.f90 index e93cc265e..974c6aae0 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -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 !! @@ -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) @@ -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(:) @@ -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 diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 4cd537632..e557cbdfa 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -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 @@ -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) @@ -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, & @@ -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 @@ -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 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 78b82d43b..26686c1d6 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -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 @@ -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 @@ -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 @@ -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. @@ -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