diff --git a/CMakeLists.txt b/CMakeLists.txt index 474451c8a..87b8789d1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -338,6 +338,7 @@ list(APPEND libopenmc_SOURCES src/cell.cpp src/chain.cpp src/cmfd_solver.cpp + src/collision_track.cpp src/cross_sections.cpp src/dagmc.cpp src/distribution.cpp diff --git a/docs/source/io_formats/collision_track.rst b/docs/source/io_formats/collision_track.rst new file mode 100644 index 000000000..a36459754 --- /dev/null +++ b/docs/source/io_formats/collision_track.rst @@ -0,0 +1,46 @@ +.. _io_collision_track: + +=========================== +Collision Track File Format +=========================== + +When collision tracking is enabled with ``mcpl=false`` (the default), OpenMC +writes binary data to an HDF5 file named ``collision_track.h5``. The same data +may also be written after each batch when multiple files are requested +(``collision_track.N.h5``) or when the run is performed in parallel. The file +contains the information needed to reconstruct each recorded collision. + +The current revision of the collision track file format is 1.0. + +**/** + +:Attributes: + - **filetype** (*char[]*) -- String indicating the type of file. + For collision-track files the value is ``"collision_track"``. + +:Datasets: + + - **collision_track_bank** (Compound type) -- Collision information + for each stored event. Each entry in the dataset corresponds to one + collision and contains the following fields: + + - ``r`` (*double[3]*) -- Position of the collision in [cm]. + - ``u`` (*double[3]*) -- Direction unit vector immediately after the collision. + - ``E`` (*double*) -- Incident particle energy before the collision in [eV]. + - ``dE`` (*double*) -- Energy loss over the collision (:math:`E_\text{before} - E_\text{after}`) in [eV]. + - ``time`` (*double*) -- Time of the collision in [s]. + - ``wgt`` (*double*) -- Particle weight at the collision. + - ``event_mt`` (*int*) -- ENDF MT number identifying the reaction. + - ``delayed_group`` (*int*) -- Delayed neutron group index (non-zero for delayed events). + - ``cell_id`` (*int*) -- ID of the cell in which the collision occurred. + - ``nuclide_id`` (*int*) -- ZA identifier of the nuclide (ZZZAAAM format). + - ``material_id`` (*int*) -- ID of the material containing the collision site. + - ``universe_id`` (*int*) -- ID of the universe containing the collision site. + - ``n_collision`` (*int*) -- Collision counter for the particle history. + - ``particle`` (*int*) -- Particle type (0=neutron, 1=photon, 2=electron, 3=positron). + - ``parent_id`` (*int64*) -- Unique ID of the parent particle. + - ``progeny_id`` (*int64*) -- Progeny ID of the particle. + +In an MPI run, OpenMC writes the combined dataset by gathering collision-track +entries from all ranks before flushing them to disk, so the final file appears +as though it were produced serially. diff --git a/docs/source/io_formats/index.rst b/docs/source/io_formats/index.rst index 4bbaa961a..5b4efea66 100644 --- a/docs/source/io_formats/index.rst +++ b/docs/source/io_formats/index.rst @@ -44,6 +44,7 @@ Output Files statepoint source + collision_track summary properties depletion_results diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 720846c85..b7874fcf2 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -20,6 +20,85 @@ source neutrons. *Default*: None +----------------------------- +```` Element +----------------------------- + +The ```` element indicates to track information about particle +collisions based on a set of criteria and store these events in a file named +``collision_track.h5``. This file records details such as the position of the +interaction, direction of the incoming particle, incident energy and deposited +energy, weight, time of the interaction, and the delayed neutron group (0 for +prompt neutrons). Additional information such as the cell ID, material ID, +universe ID, nuclide ZAID, particle type, and event MT number are also stored. +Users can specify one or more criterion to filter collisions. If no criteria are +specified, it defaults to tracking all collisions across the model. + +.. warning:: + Storing all collisions can be very memory intensive. For more targeted + tracking, users can employ a variety of parameters such as ``cell_ids``, + ``reactions``, ``universe_ids``, ``material_ids``, ``nuclides``, and + ``deposited_E_threshold`` to refine the selection of particle interactions + to be banked. + +This element can contain one or more of the following attributes or +sub-elements: + + :max_collisions: + An integer indicating the maximum number of collisions to be banked per file. + + *Default*: 1000 + + :max_collision_track_files: + An integer indicating the number of collision_track files to be used. + + *Default*: 1 + + :mcpl: + An optional boolean to enable MCPL_-format instead of the native HDF5-based + format. If activated, the output file name and type is changed to + ``collision_track.mcpl``. + + *Default*: false + + .. _MCPL: https://mctools.github.io/mcpl/mcpl.pdf + + :cell_ids: + A list of integers representing cell IDs to define specific cells in which + collisions are to be banked. + + *Default*: None + + :universe_ids: + A list of integers representing the universe IDs to define specific + universes in which collisions are to be banked. + + *Default*: None + + :material_ids: + A list of integers representing the material IDs to define specific + materials in which collisions are to be banked. + + *Default*: None + + :nuclides: + A list of strings representing the nuclide, to define specific + define specific target nuclide collisions to be banked. + + *Default*: None + + :reactions: + A list of integers representing the ENDF-6 format MT numbers or strings + (e.g. (n,fission)) to define specific reaction types to be banked. + + *Default*: None + + :deposited_E_threshold: + A float defining the minimum deposited energy per collision (in eV) to + trigger banking. + + *Default*: 0.0 + ---------------------------------- ```` Element ---------------------------------- diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 2a9d0876c..ce2f6f0f8 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -216,6 +216,9 @@ Post-processing :nosignatures: :template: myfunction.rst + openmc.read_collision_track_file + openmc.read_collision_track_hdf5 + openmc.read_collision_track_mcpl openmc.voxel_to_vtk The following classes and functions are used for functional expansion reconstruction. diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 1b2d4bc1a..f973f1465 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -756,6 +756,62 @@ instance, whereas the :meth:`openmc.Track.filter` method returns a new track_files = [f"tracks_p{rank}.h5" for rank in range(32)] openmc.Tracks.combine(track_files, "tracks.h5") +Collision Track File +--------------------- + +OpenMC can generate a collision track file that contains detailed collision +information (position, direction, energy, deposited energy, time, weight, cell +ID, material ID, universe ID, nuclide ZAID, particle type, particle delayed +group and particle ID) for each particle collision depending on user-defined +parameters. To invoke this feature, set the +:attr:`~openmc.Settings.collision_track` attribute as shown in this example:: + + settings.collision_track = { + "max_collisions": 300, + "reactions": ["(n,fission)", "(n,2n)"], + "material_ids": [1,2], + "nuclides": ["U238", "O16"], + "cell_ids": [5, 12] + } + +In this example, collision track information is written to the +collision_track.h5 file at the end of the simulation. The file contains +300 recorded collisions that occurred in materials with IDs 1 or 2, involving +fission or (n,2n) reactions on the nuclides U-238 or O-16, within cells +with IDs 5 and 12. +The file can be read using :func:`openmc.read_collision_track_file`. +The example below shows how to extract the data from the collision_track +feature and displays the fields stored in the file: + +>>> data = openmc.read_collision_track_file('collision_track.h5') +>>> data.dtype + dtype([('r', [('x', ' source_bank; extern SharedArray surf_source_bank; +extern SharedArray collision_track_bank; + extern SharedArray fission_bank; extern vector> ifp_source_delayed_group_bank; diff --git a/include/openmc/bank_io.h b/include/openmc/bank_io.h new file mode 100644 index 000000000..90ffd820f --- /dev/null +++ b/include/openmc/bank_io.h @@ -0,0 +1,103 @@ +#ifndef OPENMC_BANK_IO_H +#define OPENMC_BANK_IO_H + +#include "hdf5.h" + +#include "openmc/message_passing.h" +#include "openmc/span.h" +#include "openmc/vector.h" + +#include + +#ifdef OPENMC_MPI +#include +#endif + +namespace openmc { + +template +void write_bank_dataset(const char* dataset_name, hid_t group_id, + span bank, const vector& bank_index, hid_t banktype +#ifdef OPENMC_MPI + , + MPI_Datatype mpi_dtype +#endif +) +{ + int64_t dims_size = bank_index.back(); + int64_t count_size = bank_index[mpi::rank + 1] - bank_index[mpi::rank]; + +#ifdef PHDF5 + hsize_t dims[] {static_cast(dims_size)}; + hid_t dspace = H5Screate_simple(1, dims, nullptr); + hid_t dset = H5Dcreate(group_id, dataset_name, banktype, dspace, H5P_DEFAULT, + H5P_DEFAULT, H5P_DEFAULT); + + hsize_t count[] {static_cast(count_size)}; + hid_t memspace = H5Screate_simple(1, count, nullptr); + + hsize_t start[] {static_cast(bank_index[mpi::rank])}; + H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr); + + hid_t plist = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE); + + H5Dwrite(dset, banktype, memspace, dspace, plist, bank.data()); + + H5Sclose(dspace); + H5Sclose(memspace); + H5Dclose(dset); + H5Pclose(plist); +#else + if (mpi::master) { + hsize_t dims[] {static_cast(dims_size)}; + hid_t dspace = H5Screate_simple(1, dims, nullptr); + hid_t dset = H5Dcreate(group_id, dataset_name, banktype, dspace, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + +#ifdef OPENMC_MPI + vector temp_bank {bank.begin(), bank.end()}; +#endif + + for (int i = 0; i < mpi::n_procs; ++i) { + hsize_t count[] {static_cast(bank_index[i + 1] - bank_index[i])}; + hid_t memspace = H5Screate_simple(1, count, nullptr); + +#ifdef OPENMC_MPI + if (i > 0) { + MPI_Recv(bank.data(), count[0], mpi_dtype, i, i, mpi::intracomm, + MPI_STATUS_IGNORE); + } +#endif + + hid_t dspace_rank = H5Dget_space(dset); + hsize_t start[] {static_cast(bank_index[i])}; + H5Sselect_hyperslab( + dspace_rank, H5S_SELECT_SET, start, nullptr, count, nullptr); + + H5Dwrite(dset, banktype, memspace, dspace_rank, H5P_DEFAULT, bank.data()); + + H5Sclose(memspace); + H5Sclose(dspace_rank); + } + + H5Dclose(dset); + +#ifdef OPENMC_MPI + std::copy(temp_bank.begin(), temp_bank.end(), bank.begin()); +#endif + } +#ifdef OPENMC_MPI + else { + if (!bank.empty()) { + MPI_Send( + bank.data(), bank.size(), mpi_dtype, 0, mpi::rank, mpi::intracomm); + } + } +#endif +#endif +} + +} // namespace openmc + +#endif // OPENMC_BANK_IO_H diff --git a/include/openmc/collision_track.h b/include/openmc/collision_track.h new file mode 100644 index 000000000..208b8f6e1 --- /dev/null +++ b/include/openmc/collision_track.h @@ -0,0 +1,23 @@ +#ifndef OPENMC_COLLISION_TRACK_H +#define OPENMC_COLLISION_TRACK_H + +#include + +namespace openmc { + +class Particle; + +//! Reserve space in the collision track bank according to user settings. +void collision_track_reserve_bank(); + +//! Write collision track data to disk when the bank is full or the batch ends. +void collision_track_flush_bank(); + +//! Record the current particle as a collision-track entry when applicable. +//! +//! \param particle Particle whose collision should be recorded if eligible +void collision_track_record(Particle& particle); + +} // namespace openmc + +#endif // OPENMC_COLLISION_TRACK_H diff --git a/include/openmc/constants.h b/include/openmc/constants.h index b66193481..bba260c3a 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -34,6 +34,7 @@ constexpr array VERSION_VOXEL {2, 0}; constexpr array VERSION_MGXS_LIBRARY {1, 0}; constexpr array VERSION_PROPERTIES {1, 1}; constexpr array VERSION_WEIGHT_WINDOWS {1, 0}; +constexpr array VERSION_COLLISION_TRACK {1, 0}; // ============================================================================ // ADJUSTABLE PARAMETERS diff --git a/include/openmc/mcpl_interface.h b/include/openmc/mcpl_interface.h index a76d72e64..a9cce3e69 100644 --- a/include/openmc/mcpl_interface.h +++ b/include/openmc/mcpl_interface.h @@ -38,6 +38,21 @@ vector mcpl_source_sites(std::string path); void write_mcpl_source_point(const char* filename, span source_bank, const vector& bank_index); +//! Write an MCPL collision track file +//! +//! This function writes collision track data to an MCPL file. Additional +//! collision-specific metadata (such as energy deposition, material info, etc.) +//! is stored in the file header as blob data. +//! +//! \param[in] filename Path to MCPL file +//! \param[in] collision_track_bank Vector of CollisionTrackSites to write to +//! file for this MPI rank. +//! \param[in] bank_index Pointer to vector of site index ranges over all +//! MPI ranks. +void write_mcpl_collision_track(const char* filename, + span collision_track_bank, + const vector& bank_index); + //! Check if MCPL functionality is available bool is_mcpl_interface_available(); diff --git a/include/openmc/message_passing.h b/include/openmc/message_passing.h index a1641a906..ce993776b 100644 --- a/include/openmc/message_passing.h +++ b/include/openmc/message_passing.h @@ -18,6 +18,7 @@ extern bool master; #ifdef OPENMC_MPI extern MPI_Datatype source_site; +extern MPI_Datatype collision_track_site; extern MPI_Comm intracomm; #endif diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 5383487d4..fdacfa765 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -56,6 +56,25 @@ struct SourceSite { int64_t progeny_id; }; +struct CollisionTrackSite { + Position r; + Direction u; + double E; + double dE; + double time {0.0}; + double wgt {1.0}; + int event_mt {0}; + int delayed_group {0}; + int cell_id {0}; + int nuclide_id; + int material_id {0}; + int universe_id {0}; + int n_collision {0}; + ParticleType particle; + int64_t parent_id; + int64_t progeny_id; +}; + //! State of a particle used for particle track files struct TrackState { Position r; //!< Position in [cm] diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 78bfa088e..b369c99fe 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -32,6 +32,24 @@ enum class IFPParameter { GenerationTime, }; +struct CollisionTrackConfig { + bool mcpl_write {false}; //!< Write collision tracks using MCPL? + std::unordered_set + cell_ids; //!< Cell ids where collisions will be written + std::unordered_set + mt_numbers; //!< MT Numbers where collisions will be written + std::unordered_set + universe_ids; //!< Universe IDs where collisions will be written + std::unordered_set + material_ids; //!< Material IDs where collisions will be written + std::unordered_set + nuclides; //!< Nuclides where collisions will be written + double deposited_energy_threshold {0.0}; //!< Minimum deposited energy [eV] + int64_t max_collisions { + 1000}; //!< Maximum events recorded per collision track file + int64_t max_files {1}; //!< Maximum number of collision track files +}; + //============================================================================== // Global variable declarations //============================================================================== @@ -41,6 +59,7 @@ namespace settings { // Boolean flags extern bool assume_separate; //!< assume tallies are spatially separate? extern bool check_overlaps; //!< check overlaps in geometry? +extern bool collision_track; //!< flag to use collision track feature? extern bool confidence_intervals; //!< use confidence intervals for results? extern bool create_fission_neutrons; //!< create fission neutrons (fixed source)? @@ -145,6 +164,7 @@ extern std::unordered_set statepoint_batch; //!< Batches when state should be written extern std::unordered_set source_write_surf_id; //!< Surface ids where sources will be written +extern CollisionTrackConfig collision_track_config; extern double source_rejection_fraction; //!< Minimum fraction of source sites //!< that must be accepted extern double free_gas_threshold; //!< Threshold multiplier for free gas diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 3e4e24e1d..9a6cf1b21 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -22,6 +22,7 @@ constexpr int STATUS_EXIT_ON_TRIGGER {2}; namespace simulation { +extern int ct_current_file; //!< current collision track file index extern "C" int current_batch; //!< current batch extern "C" int current_gen; //!< current fission generation extern "C" bool initialized; //!< has simulation been initialized? diff --git a/include/openmc/urr.h b/include/openmc/urr.h index 3978c7b86..1e6037158 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -27,10 +27,10 @@ public: double heating; }; - Interpolation interp_; //!< interpolation type - int inelastic_flag_; //!< inelastic competition flag - int absorption_flag_; //!< other absorption flag - bool multiply_smooth_; //!< multiply by smooth cross section? + Interpolation interp_; //!< interpolation type + int inelastic_flag_; //!< inelastic competition flag + int absorption_flag_; //!< other absorption flag + bool multiply_smooth_; //!< multiply by smooth cross section? vector energy_; //!< incident energies auto n_energy() const { return energy_.size(); } diff --git a/openmc/settings.py b/openmc/settings.py index 2f8a2b124..43c1fe069 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -6,7 +6,7 @@ from numbers import Integral, Real from pathlib import Path import lxml.etree as ET - +import warnings import openmc import openmc.checkvalue as cv from openmc.checkvalue import PathLike @@ -47,6 +47,21 @@ class Settings: half-width of the 95% two-sided confidence interval. If False, uncertainties on tally results will be reported as the sample standard deviation. + collision_track : dict + Options for writing collision information. Acceptable keys are: + + :max_collisions: Maximum number of collisions to be banked per file. (int) + :max_collision_track_files: Maximum number of collision_track files. (int) + :mcpl: Output in the form of an MCPL-file. (bool) + :cell_ids: List of cell IDs to define cells in which collisions should be banked. (list of int) + :universe_ids: List of universe IDs to define universes in which collisions should be banked. (list of int) + :material_ids: List of material IDs to define materials in which collisions should be banked. (list of int) + :nuclides: List of nuclides to define nuclides in which collisions should be banked. + (ex: ["I135m", "U233"] ). (list of str) + :reactions: List of reaction to define specific reactions that should be banked + (ex: ["(n,fission)", 2, "(n,2n)"] ). (list of str or int) + :deposited_E_threshold: Number to define the minimum deposited energy during + per collision to trigger banking. (float) create_fission_neutrons : bool Indicate whether fission neutrons should be created or not. cutoff : dict @@ -395,6 +410,9 @@ class Settings: # Iterated Fission Probability self._ifp_n_generation = None + # Collision track feature + self._collision_track = {} + # Output options self._statepoint = {} self._sourcepoint = {} @@ -434,7 +452,8 @@ class Settings: self._max_particle_events = None self._write_initial_source = None self._weight_windows = WeightWindowsList() - self._weight_window_generators = cv.CheckedList(WeightWindowGenerator, 'weight window generators') + self._weight_window_generators = cv.CheckedList( + WeightWindowGenerator, 'weight window generators') self._weight_windows_on = None self._weight_windows_file = None self._weight_window_checkpoints = {} @@ -475,8 +494,9 @@ class Settings: @generations_per_batch.setter def generations_per_batch(self, generations_per_batch: int): - cv.check_type('generations per patch', generations_per_batch, Integral) - cv.check_greater_than('generations per batch', generations_per_batch, 0) + cv.check_type('generations per batch', generations_per_batch, Integral) + cv.check_greater_than('generations per batch', + generations_per_batch, 0) self._generations_per_batch = generations_per_batch @property @@ -506,7 +526,8 @@ class Settings: @rel_max_lost_particles.setter def rel_max_lost_particles(self, rel_max_lost_particles: float): cv.check_type('rel_max_lost_particles', rel_max_lost_particles, Real) - cv.check_greater_than('rel_max_lost_particles', rel_max_lost_particles, 0) + cv.check_greater_than('rel_max_lost_particles', + rel_max_lost_particles, 0) cv.check_less_than('rel_max_lost_particles', rel_max_lost_particles, 1) self._rel_max_lost_particles = rel_max_lost_particles @@ -516,8 +537,10 @@ class Settings: @max_write_lost_particles.setter def max_write_lost_particles(self, max_write_lost_particles: int): - cv.check_type('max_write_lost_particles', max_write_lost_particles, Integral) - cv.check_greater_than('max_write_lost_particles', max_write_lost_particles, 0) + cv.check_type('max_write_lost_particles', + max_write_lost_particles, Integral) + cv.check_greater_than('max_write_lost_particles', + max_write_lost_particles, 0) self._max_write_lost_particles = max_write_lost_particles @property @@ -538,12 +561,12 @@ class Settings: def keff_trigger(self, keff_trigger: dict): if not isinstance(keff_trigger, dict): msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \ - 'which is not a Python dictionary' + 'which is not a Python dictionary' raise ValueError(msg) elif 'type' not in keff_trigger: msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \ - 'which does not have a "type" key' + 'which does not have a "type" key' raise ValueError(msg) elif keff_trigger['type'] not in ['variance', 'std_dev', 'rel_err']: @@ -553,7 +576,7 @@ class Settings: elif 'threshold' not in keff_trigger: msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \ - 'which does not have a "threshold" key' + 'which does not have a "threshold" key' raise ValueError(msg) elif not isinstance(keff_trigger['threshold'], Real): @@ -570,7 +593,7 @@ class Settings: @energy_mode.setter def energy_mode(self, energy_mode: str): cv.check_value('energy mode', energy_mode, - ['continuous-energy', 'multi-group']) + ['continuous-energy', 'multi-group']) self._energy_mode = energy_mode @property @@ -593,7 +616,8 @@ class Settings: def source(self, source: SourceBase | Iterable[SourceBase]): if not isinstance(source, MutableSequence): source = [source] - self._source = cv.CheckedList(SourceBase, 'source distributions', source) + self._source = cv.CheckedList( + SourceBase, 'source distributions', source) @property def confidence_intervals(self) -> bool: @@ -610,7 +634,8 @@ class Settings: @electron_treatment.setter def electron_treatment(self, electron_treatment: str): - cv.check_value('electron treatment', electron_treatment, ['led', 'ttb']) + cv.check_value('electron treatment', + electron_treatment, ['led', 'ttb']) self._electron_treatment = electron_treatment @property @@ -704,7 +729,8 @@ class Settings: @trigger_max_batches.setter def trigger_max_batches(self, trigger_max_batches: int): cv.check_type('trigger maximum batches', trigger_max_batches, Integral) - cv.check_greater_than('trigger maximum batches', trigger_max_batches, 0) + cv.check_greater_than('trigger maximum batches', + trigger_max_batches, 0) self._trigger_max_batches = trigger_max_batches @property @@ -713,8 +739,10 @@ class Settings: @trigger_batch_interval.setter def trigger_batch_interval(self, trigger_batch_interval: int): - cv.check_type('trigger batch interval', trigger_batch_interval, Integral) - cv.check_greater_than('trigger batch interval', trigger_batch_interval, 0) + cv.check_type('trigger batch interval', + trigger_batch_interval, Integral) + cv.check_greater_than('trigger batch interval', + trigger_batch_interval, 0) self._trigger_batch_interval = trigger_batch_interval @property @@ -798,19 +826,22 @@ class Settings: @surf_source_write.setter def surf_source_write(self, surf_source_write: dict): - cv.check_type("surface source writing options", surf_source_write, Mapping) + cv.check_type("surface source writing options", + surf_source_write, Mapping) for key, value in surf_source_write.items(): cv.check_value( "surface source writing key", key, - ("surface_ids", "max_particles", "max_source_files", "mcpl", "cell", "cellfrom", "cellto"), + ("surface_ids", "max_particles", "max_source_files", + "mcpl", "cell", "cellfrom", "cellto"), ) if key == "surface_ids": cv.check_type( "surface ids for source banking", value, Iterable, Integral ) for surf_id in value: - cv.check_greater_than("surface id for source banking", surf_id, 0) + cv.check_greater_than( + "surface id for source banking", surf_id, 0) elif key == "mcpl": cv.check_type("write to an MCPL-format file", value, bool) @@ -827,6 +858,79 @@ class Settings: self._surf_source_write = surf_source_write + @property + def collision_track(self) -> dict: + return self._collision_track + + @collision_track.setter + def collision_track(self, collision_track: dict): + cv.check_type('Collision tracking options', collision_track, Mapping) + for key, value in collision_track.items(): + cv.check_value('collision_track key', key, + ('cell_ids', 'reactions', 'universe_ids', 'material_ids', 'nuclides', + 'deposited_E_threshold', 'max_collisions', 'max_collision_track_files', 'mcpl')) + if key == 'cell_ids': + cv.check_type('cell ids for collision tracking data banking', value, + Iterable, Integral) + for cell_id in value: + cv.check_greater_than('cell id for collision tracking data banking', + cell_id, 0) + elif key == 'reactions': + cv.check_type('MT numbers for collision tracking data banking', value, + Iterable) + for reaction in value: + if isinstance(reaction, int): + cv.check_greater_than( + 'MT number for collision tracking data banking', reaction, 0 + ) + elif isinstance(reaction, str): + # check against allowed strings? so far let C++ code handle it + pass + else: + raise TypeError( + f"MT number for collision tracking data banking must be a positive int or string, " + f"got {type(reaction).__name__}") + elif key == 'universe_ids': + cv.check_type('universe ids for collision tracking data banking', value, + Iterable, Integral) + for universe_id in value: + cv.check_greater_than('universe id for collision tracking data banking', + universe_id, 0) + elif key == 'material_ids': + cv.check_type('material ids for collision tracking data banking', value, + Iterable, Integral) + for material_id in value: + cv.check_greater_than('material id for collision tracking data banking', + material_id, 0) + elif key == 'nuclides': + cv.check_type('nuclides for collision tracking data banking', value, + Iterable, str) + for nuclide in value: + # If nuclide name doesn't look valid, give a warning + try: + openmc.data.zam(nuclide) + except ValueError: + warnings.warn(f"Nuclide {nuclide} is not valid") + elif key == 'deposited_E_threshold': + cv.check_type('Deposited Energy Threshold for collision tracking data banking', + value, Real) + cv.check_greater_than('Deposited Energy Threshold for collision tracking data banking', + value, 0) + elif key == 'max_collisions': + cv.check_type('maximum collisions banks per file', + value, Integral) + cv.check_greater_than('maximum collisions banks in collision tracking', + value, 0) + elif key == 'max_collision_track_files': + cv.check_type('maximum collisions banks', + value, Integral) + cv.check_greater_than('maximum number of collision_track files ', + value, 0) + elif key == 'mcpl': + cv.check_type('write to an MCPL-format file', value, bool) + + self._collision_track = collision_track + @property def no_reduce(self) -> bool: return self._no_reduce @@ -943,7 +1047,7 @@ class Settings: def cutoff(self, cutoff: dict): if not isinstance(cutoff, Mapping): msg = f'Unable to set cutoff from "{cutoff}" which is not a '\ - 'Python dictionary' + 'Python dictionary' raise ValueError(msg) for key in cutoff: if key == 'weight': @@ -961,7 +1065,7 @@ class Settings: cv.check_greater_than('energy cutoff', cutoff[key], 0.0) else: msg = f'Unable to set cutoff to "{key}" which is unsupported ' \ - 'by OpenMC' + 'by OpenMC' self._cutoff = cutoff @@ -1130,12 +1234,14 @@ class Settings: @weight_window_checkpoints.setter def weight_window_checkpoints(self, weight_window_checkpoints: dict): for key in weight_window_checkpoints.keys(): - cv.check_value('weight_window_checkpoints', key, ('collision', 'surface')) + cv.check_value('weight_window_checkpoints', + key, ('collision', 'surface')) self._weight_window_checkpoints = weight_window_checkpoints @property def max_splits(self): - raise AttributeError('max_splits has been deprecated. Please use max_history_splits instead') + raise AttributeError( + 'max_splits has been deprecated. Please use max_history_splits instead') @property def max_history_splits(self) -> int: @@ -1184,7 +1290,8 @@ class Settings: def weight_window_generators(self, wwgs): if not isinstance(wwgs, MutableSequence): wwgs = [wwgs] - self._weight_window_generators = cv.CheckedList(WeightWindowGenerator, 'weight window generators', wwgs) + self._weight_window_generators = cv.CheckedList( + WeightWindowGenerator, 'weight window generators', wwgs) @property def random_ray(self) -> dict: @@ -1221,7 +1328,8 @@ class Settings: for mesh, domains in value: cv.check_type('mesh', mesh, MeshBase) cv.check_type('domains', domains, Iterable) - valid_types = (openmc.Material, openmc.Cell, openmc.Universe) + valid_types = (openmc.Material, + openmc.Cell, openmc.Universe) for domain in domains: if not isinstance(domain, valid_types): raise ValueError( @@ -1255,9 +1363,12 @@ class Settings: @source_rejection_fraction.setter def source_rejection_fraction(self, source_rejection_fraction: float): - cv.check_type('source_rejection_fraction', source_rejection_fraction, Real) - cv.check_greater_than('source_rejection_fraction', source_rejection_fraction, 0) - cv.check_less_than('source_rejection_fraction', source_rejection_fraction, 1) + cv.check_type('source_rejection_fraction', + source_rejection_fraction, Real) + cv.check_greater_than('source_rejection_fraction', + source_rejection_fraction, 0) + cv.check_less_than('source_rejection_fraction', + source_rejection_fraction, 1) self._source_rejection_fraction = source_rejection_fraction @property @@ -1422,6 +1533,45 @@ class Settings: subelement = ET.SubElement(element, key) subelement.text = str(self._surf_source_write[key]) + def _create_collision_track_subelement(self, root): + if self._collision_track: + element = ET.SubElement(root, "collision_track") + if 'cell_ids' in self._collision_track: + subelement = ET.SubElement(element, "cell_ids") + subelement.text = ' '.join( + str(x) for x in self._collision_track['cell_ids']) + if 'reactions' in self._collision_track: + subelement = ET.SubElement(element, "reactions") + subelement.text = ' '.join( + str(x) for x in self._collision_track['reactions']) + if 'universe_ids' in self._collision_track: + subelement = ET.SubElement(element, "universe_ids") + subelement.text = ' '.join( + str(x) for x in self._collision_track['universe_ids']) + if 'material_ids' in self._collision_track: + subelement = ET.SubElement(element, "material_ids") + subelement.text = ' '.join( + str(x) for x in self._collision_track['material_ids']) + if 'nuclides' in self._collision_track: + subelement = ET.SubElement(element, "nuclides") + subelement.text = ' '.join( + str(x) for x in self._collision_track['nuclides']) + if 'deposited_E_threshold' in self._collision_track: + subelement = ET.SubElement(element, "deposited_E_threshold") + subelement.text = str( + self._collision_track['deposited_E_threshold']) + if 'max_collisions' in self._collision_track: + subelement = ET.SubElement(element, "max_collisions") + subelement.text = str(self._collision_track['max_collisions']) + if 'max_collision_track_files' in self._collision_track: + subelement = ET.SubElement( + element, "max_collision_track_files") + subelement.text = str( + self._collision_track['max_collision_track_files']) + if 'mcpl' in self._collision_track: + subelement = ET.SubElement(element, "mcpl") + subelement.text = str(self._collision_track['mcpl']).lower() + def _create_confidence_intervals(self, root): if self._confidence_intervals is not None: element = ET.SubElement(root, "confidence_intervals") @@ -1477,8 +1627,8 @@ class Settings: # use default heuristic for entropy mesh if not set by user if self.entropy_mesh.dimension is None: if self.particles is None: - raise RuntimeError("Number of particles must be set in order to " \ - "use entropy mesh dimension heuristic") + raise RuntimeError("Number of particles must be set in order to " + "use entropy mesh dimension heuristic") else: n = ceil((self.particles / 20.0)**(1.0 / 3.0)) d = len(self.entropy_mesh.lower_left) @@ -1568,7 +1718,8 @@ class Settings: path = f"./mesh[@id='{self.ufs_mesh.id}']" if root.find(path) is None: root.append(self.ufs_mesh.to_xml_element()) - if mesh_memo is not None: mesh_memo.add(self.ufs_mesh.id) + if mesh_memo is not None: + mesh_memo.add(self.ufs_mesh.id) def _create_use_decay_photons_subelement(self, root): if self._use_decay_photons is not None: @@ -1694,11 +1845,13 @@ class Settings: if 'collision' in self._weight_window_checkpoints: subelement = ET.SubElement(element, "collision") - subelement.text = str(self._weight_window_checkpoints['collision']).lower() + subelement.text = str( + self._weight_window_checkpoints['collision']).lower() if 'surface' in self._weight_window_checkpoints: subelement = ET.SubElement(element, "surface") - subelement.text = str(self._weight_window_checkpoints['surface']).lower() + subelement.text = str( + self._weight_window_checkpoints['surface']).lower() def _create_max_history_splits_subelement(self, root): if self._max_history_splits is not None: @@ -1730,6 +1883,9 @@ class Settings: for domain in domains: domain_elem = ET.SubElement(mesh_elem, 'domain') domain_elem.set('id', str(domain.id)) + domain_elem.set( + 'type', domain.__class__.__name__.lower()) + if mesh_memo is not None and mesh.id not in mesh_memo: domain_elem.set('type', domain.__class__.__name__.lower()) # See if a element already exists -- if not, add it path = f"./mesh[@id='{mesh.id}']" @@ -1881,6 +2037,25 @@ class Settings: value = int(value) self.surf_source_write[key] = value + def _collision_track_from_xml_element(self, root): + elem = root.find('collision_track') + if elem is not None: + for key in ('cell_ids', 'reactions', 'universe_ids', 'material_ids', 'nuclides', + 'deposited_E_threshold', 'max_collisions', "max_collision_track_files", 'mcpl'): + value = get_text(elem, key) + if value is not None: + if key in ('cell_ids', 'universe_ids', 'material_ids'): + value = [int(x) for x in value.split()] + elif key in ('reactions', 'nuclides'): + value = value.split() + elif key in ('max_collisions', 'max_collision_track_files'): + value = int(value) + elif key == 'deposited_E_threshold': + value = float(value) + elif key == 'mcpl': + value = value in ('true', '1') + self.collision_track[key] = value + def _confidence_intervals_from_xml_element(self, root): text = get_text(root, 'confidence_intervals') if text is not None: @@ -2180,7 +2355,8 @@ class Settings: elif domain_type == 'universe': domain = openmc.Universe(domain_id) domains.append(domain) - self.random_ray['source_region_meshes'].append((mesh, domains)) + self.random_ray['source_region_meshes'].append( + (mesh, domains)) def _use_decay_photons_from_xml_element(self, root): text = get_text(root, 'use_decay_photons') @@ -2223,6 +2399,7 @@ class Settings: self._create_sourcepoint_subelement(element) self._create_surf_source_read_subelement(element) self._create_surf_source_write_subelement(element) + self._create_collision_track_subelement(element) self._create_confidence_intervals(element) self._create_electron_treatment_subelement(element) self._create_energy_mode_subelement(element) @@ -2336,6 +2513,7 @@ class Settings: settings._sourcepoint_from_xml_element(elem) settings._surf_source_read_from_xml_element(elem) settings._surf_source_write_from_xml_element(elem) + settings._collision_track_from_xml_element(elem) settings._confidence_intervals_from_xml_element(elem) settings._electron_treatment_from_xml_element(elem) settings._energy_mode_from_xml_element(elem) diff --git a/openmc/source.py b/openmc/source.py index 9b730cf1d..84d8a9619 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -6,7 +6,6 @@ from numbers import Real from pathlib import Path import warnings from typing import Any -from pathlib import Path import lxml.etree as ET import numpy as np @@ -107,10 +106,12 @@ class SourceBase(ABC): cv.check_type('fissionable', value, bool) self._constraints['fissionable'] = value elif key == 'rejection_strategy': - cv.check_value('rejection strategy', value, ('resample', 'kill')) + cv.check_value('rejection strategy', + value, ('resample', 'kill')) self._constraints['rejection_strategy'] = value else: - raise ValueError(f'Unknown key in constraints dictionary: {key}') + raise ValueError( + f'Unknown key in constraints dictionary: {key}') @abstractmethod def populate_xml_element(self, element): @@ -144,13 +145,16 @@ class SourceBase(ABC): dt_elem = ET.SubElement(constraints_elem, "domain_type") dt_elem.text = constraints["domain_type"] id_elem = ET.SubElement(constraints_elem, "domain_ids") - id_elem.text = ' '.join(str(uid) for uid in constraints["domain_ids"]) + id_elem.text = ' '.join(str(uid) + for uid in constraints["domain_ids"]) if "time_bounds" in constraints: dt_elem = ET.SubElement(constraints_elem, "time_bounds") - dt_elem.text = ' '.join(str(t) for t in constraints["time_bounds"]) + dt_elem.text = ' '.join(str(t) + for t in constraints["time_bounds"]) if "energy_bounds" in constraints: dt_elem = ET.SubElement(constraints_elem, "energy_bounds") - dt_elem.text = ' '.join(str(E) for E in constraints["energy_bounds"]) + dt_elem.text = ' '.join(str(E) + for E in constraints["energy_bounds"]) if "fissionable" in constraints: dt_elem = ET.SubElement(constraints_elem, "fissionable") dt_elem.text = str(constraints["fissionable"]).lower() @@ -199,7 +203,8 @@ class SourceBase(ABC): elif source_type == 'mesh': return MeshSource.from_xml_element(elem, meshes) else: - raise ValueError(f'Source type {source_type} is not recognized') + raise ValueError( + f'Source type {source_type} is not recognized') @staticmethod def _get_constraints(elem: ET.Element) -> dict[str, Any]: @@ -316,7 +321,8 @@ class IndependentSource(SourceBase): time: openmc.stats.Univariate | None = None, strength: float = 1.0, particle: str = 'neutron', - domains: Sequence[openmc.Cell | openmc.Material | openmc.Universe] | None = None, + domains: Sequence[openmc.Cell | openmc.Material | + openmc.Universe] | None = None, constraints: dict[str, Any] | None = None ): if domains is not None: @@ -527,11 +533,12 @@ class MeshSource(SourceBase): 'fissionable', and 'rejection_strategy'. """ + def __init__( self, mesh: MeshBase, sources: Sequence[SourceBase], - constraints: dict[str, Any] | None = None, + constraints: dict[str, Any] | None = None, ): super().__init__(strength=None, constraints=constraints) self.mesh = mesh @@ -577,7 +584,8 @@ class MeshSource(SourceBase): elif isinstance(self.mesh, UnstructuredMesh): if s.ndim > 1: - raise ValueError('Sources must be a 1-D array for unstructured mesh') + raise ValueError( + 'Sources must be a 1-D array for unstructured mesh') self._sources = s for src in self._sources: @@ -646,7 +654,8 @@ class MeshSource(SourceBase): mesh_id = int(get_text(elem, 'mesh')) mesh = meshes[mesh_id] - sources = [SourceBase.from_xml_element(e) for e in elem.iterchildren('source')] + sources = [SourceBase.from_xml_element( + e) for e in elem.iterchildren('source')] constraints = cls._get_constraints(elem) return cls(mesh, sources, constraints=constraints) @@ -656,7 +665,8 @@ def Source(*args, **kwargs): A function for backward compatibility of sources. Will be removed in the future. Please update to IndependentSource. """ - warnings.warn("This class is deprecated in favor of 'IndependentSource'", FutureWarning) + warnings.warn( + "This class is deprecated in favor of 'IndependentSource'", FutureWarning) return openmc.IndependentSource(*args, **kwargs) @@ -703,6 +713,7 @@ class CompiledSource(SourceBase): 'fissionable', and 'rejection_strategy'. """ + def __init__( self, library: PathLike, @@ -914,7 +925,8 @@ class ParticleType(IntEnum): try: return cls[value.upper()] except KeyError: - raise ValueError(f"Invalid string for creation of {cls.__name__}: {value}") + raise ValueError( + f"Invalid string for creation of {cls.__name__}: {value}") @classmethod def from_pdg_number(cls, pdg_number: int) -> ParticleType: @@ -984,6 +996,7 @@ class SourceParticle: Type of the particle """ + def __init__( self, r: Iterable[float] = (0., 0., 0.), @@ -1042,7 +1055,8 @@ def write_source_file( openmc.SourceParticle """ - cv.check_iterable_type("source particles", source_particles, SourceParticle) + cv.check_iterable_type( + "source particles", source_particles, SourceParticle) pl = ParticleList(source_particles) pl.export_to_hdf5(filename, **kwargs) @@ -1104,7 +1118,8 @@ class ParticleList(list): for particle in f.particles: # Determine particle type based on the PDG number try: - particle_type = ParticleType.from_pdg_number(particle.pdgcode) + particle_type = ParticleType.from_pdg_number( + particle.pdgcode) except ValueError: particle_type = "UNKNOWN" @@ -1241,3 +1256,133 @@ def read_source_file(filename: PathLike) -> ParticleList: return ParticleList.from_hdf5(filename) else: return ParticleList.from_mcpl(filename) + + +def read_collision_track_hdf5(filename): + """Read a collision track file in HDF5 format. + + Parameters + ---------- + filename : str or path-like + Path to the HDF5 collision track file. + + Returns + ------- + numpy.ndarray + Structured array containing collision track data. + + See Also + -------- + read_collision_track_mcpl + read_collision_track_file + """ + + with h5py.File(filename, 'r') as file: + data = file['collision_track_bank'][:] + + return data + + +def read_collision_track_mcpl(file_path): + """Read a collision track file in MCPL format. + + Parameters + ---------- + file_path : str or path-like + Path to the MCPL collision track file. + + Returns + ------- + numpy.ndarray + Structured array of particle collision track information, including + position, direction, energy, weight, reaction data, and identifiers. + + See Also + -------- + read_collision_track_hdf5 + read_collision_track_file + """ + import mcpl + myfile = mcpl.MCPLFile(file_path) + data = { + 'r': [], # for position (x, y, z) + 'u': [], # for direction (ux, uy, uz) + 'E': [], 'dE': [], 'time': [], + 'wgt': [], 'event_mt': [], 'delayed_group': [], + 'cell_id': [], 'nuclide_id': [], 'material_id': [], + 'universe_id': [], 'n_collision': [], 'particle': [], + 'parent_id': [], 'progeny_id': [] + } + + # Read and collect data from the MCPL file + for i, p in enumerate(myfile.particles): + if f'blob_{i}' in myfile.blobs: + blob_data = myfile.blobs[f'blob_{i}'] + decoded_str = blob_data.decode('utf-8') + pairs = decoded_str.split(';') + values_dict = {k.strip(): v.strip() + for k, v in (pair.split(':') for pair in pairs if pair.strip())} + + data['r'].append((p.x, p.y, p.z)) # Append as tuple + data['u'].append((p.ux, p.uy, p.uz)) # Append as tuple + data['E'].append(p.ekin * 1e6) + data['dE'].append(float(values_dict.get('dE', 0))) + data['time'].append(p.time * 1e-3) + data['wgt'].append(p.weight) + data['event_mt'].append(int(values_dict.get('event_mt', 0))) + data['delayed_group'].append( + int(values_dict.get('delayed_group', 0))) + data['cell_id'].append(int(values_dict.get('cell_id', 0))) + data['nuclide_id'].append(int(values_dict.get('nuclide_id', 0))) + data['material_id'].append(int(values_dict.get('material_id', 0))) + data['universe_id'].append(int(values_dict.get('universe_id', 0))) + data['n_collision'].append(int(values_dict.get('n_collision', 0))) + data['particle'].append(ParticleType.from_pdg_number(p.pdgcode)) + data['parent_id'].append(int(values_dict.get('parent_id', 0))) + data['progeny_id'].append(int(values_dict.get('progeny_id', 0))) + + dtypes = [ + ('r', [('x', 'f8'), ('y', 'f8'), ('z', 'f8')]), + ('u', [('x', 'f8'), ('y', 'f8'), ('z', 'f8')]), + ('E', 'f8'), ('dE', 'f8'), ('time', 'f8'), ('wgt', 'f8'), + ('event_mt', 'f8'), ('delayed_group', 'i4'), ('cell_id', 'i4'), + ('nuclide_id', 'i4'), ('material_id', 'i4'), ('universe_id', 'i4'), + ('n_collision', 'i4'), ('particle', 'i4'), + ('parent_id', 'i8'), ('progeny_id', 'i8') + ] + + structured_array = np.zeros(len(data['r']), dtype=dtypes) + for key in data: + structured_array[key] = data[key] # Assign data + + return structured_array + + +def read_collision_track_file(filename): + """Read a collision track file (HDF5 or MCPL) and return its data. + + Parameters + ---------- + filename : str or path-like + Path to the collision track file to read. Must end with + ``.h5`` or ``.mcpl``. + + Returns + ------- + numpy.ndarray + Structured array containing collision track data. + + See Also + -------- + read_collision_track_hdf5 + read_collision_track_mcpl + """ + + filename = Path(filename) + if filename.suffix not in ('.h5', '.mcpl'): + raise ValueError('Collision track file must have a .h5 or .mcpl extension.') + + if filename.suffix == '.h5': + return read_collision_track_hdf5(filename) + else: + return read_collision_track_mcpl(filename) diff --git a/src/bank.cpp b/src/bank.cpp index 3e806b3c0..33790379b 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -20,6 +20,8 @@ vector source_bank; SharedArray surf_source_bank; +SharedArray collision_track_bank; + // The fission bank is allocated as a SharedArray, rather than a vector, as it // will be shared by all threads in the simulation. It will be allocated to a // fixed maximum capacity in the init_fission_bank() function. Then, Elements @@ -50,6 +52,7 @@ void free_memory_bank() { simulation::source_bank.clear(); simulation::surf_source_bank.clear(); + simulation::collision_track_bank.clear(); simulation::fission_bank.clear(); simulation::progeny_per_particle.clear(); simulation::ifp_source_delayed_group_bank.clear(); diff --git a/src/collision_track.cpp b/src/collision_track.cpp new file mode 100644 index 000000000..75e56574a --- /dev/null +++ b/src/collision_track.cpp @@ -0,0 +1,238 @@ +#include "openmc/collision_track.h" + +#include +#include + +#include + +#include "openmc/bank.h" +#include "openmc/bank_io.h" +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/file_utils.h" +#include "openmc/hdf5_interface.h" +#include "openmc/material.h" +#include "openmc/mcpl_interface.h" +#include "openmc/message_passing.h" +#include "openmc/nuclide.h" +#include "openmc/output.h" +#include "openmc/particle.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" +#include "openmc/universe.h" + +#ifdef OPENMC_MPI +#include +#endif + +namespace openmc { + +namespace { + +hid_t h5_collision_track_banktype() +{ + hid_t postype = H5Tcreate(H5T_COMPOUND, sizeof(Position)); + H5Tinsert(postype, "x", HOFFSET(Position, x), H5T_NATIVE_DOUBLE); + H5Tinsert(postype, "y", HOFFSET(Position, y), H5T_NATIVE_DOUBLE); + H5Tinsert(postype, "z", HOFFSET(Position, z), H5T_NATIVE_DOUBLE); + + hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(CollisionTrackSite)); + + H5Tinsert(banktype, "r", HOFFSET(CollisionTrackSite, r), postype); + H5Tinsert(banktype, "u", HOFFSET(CollisionTrackSite, u), postype); + H5Tinsert(banktype, "E", HOFFSET(CollisionTrackSite, E), H5T_NATIVE_DOUBLE); + H5Tinsert(banktype, "dE", HOFFSET(CollisionTrackSite, dE), H5T_NATIVE_DOUBLE); + H5Tinsert( + banktype, "time", HOFFSET(CollisionTrackSite, time), H5T_NATIVE_DOUBLE); + H5Tinsert( + banktype, "wgt", HOFFSET(CollisionTrackSite, wgt), H5T_NATIVE_DOUBLE); + H5Tinsert(banktype, "event_mt", HOFFSET(CollisionTrackSite, event_mt), + H5T_NATIVE_INT); + H5Tinsert(banktype, "delayed_group", + HOFFSET(CollisionTrackSite, delayed_group), H5T_NATIVE_INT); + H5Tinsert( + banktype, "cell_id", HOFFSET(CollisionTrackSite, cell_id), H5T_NATIVE_INT); + H5Tinsert(banktype, "nuclide_id", HOFFSET(CollisionTrackSite, nuclide_id), + H5T_NATIVE_INT); + H5Tinsert(banktype, "material_id", HOFFSET(CollisionTrackSite, material_id), + H5T_NATIVE_INT); + H5Tinsert(banktype, "universe_id", HOFFSET(CollisionTrackSite, universe_id), + H5T_NATIVE_INT); + H5Tinsert(banktype, "n_collision", HOFFSET(CollisionTrackSite, n_collision), + H5T_NATIVE_INT); + H5Tinsert(banktype, "particle", HOFFSET(CollisionTrackSite, particle), + H5T_NATIVE_INT); + H5Tinsert(banktype, "parent_id", HOFFSET(CollisionTrackSite, parent_id), + H5T_NATIVE_INT64); + H5Tinsert(banktype, "progeny_id", HOFFSET(CollisionTrackSite, progeny_id), + H5T_NATIVE_INT64); + H5Tclose(postype); + return banktype; +} + +void write_collision_track_bank(hid_t group_id, + openmc::span collision_track_bank, + const openmc::vector& bank_index) +{ + hid_t banktype = h5_collision_track_banktype(); +#ifdef OPENMC_MPI + write_bank_dataset("collision_track_bank", group_id, collision_track_bank, + bank_index, banktype, mpi::collision_track_site); +#else + write_bank_dataset("collision_track_bank", group_id, collision_track_bank, + bank_index, banktype); +#endif + + H5Tclose(banktype); +} + +void write_h5_collision_track(const char* filename, + openmc::span collision_track_bank, + const openmc::vector& bank_index) +{ +#ifdef PHDF5 + bool parallel = true; +#else + bool parallel = false; +#endif + + if (!filename) + fatal_error("write_h5_collision_track filename needs a nonempty name."); + + std::string filename_(filename); + const auto extension = get_file_extension(filename_); + if (extension.empty()) { + filename_.append(".h5"); + } else if (extension != "h5") { + warning("write_h5_collision_track was passed a file extension differing " + "from .h5, but an hdf5 file will be written."); + } + + hid_t file_id; + if (mpi::master || parallel) { + file_id = file_open(filename_.c_str(), 'w', true); + + // Write filetype and version info + write_attribute(file_id, "filetype", "collision_track"); + write_attribute(file_id, "version", VERSION_COLLISION_TRACK); + } + + write_collision_track_bank(file_id, collision_track_bank, bank_index); + + if (mpi::master || parallel) + file_close(file_id); +} + +} // namespace + +bool should_record_event(int id_cell, int mt_event, const std::string& nuclide, + int id_universe, int id_material, double energy_loss) +{ + auto matches_filter = [](const auto& filter_set, const auto& value) { + return filter_set.empty() || filter_set.count(value) > 0; + }; + + const auto& cfg = settings::collision_track_config; + return simulation::current_batch > settings::n_inactive && + !simulation::collision_track_bank.full() && + matches_filter(cfg.cell_ids, id_cell) && + matches_filter(cfg.mt_numbers, mt_event) && + matches_filter(cfg.universe_ids, id_universe) && + matches_filter(cfg.material_ids, id_material) && + matches_filter(cfg.nuclides, nuclide) && + (cfg.deposited_energy_threshold == 0 || + cfg.deposited_energy_threshold < energy_loss); +} + +void collision_track_reserve_bank() +{ + simulation::collision_track_bank.reserve( + settings::collision_track_config.max_collisions); +} + +void collision_track_flush_bank() +{ + const auto& cfg = settings::collision_track_config; + if (simulation::ct_current_file > cfg.max_files) + return; + + bool last_batch = (simulation::current_batch == settings::n_batches); + if (!simulation::collision_track_bank.full() && !last_batch) + return; + + auto size = simulation::collision_track_bank.size(); + if (size == 0 && !last_batch) + return; + + auto collision_track_work_index = mpi::calculate_parallel_index_vector(size); + openmc::span collisiontrackbankspan( + simulation::collision_track_bank.begin(), size); + + std::string ext = cfg.mcpl_write ? "mcpl" : "h5"; + auto filename = fmt::format("{}collision_track.{}.{}", settings::path_output, + simulation::ct_current_file, ext); + + if (cfg.max_files == 1 || (simulation::ct_current_file == 1 && last_batch)) { + filename = settings::path_output + "collision_track." + ext; + } + write_message("Creating {}...", filename, 4); + + if (cfg.mcpl_write) { + write_mcpl_collision_track( + filename.c_str(), collisiontrackbankspan, collision_track_work_index); + } else { + write_h5_collision_track( + filename.c_str(), collisiontrackbankspan, collision_track_work_index); + } + + simulation::collision_track_bank.clear(); + if (!last_batch && cfg.max_files >= 1) { + collision_track_reserve_bank(); + } + ++simulation::ct_current_file; +} + +void collision_track_record(Particle& particle) +{ + int cell_index = particle.lowest_coord().cell(); + if (cell_index == C_NONE) + return; + + int cell_id = model::cells[cell_index]->id_; + const auto* nuclide_ptr = data::nuclides[particle.event_nuclide()].get(); + std::string nuclide = nuclide_ptr->name_; + int universe_id = model::universes[particle.lowest_coord().universe()]->id_; + double delta_E = particle.E_last() - particle.E(); + int material_index = particle.material(); + if (material_index == C_NONE) + return; + + int material_id = model::materials[material_index]->id_; + + if (!should_record_event(cell_id, particle.event_mt(), nuclide, universe_id, + material_id, delta_E)) + return; + + CollisionTrackSite site; + site.r = particle.r(); + site.u = particle.u(); + site.E = particle.E_last(); + site.dE = delta_E; + site.time = particle.time(); + site.wgt = particle.wgt(); + site.event_mt = particle.event_mt(); + site.delayed_group = particle.delayed_group(); + site.cell_id = cell_id; + site.nuclide_id = + 10000 * nuclide_ptr->Z_ + 10 * nuclide_ptr->A_ + nuclide_ptr->metastable_; + site.material_id = material_id; + site.universe_id = universe_id; + site.n_collision = particle.n_collision(); + site.particle = particle.type(); + site.parent_id = particle.id(); + site.progeny_id = particle.n_progeny(); + simulation::collision_track_bank.thread_safe_append(site); +} + +} // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 659f390b3..9ee943409 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -3,6 +3,7 @@ #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/cmfd_solver.h" +#include "openmc/collision_track.h" #include "openmc/constants.h" #include "openmc/cross_sections.h" #include "openmc/dagmc.h" @@ -76,6 +77,7 @@ int openmc_finalize() // Reset global variables settings::assume_separate = false; settings::check_overlaps = false; + settings::collision_track_config = CollisionTrackConfig {}; settings::confidence_intervals = false; settings::create_fission_neutrons = true; settings::create_delayed_neutrons = true; @@ -177,6 +179,9 @@ int openmc_finalize() if (mpi::source_site != MPI_DATATYPE_NULL) { MPI_Type_free(&mpi::source_site); } + if (mpi::collision_track_site != MPI_DATATYPE_NULL) { + MPI_Type_free(&mpi::collision_track_site); + } #endif openmc_reset_random_ray(); diff --git a/src/initialize.cpp b/src/initialize.cpp index 36f326116..e2a5b9743 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -178,6 +178,37 @@ void initialize_mpi(MPI_Comm intracomm) MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); + + CollisionTrackSite bc; + MPI_Aint dispc[16]; + MPI_Get_address(&bc.r, &dispc[0]); // double + MPI_Get_address(&bc.u, &dispc[1]); // double + MPI_Get_address(&bc.E, &dispc[2]); // double + MPI_Get_address(&bc.dE, &dispc[3]); // double + MPI_Get_address(&bc.time, &dispc[4]); // double + MPI_Get_address(&bc.wgt, &dispc[5]); // double + MPI_Get_address(&bc.event_mt, &dispc[6]); // int + MPI_Get_address(&bc.delayed_group, &dispc[7]); // int + MPI_Get_address(&bc.cell_id, &dispc[8]); // int + MPI_Get_address(&bc.nuclide_id, &dispc[9]); // int + MPI_Get_address(&bc.material_id, &dispc[10]); // int + MPI_Get_address(&bc.universe_id, &dispc[11]); // int + MPI_Get_address(&bc.n_collision, &dispc[12]); // int + MPI_Get_address(&bc.particle, &dispc[13]); // int + MPI_Get_address(&bc.parent_id, &dispc[14]); // int64_t + MPI_Get_address(&bc.progeny_id, &dispc[15]); // int64_t + for (int i = 15; i >= 0; --i) { + dispc[i] -= dispc[0]; + } + + int blocksc[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + MPI_Datatype typesc[] = {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, + MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT, + MPI_INT, MPI_INT, MPI_INT, MPI_INT64_T, MPI_INT64_T}; + + MPI_Type_create_struct( + 16, blocksc, dispc, typesc, &mpi::collision_track_site); + MPI_Type_commit(&mpi::collision_track_site); } #endif // OPENMC_MPI diff --git a/src/mcpl_interface.cpp b/src/mcpl_interface.cpp index b8e780710..129407301 100644 --- a/src/mcpl_interface.cpp +++ b/src/mcpl_interface.cpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -61,6 +62,8 @@ using mcpl_hdr_nparticles_fpt = uint64_t (*)(mcpl_file_t* file_handle); using mcpl_read_fpt = const mcpl_particle_repr_t* (*)(mcpl_file_t* file_handle); using mcpl_close_file_fpt = void (*)(mcpl_file_t* file_handle); +using mcpl_hdr_add_data_fpt = void (*)(mcpl_outfile_t* file_handle, + const char* key, int32_t ldata, const char* data); using mcpl_create_outfile_fpt = mcpl_outfile_t* (*)(const char* filename); using mcpl_hdr_set_srcname_fpt = void (*)( mcpl_outfile_t* outfile_handle, const char* srcname); @@ -110,6 +113,7 @@ struct McplApi { mcpl_close_file_fpt close_file; mcpl_create_outfile_fpt create_outfile; mcpl_hdr_set_srcname_fpt hdr_set_srcname; + mcpl_hdr_add_data_fpt hdr_add_data; mcpl_add_particle_fpt add_particle; mcpl_close_outfile_fpt close_outfile; mcpl_hdr_add_stat_sum_fpt hdr_add_stat_sum; @@ -146,6 +150,8 @@ struct McplApi { load_symbol_platform("mcpl_create_outfile")); hdr_set_srcname = reinterpret_cast( load_symbol_platform("mcpl_hdr_set_srcname")); + hdr_add_data = reinterpret_cast( + load_symbol_platform("mcpl_hdr_add_data")); add_particle = reinterpret_cast( load_symbol_platform("mcpl_add_particle")); close_outfile = reinterpret_cast( @@ -545,4 +551,153 @@ void write_mcpl_source_point(const char* filename, span source_bank, } } +// Collision track feature with MCPL +void write_mcpl_collision_track_internal(mcpl_outfile_t* file_id, + span collision_track_bank, + const vector& bank_index_all_ranks) +{ + if (mpi::master) { + if (!file_id) { + fatal_error("MCPL: Internal error - master rank called " + "write_mcpl_source_bank_internal with null file_id."); + } + vector receive_buffer; + vector all_sites; + all_sites.reserve(static_cast(bank_index_all_ranks.back())); + vector all_blobs; + all_blobs.reserve(static_cast(bank_index_all_ranks.back())); + + for (int rank_idx = 0; rank_idx < mpi::n_procs; ++rank_idx) { + size_t num_sites_on_rank = static_cast( + bank_index_all_ranks[rank_idx + 1] - bank_index_all_ranks[rank_idx]); + if (num_sites_on_rank == 0) + continue; + + span sites_to_process; +#ifdef OPENMC_MPI + if (rank_idx == mpi::rank) { + sites_to_process = openmc::span( + collision_track_bank.data(), num_sites_on_rank); + } else { + receive_buffer.resize(num_sites_on_rank); + MPI_Recv(receive_buffer.data(), num_sites_on_rank, + mpi::collision_track_site, rank_idx, rank_idx, mpi::intracomm, + MPI_STATUS_IGNORE); + sites_to_process = openmc::span( + receive_buffer.data(), num_sites_on_rank); + } +#else + sites_to_process = openmc::span( + collision_track_bank.data(), num_sites_on_rank); +#endif + + for (const auto& site : sites_to_process) { + std::ostringstream custom_data_stream; + custom_data_stream << " dE : " << site.dE + << " ; event_mt : " << site.event_mt + << " ; delayed_group : " << site.delayed_group + << " ; cell_id : " << site.cell_id + << " ; nuclide_id : " << site.nuclide_id + << " ; material_id : " << site.material_id + << " ; universe_id : " << site.universe_id + << " ; n_collision : " << site.n_collision + << " ; parent_id : " << site.parent_id + << " ; progeny_id : " << site.progeny_id; + + all_blobs.push_back(custom_data_stream.str()); + all_sites.push_back(site); + } + } + + for (size_t idx = 0; idx < all_blobs.size(); ++idx) { + const auto& blob = all_blobs[idx]; + std::string key = "blob_" + std::to_string(idx); + g_mcpl_api->hdr_add_data(file_id, key.c_str(), blob.size(), blob.c_str()); + } + + for (const auto& site : all_sites) { + mcpl_particle_repr_t p_repr {}; + p_repr.position[0] = site.r.x; + p_repr.position[1] = site.r.y; + p_repr.position[2] = site.r.z; + p_repr.direction[0] = site.u.x; + p_repr.direction[1] = site.u.y; + p_repr.direction[2] = site.u.z; + p_repr.ekin = site.E * 1e-6; + p_repr.time = site.time * 1e3; + p_repr.weight = site.wgt; + switch (site.particle) { + case ParticleType::neutron: + p_repr.pdgcode = 2112; + break; + case ParticleType::photon: + p_repr.pdgcode = 22; + break; + case ParticleType::electron: + p_repr.pdgcode = 11; + break; + case ParticleType::positron: + p_repr.pdgcode = -11; + break; + default: + continue; + } + g_mcpl_api->add_particle(file_id, &p_repr); + } + } else { +#ifdef OPENMC_MPI + if (!collision_track_bank.empty()) { + MPI_Send(collision_track_bank.data(), collision_track_bank.size(), + mpi::collision_track_site, 0, mpi::rank, mpi::intracomm); + } +#endif + } +} + +void write_mcpl_collision_track(const char* filename, + span collision_track_bank, + const vector& bank_index) +{ + ensure_mcpl_ready_or_fatal(); + + std::string filename_(filename); + const auto extension = get_file_extension(filename_); + if (extension.empty()) { + filename_.append(".mcpl"); + } else if (extension != "mcpl") { + warning(fmt::format("Specified filename '{}' has an extension '.{}', but " + "an MCPL file (.mcpl) will be written using this name.", + filename, extension)); + } + + mcpl_outfile_t* file_id = nullptr; + + if (mpi::master) { + file_id = g_mcpl_api->create_outfile(filename_.c_str()); + if (!file_id) { + fatal_error(fmt::format( + "MCPL: Failed to create output file '{}'. Check permissions and path.", + filename_)); + } + std::string src_line; + if (VERSION_DEV) { + src_line = fmt::format("OpenMC {}.{}.{}-dev{}", VERSION_MAJOR, + VERSION_MINOR, VERSION_RELEASE, VERSION_COMMIT_COUNT); + } else { + src_line = fmt::format( + "OpenMC {}.{}.{}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE); + } + + g_mcpl_api->hdr_set_srcname(file_id, src_line.c_str()); + } + write_mcpl_collision_track_internal( + file_id, collision_track_bank, bank_index); + + if (mpi::master) { + if (file_id) { + g_mcpl_api->close_outfile(file_id); + } + } +} + } // namespace openmc diff --git a/src/message_passing.cpp b/src/message_passing.cpp index 374c1aa72..a160f6d73 100644 --- a/src/message_passing.cpp +++ b/src/message_passing.cpp @@ -10,6 +10,7 @@ bool master {true}; #ifdef OPENMC_MPI MPI_Comm intracomm {MPI_COMM_NULL}; MPI_Datatype source_site {MPI_DATATYPE_NULL}; +MPI_Datatype collision_track_site {MPI_DATATYPE_NULL}; #endif extern "C" bool openmc_master() diff --git a/src/particle.cpp b/src/particle.cpp index 6ba8ebf12..2d70d715e 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -8,6 +8,7 @@ #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/cell.h" +#include "openmc/collision_track.h" #include "openmc/constants.h" #include "openmc/dagmc.h" #include "openmc/error.h" @@ -351,6 +352,11 @@ void Particle::event_collide() collision_mg(*this); } + // Collision track feature to recording particle interaction + if (settings::collision_track) { + collision_track_record(*this); + } + // Score collision estimator tallies -- this is done after a collision // has occurred rather than before because we need information on the // outgoing energy for any tallies with an outgoing energy filter diff --git a/src/physics.cpp b/src/physics.cpp index 3a17077b3..41509af97 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -122,6 +122,7 @@ void sample_neutron_reaction(Particle& p) "with k-effective close to or greater than one."); } } + p.event_mt() = rx.mt_; } // Create secondary photons @@ -663,7 +664,9 @@ void absorption(Particle& p, int i_nuclide) p.wgt() = 0.0; p.event() = TallyEvent::ABSORB; - p.event_mt() = N_DISAPPEAR; + if (!p.fission()) { + p.event_mt() = N_DISAPPEAR; + } } } } diff --git a/src/settings.cpp b/src/settings.cpp index 13b91b0e4..9dcf7c8db 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -11,6 +11,7 @@ #endif #include "openmc/capi.h" +#include "openmc/collision_track.h" #include "openmc/constants.h" #include "openmc/container_util.h" #include "openmc/distribution.h" @@ -26,6 +27,7 @@ #include "openmc/plot.h" #include "openmc/random_lcg.h" #include "openmc/random_ray/random_ray.h" +#include "openmc/reaction.h" #include "openmc/simulation.h" #include "openmc/source.h" #include "openmc/string_utils.h" @@ -45,6 +47,7 @@ namespace settings { // Default values for boolean flags bool assume_separate {false}; bool check_overlaps {false}; +bool collision_track {false}; bool cmfd_run {false}; bool confidence_intervals {false}; bool create_delayed_neutrons {true}; @@ -128,6 +131,7 @@ std::unordered_set statepoint_batch; double source_rejection_fraction {0.05}; double free_gas_threshold {400.0}; std::unordered_set source_write_surf_id; +CollisionTrackConfig collision_track_config {}; int64_t ssw_max_particles; int64_t ssw_max_files; int64_t ssw_cell_id {C_NONE}; @@ -925,8 +929,72 @@ void read_settings_xml(pugi::xml_node root) } } - // If source is not separate and is to be written out in the statepoint file, - // make sure that the sourcepoint batch numbers are contained in the + // Check if the user has specified to write specific collisions + if (check_for_node(root, "collision_track")) { + settings::collision_track = true; + // Get collision track node + xml_node node_ct = root.child("collision_track"); + collision_track_config = CollisionTrackConfig {}; + + // Determine cell ids at which crossing particles are to be banked + if (check_for_node(node_ct, "cell_ids")) { + auto temp = get_node_array(node_ct, "cell_ids"); + for (const auto& b : temp) { + collision_track_config.cell_ids.insert(b); + } + } + if (check_for_node(node_ct, "reactions")) { + auto temp = get_node_array(node_ct, "reactions"); + for (const auto& b : temp) { + int reaction_int = reaction_type(b); + if (reaction_int > 0) { + collision_track_config.mt_numbers.insert(reaction_int); + } + } + } + if (check_for_node(node_ct, "universe_ids")) { + auto temp = get_node_array(node_ct, "universe_ids"); + for (const auto& b : temp) { + collision_track_config.universe_ids.insert(b); + } + } + if (check_for_node(node_ct, "material_ids")) { + auto temp = get_node_array(node_ct, "material_ids"); + for (const auto& b : temp) { + collision_track_config.material_ids.insert(b); + } + } + if (check_for_node(node_ct, "nuclides")) { + auto temp = get_node_array(node_ct, "nuclides"); + for (const auto& b : temp) { + collision_track_config.nuclides.insert(b); + } + } + if (check_for_node(node_ct, "deposited_E_threshold")) { + collision_track_config.deposited_energy_threshold = + std::stod(get_node_value(node_ct, "deposited_E_threshold")); + } + // Get maximum number of particles to be banked per collision + if (check_for_node(node_ct, "max_collisions")) { + collision_track_config.max_collisions = + std::stoll(get_node_value(node_ct, "max_collisions")); + } else { + warning("A maximum number of collisions needs to be specified. " + "By default the code sets 'max_collisions' parameter equals to " + "1000."); + } + // Get maximum number of collision_track files to be created + if (check_for_node(node_ct, "max_collision_track_files")) { + collision_track_config.max_files = + std::stoll(get_node_value(node_ct, "max_collision_track_files")); + } + if (check_for_node(node_ct, "mcpl")) { + collision_track_config.mcpl_write = get_node_value_bool(node_ct, "mcpl"); + } + } + + // If source is not separate and is to be written out in the statepoint + // file, make sure that the sourcepoint batch numbers are contained in the // statepoint list if (!source_separate) { for (const auto& b : sourcepoint_batch) { @@ -1171,8 +1239,8 @@ void read_settings_xml(pugi::xml_node root) variance_reduction::weight_windows_generators.emplace_back( std::make_unique(node_wwg)); } - // if any of the weight windows are intended to be generated otf, make sure - // they're applied + // if any of the weight windows are intended to be generated otf, make + // sure they're applied for (const auto& wwg : variance_reduction::weight_windows_generators) { if (wwg->on_the_fly_) { settings::weight_windows_on = true; diff --git a/src/simulation.cpp b/src/simulation.cpp index 05d554526..b536ae588 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -2,6 +2,7 @@ #include "openmc/bank.h" #include "openmc/capi.h" +#include "openmc/collision_track.h" #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" @@ -9,7 +10,6 @@ #include "openmc/geometry_aux.h" #include "openmc/ifp.h" #include "openmc/material.h" -#include "openmc/mcpl_interface.h" #include "openmc/message_passing.h" #include "openmc/nuclide.h" #include "openmc/output.h" @@ -118,6 +118,7 @@ int openmc_simulation_init() // Reset global variables -- this is done before loading state point (as that // will potentially populate k_generation and entropy) simulation::current_batch = 0; + simulation::ct_current_file = 1; simulation::ssw_current_file = 1; simulation::k_generation.clear(); simulation::entropy.clear(); @@ -297,6 +298,7 @@ namespace openmc { namespace simulation { +int ct_current_file; int current_batch; int current_gen; bool initialized {false}; @@ -347,6 +349,11 @@ void allocate_banks() // Allocate surface source bank simulation::surf_source_bank.reserve(settings::ssw_max_particles); } + + if (settings::collision_track) { + // Allocate collision track bank + collision_track_reserve_bank(); + } } void initialize_batch() @@ -490,6 +497,10 @@ void finalize_batch() ++simulation::ssw_current_file; } } + // Write collision track file if requested + if (settings::collision_track) { + collision_track_flush_bank(); + } } void initialize_generation() diff --git a/src/state_point.cpp b/src/state_point.cpp index 0b0fed132..47296da3a 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -9,6 +9,7 @@ #include #include "openmc/bank.h" +#include "openmc/bank_io.h" #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/eigenvalue.h" @@ -642,91 +643,12 @@ void write_source_bank(hid_t group_id, span source_bank, { hid_t banktype = h5banktype(); - // Set total and individual process dataspace sizes for source bank - int64_t dims_size = bank_index.back(); - int64_t count_size = bank_index[mpi::rank + 1] - bank_index[mpi::rank]; - -#ifdef PHDF5 - // Set size of total dataspace for all procs and rank - hsize_t dims[] {static_cast(dims_size)}; - hid_t dspace = H5Screate_simple(1, dims, nullptr); - hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace, H5P_DEFAULT, - H5P_DEFAULT, H5P_DEFAULT); - - // Create another data space but for each proc individually - hsize_t count[] {static_cast(count_size)}; - hid_t memspace = H5Screate_simple(1, count, nullptr); - - // Select hyperslab for this dataspace - hsize_t start[] {static_cast(bank_index[mpi::rank])}; - H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr); - - // Set up the property list for parallel writing - hid_t plist = H5Pcreate(H5P_DATASET_XFER); - H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE); - - // Write data to file in parallel - H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank.data()); - - // Free resources - H5Sclose(dspace); - H5Sclose(memspace); - H5Dclose(dset); - H5Pclose(plist); - +#ifdef OPENMC_MPI + write_bank_dataset("source_bank", group_id, source_bank, bank_index, banktype, + mpi::source_site); #else - - if (mpi::master) { - // Create dataset big enough to hold all source sites - hsize_t dims[] {static_cast(dims_size)}; - hid_t dspace = H5Screate_simple(1, dims, nullptr); - hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - - // Save source bank sites since the array is overwritten below -#ifdef OPENMC_MPI - vector temp_source {source_bank.begin(), source_bank.end()}; -#endif - - for (int i = 0; i < mpi::n_procs; ++i) { - // Create memory space - hsize_t count[] {static_cast(bank_index[i + 1] - bank_index[i])}; - hid_t memspace = H5Screate_simple(1, count, nullptr); - -#ifdef OPENMC_MPI - // Receive source sites from other processes - if (i > 0) - MPI_Recv(source_bank.data(), count[0], mpi::source_site, i, i, - mpi::intracomm, MPI_STATUS_IGNORE); -#endif - - // Select hyperslab for this dataspace - dspace = H5Dget_space(dset); - hsize_t start[] {static_cast(bank_index[i])}; - H5Sselect_hyperslab( - dspace, H5S_SELECT_SET, start, nullptr, count, nullptr); - - // Write data to hyperslab - H5Dwrite( - dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank.data()); - - H5Sclose(memspace); - H5Sclose(dspace); - } - - // Close all ids - H5Dclose(dset); - -#ifdef OPENMC_MPI - // Restore state of source bank - std::copy(temp_source.begin(), temp_source.end(), source_bank.begin()); -#endif - } else { -#ifdef OPENMC_MPI - MPI_Send(source_bank.data(), count_size, mpi::source_site, 0, mpi::rank, - mpi::intracomm); -#endif - } + write_bank_dataset( + "source_bank", group_id, source_bank, bank_index, banktype); #endif H5Tclose(banktype); diff --git a/tests/regression_tests/collision_track/__init__.py b/tests/regression_tests/collision_track/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat b/tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat new file mode 100644 index 000000000..7533616c0 --- /dev/null +++ b/tests/regression_tests/collision_track/case_1_Reactions/inputs_true.dat @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 5 + 1 + + + -2.0 -2.0 -2.0 2.0 2.0 2.0 + + + true + + + + (n,fission) 101 + 300 + + 1 + + diff --git a/tests/regression_tests/collision_track/case_1_Reactions/results_true.dat b/tests/regression_tests/collision_track/case_1_Reactions/results_true.dat new file mode 100644 index 000000000..d4d1d1e5a --- /dev/null +++ b/tests/regression_tests/collision_track/case_1_Reactions/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat b/tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat new file mode 100644 index 000000000..55fb835de --- /dev/null +++ b/tests/regression_tests/collision_track/case_2_Cell_ID/inputs_true.dat @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 5 + 1 + + + -2.0 -2.0 -2.0 2.0 2.0 2.0 + + + true + + + + 22 + 300 + + 1 + + diff --git a/tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat b/tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat new file mode 100644 index 000000000..d4d1d1e5a --- /dev/null +++ b/tests/regression_tests/collision_track/case_2_Cell_ID/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat b/tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat new file mode 100644 index 000000000..61890414b --- /dev/null +++ b/tests/regression_tests/collision_track/case_3_Material_ID/inputs_true.dat @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 5 + 1 + + + -2.0 -2.0 -2.0 2.0 2.0 2.0 + + + true + + + + 1 + 300 + + 1 + + diff --git a/tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat b/tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat new file mode 100644 index 000000000..d4d1d1e5a --- /dev/null +++ b/tests/regression_tests/collision_track/case_3_Material_ID/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat b/tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat new file mode 100644 index 000000000..8960dde5c --- /dev/null +++ b/tests/regression_tests/collision_track/case_4_Nuclide_ID/inputs_true.dat @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 5 + 1 + + + -2.0 -2.0 -2.0 2.0 2.0 2.0 + + + true + + + + O16 U235 + 300 + + 1 + + diff --git a/tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat b/tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat new file mode 100644 index 000000000..d4d1d1e5a --- /dev/null +++ b/tests/regression_tests/collision_track/case_4_Nuclide_ID/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat b/tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat new file mode 100644 index 000000000..8c0d7aa8e --- /dev/null +++ b/tests/regression_tests/collision_track/case_5_Universe_ID/inputs_true.dat @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 5 + 1 + + + -2.0 -2.0 -2.0 2.0 2.0 2.0 + + + true + + + + 22 + 77 + 300 + + 1 + + diff --git a/tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat b/tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat new file mode 100644 index 000000000..d4d1d1e5a --- /dev/null +++ b/tests/regression_tests/collision_track/case_5_Universe_ID/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat new file mode 100644 index 000000000..5173dc35c --- /dev/null +++ b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/inputs_true.dat @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 5 + 1 + + + -2.0 -2.0 -2.0 2.0 2.0 2.0 + + + true + + + + 550000.0 + 300 + + 1 + + diff --git a/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat new file mode 100644 index 000000000..d4d1d1e5a --- /dev/null +++ b/tests/regression_tests/collision_track/case_6_deposited_energy_threshold/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat new file mode 100644 index 000000000..005d9feb2 --- /dev/null +++ b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/inputs_true.dat @@ -0,0 +1,63 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 5 + 1 + + + -2.0 -2.0 -2.0 2.0 2.0 2.0 + + + true + + + + 22 33 + elastic 18 (n,disappear) + 77 + 1 11 + U238 U235 H1 U234 + 100000.0 + 300 + + 1 + + diff --git a/tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat new file mode 100644 index 000000000..d4d1d1e5a --- /dev/null +++ b/tests/regression_tests/collision_track/case_7_all_parameters_used_together/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/case_8_2threads/inputs_true.dat b/tests/regression_tests/collision_track/case_8_2threads/inputs_true.dat new file mode 100644 index 000000000..514932c1a --- /dev/null +++ b/tests/regression_tests/collision_track/case_8_2threads/inputs_true.dat @@ -0,0 +1,57 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 5 + 1 + + + -2.0 -2.0 -2.0 2.0 2.0 2.0 + + + true + + + + 200 + + 1 + + diff --git a/tests/regression_tests/collision_track/case_8_2threads/results_true.dat b/tests/regression_tests/collision_track/case_8_2threads/results_true.dat new file mode 100644 index 000000000..d4d1d1e5a --- /dev/null +++ b/tests/regression_tests/collision_track/case_8_2threads/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +5.642735E-02 1.494035E-02 diff --git a/tests/regression_tests/collision_track/test.py b/tests/regression_tests/collision_track/test.py new file mode 100644 index 000000000..00e1e3de4 --- /dev/null +++ b/tests/regression_tests/collision_track/test.py @@ -0,0 +1,255 @@ +"""Test the 'collision_track' setting. + +Results +------- + +All results are generated using only 1 MPI process. + +All results are generated using 1 thread except for "test_consistency_low_realization_number". +This specific test verifies that when the number of realization (i.e., point being candidate +to be stored) is lower than the capacity, results are reproducible even with multiple +threads (i.e., there is no potential thread competition that would produce different +results in that case). + +All results are generated using the history-based mode except for cases e01 to e03. + +All results are visually verified using the '_visualize.py' script in the regression test folder. + +OpenMC models +------------- + +Four OpenMC models with CSG-only geometries are used to cover the transmission, vacuum, +reflective and periodic Boundary Conditions (BC): + +- model_1: cylindrical core in 2 boxes (vacuum and transmission BC), + +# Test cases for simulation parameters using CSG-only geometries +# ============================================================ +# Each test case is defined by a combination of folder name, model name, and specific parameters. +# Below is a summary of the parameters used in the test cases: +# +# - max_collisions: Maximum number of particles to track in the simulation. +# - reactions: List of MT numbers (reaction types- 2 for scattering, 18 for fission, 101 for absorbtion). +# - cell_ids: IDs of specific cells in the model. +# - mat_ids: Material IDs for filtering particles. +# - nuclides: Nuclides for filtering particles. +# - univ_ids: Universe IDs for filtering particles. +# - E_threshold: Energy threshold for filtering particles (optional). +# +# The test cases are designed to validate the behavior of the simulation under various configurations. + +*: BC stands for Boundary Conditions, T for Transmission, R for Reflective, and V for Vacuum. + +An additional case, called 'case-a01', is used to check that the results are comparable when +the number of threads is set to 2 if the number of realization is lower than the capacity. + + +*: BC stands for Boundary Conditions, T for Transmission, and V for Vacuum. + +Notes: + +- The test cases list is non-exhaustive compared to the number of possible combinations. + Test cases have been selected based on use and internal code logic. + + + +TODO: + +- Test with a lattice. + +""" + +import os + +import openmc +import openmc.lib +import pytest + +from tests.testing_harness import CollisionTrackTestHarness +from tests.regression_tests import config + + +@pytest.fixture(scope="function") +def two_threads(monkeypatch): + """Set the number of OMP threads to 2 for the test.""" + monkeypatch.setenv("OMP_NUM_THREADS", "2") + + +@pytest.fixture(scope="function") +def single_process(monkeypatch): + """Set the number of MPI process to 1 for the test.""" + monkeypatch.setitem(config, "mpi_np", "1") + + +@pytest.fixture(scope="module") +def model_1(): + """Cylindrical core contained in a first box which is contained in a larger box. + A lower universe is used to describe the interior of the first box which + contains the core and its surrounding space. + + """ + openmc.reset_auto_ids() + model = openmc.Model() + + # ============================================================================= + # Materials + # ============================================================================= + + fuel = openmc.Material(material_id=1) + fuel.add_nuclide("U234", 0.0004524) + fuel.add_nuclide("U235", 0.0506068) + fuel.add_nuclide("U238", 0.9487090) + fuel.add_nuclide("U236", 0.0002318) + fuel.add_nuclide("O16", 2.0) + fuel.set_density("g/cm3", 11.0) + + water = openmc.Material(material_id=11) + water.add_nuclide("H1", 2.0) + water.add_nuclide("O16", 1.0) + water.set_density("g/cm3", 1.0) + + # ============================================================================= + # Geometry + # ============================================================================= + + # ----------------------------------------------------------------------------- + # Cylindrical core + # ----------------------------------------------------------------------------- + + # Parameters + core_radius = 2.0 + core_height = 4.0 + + # Surfaces + core_cylinder = openmc.ZCylinder(r=core_radius) + core_lower_plane = openmc.ZPlane(-core_height / 2.0) + core_upper_plane = openmc.ZPlane(core_height / 2.0) + + # Region + core_region = -core_cylinder & +core_lower_plane & -core_upper_plane + + # Cells + core = openmc.Cell(fill=fuel, region=core_region, cell_id=22) + outside_core_region = +core_cylinder | -core_lower_plane | +core_upper_plane + outside_core = openmc.Cell( + fill=water, region=outside_core_region, cell_id=33) + + # Universe + inside_box1_universe = openmc.Universe( + cells=[core, outside_core], universe_id=77) + + # ----------------------------------------------------------------------------- + # Box 1 + # ----------------------------------------------------------------------------- + + # Parameters + box1_size = 6.0 + + # Surfaces + box1_rpp = openmc.model.RectangularParallelepiped( + -box1_size / 2.0, box1_size / 2.0, + -box1_size / 2.0, box1_size / 2.0, + -box1_size / 2.0, box1_size / 2.0, + ) + + # Cell + box1 = openmc.Cell(fill=inside_box1_universe, region=-box1_rpp, cell_id=5) + + # ----------------------------------------------------------------------------- + # Box 2 + # ----------------------------------------------------------------------------- + + # Parameters + box2_size = 8 + + # Surfaces + box2_rpp = openmc.model.RectangularParallelepiped( + -box2_size / 2.0, box2_size / 2.0, + -box2_size / 2.0, box2_size / 2.0, + -box2_size / 2.0, box2_size / 2.0, + boundary_type="vacuum" + ) + + # Cell + box2 = openmc.Cell(fill=water, region=-box2_rpp & +box1_rpp, cell_id=8) + + # Register geometry + model.geometry = openmc.Geometry([box1, box2]) + + # ============================================================================= + # Settings + # ============================================================================= + + model.settings = openmc.Settings() + model.settings.particles = 100 + model.settings.batches = 5 + model.settings.inactive = 1 + model.settings.seed = 1 + + bounds = [ + -core_radius, + -core_radius, + -core_height / 2.0, + core_radius, + core_radius, + core_height / 2.0, + ] + distribution = openmc.stats.Box(bounds[:3], bounds[3:]) + model.settings.source = openmc.IndependentSource( + space=distribution, constraints={'fissionable': True}) + + return model + + +@pytest.mark.parametrize( + "folder, model_name, parameter", + [("case_1_Reactions", "model_1", {"max_collisions": 300, "reactions": ["(n,fission)", 101]}), + ("case_2_Cell_ID", "model_1", { + "max_collisions": 300, "cell_ids": [22]}), + ("case_3_Material_ID", "model_1", { + "max_collisions": 300, "material_ids": [1]}), + ("case_4_Nuclide_ID", "model_1", { + "max_collisions": 300, "nuclides": ["O16", "U235"]}), + ("case_5_Universe_ID", "model_1", { + "max_collisions": 300, "cell_ids": [22], "universe_ids": [77]}), + ("case_6_deposited_energy_threshold", "model_1", { + "max_collisions": 300, "deposited_E_threshold": 5.5e5}), + ("case_7_all_parameters_used_together", "model_1", { + "max_collisions": 300, + "reactions": ["elastic", 18, "(n,disappear)"], + "material_ids": [1, 11], + "universe_ids": [77], + "nuclides": ["U238", "U235", "H1", "U234"], + "cell_ids": [22, 33], + "deposited_E_threshold": 1e5}) + ], +) +def test_collision_track_several_cases( + folder, model_name, parameter, request +): + # Since for these tests the actual number of collisions recorded is < max_collisions, + # we can run them with 1 or 2 threads, and in history or event mode. + model = request.getfixturevalue(model_name) + model.settings.collision_track = parameter + harness = CollisionTrackTestHarness( + "statepoint.5.h5", model=model, workdir=folder + ) + harness.main() + + +@pytest.mark.skipif(config["event"], reason="Results from history-based mode.") +def test_collision_track_2threads(model_1, two_threads, single_process): + # This test checks that the `max_collisions` setting is honored: + # no collisions beyond the specified limit should be recorded. + # + # For the result to be reproducible, the number of threads and + # the transport mode (history vs. event) must remain fixed. + assert os.environ["OMP_NUM_THREADS"] == "2" + assert config["mpi_np"] == "1" + model_1.settings.collision_track = { + "max_collisions": 200 + } + harness = CollisionTrackTestHarness( + "statepoint.5.h5", model=model_1, workdir="case_8_2threads" + ) + harness.main() diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 11ced5b38..1ad91b7a8 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -68,7 +68,7 @@ class TestHarness: if config['mpi']: mpi_args = [config['mpiexec'], '-n', config['mpi_np']] openmc.run(openmc_exec=config['exe'], mpi_args=mpi_args, - event_based=config['event']) + event_based=config['event']) else: openmc.run(openmc_exec=config['exe'], event_based=config['event']) @@ -305,9 +305,12 @@ class PyAPITestHarness(TestHarness): else: self.execute_test() - def execute_test(self): + def execute_test(self, change_dir=False): """Build input XMLs, run OpenMC, and verify correct results.""" + base_dir = os.getcwd() if change_dir else None try: + if change_dir: + os.chdir(self.workdir) self._build_inputs() inputs = self._get_inputs() self._write_inputs(inputs) @@ -319,10 +322,15 @@ class PyAPITestHarness(TestHarness): self._compare_results() finally: self._cleanup() + if base_dir: + os.chdir(base_dir) - def update_results(self): + def update_results(self, change_dir=False): """Update results_true.dat and inputs_true.dat""" + base_dir = os.getcwd() if change_dir else None try: + if change_dir: + os.chdir(self.workdir) self._build_inputs() inputs = self._get_inputs() self._write_inputs(inputs) @@ -334,6 +342,8 @@ class PyAPITestHarness(TestHarness): self._overwrite_results() finally: self._cleanup() + if base_dir: + os.chdir(base_dir) def _build_inputs(self): """Write input XML files.""" @@ -371,7 +381,8 @@ class PyAPITestHarness(TestHarness): """Delete XMLs, statepoints, tally, and test files.""" super()._cleanup() output = ['materials.xml', 'geometry.xml', 'settings.xml', - 'tallies.xml', 'plots.xml', 'inputs_test.dat', 'model.xml'] + 'tallies.xml', 'plots.xml', 'inputs_test.dat', 'model.xml', + 'collision_track.h5', 'collision_track.mcpl'] for f in output: if os.path.exists(f): os.remove(f) @@ -389,6 +400,7 @@ class TolerantPyAPITestHarness(PyAPITestHarness): due to single precision usage (e.g., as in the random ray solver). """ + def _are_files_equal(self, actual_path, expected_path, tolerance): def isfloat(value): try: @@ -428,7 +440,8 @@ class TolerantPyAPITestHarness(PyAPITestHarness): def _compare_results(self): """Make sure the current results agree with the reference.""" - compare = self._are_files_equal('results_test.dat', 'results_true.dat', 1e-6) + compare = self._are_files_equal( + 'results_test.dat', 'results_true.dat', 1e-6) if not compare: expected = open('results_true.dat').readlines() actual = open('results_test.dat').readlines() @@ -476,6 +489,7 @@ class WeightWindowPyAPITestHarness(PyAPITestHarness): class PlotTestHarness(TestHarness): """Specialized TestHarness for running OpenMC plotting tests.""" + def __init__(self, plot_names, voxel_convert_checks=[]): super().__init__(None) self._plot_names = plot_names @@ -523,3 +537,105 @@ class PlotTestHarness(TestHarness): outstr = sha512.hexdigest() return outstr + + +class CollisionTrackTestHarness(PyAPITestHarness): + def __init__(self, statepoint_name, model=None, inputs_true=None, workdir=None): + super().__init__(statepoint_name, model, inputs_true) + self.workdir = workdir + + def _test_output_created(self): + """Make sure collision_track.h5 has also been created.""" + super()._test_output_created() + if self._model.settings.collision_track: + assert os.path.exists( + "collision_track.h5" + ), "collision_track file has not been created." + + def _compare_output(self): + """Compare collision_track.h5 files.""" + if self._model.settings.collision_track: + collision_track_true = self._return_collision_track_data( + "collision_track_true.h5") + collision_track_test = self._return_collision_track_data( + "collision_track.h5") + np.testing.assert_allclose( + collision_track_true, collision_track_test, rtol=1e-07) + + def main(self): + """Accept commandline arguments and either run or update tests.""" + if config["build_inputs"]: + self.build_inputs() + elif config["update"]: + self.update_results(change_dir=True) + else: + self.execute_test(change_dir=True) + + def build_inputs(self): + """Build inputs.""" + base_dir = os.getcwd() + try: + os.chdir(self.workdir) + self._build_inputs() + finally: + os.chdir(base_dir) + + def _overwrite_results(self): + """Also add the 'collision_track.h5' file during overwriting.""" + super()._overwrite_results() + if os.path.exists("collision_track.h5"): + shutil.copyfile("collision_track.h5", "collision_track_true.h5") + + @staticmethod + def _return_collision_track_data(filepath): + """ + Read a collision_track file and return a sorted array composed + of flatten arrays of collision information. + + Parameters + ---------- + filepath : str + Path to the collision_track file + + Returns + ------- + data : np.array + Sorted array composed of flatten arrays of collision_track data for + each collision information + """ + data = [] + keys = [] + + # Read source file + source = openmc.read_collision_track_file(filepath) + for src in source: + r = src['r'] + u = src['u'] + e = src['E'] + de = src['dE'] + time = src['time'] + wgt = src['wgt'] + delayed_group = src['delayed_group'] + cell_id = src['cell_id'] + nuclide_id = src['nuclide_id'] + material_id = src['material_id'] + universe_id = src['universe_id'] + n_collision = src['n_collision'] + event_mt = src['event_mt'] + key = ( + f"{r[0]:.10e} {r[1]:.10e} {r[2]:.10e} {u[0]:.10e} {u[1]:.10e} {u[2]:.10e}" + f"{e:.10e} {de:.10e} {time:.10e} {wgt:.10e} {event_mt} {delayed_group} {cell_id}" + f"{nuclide_id} {material_id} {universe_id} {n_collision} " + ) + keys.append(key) + values = [*r, *u, e, de, time, wgt, event_mt, + delayed_group, cell_id, nuclide_id, material_id, + universe_id, n_collision] + assert len(values) == 17 + data.append(values) + + data = np.array(data) + keys = np.array(keys) + sorted_idx = np.argsort(keys, kind='stable') + + return data[sorted_idx] diff --git a/tests/unit_tests/test_collision_track.py b/tests/unit_tests/test_collision_track.py new file mode 100644 index 000000000..9bc6a8c15 --- /dev/null +++ b/tests/unit_tests/test_collision_track.py @@ -0,0 +1,127 @@ +"""Test the 'collision_track' setting used to store particle information +during specified collision conditions in a file for a given simulation.""" + +import openmc +import pytest +import h5py +import numpy as np + +from tests.testing_harness import CollisionTrackTestHarness as ctt + + +@pytest.fixture(scope="module") +def geometry(): + """Simple hydrogen sphere geometry""" + openmc.reset_auto_ids() + material = openmc.Material(name="H1") + material.add_element("H", 1.0) + sphere = openmc.Sphere(r=1.0, boundary_type="vacuum") + cell = openmc.Cell(region=-sphere, fill=material) + return openmc.Geometry([cell]) + + +@pytest.mark.parametrize( + "parameter", + [ + {"max_collisions": 200}, + {"max_collisions": 200, "reactions": ["(n,disappear)"]}, + {"max_collisions": 200, "cell_ids": [1]}, + {"max_collisions": 200, "material_ids": [1]}, + {"max_collisions": 200, "universe_ids": [1]}, + {"max_collisions": 200, "nuclides": ["H1"]}, + {"max_collisions": 200, "deposited_E_threshold": 200000.0}, + {"max_collisions": 200, "mcpl": True} + + ], +) +def test_xml_serialization(parameter, run_in_tmpdir): + """Check that the different use cases can be written and read in XML.""" + settings = openmc.Settings() + settings.collision_track = parameter + settings.export_to_xml() + + read_settings = openmc.Settings.from_xml() + assert read_settings.collision_track == parameter + + +@pytest.fixture(scope="module") +def model(): + """Simple hydrogen sphere divided in two hemispheres + by a z-plane to form 2 cells.""" + openmc.reset_auto_ids() + model = openmc.Model() + + # Material + material = openmc.Material(name="H1") + material.add_element("H", 1.0) + + # Geometry + radius = 1.0 + sphere = openmc.Sphere(r=radius, boundary_type="reflective") + plane = openmc.ZPlane(0.0) + cell_1 = openmc.Cell(region=-sphere & -plane, fill=material, cell_id=1) + cell_2 = openmc.Cell(region=-sphere & +plane, fill=material, cell_id=2) + root = openmc.Universe(cells=[cell_1, cell_2]) + model.geometry = openmc.Geometry(root) + + # Settings + model.settings = openmc.Settings() + model.settings.run_mode = "fixed source" + model.settings.particles = 1 + model.settings.batches = 1 + model.settings.seed = 2 + + bounds = [-radius, -radius, -radius, radius, radius, radius] + distribution = openmc.stats.Box(bounds[:3], bounds[3:]) + model.settings.source = openmc.IndependentSource(space=distribution) + + return model + + +def test_particle_location(run_in_tmpdir, model): + """Test the location of particles with respected to the "cell_ids" + and the location x, y, z of the particle itself. the upper sphere will + have positive z component and the bottom sphere a negative z compnent. + + """ + model.settings.collision_track = { + "max_collisions": 200, + "reactions": ["elastic"], + "cell_ids": [1, 2] + } + model.run() + + with h5py.File("collision_track.h5", "r") as f: + source = f["collision_track_bank"] + + assert len(source) == 60 + + # We want to verify that the collisions happenening are in the right cells + # and the position of the particle is either positive or negative relative + # to the z plane. In this case, we track the position of the particle + # relative to the cell_id already set. + for point in source: + if point['cell_id'] == 1: + assert point['r'][2] < 0.0 # z component negative + elif point['cell_id'] == 2: + assert point["r"][2] > 0.0 # z component positive + else: + assert False + + +def test_format_similarity(run_in_tmpdir, model): + model.settings.collision_track = {"max_collisions": 200, "reactions": ['elastic'], + "cell_ids": [1, 2], "mcpl": False} + model.run() + data_h5 = ctt._return_collision_track_data('collision_track.h5') + + model.settings.collision_track["mcpl"] = True + model.run() + data_mcpl = ctt._return_collision_track_data('collision_track.mcpl') + + assert len(data_h5) == 60 + assert len(data_mcpl) == 60 + + np.testing.assert_allclose(data_h5, data_mcpl, rtol=1e-05) + # tolerance not that low due to the strings that is saved in MCPL, + # not enough precision!