diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index a664ef49e..914f8ed11 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -203,7 +203,7 @@ def data_stream(self, frame=0): minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades sim.set_parameter(fragmentation=True, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) - sim.run(dt=1e-4, tstop=1.0e-3, istep_out=1, dump_cadence=0) + sim.run(dt=1e-4, tstop=1.0e-3, istep_out=1, dump_cadence=4) print("Generating animation") anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) \ No newline at end of file diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 67020fb13..86a0c60ad 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -23,7 +23,6 @@ module subroutine encounter_io_dump(self, param) ! Internals integer(I4B) :: i - ! Most of this is just temporary test code just to get something working. Eventually this should get cleaned up. do i = 1, self%nframes if (allocated(self%frame(i)%item)) then diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 56ceafd9f..bfc0b38c6 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -18,7 +18,7 @@ program swiftest_driver use swiftest implicit none - class(swiftest_nbody_system), allocatable :: nbody_system !! Polymorphic object containing the nbody system to be integrated + class(swiftest_nbody_system), allocatable :: 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 swiftest_globals for symbolic names) character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters @@ -52,33 +52,21 @@ program swiftest_driver param%integrator = trim(adjustl(integrator)) call param%set_display(display_style) call param%read_in(param_file_name) - call setup_construct_system(nbody_system, param) - - !> Define the maximum number of threads - nthreads = 1 ! In the *serial* case - !$ nthreads = omp_get_max_threads() ! In the *parallel* case - !$ write(param%display_unit,'(a)') ' OpenMP parameters:' - !$ write(param%display_unit,'(a)') ' ------------------' - !$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads - !$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads - - - associate(t => nbody_system%t, & - t0 => param%t0, & - tstart => param%tstart, & - dt => param%dt, & - tstop => param%tstop, & - iloop => param%iloop, & - istep_out => param%istep_out, & - dump_cadence => param%dump_cadence, & - ioutput => param%ioutput, & - display_style => param%display_style, & - display_unit => param%display_unit) - - call nbody_system%initialize(param) + + + associate(t0 => param%t0, & + tstart => param%tstart, & + dt => param%dt, & + tstop => param%tstop, & + iloop => param%iloop, & + istep_out => param%istep_out, & + dump_cadence => param%dump_cadence, & + ioutput => param%ioutput, & + display_style => param%display_style, & + display_unit => param%display_unit) + ! Set up loop and output cadence variables - t = tstart nloops = ceiling((tstop - t0) / dt, kind=I8B) istart = ceiling((tstart - t0) / dt + 1.0_DP, kind=I8B) ioutput = max(int(istart / istep_out, kind=I4B),1) @@ -87,13 +75,27 @@ program swiftest_driver if (dump_cadence == 0) dump_cadence = ceiling(nloops / (1.0_DP * istep_out), kind=I8B) allocate(swiftest_storage(dump_cadence) :: system_history) + ! Construct the main n-body system using the user-input integrator to choose the type of system + call setup_construct_system(system, param) + + !> Define the maximum number of threads + nthreads = 1 ! In the *serial* case + !$ nthreads = omp_get_max_threads() ! In the *parallel* case + !$ write(param%display_unit,'(a)') ' OpenMP parameters:' + !$ write(param%display_unit,'(a)') ' ------------------' + !$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads + !$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads + + call system%initialize(param) + + ! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file. if (param%lrestart) then - if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) + if (param%lenergy) call system%conservation_report(param, lterminal=.true.) else - if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum - call nbody_system%write_frame(param) - call nbody_system%dump(param) + if (param%lenergy) call system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum + call system%write_frame(param) + call system%dump(param) end if write(display_unit, *) " *************** Main Loop *************** " @@ -104,21 +106,22 @@ program swiftest_driver 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) + call system%compact_output(param,integration_timer) end if iout = 0 idump = 0 + system%t = tstart do iloop = istart, nloops !> Step the system forward in time call integration_timer%start() - call nbody_system%step(param, t, dt) + call system%step(param, system%t, dt) call integration_timer%stop() - t = t0 + iloop * dt + system%t = t0 + iloop * dt !> Evaluate any discards or collisional outcomes - call nbody_system%discard(param) + call system%discard(param) if (display_style == "PROGRESS") call pbar%update(iloop) !> If the loop counter is at the output cadence value, append the data file with a single frame @@ -127,30 +130,30 @@ program swiftest_driver if (iout == istep_out) then iout = 0 idump = idump + 1 - system_history%frame(idump) = nbody_system ! Store a snapshot of the system for posterity + system_history%frame(idump) = system ! Store a snapshot of the system for posterity if (idump == dump_cadence) then idump = 0 - call nbody_system%dump(param) + call system%dump(param) call system_history%dump(param) end if - tfrac = (t - t0) / (tstop - t0) + tfrac = (system%t - t0) / (tstop - t0) - select type(pl => nbody_system%pl) + select type(pl => system%pl) class is (symba_pl) - write(display_unit, symbastatfmt) t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody + write(display_unit, symbastatfmt) system%t, tfrac, pl%nplm, pl%nbody, system%tp%nbody class default - write(display_unit, statusfmt) t, tfrac, pl%nbody, nbody_system%tp%nbody + write(display_unit, statusfmt) system%t, tfrac, pl%nbody, system%tp%nbody end select - if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) + if (param%lenergy) call 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) t, tstop + write(pbarmessage,fmt=pbarfmt) system%t, tstop call pbar%update(1,message=pbarmessage) else if (display_style == "COMPACT") then - call nbody_system%compact_output(param,integration_timer) + call system%compact_output(param,integration_timer) end if call integration_timer%reset() @@ -160,6 +163,7 @@ program swiftest_driver end do ! Dump any remaining history if it exists + call system%dump(param) call system_history%dump(param) if (display_style == "COMPACT") write(*,*) "SWIFTEST STOP" // param%integrator end associate diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index ad73dc4ad..3329fde02 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -60,6 +60,7 @@ module encounter_classes integer(I4B) :: loop_varid !! ID for the recursion level variable integer(I4B) :: time_dimsize = 0 !! Number of time values in snapshot integer(I4B) :: id_dimsize = 0 !! Number of potential id values in snapshot + integer(I4B) :: file_number = 1 !! The number to append on the output file contains procedure :: initialize => encounter_io_initialize !! Initialize a set of parameters used to identify a NetCDF output object end type encounter_io_parameters diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 4cd5d130a..13e3ec9a1 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -74,6 +74,7 @@ module subroutine setup_construct_system(system, param) if (param%lencounter_save) then if (.not. allocated(system%encounter_history)) allocate(encounter_storage :: system%encounter_history) call system%encounter_history%reset() + system%encounter_history%nc%file_number = param%iloop / param%dump_cadence end if end select end select diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 3dbf4ebbd..606d359cf 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -22,10 +22,10 @@ module subroutine symba_io_dump_encounter(self, param) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters if (self%encounter_history%iframe == 0) return ! No enounters in this interval - + self%encounter_history%nc%file_number = self%encounter_history%nc%file_number + 1 ! Create and save the output file for this encounter self%encounter_history%nc%time_dimsize = maxval(self%encounter_history%tslot(:)) - write(self%encounter_history%nc%enc_file, '("encounter_",I0.6,".nc")') param%iloop / param%dump_cadence + write(self%encounter_history%nc%enc_file, '("encounter_",I0.6,".nc")') self%encounter_history%nc%file_number call self%encounter_history%nc%initialize(param) call self%encounter_history%dump(param) call self%encounter_history%nc%close()