diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ccc457f7..021b274c2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -253,6 +253,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_cellborn.cpp src/tallies/filter_cellfrom.cpp src/tallies/filter_cell.cpp + src/tallies/filter_cell_instance.cpp src/tallies/filter_delayedgroup.cpp src/tallies/filter_distribcell.cpp src/tallies/filter_energyfunc.cpp diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 3a5e347f9..6ec067a12 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -195,6 +195,18 @@ based on the recommended value in LA-UR-14-24530_. .. _LA-UR-14-24530: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf +--------------------------- +```` +--------------------------- + +By default, OpenMC will count the number of instances of each cell filled with a +material and generate "offset tables" that are used for cell instance tallies. +The ```` element allows a user to override this default +setting and turn off the generation of offset tables, if desired, by setting it +to false. + + *Default*: true + --------------------------- ```` Element --------------------------- @@ -423,16 +435,16 @@ attributes/sub-elements: :type: The type of spatial distribution. Valid options are "box", "fission", - "point", "cartesian", "cylindrical", and "spherical". A "box" spatial - distribution has coordinates sampled uniformly in a parallelepiped. A - "fission" spatial distribution samples locations from a "box" - distribution but only locations in fissionable materials are accepted. - A "point" spatial distribution has coordinates specified by a triplet. - A "cartesian" spatial distribution specifies independent distributions of + "point", "cartesian", "cylindrical", and "spherical". A "box" spatial + distribution has coordinates sampled uniformly in a parallelepiped. A + "fission" spatial distribution samples locations from a "box" + distribution but only locations in fissionable materials are accepted. + A "point" spatial distribution has coordinates specified by a triplet. + A "cartesian" spatial distribution specifies independent distributions of x-, y-, and z-coordinates. A "cylindrical" spatial distribution specifies independent distributions of r-, phi-, and z-coordinates where phi is the azimuthal angle and the origin for the cylindrical coordinate system is - specified by origin. A "spherical" spatial distribution specifies + specified by origin. A "spherical" spatial distribution specifies independent distributions of r-, theta-, and phi-coordinates where theta is the angle with respect to the z-axis, phi is the azimuthal angle, and the sphere is centered on the coordinate (x0,y0,z0). @@ -452,7 +464,7 @@ attributes/sub-elements: For an "cartesian" distribution, no parameters are specified. Instead, the ``x``, ``y``, and ``z`` elements must be specified. - + For a "cylindrical" distribution, no parameters are specified. Instead, the ``r``, ``phi``, ``z``, and ``origin`` elements must be specified. @@ -474,15 +486,15 @@ attributes/sub-elements: :ref:`univariate`). :z: - For both "cartesian" and "cylindrical" distributions, this element - specifies the distribution of z-coordinates. The necessary - sub-elements/attributes are those of a univariate probability + For both "cartesian" and "cylindrical" distributions, this element + specifies the distribution of z-coordinates. The necessary + sub-elements/attributes are those of a univariate probability distribution (see the description in :ref:`univariate`). :r: For "cylindrical" and "spherical" distributions, this element specifies - the distribution of r-coordinates (cylindrical radius and spherical - radius, respectively). The necessary sub-elements/attributes are those + the distribution of r-coordinates (cylindrical radius and spherical + radius, respectively). The necessary sub-elements/attributes are those of a univariate probability distribution (see the description in :ref:`univariate`). @@ -493,13 +505,13 @@ attributes/sub-elements: :ref:`univariate`). :phi: - For "cylindrical" and "spherical" distributions, this element specifies - the distribution of phi-coordinates. The necessary - sub-elements/attributes are those of a univariate probability + For "cylindrical" and "spherical" distributions, this element specifies + the distribution of phi-coordinates. The necessary + sub-elements/attributes are those of a univariate probability distribution (see the description in :ref:`univariate`). :origin: - For "cylindrical and "spherical" distributions, this element specifies + For "cylindrical and "spherical" distributions, this element specifies the coordinates for the origin of the coordinate system. :angle: diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 2c0ac8797..e7d139938 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -107,6 +107,7 @@ Constructing Tallies openmc.CellFilter openmc.CellFromFilter openmc.CellbornFilter + openmc.CellInstanceFilter openmc.SurfaceFilter openmc.MeshFilter openmc.MeshSurfaceFilter diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 8a330060c..3df2f40d5 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -2,12 +2,14 @@ #define OPENMC_CELL_H #include +#include // for hash #include #include // for unique_ptr #include #include #include +#include #include "hdf5.h" #include "pugixml.hpp" #include "dagmc.h" @@ -286,6 +288,26 @@ private: std::vector> partitions_; }; +//============================================================================== +//! Define an instance of a particular cell +//============================================================================== + +struct CellInstance { + //! Check for equality + bool operator==(const CellInstance& other) const + { return index_cell == other.index_cell && instance == other.instance; } + + gsl::index index_cell; + gsl::index instance; +}; + +struct CellInstanceHash { + std::size_t operator()(const CellInstance& k) const + { + return 4096*k.index_cell + k.instance; + } +}; + //============================================================================== // Non-member functions //============================================================================== diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 45cb7394c..6570a25f7 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -31,6 +31,7 @@ extern "C" bool cmfd_run; //!< is a CMFD run? extern "C" bool dagmc; //!< indicator of DAGMC geometry extern "C" bool entropy_on; //!< calculate Shannon entropy? extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular? +extern bool material_cell_offsets; //!< create material cells offsets? extern "C" bool output_summary; //!< write summary.h5? extern bool output_tallies; //!< write tallies.out? extern bool particle_restart_run; //!< particle restart run? diff --git a/include/openmc/tallies/filter_cell_instance.h b/include/openmc/tallies/filter_cell_instance.h new file mode 100644 index 000000000..1b2cdd425 --- /dev/null +++ b/include/openmc/tallies/filter_cell_instance.h @@ -0,0 +1,63 @@ +#ifndef OPENMC_TALLIES_FILTER_CELL_INSTANCE_H +#define OPENMC_TALLIES_FILTER_CELL_INSTANCE_H + +#include +#include +#include + +#include + +#include "openmc/cell.h" +#include "openmc/tallies/filter.h" + + +namespace openmc { + +//============================================================================== +//! Specifies cell instances that tally events reside in. +//============================================================================== + +class CellInstanceFilter : public Filter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + CellInstanceFilter() = default; + CellInstanceFilter(gsl::span instances); + ~CellInstanceFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type() const override {return "cellinstance";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //---------------------------------------------------------------------------- + // Accessors + + const std::vector& cell_instances() const { return cell_instances_; } + + void set_cell_instances(gsl::span instances); + +private: + //---------------------------------------------------------------------------- + // Data members + + //! The indices of the cells binned by this filter. + std::vector cell_instances_; + + //! A map from cell/instance indices to filter bin indices. + std::unordered_map map_; +}; + +} // namespace openmc + +#endif // OPENMC_TALLIES_FILTER_CELL_INSTANCE_H diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 71c83b5e8..aebade144 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -62,6 +62,8 @@ public: //---------------------------------------------------------------------------- // Other methods. + void add_filter(Filter* filter) { set_filters({&filter, 1}); } + void init_triggers(pugi::xml_node node); void init_results(); diff --git a/openmc/filter.py b/openmc/filter.py index 8527ca9b9..1988aee81 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -23,7 +23,7 @@ _FILTER_TYPES = ( 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre', - 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle' + 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance' ) _CURRENT_NAMES = ( @@ -266,7 +266,7 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): # Merge unique filter bins merged_bins = np.concatenate((self.bins, other.bins)) - merged_bins = np.unique(merged_bins) + merged_bins = np.unique(merged_bins, axis=0) # Create a new filter with these bins and a new auto-generated ID return type(self)(merged_bins) @@ -518,6 +518,105 @@ class CellbornFilter(WithIDFilter): expected_type = Cell +class CellInstanceFilter(Filter): + """Bins tally events based on which cell instance a particle is in. + + This filter is similar to :class:`DistribcellFilter` but allows one to + select particular instances to be tallied (instead of obtaining *all* + instances by default) and allows instances from different cells to be + specified in a single filter. + + Parameters + ---------- + bins : iterable of 2-tuples or numpy.ndarray + The cell instances to tally, given as 2-tuples. For the first value in + the tuple, either openmc.Cell objects or their integral ID numbers can + be used. + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + bins : numpy.ndarray + 2D numpy array of cell IDs and instances + id : int + Unique identifier for the filter + num_bins : Integral + The number of filter bins + + See Also + -------- + DistribcellFilter + + """ + def __init__(self, bins, filter_id=None): + self.bins = bins + self.id = filter_id + + @Filter.bins.setter + def bins(self, bins): + pairs = np.empty((len(bins), 2), dtype=int) + for i, (cell, instance) in enumerate(bins): + cv.check_type('cell', cell, (openmc.Cell, Integral)) + cv.check_type('instance', instance, Integral) + pairs[i, 0] = cell if isinstance(cell, Integral) else cell.id + pairs[i, 1] = instance + self._bins = pairs + + def get_pandas_dataframe(self, data_size, stride, **kwargs): + """Builds a Pandas DataFrame for the Filter's bins. + + This method constructs a Pandas DataFrame object for the filter with + columns annotated by filter bin information. This is a helper method for + :meth:`Tally.get_pandas_dataframe`. + + Parameters + ---------- + data_size : int + The total number of bins in the tally corresponding to this filter + stride : int + Stride in memory for the filter + + Returns + ------- + pandas.DataFrame + A Pandas DataFrame with a multi-index column for the cell instance. + The number of rows in the DataFrame is the same as the total number + of bins in the corresponding tally, with the filter bin appropriately + tiled to map to the corresponding tally bins. + + See also + -------- + Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe() + + """ + # Repeat and tile bins as necessary to account for other filters. + bins = np.repeat(self.bins, stride, axis=0) + tile_factor = data_size // len(bins) + bins = np.tile(bins, (tile_factor, 1)) + + columns = pd.MultiIndex.from_product([[self.short_name.lower()], + ['cell', 'instance']]) + return pd.DataFrame(bins, columns=columns) + + def to_xml_element(self): + """Return XML Element representing the Filter. + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing filter data + + """ + element = ET.Element('filter') + element.set('id', str(self.id)) + element.set('type', self.short_name.lower()) + + subelement = ET.SubElement(element, 'bins') + subelement.text = ' '.join(str(i) for i in self.bins.ravel()) + return element + + class SurfaceFilter(WithIDFilter): """Filters particles by surface crossing @@ -564,11 +663,7 @@ class ParticleFilter(Filter): The number of filter bins """ - @property - def bins(self): - return self._bins - - @bins.setter + @Filter.bins.setter def bins(self, bins): bins = np.atleast_1d(bins) cv.check_iterable_type('filter bins', bins, str) @@ -1174,6 +1269,11 @@ def _path_to_levels(path): class DistribcellFilter(Filter): """Bins tally event locations on instances of repeated cells. + This filter provides a separate score for each unique instance of a repeated + cell in a geometry. Note that only one cell can be specified in this filter. + The related :class:`CellInstanceFilter` allows one to obtain scores for + particular cell instances as well as instances from different cells. + Parameters ---------- cell : openmc.Cell or Integral @@ -1194,6 +1294,10 @@ class DistribcellFilter(Filter): The paths traversed through the CSG tree to reach each distribcell instance (for 'distribcell' filters only) + See Also + -------- + CellInstanceFilter + """ def __init__(self, cell, filter_id=None): diff --git a/openmc/settings.py b/openmc/settings.py index 784c2ac6e..f25fd681d 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -61,6 +61,9 @@ class Settings(object): relative error used. log_grid_bins : int Number of bins for logarithmic energy grid search + material_cell_offsets : bool + Generate an "offset table" for material cells by default. These tables + are necessary when a particular instance of a cell needs to be tallied. max_order : None or int Maximum scattering order to apply globally when in multi-group mode. no_reduce : bool @@ -216,6 +219,7 @@ class Settings(object): VolumeCalculation, 'volume calculations') self._create_fission_neutrons = None + self._material_cell_offsets = None self._log_grid_bins = None self._dagmc = False @@ -352,6 +356,10 @@ class Settings(object): def create_fission_neutrons(self): return self._create_fission_neutrons + @property + def material_cell_offsets(self): + return self._material_cell_offsets + @property def log_grid_bins(self): return self._log_grid_bins @@ -683,6 +691,11 @@ class Settings(object): create_fission_neutrons, bool) self._create_fission_neutrons = create_fission_neutrons + @material_cell_offsets.setter + def material_cell_offsets(self, value): + cv.check_type('material cell offsets', value, bool) + self._material_cell_offsets = value + @log_grid_bins.setter def log_grid_bins(self, log_grid_bins): cv.check_type('log grid bins', log_grid_bins, Real) @@ -917,6 +930,11 @@ class Settings(object): elem = ET.SubElement(root, "create_fission_neutrons") elem.text = str(self._create_fission_neutrons).lower() + def _create_material_cell_offsets_subelement(self, root): + if self._material_cell_offsets is not None: + elem = ET.SubElement(root, "material_cell_offsets") + elem.text = str(self._material_cell_offsets).lower() + def _create_log_grid_bins_subelement(self, root): if self._log_grid_bins is not None: elem = ET.SubElement(root, "log_grid_bins") @@ -1148,6 +1166,11 @@ class Settings(object): if text is not None: self.create_fission_neutrons = text in ('true', '1') + def _material_cell_offsets_from_xml_element(self, root): + text = get_text(root, 'material_cell_offsets') + if text is not None: + self.material_cell_offsets = text in ('true', '1') + def _log_grid_bins_from_xml_element(self, root): text = get_text(root, 'log_grid_bins') if text is not None: @@ -1202,6 +1225,7 @@ class Settings(object): self._create_resonance_scattering_subelement(root_element) self._create_volume_calcs_subelement(root_element) self._create_create_fission_neutrons_subelement(root_element) + self._create_material_cell_offsets_subelement(root_element) self._create_log_grid_bins_subelement(root_element) self._create_dagmc_subelement(root_element) @@ -1267,6 +1291,7 @@ class Settings(object): settings._ufs_mesh_from_xml_element(root) settings._resonance_scattering_from_xml_element(root) settings._create_fission_neutrons_from_xml_element(root) + settings._material_cell_offsets_from_xml_element(root) settings._log_grid_bins_from_xml_element(root) settings._dagmc_from_xml_element(root) diff --git a/src/cell.cpp b/src/cell.cpp index 94bae32eb..9496dbba2 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -258,8 +258,13 @@ Cell::set_temperature(double T, int32_t instance) } if (instance >= 0) { + // If temperature vector is not big enough, resize it first + if (sqrtkT_.size() != n_instances_) sqrtkT_.resize(n_instances_, sqrtkT_[0]); + + // Set temperature for the corresponding instance sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T); } else { + // Set temperature for all instances for (auto& T_ : sqrtkT_) { T_ = std::sqrt(K_BOLTZMANN * T); } diff --git a/src/finalize.cpp b/src/finalize.cpp index 6169e2819..d359f3885 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -75,6 +75,7 @@ int openmc_finalize() settings::gen_per_batch = 1; settings::legendre_to_tabular = true; settings::legendre_to_tabular_points = -1; + settings::material_cell_offsets = true; settings::n_particles = -1; settings::output_summary = true; settings::output_tallies = true; diff --git a/src/geometry.cpp b/src/geometry.cpp index 88b48c866..6bcfc9c58 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -132,8 +132,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) //! Found a material cell which means this is the lowest coord level. // Find the distribcell instance number. - if (c.material_.size() > 1 || c.sqrtkT_.size() > 1) { - int offset = 0; + int offset = 0; + if (c.distribcell_index_ >= 0) { for (int i = 0; i < p->n_coord_; i++) { const auto& c_i {*model::cells[p->coord_[i].cell]}; if (c_i.type_ == FILL_UNIVERSE) { @@ -148,10 +148,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) } } } - p->cell_instance_ = offset; - } else { - p->cell_instance_ = 0; } + p->cell_instance_ = offset; // Set the material and temperature. p->material_last_ = p->material_; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 5ebee4aca..011dc2075 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -18,6 +18,7 @@ #include "openmc/settings.h" #include "openmc/surface.h" #include "openmc/tallies/filter.h" +#include "openmc/tallies/filter_cell_instance.h" #include "openmc/tallies/filter_distribcell.h" @@ -319,7 +320,9 @@ find_root_universe() void prepare_distribcell() { - // Find all cells listed in a DistribcellFilter. + write_message("Preparing distributed cell instances...", 5); + + // Find all cells listed in a DistribcellFilter or CellInstanceFilter std::unordered_set distribcells; for (auto& filt : model::tally_filters) { auto* distrib_filt = dynamic_cast(filt.get()); @@ -328,8 +331,15 @@ prepare_distribcell() } } - // Find all cells with distributed materials or temperatures. Make sure that - // the number of materials/temperatures matches the number of cell instances. + // By default, add material cells to the list of distributed cells + if (settings::material_cell_offsets) { + for (gsl::index i = 0; i < model::cells.size(); ++i) { + if (model::cells[i]->type_ == FILL_MATERIAL) distribcells.insert(i); + } + } + + // Make sure that the number of materials/temperatures matches the number of + // cell instances. for (int i = 0; i < model::cells.size(); i++) { Cell& c {*model::cells[i]}; @@ -342,7 +352,6 @@ prepare_distribcell() "one or the number of instances."; fatal_error(err_msg); } - distribcells.insert(i); } if (c.sqrtkT_.size() > 1) { @@ -354,20 +363,18 @@ prepare_distribcell() "one or the number of instances."; fatal_error(err_msg); } - distribcells.insert(i); } } - // Search through universes for distributed cells and assign each one a + // Search through universes for material cells and assign each one a // unique distribcell array index. int distribcell_index = 0; std::vector target_univ_ids; for (const auto& u : model::universes) { - for (auto cell_indx : u->cells_) { - if (distribcells.find(cell_indx) != distribcells.end()) { - model::cells[cell_indx]->distribcell_index_ = distribcell_index; + for (auto idx : u->cells_) { + if (distribcells.find(idx) != distribcells.end()) { + model::cells[idx]->distribcell_index_ = distribcell_index++; target_univ_ids.push_back(u->id_); - ++distribcell_index; } } } @@ -387,7 +394,7 @@ prepare_distribcell() for (int map = 0; map < target_univ_ids.size(); map++) { auto target_univ_id = target_univ_ids[map]; for (const auto& univ : model::universes) { - int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation. + int32_t offset = 0; for (int32_t cell_indx : univ->cells_) { Cell& c = *model::cells[cell_indx]; diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index eb27d8591..c69624c53 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -7,9 +7,15 @@ element settings { element cutoff { (element weight { xsd:double } | attribute weight { xsd:double })? & - (element weight_avg { xsd:double } | attribute weight_avg { xsd:double })? + (element weight_avg { xsd:double } | attribute weight_avg { xsd:double })? & + (element energy_neutron { xsd:double } | attribute energy_neutron { xsd:double })? & + (element energy_photon { xsd:double } | attribute energy_photon { xsd:double })? & + (element energy_electron { xsd:double } | attribute energy_electron { xsd:double })? & + (element energy_positron { xsd:double } | attribute energy_positron { xsd:double })? }? & + element electron_treatment { ( "led" | "ttb" ) }? & + element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" | "material-union" | "union" ) }? & element energy_mode { ( "continuous-energy" | "ce" | "CE" | "multi-group" | "mg" | "MG" ) }? & @@ -27,6 +33,8 @@ element settings { element log_grid_bins { xsd:positiveInteger }? & + element material_cell_offsets { xsd:boolean }? & + element max_order { xsd:nonNegativeInteger }? & element mesh { @@ -55,6 +63,8 @@ element settings { element particles { xsd:positiveInteger }? & + element photon_transport { xsd:boolean }? & + element ptables { xsd:boolean }? & element dagmc { xsd:boolean }? & @@ -78,8 +88,8 @@ element settings { element z { distribution }? & element r { distribution }? & element theta { distribution }? & - element phi { distribution }? & - element origin { list { xsd:double, xsd:double, xsd:double } }? + element phi { distribution }? & + element origin { list { xsd:double, xsd:double, xsd:double } }? }? & element angle { (element type { xsd:string } | attribute type { xsd:string }) & @@ -88,8 +98,7 @@ element settings { element mu { distribution }? & element phi { distribution }? }? & - element energy { distribution }? & - (element write_initial { xsd:boolean } | attribute write_initial { xsd:boolean })? + element energy { distribution }? distribution = (element type { xsd:string { maxLength = "16" } } | attribute type { xsd:string { maxLength = "16" } }) & @@ -165,6 +174,8 @@ element settings { attribute upper_right { list { xsd:double+ } }) }* & + element write_initial_source { xsd:boolean }? & + element resonance_scattering { (element enable { xsd:boolean } | attribute enable { xsd:boolean })? & (element method { xsd:string } | attribute method { xsd:string })? & diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 0c825d0d8..46845b094 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -39,9 +39,57 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + led + ttb + + + @@ -108,6 +156,11 @@ + + + + + @@ -249,6 +302,11 @@ + + + + + @@ -422,16 +480,6 @@ - - - - - - - - - - @@ -743,6 +791,11 @@ + + + + + diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index ec511481f..3d504a66c 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -51,11 +51,11 @@ element tallies { ( (element type { ( "cell" | "cellfrom" | "cellborn" | "material" | "universe" | "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" | - "energyfunction" | "meshsurface") } | + "energyfunction" | "meshsurface" | "cellinstance") } | attribute type { ( "cell" | "cellfrom" | "cellborn" | "material" | "universe" | "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" | - "energyfunction" | "meshsurface") }) & + "energyfunction" | "meshsurface" | "cellinstance") }) & (element bins { list { xsd:double+ } } | attribute bins { list { xsd:double+ } }) ) | diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index 98a48eeb8..350c3cb4a 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -245,6 +245,7 @@ delayedgroup energyfunction meshsurface + cellinstance @@ -265,6 +266,7 @@ delayedgroup energyfunction meshsurface + cellinstance diff --git a/src/settings.cpp b/src/settings.cpp index b647d32e8..59d2479be 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -46,6 +46,7 @@ bool create_fission_neutrons {true}; bool dagmc {false}; bool entropy_on {false}; bool legendre_to_tabular {true}; +bool material_cell_offsets {true}; bool output_summary {true}; bool output_tallies {true}; bool particle_restart_run {false}; @@ -776,6 +777,11 @@ void read_settings_xml() create_fission_neutrons = get_node_value_bool(root, "create_fission_neutrons"); } } + + // Check whether material cell offsets should be generated + if (check_for_node(root, "material_cell_offsets")) { + material_cell_offsets = get_node_value_bool(root, "material_cell_offsets"); + } } void free_memory_settings() { diff --git a/src/simulation.cpp b/src/simulation.cpp index 2f3230c6e..eaed2badf 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -5,6 +5,7 @@ #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/geometry_aux.h" #include "openmc/material.h" #include "openmc/message_passing.h" #include "openmc/nuclide.h" diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index d8fe58cc7..b1b6f4d61 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -12,6 +12,7 @@ #include "openmc/tallies/filter_cell.h" #include "openmc/tallies/filter_cellborn.h" #include "openmc/tallies/filter_cellfrom.h" +#include "openmc/tallies/filter_cell_instance.h" #include "openmc/tallies/filter_delayedgroup.h" #include "openmc/tallies/filter_distribcell.h" #include "openmc/tallies/filter_energyfunc.h" @@ -100,6 +101,8 @@ Filter* Filter::create(const std::string& type, int32_t id) model::tally_filters.push_back(std::make_unique()); } else if (type == "cellfrom") { model::tally_filters.push_back(std::make_unique()); + } else if (type == "cellinstance") { + model::tally_filters.push_back(std::make_unique()); } else if (type == "distribcell") { model::tally_filters.push_back(std::make_unique()); } else if (type == "delayedgroup") { diff --git a/src/tallies/filter_cell_instance.cpp b/src/tallies/filter_cell_instance.cpp new file mode 100644 index 000000000..a79e3017c --- /dev/null +++ b/src/tallies/filter_cell_instance.cpp @@ -0,0 +1,106 @@ +#include "openmc/tallies/filter_cell_instance.h" + +#include +#include + +#include "openmc/capi.h" +#include "openmc/cell.h" +#include "openmc/error.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +CellInstanceFilter::CellInstanceFilter(gsl::span instances) +{ + this->set_cell_instances(instances); +} + +void +CellInstanceFilter::from_xml(pugi::xml_node node) +{ + // Get cell IDs/instances + auto cells = get_node_array(node, "bins"); + Expects(cells.size() % 2 == 0); + + // Convert into vector of CellInstance + std::vector instances; + for (gsl::index i = 0; i < cells.size() / 2; ++i) { + int32_t cell_id = cells[2*i]; + gsl::index instance = cells[2*i + 1]; + auto search = model::cell_map.find(cell_id); + if (search == model::cell_map.end()) { + std::stringstream err_msg; + err_msg << "Could not find cell " << cell_id + << " specified on tally filter."; + throw std::runtime_error{err_msg.str()}; + } + gsl::index index = search->second; + instances.push_back({index, instance}); + } + + this->set_cell_instances(instances); +} + +void +CellInstanceFilter::set_cell_instances(gsl::span instances) +{ + // Clear existing cells + cell_instances_.clear(); + cell_instances_.reserve(instances.size()); + map_.clear(); + + // Update cells and mapping + for (auto& x : instances) { + Expects(x.index_cell >= 0); + Expects(x.index_cell < model::cells.size()); + const auto& c {model::cells[x.index_cell]}; + if (c->type_ != FILL_MATERIAL) { + throw std::invalid_argument{"Cell " + std::to_string(c->id_) + " is not " + "filled with a material. Only material cells can be used in a cell " + "instance filter."}; + } + cell_instances_.push_back(x); + map_[x] = cell_instances_.size() - 1; + } + + n_bins_ = cell_instances_.size(); +} + +void +CellInstanceFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + gsl::index index_cell = p->coord_[p->n_coord_ - 1].cell; + gsl::index instance = p->cell_instance_; + auto search = map_.find({index_cell, instance}); + if (search != map_.end()) { + int index_bin = search->second; + match.bins_.push_back(index_bin); + match.weights_.push_back(1.0); + } +} + +void +CellInstanceFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + size_t n = cell_instances_.size(); + xt::xtensor data({n, 2}); + for (gsl::index i = 0; i < n; ++i) { + const auto& x = cell_instances_[i]; + data(i, 0) = model::cells[x.index_cell]->id_; + data(i, 1) = x.instance; + } + write_dataset(filter_group, "bins", data); +} + +std::string +CellInstanceFilter::text_label(int bin) const +{ + const auto& x = cell_instances_[bin]; + auto cell_id = model::cells[x.index_cell]->id_; + return "Cell " + std::to_string(cell_id) + ", Instance " + + std::to_string(x.instance); +} + +} // namespace openmc diff --git a/tests/regression_tests/filter_cellinstance/__init__.py b/tests/regression_tests/filter_cellinstance/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/filter_cellinstance/inputs_true.dat b/tests/regression_tests/filter_cellinstance/inputs_true.dat new file mode 100644 index 000000000..5abf73433 --- /dev/null +++ b/tests/regression_tests/filter_cellinstance/inputs_true.dat @@ -0,0 +1,64 @@ + + + + + + + + + 2 2 + 4 4 + -4 -4 + +1 2 2 2 +2 1 2 2 +2 2 1 2 +2 2 2 1 + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 5 + 0 + + + 0.0 0.0 0.0 + + + + + + + 3 0 3 1 3 2 3 3 3 4 3 5 3 6 3 7 3 8 3 9 3 10 3 11 2 0 2 1 2 2 2 3 + + + 2 3 2 2 2 1 2 0 3 11 3 10 3 9 3 8 3 7 3 6 3 5 3 4 3 3 3 2 3 1 3 0 + + + 1 + total + + + 2 + total + + diff --git a/tests/regression_tests/filter_cellinstance/results_true.dat b/tests/regression_tests/filter_cellinstance/results_true.dat new file mode 100644 index 000000000..c1b48be8a --- /dev/null +++ b/tests/regression_tests/filter_cellinstance/results_true.dat @@ -0,0 +1,68 @@ +k-combined: +1.060380E+00 4.511966E-03 +tally 1: +8.636855E-02 +1.795640E-03 +1.478489E-01 +5.007770E-03 +1.687036E-01 +6.223488E-03 +1.035433E-01 +2.767195E-03 +2.789537E-01 +1.720816E-02 +1.420504E-01 +4.539772E-03 +1.325589E-01 +3.998010E-03 +2.934613E-01 +2.064779E-02 +1.231465E-01 +3.847169E-03 +1.282909E-01 +3.680695E-03 +1.545480E-01 +5.588411E-03 +5.110026E-02 +5.947113E-04 +1.119824E+01 +3.030467E+01 +2.847726E+01 +1.920547E+02 +2.819853E+01 +1.906807E+02 +8.816481E+00 +2.016785E+01 +tally 2: +8.816481E+00 +2.016785E+01 +2.819853E+01 +1.906807E+02 +2.847726E+01 +1.920547E+02 +1.119824E+01 +3.030467E+01 +5.110026E-02 +5.947113E-04 +1.545480E-01 +5.588411E-03 +1.282909E-01 +3.680695E-03 +1.231465E-01 +3.847169E-03 +2.934613E-01 +2.064779E-02 +1.325589E-01 +3.998010E-03 +1.420504E-01 +4.539772E-03 +2.789537E-01 +1.720816E-02 +1.035433E-01 +2.767195E-03 +1.687036E-01 +6.223488E-03 +1.478489E-01 +5.007770E-03 +8.636855E-02 +1.795640E-03 diff --git a/tests/regression_tests/filter_cellinstance/test.py b/tests/regression_tests/filter_cellinstance/test.py new file mode 100644 index 000000000..dc00394fc --- /dev/null +++ b/tests/regression_tests/filter_cellinstance/test.py @@ -0,0 +1,69 @@ +import openmc +import openmc.model +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + model = openmc.model.Model() + + # Materials + m1 = openmc.Material() + m1.set_density('g/cc', 4.5) + m1.add_nuclide('U235', 1.0) + m2 = openmc.Material() + m2.set_density('g/cc', 1.0) + m2.add_nuclide('H1', 1.0) + model.materials += [m1, m2] + + # Geometry + cyl1 = openmc.ZCylinder(r=0.7) + c1 = openmc.Cell(fill=m1, region=-cyl1) + c2 = openmc.Cell(fill=m2, region=+cyl1) + u1 = openmc.Universe(cells=[c1, c2]) + + cyl2 = openmc.ZCylinder(r=0.5) + c3 = openmc.Cell(fill=m1, region=-cyl2) + c4 = openmc.Cell(fill=m2, region=+cyl2) + u2 = openmc.Universe(cells=[c3, c4]) + + lat = openmc.RectLattice() + lat.lower_left = (-4, -4) + lat.pitch = (2, 2) + lat.universes = [ + [u1, u2, u2, u2], + [u2, u1, u2, u2], + [u2, u2, u1, u2], + [u2, u2, u2, u1] + ] + box = openmc.model.rectangular_prism(8.0, 8.0, boundary_type='reflective') + main_cell = openmc.Cell(fill=lat, region=box) + model.geometry.root_universe = openmc.Universe(cells=[main_cell]) + model.geometry.determine_paths() + + # Settings + model.settings.batches = 5 + model.settings.inactive = 0 + model.settings.particles = 1000 + model.settings.source = openmc.Source(space=openmc.stats.Point()) + + instances = ([(c3, i) for i in range(c3.num_instances)] + + [(c2, i) for i in range(c2.num_instances)]) + f1 = openmc.CellInstanceFilter(instances) + f2 = openmc.CellInstanceFilter(instances[::-1]) + t1 = openmc.Tally() + t1.filters = [f1] + t1.scores = ['total'] + t2 = openmc.Tally() + t2.filters = [f2] + t2.scores = ['total'] + model.tallies += [t1, t2] + + return model + + +def test_cell_instance(model): + harness = PyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py index 676ee452e..10587a812 100644 --- a/tests/unit_tests/test_filters.py +++ b/tests/unit_tests/test_filters.py @@ -22,6 +22,30 @@ def box_model(): return model +def test_cell_instance(): + c1 = openmc.Cell() + c2 = openmc.Cell() + f = openmc.CellInstanceFilter([(c1, 0), (c1, 1), (c1, 2), (c2, 0), (c2, 1)]) + + # Make sure __repr__ works + repr(f) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'cellinstance' + bins = [int(x) for x in elem.find('bins').text.split()] + assert all(x == c1.id for x in bins[:6:2]) + assert all(x == c2.id for x in bins[6::2]) + + # get_pandas_dataframe() + df = f.get_pandas_dataframe(f.num_bins, 1) + cells = df['cellinstance', 'cell'] + instances = df['cellinstance', 'instance'] + assert cells.apply(lambda x: x in (c1.id, c2.id)).all() + assert instances.apply(lambda x: x in (0, 1, 2)).all() + + def test_legendre(): n = 5 f = openmc.LegendreFilter(n) @@ -69,7 +93,7 @@ def test_spherical_harmonics(): f.cosine = 'particle' assert f.order == n assert f.bins[0] == 'Y0,0' - assert f.bins[-1] == 'Y{0},{0}'.format(n, n) + assert f.bins[-1] == 'Y{0},{0}'.format(n) assert len(f.bins) == (n + 1)**2 # Make sure __repr__ works