diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index a4f370773..f33535327 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -64,7 +64,7 @@ # handles strings differently than Python's Xarray. string_varnames = ["name", "particle_type", "status", "origin_type"] char_varnames = ["space"] -int_varnames = ["id", "ntp", "npl", "nplm", "discard_body_id", "collision_id"] +int_varnames = ["id", "ntp", "npl", "nplm", "discard_body_id", "collision_id", "loopnum"] def bool2yesno(boolval): """ @@ -816,8 +816,8 @@ def process_netcdf_input(ds, param): ------- ds : xarray dataset """ - - ds = ds.where(~np.isnan(ds.id) ,drop=True) + # + ds = ds.where(ds.id >=0,drop=True) if param['OUT_TYPE'] == "NETCDF_DOUBLE": ds = fix_types(ds,ftype=np.float64) elif param['OUT_TYPE'] == "NETCDF_FLOAT": diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 9b9bd02dc..8fb77ab4a 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -2780,6 +2780,11 @@ def _preprocess(ds, param): tgood,tid = np.unique(self.enc.time,return_index=True) self.enc = self.enc.isel(time=tid) + # Reduce the dimensionality of variables that got expanded in the combine process + self.enc['loopnum'] = self.enc['loopnum'].max(dim="name") + self.enc['id'] = self.enc['id'].max(dim="time") + self.enc['particle_type'] = self.enc['particle_type'].max(dim="time") + return diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 0cc07b009..67020fb13 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -53,7 +53,7 @@ module subroutine encounter_io_initialize(self, param) class(encounter_io_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: nvar, varid, vartype + integer(I4B) :: i, nvar, varid, vartype real(DP) :: dfill real(SP) :: sfill logical :: fileExists @@ -97,7 +97,7 @@ module subroutine encounter_io_initialize(self, param) call check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%rh_varid), "encounter_io_initialize nf90_def_var rh_varid" ) call check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type, [nc%space_dimid, nc%id_dimid, nc%time_dimid], nc%vh_varid), "encounter_io_initialize nf90_def_var vh_varid" ) call check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%Gmass_varid), "encounter_io_initialize nf90_def_var Gmass_varid" ) - call check( nf90_def_var(nc%id, nc%level_varname, NF90_INT, [nc%id_dimid, nc%time_dimid], nc%level_varid), "encounter_io_initialize nf90_def_var level_varid" ) + call check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, [nc%time_dimid], nc%loop_varid), "encounter_io_initialize nf90_def_var loop_varid" ) if (param%lclose) then call check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%id_dimid, nc%time_dimid], nc%radius_varid), "encounter_io_initialize nf90_def_var radius_varid" ) end if @@ -126,6 +126,10 @@ module subroutine encounter_io_initialize(self, param) ! Add in the space dimension coordinates call check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), "encounter_io_initialize nf90_put_var space" ) + + ! Pre-fill id slots with ids + + call check( nf90_put_var(nc%id, nc%id_varid, [(-1,i=1,param%maxid)], start=[1], count=[param%maxid]), "encounter_io_initialize nf90_put_var pl id_varid" ) end associate return @@ -145,7 +149,7 @@ module subroutine encounter_io_write_frame(self, nc, param) ! Arguments class(encounter_snapshot), intent(in) :: self !! Swiftest encounter structure class(encounter_io_parameters), intent(inout) :: nc !! Parameters used to identify a particular encounter io NetCDF dataset - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, tslot, idslot, old_mode, npl, ntp character(len=NAMELEN) :: charstring @@ -155,19 +159,16 @@ module subroutine encounter_io_write_frame(self, nc, param) call check( nf90_set_fill(nc%id, nf90_nofill, old_mode), "encounter_io_write_frame nf90_set_fill" ) call check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), "encounter_io_write_frame nf90_put_var time_varid" ) + call check( nf90_put_var(nc%id, nc%loop_varid, int(self%iloop,kind=I4B), start=[tslot]), "encounter_io_write_frame nf90_put_var pl loop_varid" ) associate(pl => self%pl, tp => self%tp) npl = pl%nbody do i = 1, npl idslot = pl%id(i) - call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var pl id_varid" ) + call check( nf90_put_var(nc%id, nc%id_varid, pl%id(i), start=[idslot]), "encounter_io_write_frame nf90_put_var pl id_varid" ) call check( nf90_put_var(nc%id, nc%rh_varid, pl%rh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl rh_varid" ) call check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1,idslot,tslot], count=[NDIM,1,1]), "encounter_io_write_frame nf90_put_var pl vh_varid" ) call check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl Gmass_varid" ) - select type(pl) - class is (symba_pl) - call check( nf90_put_var(nc%id, nc%level_varid, pl%levelg(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl level_varid" ) - end select if (param%lclose) call check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[idslot, tslot]), "encounter_io_write_frame nf90_put_var pl radius_varid" ) diff --git a/src/io/io.f90 b/src/io/io.f90 index d14f0a694..f0c1cb994 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -264,6 +264,17 @@ module subroutine io_dump_system(self, param) idx = idx + 1 if (idx > NDUMPFILES) idx = 1 + ! Dump the encounter history if necessary + select type(param) + class is (symba_parameters) + if (param%lencounter_save) then + select type(self) + class is (symba_nbody_system) + call self%dump_encounter(param) + end select + end if + end select + return end subroutine io_dump_system diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 6d5abce79..56ceafd9f 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -51,6 +51,8 @@ program swiftest_driver end select 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 @@ -60,8 +62,6 @@ program swiftest_driver !$ write(param%display_unit,'(a,i3,/)') ' Number of threads = ', nthreads !$ if (param%log_output) write(*,'(a,i3)') ' OpenMP: Number of threads = ',nthreads - call setup_construct_system(nbody_system, param) - call param%read_in(param_file_name) associate(t => nbody_system%t, & t0 => param%t0, & diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index c7edc5e8a..ad73dc4ad 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -43,9 +43,10 @@ module encounter_classes type :: encounter_snapshot !! A simplified version of a SyMBA nbody system object for storing minimal snapshots of the system state during encounters - class(swiftest_pl), allocatable :: pl !! Massive body data structure - class(swiftest_tp), allocatable :: tp !! Test particle data structure - real(DP) :: t !! Simulation time when snapshot was taken + class(swiftest_pl), allocatable :: pl !! Massive body data structure + class(swiftest_tp), allocatable :: tp !! Test particle data structure + real(DP) :: t !! Simulation time when snapshot was taken + integer(I8B) :: iloop !! Loop number at time of snapshot contains procedure :: write_frame => encounter_io_write_frame !! Writes a frame of encounter data to file final :: encounter_util_final_snapshot @@ -53,19 +54,19 @@ module encounter_classes !> NetCDF dimension and variable names for the enounter save object type, extends(netcdf_parameters) :: encounter_io_parameters - integer(I4B) :: ienc_frame = 1 !! Current frame number for the encounter history - character(STRMAX) :: enc_file !! Encounter output file name - character(NAMELEN) :: level_varname = "level" !! Recursion depth - integer(I4B) :: level_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) :: ienc_frame = 1 !! Current frame number for the encounter history + character(STRMAX) :: enc_file !! Encounter output file name + character(NAMELEN) :: loop_varname = "loopnum" !! Loop number for encounter + 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 contains procedure :: initialize => encounter_io_initialize !! Initialize a set of parameters used to identify a NetCDF output object end type encounter_io_parameters !> A class that that is used to store simulation history data between file output type, extends(swiftest_storage) :: encounter_storage - type(encounter_io_parameters) :: nc !! NetCDF parameter object containing the details about the file attached to this storage object + type(encounter_io_parameters) :: nc !! NetCDF parameter object containing the details about the file attached to this storage object contains procedure :: dump => encounter_io_dump !! Dumps contents of encounter history to file final :: encounter_util_final_storage diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index faa1e2e80..4c56ff6f7 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -199,8 +199,7 @@ module symba_classes procedure :: recursive_step => symba_step_recur_system !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current recursion level, if applicable, and descend to the next deeper level if necessary procedure :: reset => symba_step_reset_system !! Resets pl, tp,and encounter structures at the start of a new step procedure :: snapshot => symba_util_take_encounter_snapshot !! Take a minimal snapshot of the system through an encounter - procedure :: start_encounter => symba_io_start_encounter !! Initializes the new encounter history - procedure :: stop_encounter => symba_io_stop_encounter !! Saves the encounter and/or fragmentation data to file(s) + procedure :: dump_encounter => symba_io_dump_encounter !! Saves the encounter and/or fragmentation data to file(s) final :: symba_util_final_system !! Finalizes the SyMBA nbody system object - deallocates all allocatables end type symba_nbody_system @@ -409,19 +408,11 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 end subroutine symba_io_param_writer - module subroutine symba_io_start_encounter(self, param, t) + module subroutine symba_io_dump_encounter(self, param) implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(symba_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - end subroutine symba_io_start_encounter - - module subroutine symba_io_stop_encounter(self, param, t) - implicit none - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - end subroutine symba_io_stop_encounter + end subroutine symba_io_dump_encounter module subroutine symba_io_write_discard(self, param) use swiftest_classes, only : swiftest_parameters diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index e95505b9b..4cd5d130a 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -68,6 +68,14 @@ module subroutine setup_construct_system(system, param) allocate(symba_pltpenc :: system%pltpenc_list) allocate(symba_plplenc :: system%plplenc_list) allocate(symba_plplenc :: system%plplcollision_list) + + select type(param) + class is (symba_parameters) + if (param%lencounter_save) then + if (.not. allocated(system%encounter_history)) allocate(encounter_storage :: system%encounter_history) + call system%encounter_history%reset() + end if + end select end select case (RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index e790e1685..3dbf4ebbd 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -11,6 +11,30 @@ use swiftest contains + + module subroutine symba_io_dump_encounter(self, param) + !! author: David A. Minton + !! + !! Saves the encounter and/or fragmentation data to file with the name of the current output interval number attached + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters + + if (self%encounter_history%iframe == 0) return ! No enounters in this interval + + ! 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 + call self%encounter_history%nc%initialize(param) + call self%encounter_history%dump(param) + call self%encounter_history%nc%close() + call self%encounter_history%reset() + + return + end subroutine symba_io_dump_encounter + + module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! @@ -188,55 +212,6 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms end subroutine symba_io_param_writer - module subroutine symba_io_start_encounter(self, param, t) - !! author: David A. Minton - !! - !! Initializes the new encounter and/or fragmentation history - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - - if (.not. allocated(self%encounter_history)) then - allocate(encounter_storage :: self%encounter_history) - end if - call self%encounter_history%reset() - - ! Take the snapshot at the start of the encounter - call self%snapshot(param, t) - - return - end subroutine symba_io_start_encounter - - - module subroutine symba_io_stop_encounter(self, param, t) - !! author: David A. Minton - !! - !! Saves the encounter and/or fragmentation data to file(s) - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - ! Internals - integer(I4B) :: i - - ! 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 - - call self%encounter_history%nc%initialize(param) - call self%encounter_history%dump(param) - call self%encounter_history%nc%close() - call self%encounter_history%reset() - - return - end subroutine symba_io_stop_encounter - - module subroutine symba_io_write_discard(self, param) !! author: David A. Minton !! diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index e727ed9f3..0b8879464 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -37,9 +37,9 @@ module subroutine symba_step_system(self, param, t, dt) call self%reset(param) lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0) if (lencounter) then - if (param%lencounter_save) call self%start_encounter(param, t) + if (param%lencounter_save) call self%snapshot(param, t) call self%interp(param, t, dt) - if (param%lencounter_save) call self%stop_encounter(param, t+dt) + if (param%lencounter_save) call self%snapshot(param, t+dt) else self%irec = -1 call helio_step_system(self, param, t, dt) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index d53bd442b..c2a837599 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -1302,6 +1302,7 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t) associate(npl => self%pl%nbody, ntp => self%tp%nbody) snapshot%t = t + snapshot%iloop = param%iloop if (npl + ntp == 0) return npl_snap = npl