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

Commit

Permalink
Improved the consistency of the output stream. Still having a repeata…
Browse files Browse the repository at this point in the history 
…bility problem with restarts
  • Loading branch information
@daminton
daminton committed Jan 6, 2023
1 parent 197347d commit 2ed62ee
Show file tree
Hide file tree
Showing 4 changed files with 117 additions and 79 deletions.
1 change: 1 addition & 0 deletions src/base/base_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ module base
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) :: nloops = 0_I8B !! Total number of loops to execute
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
character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles
Expand Down
54 changes: 8 additions & 46 deletions src/swiftest/swiftest_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,30 +16,17 @@ program swiftest_driver
!! Adapted from Swifter by David E. Kaufmann's Swifter driver programs swifter_[bs,helio,ra15,rmvs,symba,tu4,whm].f90
!! Adapted from Hal Levison and Martin Duncan's Swift driver programs
use swiftest
use symba
implicit none

class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated
class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated
class(swiftest_parameters), allocatable :: param !! Run configuration parameters
character(len=:), allocatable :: integrator !! Integrator type code (see globals for symbolic names)
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(I8B) :: istart !! Starting index for loop counter
integer(I8B) :: nloops !! Number of steps to take in the simulation
integer(I4B) :: iout !! Output cadence counter
integer(I4B) :: idump !! Dump cadence counter
type(walltimer) :: integration_timer !! Object used for computing elapsed wall time
real(DP) :: tfrac !! Fraction of total simulation time completed
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 = ", 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,$)'
!type(base_storage(nframes=:)), allocatable :: system_history

call swiftest_io_get_args(integrator, param_file_name, display_style)

Expand All @@ -54,9 +41,9 @@ program swiftest_driver
dt => param%dt, &
tstop => param%tstop, &
iloop => param%iloop, &
nloops => param%nloops, &
istep_out => param%istep_out, &
dump_cadence => param%dump_cadence, &
display_style => param%display_style, &
display_unit => param%display_unit)

! Set up loop and output cadence variables
Expand Down Expand Up @@ -84,29 +71,18 @@ program swiftest_driver

associate (system_history => param%system_history)
! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file.
call nbody_system%display_run_information(param, integration_timer, phase="FIRST")
if (param%lenergy) then
if (param%lrestart) then
call nbody_system%get_t0_values(param)
call nbody_system%conservation_report(param, lterminal=.true.)
else
call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum
end if
end if
call nbody_system%conservation_report(param, lterminal=.true.)
call system_history%take_snapshot(param,nbody_system)
call nbody_system%dump(param)

write(display_unit, *) " *************** Main Loop *************** "

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
write(*,*) "SWIFTEST START " // param%integrator
call nbody_system%compact_output(param,integration_timer)
end if


do iloop = istart, nloops
!> Step the nbody_system forward in time
call integration_timer%start()
Expand All @@ -132,25 +108,10 @@ program swiftest_driver

end if

tfrac = (nbody_system%t - t0) / (tstop - t0)

select type(pl => nbody_system%pl)
class is (symba_pl)
write(display_unit, symbastatfmt) nbody_system%t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody
class default
write(display_unit, statusfmt) nbody_system%t, tfrac, pl%nbody, nbody_system%tp%nbody
end select
if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.)
call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out)

if (display_style == "PROGRESS") then
write(pbarmessage,fmt=pbarfmt) nbody_system%t, tstop
call pbar%update(1,message=pbarmessage)
else if (display_style == "COMPACT") then
call nbody_system%compact_output(param,integration_timer)
end if

call nbody_system%display_run_information(param, integration_timer)
call integration_timer%reset()
if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.)

end if
end if
Expand All @@ -159,7 +120,8 @@ program swiftest_driver
! Dump any remaining history if it exists
call nbody_system%dump(param)
call system_history%dump(param)
if (display_style == "COMPACT") write(*,*) "SWIFTEST STOP" // param%integrator
call nbody_system%display_run_information(param, integration_timer, phase="LAST")

end associate
end associate

Expand Down
106 changes: 86 additions & 20 deletions src/swiftest/swiftest_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,74 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal)
end subroutine swiftest_io_conservation_report


module subroutine swiftest_io_display_run_information(self, param, integration_timer, phase)
!! author: David A. Minton
!!
!! Displays helpful information while a run is executing. The format of the output depends on user parameters
implicit none
! Arguments
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
type(walltimer), intent(inout) :: integration_timer !! Object used for computing elapsed wall time
character(len=*), optional, intent(in) :: phase !! One of "FIRST" or "LAST"
! Internals
integer(I4B) :: phase_val
real(DP) :: tfrac !! Fraction of total simulation time completed
type(progress_bar), save :: pbar !! Object used to print out a progress bar
character(len=64) :: pbarmessage
character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)'
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)'

phase_val = 1
if (present(phase)) then
if (phase == "FIRST") then
phase_val = 0
else if (phase == "LAST") then
phase_val = -1
end if
end if

tfrac = (self%t - param%t0) / (param%tstop - param%t0)

if (phase_val == 0) then
if (param%lrestart) then
write(param%display_unit, *) " *************** Swiftest restart " // param%integrator // " *************** "
else
write(param%display_unit, *) " *************** Swiftest start " // param%integrator // " *************** "
end if
if (param%display_style == "PROGRESS") then
call pbar%reset(param%nloops)
else if (param%display_style == "COMPACT") then
write(*,*) "SWIFTEST START " // param%integrator
end if
end if

if (param%display_style == "PROGRESS") then
write(pbarmessage,fmt=pbarfmt) self%t, param%tstop
call pbar%update(1,message=pbarmessage)
else if (param%display_style == "COMPACT") then
call self%compact_output(param,integration_timer)
end if

if (self%pl%nplm > 0) then
write(param%display_unit, symbastatfmt) self%t, tfrac, self%pl%nplm, self%pl%nbody, self%tp%nbody
else
write(param%display_unit, statusfmt) self%t, tfrac, self%pl%nbody, self%tp%nbody
end if

if (phase_val == -1) then
write(param%display_unit, *)" *************** Swiftest stop " // param%integrator // " *************** "
if (param%display_style == "COMPACT") write(*,*) "SWIFTEST STOP" // param%integrator
end if

return
end subroutine swiftest_io_display_run_information


module subroutine swiftest_io_dump_param(self, param_file_name)
!! author: David A. Minton
!!
Expand Down Expand Up @@ -504,56 +572,54 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param)
real(DP), dimension(NDIM) :: rot0, Ip0, Lnow
real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, BE_orig

associate (nc => param%system_history%nc, cb => self%cb)
associate (nc => param%system_history%nc, cb => self%cb, tslot => param%system_history%nc%tslot)
call nc%open(param, readonly=.true.)
call nc%find_tslot(param%t0)
call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=itmax), "netcdf_io_get_t0_values_system time_dimid" )
call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "netcdf_io_get_t0_values_system name_dimid" )
allocate(vals(idmax))
call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_t0_values_system time_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system time_varid" )

if (param%lenergy) then
call netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_t0_values_system KE_orb_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system KE_orb_varid" )
KE_orb_orig = rtemp(1)

call netcdf_io_check( nf90_get_var(nc%id, nc%KE_spin_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_t0_values_system KE_spin_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%KE_spin_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system KE_spin_varid" )
KE_spin_orig = rtemp(1)

call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_t0_values_system PE_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system PE_varid" )
PE_orig = rtemp(1)

call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_t0_values_system BE_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system BE_varid" )
BE_orig = rtemp(1)

call netcdf_io_check( nf90_get_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[1]), "netcdf_io_get_t0_values_system E_collisions_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[1]), "netcdf_io_get_t0_values_system E_untracked_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[tslot]), "netcdf_io_get_t0_values_system E_collisions_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[tslot]), "netcdf_io_get_t0_values_system E_untracked_varid" )

self%E_orbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + BE_orig + self%E_collisions + self%E_untracked

call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_orbit_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%L_spin_varid, self%L_spin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_spin_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_escape_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit_orig(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_orbit_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%L_spin_varid, self%L_spin_orig(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_spin_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_escape_varid" )

self%L_total_orig(:) = self%L_orbit_orig(:) + self%L_spin_orig(:) + self%L_escape(:)

call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "netcdf_io_get_t0_values_system Gmass_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[1]), "netcdf_io_get_t0_values_system GMescape_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,tslot], count=[idmax,1]), "netcdf_io_get_t0_values_system Gmass_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_get_t0_values_system GMescape_varid" )
self%GMtot_orig = vals(1) + sum(vals(2:idmax), vals(2:idmax) == vals(2:idmax)) + self%GMescape

cb%GM0 = vals(1)
cb%dGM = cb%Gmass - cb%GM0

call netcdf_io_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1,1], count=[1,1]), "netcdf_io_get_t0_values_system radius_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1,tslot], count=[1,1]), "netcdf_io_get_t0_values_system radius_varid" )
cb%R0 = rtemp(1)

if (param%lrotation) then

call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system rot_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,1], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system Ip_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system rot_varid" )
call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system Ip_varid" )

cb%L0(:) = Ip0(3) * cb%GM0 * cb%R0**2 * rot0(:)

Lnow(:) = cb%Ip(3) * cb%Gmass * cb%radius**2 * cb%rot(:)
cb%dL(:) = Lnow(:) - cb%L0(:)
end if

end if
Expand Down Expand Up @@ -2198,7 +2264,7 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i
call param%set_display(param%display_style)

! Print the contents of the parameter file to standard output
call param%writer(unit = param%display_unit, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg)
if (.not.param%lrestart) call param%writer(unit = param%display_unit, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg)

end associate

Expand Down
Loading

0 comments on commit 2ed62ee

Please sign in to comment.