diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 9401156a64f..8edd99c0785 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -29,6 +29,9 @@ int openmc_cell_set_temperature( int32_t index, double T, const int32_t* instance, bool set_contained = false); int openmc_cell_set_translation(int32_t index, const double xyz[]); int openmc_cell_set_rotation(int32_t index, const double rot[], size_t rot_len); +int openmc_dagmc_universe_get_cell_ids( + int32_t univ_id, int32_t* ids, size_t* n); +int openmc_dagmc_universe_get_num_cells(int32_t univ_id, size_t* n); int openmc_energy_filter_get_bins( int32_t index, const double** energies, size_t* n); int openmc_energy_filter_set_bins( diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 70843140bad..032475ce982 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -320,7 +320,6 @@ class Cell { int32_t universe_; //!< Universe # this cell is in int32_t fill_; //!< Universe # filling this cell int32_t n_instances_ {0}; //!< Number of instances of this cell - GeometryType geom_type_; //!< Geometric representation type (CSG, DAGMC) //! \brief Index corresponding to this cell in distribcell arrays int distribcell_index_ {C_NONE}; @@ -350,6 +349,13 @@ class Cell { vector rotation_; vector offset_; //!< Distribcell offset table + + // Accessors + const GeometryType& geom_type() const { return geom_type_; } + GeometryType& geom_type() { return geom_type_; } + +private: + GeometryType geom_type_; //!< Geometric representation type (CSG, DAGMC) }; struct CellInstanceItem { diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 2facf4fc05e..3a636454fb5 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -133,6 +133,12 @@ class DAGUniverse : public Universe { void legacy_assign_material( std::string mat_string, std::unique_ptr& c) const; + //! Assign a material overriding normal assignement to a cell + //! \param[in] key The material key to override + //! \param[in] c The OpenMC cell to which the material is assigned + void override_assign_material(std::string key, moab::EntityHandle vol_handle, + std::unique_ptr& c) const; + //! Return the index into the model cells vector for a given DAGMC volume //! handle in the universe //! \param[in] vol MOAB handle to the DAGMC volume set @@ -184,6 +190,12 @@ class DAGUniverse : public Universe { //!< generate new material IDs for the universe bool has_graveyard_; //!< Indicates if the DAGMC geometry has a "graveyard" //!< volume + std::map> + instance_material_overrides; ///!< Map of material overrides + ///!< keys correspond to the material name + ///!< or id + ///!< values are a list of materials used + ///!< git for the override }; //============================================================================== diff --git a/include/openmc/string_utils.h b/include/openmc/string_utils.h index 2e8b0d14f39..e7e613534a4 100644 --- a/include/openmc/string_utils.h +++ b/include/openmc/string_utils.h @@ -19,6 +19,7 @@ void to_lower(std::string& str); int word_count(const std::string& str); vector split(const std::string& in); +vector split(const std::string& in, char delim); bool ends_with(const std::string& value, const std::string& ending); diff --git a/include/openmc/surface.h b/include/openmc/surface.h index af235301c14..498f71d4f9b 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -38,7 +38,6 @@ class Surface { int id_; //!< Unique ID std::string name_; //!< User-defined name unique_ptr bc_; //!< Boundary condition - GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC) bool surf_source_ {false}; //!< Activate source banking for the surface? explicit Surface(pugi::xml_node surf_node); @@ -91,6 +90,13 @@ class Surface { //! Get the BoundingBox for this surface. virtual BoundingBox bounding_box(bool /*pos_side*/) const { return {}; } + // Accessors + const GeometryType& geom_type() const { return geom_type_; } + GeometryType& geom_type() { return geom_type_; } + +private: + GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC) + protected: virtual void to_hdf5_inner(hid_t group_id) const = 0; }; diff --git a/openmc/__init__.py b/openmc/__init__.py index 566d287068f..bb972b4e6ad 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -15,6 +15,7 @@ from openmc.weight_windows import * from openmc.surface import * from openmc.universe import * +from openmc.dagmc import * from openmc.source import * from openmc.settings import * from openmc.lattice import * diff --git a/openmc/cell.py b/openmc/cell.py index fe3939bbe39..e5c836c46a7 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -780,3 +780,4 @@ def from_xml_element(cls, elem, surfaces, materials, get_universe): univ_id = int(get_text(elem, 'universe', 0)) get_universe(univ_id).add_cell(c) return c + diff --git a/openmc/dagmc.py b/openmc/dagmc.py new file mode 100644 index 00000000000..18d6500b5b5 --- /dev/null +++ b/openmc/dagmc.py @@ -0,0 +1,617 @@ +from collections.abc import Iterable +from numbers import Integral +from pathlib import Path + +import h5py +import lxml.etree as ET +import numpy as np +import warnings + +import openmc +import openmc.checkvalue as cv +from ._xml import get_text +from .checkvalue import check_type, check_value +from .surface import _BOUNDARY_TYPES +from .bounding_box import BoundingBox +from .utility_funcs import input_path + + +class DAGMCUniverse(openmc.UniverseBase): + """A reference to a DAGMC file to be used in the model. + + .. versionadded:: 0.13.0 + + .. versionadded:: 0.15.1-dev + Moved this classe from openmc.universe to openmc.dagmc + + Parameters + ---------- + filename : str + Path to the DAGMC file used to represent this universe. + universe_id : int, optional + Unique identifier of the universe. If not specified, an identifier will + automatically be assigned. + name : str, optional + Name of the universe. If not specified, the name is the empty string. + auto_geom_ids : bool + Set IDs automatically on initialization (True) or report overlaps in ID + space between CSG and DAGMC (False) + auto_mat_ids : bool + Set IDs automatically on initialization (True) or report overlaps in ID + space between OpenMC and UWUW materials (False) + material_overrides : dict + A dictionary of material overrides. The keys are material name + strings and the values are Iterables of openmc.Material objects. If a + material name is found in the DAGMC file, the material will be replaced + with the openmc.Material object in the value. + + Attributes + ---------- + id : int + Unique identifier of the universe + name : str + Name of the universe + filename : str + Path to the DAGMC file used to represent this universe. + auto_geom_ids : bool + Set IDs automatically on initialization (True) or report overlaps in ID + space between CSG and DAGMC (False) + auto_mat_ids : bool + Set IDs automatically on initialization (True) or report overlaps in ID + space between OpenMC and UWUW materials (False) + bounding_box : openmc.BoundingBox + Lower-left and upper-right coordinates of an axis-aligned bounding box + of the universe. + + .. versionadded:: 0.13.1 + material_names : list of str + Return a sorted list of materials names that are contained within the + DAGMC h5m file. This is useful when naming openmc.Material() objects + as each material name present in the DAGMC h5m file must have a + matching openmc.Material() with the same name. + + .. versionadded:: 0.13.2 + n_cells : int + The number of cells in the DAGMC model. This is the number of cells at + runtime and accounts for the implicit complement whether or not is it + present in the DAGMC file. + + .. versionadded:: 0.13.2 + n_surfaces : int + The number of surfaces in the model. + + .. versionadded:: 0.15 + material_overrides : dict + A dictionary of material overrides. The keys are material name + strings and the values are Iterables of openmc.Material objects. If a + material name is found in the DAGMC file, the material will be replaced + with the openmc.Material object in the value. + """ + + def __init__(self, + filename: cv.PathLike, + universe_id=None, + name='', + auto_geom_ids=False, + auto_mat_ids=False, + mat_overrides=None): + super().__init__(universe_id, name) + # Initialize class attributes + self.filename = filename + self.auto_geom_ids = auto_geom_ids + self.auto_mat_ids = auto_mat_ids + self.material_overrides = mat_overrides + self._dagmc_cells = [] + + def __repr__(self): + string = super().__repr__() + string += '{: <16}=\t{}\n'.format('\tGeom', 'DAGMC') + string += '{: <16}=\t{}\n'.format('\tFile', self.filename) + return string + + @property + def bounding_box(self): + with h5py.File(self.filename) as dagmc_file: + coords = dagmc_file['tstt']['nodes']['coordinates'][()] + lower_left_corner = coords.min(axis=0) + upper_right_corner = coords.max(axis=0) + return openmc.BoundingBox(lower_left_corner, upper_right_corner) + + @property + def filename(self): + return self._filename + + @filename.setter + def filename(self, val: cv.PathLike): + cv.check_type('DAGMC filename', val, cv.PathLike) + self._filename = input_path(val) + + @property + def material_overrides(self): + return self._material_overrides + + @material_overrides.setter + def material_overrides(self, val): + if val is None: + self._material_overrides = val + return + else: + cv.check_type('material overrides', val, dict) + for key, value in val.items(): + # ensuring key is a string and exists in the DAGMC file + cv.check_type('material name', key, str) + if key not in self.material_names: + raise ValueError( + f"Material name '{key}' not found in DAGMC file") + # ensuring overrides is an iterable of material name (strings) + cv.check_iterable_type('material objects', value, str) + + self._material_overrides = val + + def add_material_override(self, mat_name=None, cell_id=None, overrides=None): + """Add a material override to the universe. + + .. versionadded:: 0.15 + + Parameters + ---------- + key : str + Material name to override + value : Iterable of str + Material names to replace the key with + + """ + key = "" + if mat_name and cell_id: + raise ValueError("Only one of 'mat_name' or 'cell_id' can be set") + elif cell_id: + cv.check_type('cell id', cell_id, int) + if cell_id not in self.cells: + raise ValueError( + f"Cell ID '{cell_id}' not found in DAGMC universe") + else: + key = str(cell_id) + elif mat_name: + cv.check_type('material name', mat_name, str) + if mat_name not in self.material_names: + raise ValueError( + f"Material name '{mat_name}' not found in DAGMC file") + else: + key = mat_name + else: + raise ValueError("Either 'mat_name' or 'cell_id' must be set") + + cv.check_iterable_type('material objects', overrides, str) + self.material_overrides[mat_name] = overrides + + @property + def auto_geom_ids(self): + return self._auto_geom_ids + + @auto_geom_ids.setter + def auto_geom_ids(self, val): + cv.check_type('DAGMC automatic geometry ids', val, bool) + self._auto_geom_ids = val + + @property + def auto_mat_ids(self): + return self._auto_mat_ids + + @auto_mat_ids.setter + def auto_mat_ids(self, val): + cv.check_type('DAGMC automatic material ids', val, bool) + self._auto_mat_ids = val + + @property + def material_names(self): + dagmc_file_contents = h5py.File(self.filename) + material_tags_hex = dagmc_file_contents['/tstt/tags/NAME'].get( + 'values') + material_tags_ascii = [] + for tag in material_tags_hex: + candidate_tag = tag.tobytes().decode().replace('\x00', '') + # tags might be for temperature or reflective surfaces + if candidate_tag.startswith('mat:'): + # if name ends with _comp remove it, it is not parsed + if candidate_tag.endswith('_comp'): + candidate_tag = candidate_tag[:-5] + # removes first 4 characters as openmc.Material name should be + # set without the 'mat:' part of the tag + material_tags_ascii.append(candidate_tag[4:]) + + return sorted(set(material_tags_ascii)) + + def _n_geom_elements(self, geom_type): + """ + Helper function for retrieving the number geometric entities in a DAGMC + file + + Parameters + ---------- + geom_type : str + The type of geometric entity to count. One of {'Volume', 'Surface'}. Returns + the runtime number of voumes in the DAGMC model (includes implicit complement). + + Returns + ------- + int + Number of geometry elements of the specified type + """ + cv.check_value('geometry type', geom_type, ('volume', 'surface')) + + def decode_str_tag(tag_val): + return tag_val.tobytes().decode().replace('\x00', '') + + with h5py.File(self.filename) as dagmc_file: + category_data = dagmc_file['tstt/tags/CATEGORY/values'] + category_strs = map(decode_str_tag, category_data) + n = sum([v == geom_type.capitalize() for v in category_strs]) + + # check for presence of an implicit complement in the file and + # increment the number of cells if it doesn't exist + if geom_type == 'volume': + name_data = dagmc_file['tstt/tags/NAME/values'] + name_strs = map(decode_str_tag, name_data) + if not sum(['impl_complement' in n for n in name_strs]): + n += 1 + return n + + @property + def n_cells(self): + return self._n_geom_elements('volume') + + @property + def n_surfaces(self): + return self._n_geom_elements('surface') + + def create_xml_subelement(self, xml_element, memo=None): + if memo is None: + memo = set() + + if self in memo: + return + + memo.add(self) + + # Ensure that the material overrides are up-t-date + self.build_overide_mat_from_cells() + + # Set xml element values + dagmc_element = ET.Element('dagmc_universe') + dagmc_element.set('id', str(self.id)) + + if self.auto_geom_ids: + dagmc_element.set('auto_geom_ids', 'true') + if self.auto_mat_ids: + dagmc_element.set('auto_mat_ids', 'true') + dagmc_element.set('filename', str(self.filename)) + if self.material_overrides: + mat_element = ET.Element('material_overrides') + for key in self.material_overrides: + print(key, self.material_overrides[key]) + mat_element.set('id_{}'.format(key), ';'.join( + t for t in self.material_overrides[key])) + dagmc_element.append(mat_element) + xml_element.append(dagmc_element) + + def build_overide_mat_from_cells(self): + """ + Builds the material override dictionary for cells with multiple instances. + + Returns: + None + """ + self.material_overrides = {} + for cell in self.cells.values(): + if isinstance(cell.fill, Iterable): + self.material_overrides[str(cell.id)] = [mat.name for mat in cell.fill if mat] + + def bounding_region( + self, + bounded_type: str = 'box', + boundary_type: str = 'vacuum', + starting_id: int = 10000, + padding_distance: float = 0. + ): + """Creates a either a spherical or box shaped bounding region around + the DAGMC geometry. + + .. versionadded:: 0.13.1 + + Parameters + ---------- + bounded_type : str + The type of bounding surface(s) to use when constructing the region. + Options include a single spherical surface (sphere) or a rectangle + made from six planes (box). + boundary_type : str + Boundary condition that defines the behavior for particles hitting + the surface. Defaults to vacuum boundary condition. Passed into the + surface construction. + starting_id : int + Starting ID of the surface(s) used in the region. For bounded_type + 'box', the next 5 IDs will also be used. Defaults to 10000 to reduce + the chance of an overlap of surface IDs with the DAGMC geometry. + padding_distance : float + Distance between the bounding region surfaces and the minimal + bounding box. Allows for the region to be larger than the DAGMC + geometry. + + Returns + ------- + openmc.Region + Region instance + """ + + check_type('boundary type', boundary_type, str) + check_value('boundary type', boundary_type, _BOUNDARY_TYPES) + check_type('starting surface id', starting_id, Integral) + check_type('bounded type', bounded_type, str) + check_value('bounded type', bounded_type, ('box', 'sphere')) + + bbox = self.bounding_box.expand(padding_distance, True) + + if bounded_type == 'sphere': + radius = np.linalg.norm(bbox.upper_right - bbox.center) + bounding_surface = openmc.Sphere( + surface_id=starting_id, + x0=bbox.center[0], + y0=bbox.center[1], + z0=bbox.center[2], + boundary_type=boundary_type, + r=radius, + ) + + return -bounding_surface + + if bounded_type == 'box': + # defines plane surfaces for all six faces of the bounding box + lower_x = openmc.XPlane(bbox[0][0], surface_id=starting_id) + upper_x = openmc.XPlane(bbox[1][0], surface_id=starting_id+1) + lower_y = openmc.YPlane(bbox[0][1], surface_id=starting_id+2) + upper_y = openmc.YPlane(bbox[1][1], surface_id=starting_id+3) + lower_z = openmc.ZPlane(bbox[0][2], surface_id=starting_id+4) + upper_z = openmc.ZPlane(bbox[1][2], surface_id=starting_id+5) + + region = +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z + + for surface in region.get_surfaces().values(): + surface.boundary_type = boundary_type + + return region + + def bounded_universe(self, bounding_cell_id=10000, **kwargs): + """Returns an openmc.Universe filled with this DAGMCUniverse and bounded + with a cell. Defaults to a box cell with a vacuum surface however this + can be changed using the kwargs which are passed directly to + DAGMCUniverse.bounding_region(). + + Parameters + ---------- + bounding_cell_id : int + The cell ID number to use for the bounding cell, defaults to 10000 to reduce + the chance of overlapping ID numbers with the DAGMC geometry. + + Returns + ------- + openmc.Universe + Universe instance + """ + bounding_cell = openmc.Cell( + fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs)) + return openmc.Universe(cells=[bounding_cell]) + + @classmethod + def from_hdf5(cls, group): + """Create DAGMC universe from HDF5 group + + Parameters + ---------- + group : h5py.Group + Group in HDF5 file + + Returns + ------- + openmc.DAGMCUniverse + DAGMCUniverse instance + + """ + id = int(group.name.split('/')[-1].lstrip('universe ')) + fname = group['filename'][()].decode() + name = group['name'][()].decode() if 'name' in group else None + + out = cls(fname, universe_id=id, name=name) + + out.auto_geom_ids = bool(group.attrs['auto_geom_ids']) + out.auto_mat_ids = bool(group.attrs['auto_mat_ids']) + + return out + + @classmethod + def from_xml_element(cls, elem): + """Generate DAGMC universe from XML element + + Parameters + ---------- + elem : lxml.etree._Element + `` element + + Returns + ------- + openmc.DAGMCUniverse + DAGMCUniverse instance + + """ + id = int(get_text(elem, 'id')) + fname = get_text(elem, 'filename') + + out = cls(fname, universe_id=id) + + name = get_text(elem, 'name') + if name is not None: + out.name = name + + out.auto_geom_ids = bool(elem.get('auto_geom_ids')) + out.auto_mat_ids = bool(elem.get('auto_mat_ids')) + + for item in elem.find('material_overrides').attrib: + origin_mat, overwrite = item + for mat_name in overwrite.split(): + out.material_overrides.setdefault( + origin_mat.lower(), []).append(mat_name) + + return out + + def _partial_deepcopy(self): + """Clone all of the openmc.DAGMCUniverse object's attributes except for + its cells, as they are copied within the clone function. This should + only to be used within the openmc.UniverseBase.clone() context. + """ + clone = openmc.DAGMCUniverse(name=self.name, filename=self.filename) + clone.volume = self.volume + clone.auto_geom_ids = self.auto_geom_ids + clone.auto_mat_ids = self.auto_mat_ids + return clone + + def add_cell(self, cell): + """Add a cell to the universe. + + Parameters + ---------- + cell : openmc.DAGMCCell + Cell to add + + """ + if not isinstance(cell, openmc.DAGMCCell): + msg = f'Unable to add a DAGMCCell to DAGMCUniverse ' \ + f'ID="{self._id}" since "{cell}" is not a DAGMCCell' + raise TypeError(msg) + + cell_id = cell.id + + if cell_id not in self._cells: + self._cells[cell_id] = cell + + def remove_cell(self, cell): + """Remove a cell from the universe. + + Parameters + ---------- + cell : openmc.Cell + Cell to remove + + """ + + if not isinstance(cell, openmc.DAGMCCell): + msg = f'Unable to remove a Cell from Universe ID="{self._id}" ' \ + f'since "{cell}" is not a Cell' + raise TypeError(msg) + + # If the Cell is in the Universe's list of Cells, delete it + self._cells.pop(cell.id, None) + + def sync_dagmc_cells(self, mats): + """Synchronize DAGMC cell information between Python and C API + + .. versionadded:: 0.15.1-dev + + """ + import openmc.lib + if not openmc.lib.is_initialized: + raise RuntimeError("This universe must be part of an openmc.Model " + "initialized via Model.init_lib before calling " + "this method.") + + dagmc_cell_ids = openmc.lib.dagmc.get_dagmc_cell_ids(self.id) + if len(dagmc_cell_ids) != self.n_cells: + raise ValueError( + f"Number of cells in DAGMC universe {self.id} does not match " + f"the number of cells in the Python universe." + ) + + mats_per_id = {mat.id: mat for mat in mats} + for dag_cell_id in dagmc_cell_ids: + dag_cell = openmc.lib.cells[dag_cell_id] + if isinstance(dag_cell.fill, Iterable): + fill = [mats_per_id[mat.id] + for mat in dag_cell.fill if mat] + else: + if dag_cell.fill: + fill = mats_per_id[dag_cell.fill.id] + else: + fill = None + dag_pseudo_cell = openmc.DAGMCCell( + cell_id=dag_cell_id, fill=fill + ) + self.add_cell(dag_pseudo_cell) + + +class DAGMCCell(openmc.Cell): + """ + .. versionadded:: 0.15.1-dev + A cell class for DAGMC-based geometries. + + Parameters + ---------- + cell_id : int or None, optional + Unique identifier for the cell. If None, an identifier will be + automatically assigned. + name : str, optional + Name of the cell. + fill : openmc.Material or None, optional + Material filling the cell. If None, the cell is filled with vacuum. + + Attributes + ---------- + DAG_parent_universe : int + The parent universe of the cell. + + """ + def __init__(self, cell_id=None, name='', fill=None): + super().__init__(cell_id, name, fill, None) + + @property + def DAG_parent_universe(self): + """Get the parent universe of the cell.""" + return self._parent_universe + + @DAG_parent_universe.setter + def DAG_parent_universe(self, universe): + """Set the parent universe of the cell.""" + self._parent_universe = universe.id + + def boundingbox(self): + warnings.warn("Bounding box is not available for cells in a DAGMC " + "universe", Warning) + return BoundingBox.infinite() + + def get_all_cells(self, memo=None): + warnings.warn("get_all_cells is not available for cells in a DAGMC " + "universe", Warning) + return {} + + def get_all_universes(self, memo=None): + warnings.warn("get_all_universes is not available for cells in a " + "DAGMC universe", Warning) + return {} + + def clone(self, clone_materials=True, clone_regions=True, memo=None): + warnings.warn("clone is not available for cells in a DAGMC universe", + Warning) + return None + + def plot(self, *args, **kwargs): + warnings.warn("plot is not available for cells in a DAGMC universe", + Warning) + return None + + def create_xml_subelement(self, xml_element, memo=None): + warnings.warn("create_xml_subelement is not available for cells in a " + "DAGMC universe", Warning) + return None + + @classmethod + def from_xml_element(cls, elem, surfaces, materials, get_universe): + warnings.warn("from_xml_element is not available for cells in a DAGMC " + "universe", Warning) + return None diff --git a/openmc/geometry.py b/openmc/geometry.py index a52dc4418d2..ef104b13123 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -402,6 +402,23 @@ def get_all_nuclides(self) -> list[str]: for material in self.get_all_materials().values(): all_nuclides |= set(material.get_nuclides()) return sorted(all_nuclides) + + def get_all_dagmc_universes(self) -> typing.Dict[int, openmc.DAGMCUniverse]: + """Return all universes in the geometry. + + Returns + ------- + dict + Dictionary mapping universe IDs to :class:`openmc.DAGMCUniverse` + instances + + """ + universes = self.get_all_universes() + dag_universes = {} + for id, uni in universes.items(): + if isinstance(uni, openmc.DAGMCUniverse): + dag_universes[id] = uni + return dag_universes def get_all_materials(self) -> dict[int, openmc.Material]: """Return all materials within the geometry. diff --git a/openmc/lib/__init__.py b/openmc/lib/__init__.py index 9bb2efb38af..4cdd53b6a70 100644 --- a/openmc/lib/__init__.py +++ b/openmc/lib/__init__.py @@ -68,6 +68,7 @@ def _uwuw_enabled(): from .math import * from .plot import * from .weight_windows import * +from .dagmc import get_dagmc_cell_ids, get_dagmc_universe_num_cells # Flag to denote whether or not openmc.lib.init has been called # TODO: Establish and use a flag in the C++ code to represent the status of the diff --git a/openmc/lib/dagmc.py b/openmc/lib/dagmc.py new file mode 100644 index 00000000000..d42b5ea857f --- /dev/null +++ b/openmc/lib/dagmc.py @@ -0,0 +1,63 @@ +import sys +from ctypes import c_int, c_int32, POINTER, c_size_t + +import numpy as np + +from . import _dll +from .error import _error_handler + + +# DAGMC functions +_dll.openmc_dagmc_universe_get_cell_ids.argtypes = [c_int32, POINTER(c_int32), POINTER(c_size_t)] +_dll.openmc_dagmc_universe_get_cell_ids.restype = c_int +_dll.openmc_dagmc_universe_get_cell_ids.errcheck = _error_handler +_dll.openmc_dagmc_universe_get_num_cells.argtypes = [c_int32, POINTER(c_size_t)] +_dll.openmc_dagmc_universe_get_num_cells.restype = c_int +_dll.openmc_dagmc_universe_get_num_cells.errcheck = _error_handler + + +def get_dagmc_cell_ids(dagmc_id): + """Get the DAGMC cell IDs for a volume. + + Parameters + ---------- + dagmc_id : int + ID of the DAGMC Universe to get cell IDs from. + n_cells : int + Number of cells in the DAGMC Universe. + + Returns + ------- + numpy.ndarray + DAGMC cell IDs for the volume. + + """ + n = c_size_t() + _dll.openmc_dagmc_universe_get_num_cells(dagmc_id, n) + cell_ids = np.empty(n.value, dtype=np.int32) + + _dll.openmc_dagmc_universe_get_cell_ids( + dagmc_id, + cell_ids.ctypes.data_as(POINTER(c_int32)), + n + ) + return cell_ids + + +def get_dagmc_universe_num_cells(dagmc_id): + """Get the number of cells in a DAGMC universe. + + Parameters + ---------- + dagmc_id : int + ID of the DAGMC Universe to get the number of cell from. + + Returns + ------- + int + Number of cells in the DAGMC Universe. + + """ + n = c_size_t() + _dll.openmc_dagmc_universe_get_num_cells(dagmc_id, n) + return n.value diff --git a/openmc/model/model.py b/openmc/model/model.py index 0ea6db8db88..188be487598 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -17,6 +17,7 @@ from openmc.executor import _process_CLI_arguments from openmc.checkvalue import check_type, check_value, PathLike from openmc.exceptions import InvalidIDError +import openmc.lib from openmc.utility_funcs import change_directory @@ -324,6 +325,26 @@ def init_lib(self, threads=None, geometry_debug=False, restart_file=None, # communicator openmc.lib.init(args=args, intracomm=intracomm, output=output) + def sync_dagmc_universe(self): + """ + Synchronize all DAGMC universes in the current geometry. + This method iterates over all DAGMC universes in the geometry and + synchronizes their cells with the current material assignments. + + .. versionadded:: 0.15.1-dev + + """ + if self.is_initialized: + if self.materials: + materials = self.materials + else: + materials = self.geometry.get_all_materials().values() + for dagmc_universe in self.geometry.get_all_dagmc_universes().values(): + dagmc_universe.sync_dagmc_cells(materials) + else: + raise ValueError("The model must be initialized before calling " + "this method") + def finalize_lib(self): """Finalize simulation and free memory allocated for the C API @@ -1154,32 +1175,63 @@ def update_material_volumes(self, names_or_ids, volume): self._change_py_lib_attribs(names_or_ids, volume, 'material', 'volume') - def differentiate_depletable_mats(self, diff_volume_method: str): + def differentiate_depletable_mats(self, diff_volume_method : str = None): """Assign distribmats for each depletable material .. versionadded:: 0.14.0 + .. version added:: 0.15.1-dev + diff_volume_method default is None, do not apply volume to the new + materials. Is now a convenience method for + differentiate_mats(diff_volume_method, depletable_only=True) + Parameters ---------- diff_volume_method : str Specifies how the volumes of the new materials should be found. - Default is to 'divide equally' which divides the original material - volume equally between the new materials, 'match cell' sets the - volume of the material to volume of the cell they fill. + Default is to 'None', do not apply volume to the new materials, + 'divide equally' which divides the original material + volume equally between the new materials, + 'match cell' sets the volume of the material to volume of the cell + they fill. """ + self.differentiate_mats(diff_volume_method, depletable_only=True) + + def differentiate_mats(self, diff_volume_method: str = None, depletable_only: bool = True): + """Assign distribmats for each material + + .. versionadded:: 0.15.1-dev + + Parameters + ---------- + diff_volume_method : str + Specifies how the volumes of the new materials should be found. + Default is to 'None', do not apply volume to the new materials, + 'divide equally' which divides the original material + volume equally between the new materials, + 'match cell' sets the volume of the material to volume of the cell + they fill. + depletable_only : bool + Default is True, only depletable materials will be differentiated all materials will be + differentiated otherwise. + """ + check_value('volume differentiation method', diff_volume_method, ["divide equally", "match cell", None]) + # Count the number of instances for each cell and material self.geometry.determine_paths(instances_only=True) - # Extract all depletable materials which have multiple instances + # Extract all or depletable_only materials which have multiple instance distribmats = set( [mat for mat in self.materials - if mat.depletable and mat.num_instances > 1]) + if (mat.depletable or not depletable_only) and mat.num_instances > 1]) - if diff_volume_method == 'divide equally': + if diff_volume_method == "divide equally": for mat in distribmats: if mat.volume is None: - raise RuntimeError("Volume not specified for depletable " - f"material with ID={mat.id}.") + raise RuntimeError( + "Volume not specified for depletable " + f"material with ID={mat.id}." + ) mat.volume /= mat.num_instances if distribmats: @@ -1187,11 +1239,11 @@ def differentiate_depletable_mats(self, diff_volume_method: str): for cell in self.geometry.get_all_material_cells().values(): if cell.fill in distribmats: mat = cell.fill - if diff_volume_method == 'divide equally': + if diff_volume_method != 'match cell': cell.fill = [mat.clone() for _ in range(cell.num_instances)] elif diff_volume_method == 'match cell': - for _ in range(cell.num_instances): - cell.fill = mat.clone() + cell.fill = mat.clone() + for i in range(cell.num_instances): if not cell.volume: raise ValueError( f"Volume of cell ID={cell.id} not specified. " @@ -1199,6 +1251,9 @@ def differentiate_depletable_mats(self, diff_volume_method: str): "diff_volume_method='match cell'." ) cell.fill.volume = cell.volume + if isinstance(cell, openmc.DAGMCCell): + for i in range(cell.num_instances): + cell.fill[i].name = f"{cell.fill[i].name}_{cell.id}_{i}" if self.materials is not None: self.materials = openmc.Materials( diff --git a/openmc/universe.py b/openmc/universe.py index 648b773df6f..5409222229d 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -55,6 +55,10 @@ def __repr__(self): def name(self): return self._name + @property + def cells(self): + return self._cells + @name.setter def name(self, name): if name is not None: @@ -135,6 +139,130 @@ def create_xml_subelement(self, xml_element, memo=None): """ + def _determine_paths(self, path='', instances_only=False): + """Count the number of instances for each cell in the universe, and + record the count in the :attr:`Cell.num_instances` properties.""" + + univ_path = path + f'u{self.id}' + + for cell in self.cells.values(): + cell_path = f'{univ_path}->c{cell.id}' + fill = cell._fill + fill_type = cell.fill_type + + # If universe-filled, recursively count cells in filling universe + if fill_type == 'universe': + fill._determine_paths(cell_path + '->', instances_only) + # If lattice-filled, recursively call for all universes in lattice + elif fill_type == 'lattice': + latt = fill + + # Count instances in each universe in the lattice + for index in latt._natural_indices: + latt_path = '{}->l{}({})->'.format( + cell_path, latt.id, ",".join(str(x) for x in index)) + univ = latt.get_universe(index) + univ._determine_paths(latt_path, instances_only) + + else: + if fill_type == 'material': + mat = fill + elif fill_type == 'distribmat': + mat = fill[cell._num_instances] + else: + mat = None + + if mat is not None: + mat._num_instances += 1 + if not instances_only: + mat._paths.append(f'{cell_path}->m{mat.id}') + + # Append current path + cell._num_instances += 1 + if not instances_only: + cell._paths.append(cell_path) + + def add_cells(self, cells): + """Add multiple cells to the universe. + + Parameters + ---------- + cells : Iterable of openmc.Cell + Cells to add + + """ + + if not isinstance(cells, Iterable): + msg = f'Unable to add Cells to Universe ID="{self._id}" since ' \ + f'"{cells}" is not iterable' + raise TypeError(msg) + + for cell in cells: + self.add_cell(cell) + + @abstractmethod + def add_cell(self, cell): + pass + + @abstractmethod + def remove_cell(self, cell): + pass + + def clear_cells(self): + """Remove all cells from the universe.""" + + self._cells.clear() + + def get_all_cells(self, memo=None): + """Return all cells that are contained within the universe + + Returns + ------- + cells : dict + Dictionary whose keys are cell IDs and values are :class:`Cell` + instances + + """ + + if memo is None: + memo = set() + elif self in memo: + return {} + memo.add(self) + + # Add this Universe's cells to the dictionary + cells = {} + cells.update(self._cells) + + # Append all Cells in each Cell in the Universe to the dictionary + for cell in self._cells.values(): + cells.update(cell.get_all_cells(memo)) + + return cells + + def get_all_materials(self, memo=None): + """Return all materials that are contained within the universe + + Returns + ------- + materials : dict + Dictionary whose keys are material IDs and values are + :class:`Material` instances + + """ + + if memo is None: + memo = set() + + materials = {} + + # Append all Cells in each Cell in the Universe to the dictionary + cells = self.get_all_cells(memo) + for cell in cells.values(): + materials.update(cell.get_all_materials(memo)) + + return materials + @abstractmethod def _partial_deepcopy(self): """Deepcopy all parameters of an openmc.UniverseBase object except its cells. @@ -182,93 +310,6 @@ def clone(self, clone_materials=True, clone_regions=True, memo=None): return memo[self] - -class Universe(UniverseBase): - """A collection of cells that can be repeated. - - Parameters - ---------- - universe_id : int, optional - Unique identifier of the universe. If not specified, an identifier will - automatically be assigned - name : str, optional - Name of the universe. If not specified, the name is the empty string. - cells : Iterable of openmc.Cell, optional - Cells to add to the universe. By default no cells are added. - - Attributes - ---------- - id : int - Unique identifier of the universe - name : str - Name of the universe - cells : dict - Dictionary whose keys are cell IDs and values are :class:`Cell` - instances - volume : float - Volume of the universe in cm^3. This can either be set manually or - calculated in a stochastic volume calculation and added via the - :meth:`Universe.add_volume_information` method. - bounding_box : openmc.BoundingBox - Lower-left and upper-right coordinates of an axis-aligned bounding box - of the universe. - - """ - - def __init__(self, universe_id=None, name='', cells=None): - super().__init__(universe_id, name) - - if cells is not None: - self.add_cells(cells) - - def __repr__(self): - string = super().__repr__() - string += '{: <16}=\t{}\n'.format('\tGeom', 'CSG') - string += '{: <16}=\t{}\n'.format('\tCells', list(self._cells.keys())) - return string - - @property - def cells(self): - return self._cells - - @property - def bounding_box(self) -> openmc.BoundingBox: - regions = [c.region for c in self.cells.values() - if c.region is not None] - if regions: - return openmc.Union(regions).bounding_box - else: - return openmc.BoundingBox.infinite() - - @classmethod - def from_hdf5(cls, group, cells): - """Create universe from HDF5 group - - Parameters - ---------- - group : h5py.Group - Group in HDF5 file - cells : dict - Dictionary mapping cell IDs to instances of :class:`openmc.Cell`. - - Returns - ------- - openmc.Universe - Universe instance - - """ - universe_id = int(group.name.split('/')[-1].lstrip('universe ')) - cell_ids = group['cells'][()] - - # Create this Universe - universe = cls(universe_id) - - # Add each Cell to the Universe - for cell_id in cell_ids: - universe.add_cell(cells[cell_id]) - - return universe - def find(self, point): """Find cells/universes/lattices which contain a given point @@ -528,98 +569,37 @@ def plot(self, origin=None, width=None, pixels=40000, axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs) return axes - def add_cell(self, cell): - """Add a cell to the universe. + def get_nuclides(self): + """Returns all nuclides in the universe - Parameters - ---------- - cell : openmc.Cell - Cell to add + Returns + ------- + nuclides : list of str + List of nuclide names """ - if not isinstance(cell, openmc.Cell): - msg = f'Unable to add a Cell to Universe ID="{self._id}" since ' \ - f'"{cell}" is not a Cell' - raise TypeError(msg) + nuclides = [] - cell_id = cell.id + # Append all Nuclides in each Cell in the Universe to the dictionary + for cell in self.cells.values(): + for nuclide in cell.get_nuclides(): + if nuclide not in nuclides: + nuclides.append(nuclide) - if cell_id not in self._cells: - self._cells[cell_id] = cell + return nuclides - def add_cells(self, cells): - """Add multiple cells to the universe. + def get_nuclide_densities(self): + """Return all nuclides contained in the universe - Parameters - ---------- - cells : Iterable of openmc.Cell - Cells to add + Returns + ------- + nuclides : dict + Dictionary whose keys are nuclide names and values are 2-tuples of + (nuclide, density) """ - - if not isinstance(cells, Iterable): - msg = f'Unable to add Cells to Universe ID="{self._id}" since ' \ - f'"{cells}" is not iterable' - raise TypeError(msg) - - for cell in cells: - self.add_cell(cell) - - def remove_cell(self, cell): - """Remove a cell from the universe. - - Parameters - ---------- - cell : openmc.Cell - Cell to remove - - """ - - if not isinstance(cell, openmc.Cell): - msg = f'Unable to remove a Cell from Universe ID="{self._id}" ' \ - f'since "{cell}" is not a Cell' - raise TypeError(msg) - - # If the Cell is in the Universe's list of Cells, delete it - self._cells.pop(cell.id, None) - - def clear_cells(self): - """Remove all cells from the universe.""" - - self._cells.clear() - - def get_nuclides(self): - """Returns all nuclides in the universe - - Returns - ------- - nuclides : list of str - List of nuclide names - - """ - - nuclides = [] - - # Append all Nuclides in each Cell in the Universe to the dictionary - for cell in self.cells.values(): - for nuclide in cell.get_nuclides(): - if nuclide not in nuclides: - nuclides.append(nuclide) - - return nuclides - - def get_nuclide_densities(self): - """Return all nuclides contained in the universe - - Returns - ------- - nuclides : dict - Dictionary whose keys are nuclide names and values are 2-tuples of - (nuclide, density) - - """ - nuclides = {} + nuclides = {} if self._atoms: volume = self.volume @@ -636,150 +616,20 @@ def get_nuclide_densities(self): return nuclides - def get_all_cells(self, memo=None): - """Return all cells that are contained within the universe - - Returns - ------- - cells : dict - Dictionary whose keys are cell IDs and values are :class:`Cell` - instances - - """ - - if memo is None: - memo = set() - elif self in memo: - return {} - memo.add(self) - - # Add this Universe's cells to the dictionary - cells = {} - cells.update(self._cells) - - # Append all Cells in each Cell in the Universe to the dictionary - for cell in self._cells.values(): - cells.update(cell.get_all_cells(memo)) - - return cells - - def get_all_materials(self, memo=None): - """Return all materials that are contained within the universe - - Returns - ------- - materials : dict - Dictionary whose keys are material IDs and values are - :class:`Material` instances - - """ - - if memo is None: - memo = set() - - materials = {} - - # Append all Cells in each Cell in the Universe to the dictionary - cells = self.get_all_cells(memo) - for cell in cells.values(): - materials.update(cell.get_all_materials(memo)) - - return materials - - def create_xml_subelement(self, xml_element, memo=None): - if memo is None: - memo = set() - - # Iterate over all Cells - for cell in self._cells.values(): - - # If the cell was already written, move on - if cell in memo: - continue - - memo.add(cell) - - # Create XML subelement for this Cell - cell_element = cell.create_xml_subelement(xml_element, memo) - - # Append the Universe ID to the subelement and add to Element - cell_element.set("universe", str(self._id)) - xml_element.append(cell_element) - - def _determine_paths(self, path='', instances_only=False): - """Count the number of instances for each cell in the universe, and - record the count in the :attr:`Cell.num_instances` properties.""" - - univ_path = path + f'u{self.id}' - - for cell in self.cells.values(): - cell_path = f'{univ_path}->c{cell.id}' - fill = cell._fill - fill_type = cell.fill_type - - # If universe-filled, recursively count cells in filling universe - if fill_type == 'universe': - fill._determine_paths(cell_path + '->', instances_only) - - # If lattice-filled, recursively call for all universes in lattice - elif fill_type == 'lattice': - latt = fill - # Count instances in each universe in the lattice - for index in latt._natural_indices: - latt_path = '{}->l{}({})->'.format( - cell_path, latt.id, ",".join(str(x) for x in index)) - univ = latt.get_universe(index) - univ._determine_paths(latt_path, instances_only) - else: - if fill_type == 'material': - mat = fill - elif fill_type == 'distribmat': - mat = fill[cell._num_instances] - else: - mat = None - - if mat is not None: - mat._num_instances += 1 - if not instances_only: - mat._paths.append(f'{cell_path}->m{mat.id}') - - # Append current path - cell._num_instances += 1 - if not instances_only: - cell._paths.append(cell_path) - - def _partial_deepcopy(self): - """Clone all of the openmc.Universe object's attributes except for its cells, - as they are copied within the clone function. This should only to be - used within the openmc.UniverseBase.clone() context. - """ - clone = openmc.Universe(name=self.name) - clone.volume = self.volume - return clone - - -class DAGMCUniverse(UniverseBase): - """A reference to a DAGMC file to be used in the model. - - .. versionadded:: 0.13.0 +class Universe(UniverseBase): + """A collection of cells that can be repeated. Parameters ---------- - filename : path-like - Path to the DAGMC file used to represent this universe. universe_id : int, optional Unique identifier of the universe. If not specified, an identifier will - automatically be assigned. + automatically be assigned name : str, optional Name of the universe. If not specified, the name is the empty string. - auto_geom_ids : bool - Set IDs automatically on initialization (True) or report overlaps in ID - space between CSG and DAGMC (False) - auto_mat_ids : bool - Set IDs automatically on initialization (True) or report overlaps in ID - space between OpenMC and UWUW materials (False) + cells : Iterable of openmc.Cell, optional + Cells to add to the universe. By default no cells are added. Attributes ---------- @@ -787,334 +637,133 @@ class DAGMCUniverse(UniverseBase): Unique identifier of the universe name : str Name of the universe - filename : str - Path to the DAGMC file used to represent this universe. - auto_geom_ids : bool - Set IDs automatically on initialization (True) or report overlaps in ID - space between CSG and DAGMC (False) - auto_mat_ids : bool - Set IDs automatically on initialization (True) or report overlaps in ID - space between OpenMC and UWUW materials (False) + cells : dict + Dictionary whose keys are cell IDs and values are :class:`Cell` + instances + volume : float + Volume of the universe in cm^3. This can either be set manually or + calculated in a stochastic volume calculation and added via the + :meth:`Universe.add_volume_information` method. bounding_box : openmc.BoundingBox Lower-left and upper-right coordinates of an axis-aligned bounding box of the universe. - .. versionadded:: 0.13.1 - material_names : list of str - Return a sorted list of materials names that are contained within the - DAGMC h5m file. This is useful when naming openmc.Material() objects - as each material name present in the DAGMC h5m file must have a - matching openmc.Material() with the same name. - - .. versionadded:: 0.13.2 - n_cells : int - The number of cells in the DAGMC model. This is the number of cells at - runtime and accounts for the implicit complement whether or not is it - present in the DAGMC file. - - .. versionadded:: 0.13.2 - n_surfaces : int - The number of surfaces in the model. - - .. versionadded:: 0.13.2 - """ - def __init__(self, - filename: cv.PathLike, - universe_id=None, - name='', - auto_geom_ids=False, - auto_mat_ids=False): + def __init__(self, universe_id=None, name='', cells=None): super().__init__(universe_id, name) - # Initialize class attributes - self.filename = filename - self.auto_geom_ids = auto_geom_ids - self.auto_mat_ids = auto_mat_ids + + if cells is not None: + self.add_cells(cells) def __repr__(self): string = super().__repr__() - string += '{: <16}=\t{}\n'.format('\tGeom', 'DAGMC') - string += '{: <16}=\t{}\n'.format('\tFile', self.filename) + string += '{: <16}=\t{}\n'.format('\tGeom', 'CSG') + string += '{: <16}=\t{}\n'.format('\tCells', list(self._cells.keys())) return string @property - def bounding_box(self): - with h5py.File(self.filename) as dagmc_file: - coords = dagmc_file['tstt']['nodes']['coordinates'][()] - lower_left_corner = coords.min(axis=0) - upper_right_corner = coords.max(axis=0) - return openmc.BoundingBox(lower_left_corner, upper_right_corner) - - @property - def filename(self): - return self._filename - - @filename.setter - def filename(self, val: cv.PathLike): - cv.check_type('DAGMC filename', val, cv.PathLike) - self._filename = input_path(val) - - @property - def auto_geom_ids(self): - return self._auto_geom_ids - - @auto_geom_ids.setter - def auto_geom_ids(self, val): - cv.check_type('DAGMC automatic geometry ids', val, bool) - self._auto_geom_ids = val - - @property - def auto_mat_ids(self): - return self._auto_mat_ids - - @auto_mat_ids.setter - def auto_mat_ids(self, val): - cv.check_type('DAGMC automatic material ids', val, bool) - self._auto_mat_ids = val - - @property - def material_names(self): - dagmc_file_contents = h5py.File(self.filename) - material_tags_hex = dagmc_file_contents['/tstt/tags/NAME'].get( - 'values') - material_tags_ascii = [] - for tag in material_tags_hex: - candidate_tag = tag.tobytes().decode().replace('\x00', '') - # tags might be for temperature or reflective surfaces - if candidate_tag.startswith('mat:'): - # removes first 4 characters as openmc.Material name should be - # set without the 'mat:' part of the tag - material_tags_ascii.append(candidate_tag[4:]) - - return sorted(set(material_tags_ascii)) - - def get_all_cells(self, memo=None): - return {} - - def get_all_materials(self, memo=None): - return {} + def bounding_box(self) -> openmc.BoundingBox: + regions = [c.region for c in self.cells.values() + if c.region is not None] + if regions: + return openmc.Union(regions).bounding_box + else: + return openmc.BoundingBox.infinite() - def _n_geom_elements(self, geom_type): - """ - Helper function for retrieving the number geometric entities in a DAGMC - file + @classmethod + def from_hdf5(cls, group, cells): + """Create universe from HDF5 group Parameters ---------- - geom_type : str - The type of geometric entity to count. One of {'Volume', 'Surface'}. Returns - the runtime number of voumes in the DAGMC model (includes implicit complement). + group : h5py.Group + Group in HDF5 file + cells : dict + Dictionary mapping cell IDs to instances of :class:`openmc.Cell`. Returns ------- - int - Number of geometry elements of the specified type - """ - cv.check_value('geometry type', geom_type, ('volume', 'surface')) - - def decode_str_tag(tag_val): - return tag_val.tobytes().decode().replace('\x00', '') - - with h5py.File(self.filename) as dagmc_file: - category_data = dagmc_file['tstt/tags/CATEGORY/values'] - category_strs = map(decode_str_tag, category_data) - n = sum([v == geom_type.capitalize() for v in category_strs]) - - # check for presence of an implicit complement in the file and - # increment the number of cells if it doesn't exist - if geom_type == 'volume': - name_data = dagmc_file['tstt/tags/NAME/values'] - name_strs = map(decode_str_tag, name_data) - if not sum(['impl_complement' in n for n in name_strs]): - n += 1 - return n + openmc.Universe + Universe instance - @property - def n_cells(self): - return self._n_geom_elements('volume') + """ + universe_id = int(group.name.split('/')[-1].lstrip('universe ')) + cell_ids = group['cells'][()] - @property - def n_surfaces(self): - return self._n_geom_elements('surface') + # Create this Universe + universe = cls(universe_id) - def create_xml_subelement(self, xml_element, memo=None): - if memo is None: - memo = set() + # Add each Cell to the Universe + for cell_id in cell_ids: + universe.add_cell(cells[cell_id]) - if self in memo: - return + return universe - memo.add(self) - # Set xml element values - dagmc_element = ET.Element('dagmc_universe') - dagmc_element.set('id', str(self.id)) - - if self.auto_geom_ids: - dagmc_element.set('auto_geom_ids', 'true') - if self.auto_mat_ids: - dagmc_element.set('auto_mat_ids', 'true') - dagmc_element.set('filename', str(self.filename)) - xml_element.append(dagmc_element) - - def bounding_region( - self, - bounded_type: str = 'box', - boundary_type: str = 'vacuum', - starting_id: int = 10000, - padding_distance: float = 0. - ): - """Creates a either a spherical or box shaped bounding region around - the DAGMC geometry. - - .. versionadded:: 0.13.1 + def add_cell(self, cell): + """Add a cell to the universe. Parameters ---------- - bounded_type : str - The type of bounding surface(s) to use when constructing the region. - Options include a single spherical surface (sphere) or a rectangle - made from six planes (box). - boundary_type : str - Boundary condition that defines the behavior for particles hitting - the surface. Defaults to vacuum boundary condition. Passed into the - surface construction. - starting_id : int - Starting ID of the surface(s) used in the region. For bounded_type - 'box', the next 5 IDs will also be used. Defaults to 10000 to reduce - the chance of an overlap of surface IDs with the DAGMC geometry. - padding_distance : float - Distance between the bounding region surfaces and the minimal - bounding box. Allows for the region to be larger than the DAGMC - geometry. + cell : openmc.Cell + Cell to add - Returns - ------- - openmc.Region - Region instance """ - check_type('boundary type', boundary_type, str) - check_value('boundary type', boundary_type, _BOUNDARY_TYPES) - check_type('starting surface id', starting_id, Integral) - check_type('bounded type', bounded_type, str) - check_value('bounded type', bounded_type, ('box', 'sphere')) - - bbox = self.bounding_box.expand(padding_distance, True) - - if bounded_type == 'sphere': - radius = np.linalg.norm(bbox.upper_right - bbox.center) - bounding_surface = openmc.Sphere( - surface_id=starting_id, - x0=bbox.center[0], - y0=bbox.center[1], - z0=bbox.center[2], - boundary_type=boundary_type, - r=radius, - ) - - return -bounding_surface - - if bounded_type == 'box': - # defines plane surfaces for all six faces of the bounding box - lower_x = openmc.XPlane(bbox[0][0], surface_id=starting_id) - upper_x = openmc.XPlane(bbox[1][0], surface_id=starting_id+1) - lower_y = openmc.YPlane(bbox[0][1], surface_id=starting_id+2) - upper_y = openmc.YPlane(bbox[1][1], surface_id=starting_id+3) - lower_z = openmc.ZPlane(bbox[0][2], surface_id=starting_id+4) - upper_z = openmc.ZPlane(bbox[1][2], surface_id=starting_id+5) - - region = +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z - - for surface in region.get_surfaces().values(): - surface.boundary_type = boundary_type - - return region - - def bounded_universe(self, bounding_cell_id=10000, **kwargs): - """Returns an openmc.Universe filled with this DAGMCUniverse and bounded - with a cell. Defaults to a box cell with a vacuum surface however this - can be changed using the kwargs which are passed directly to - DAGMCUniverse.bounding_region(). + if not isinstance(cell, openmc.Cell): + msg = f'Unable to add a Cell to Universe ID="{self._id}" since ' \ + f'"{cell}" is not a Cell' + raise TypeError(msg) - Parameters - ---------- - bounding_cell_id : int - The cell ID number to use for the bounding cell, defaults to 10000 to reduce - the chance of overlapping ID numbers with the DAGMC geometry. + cell_id = cell.id - Returns - ------- - openmc.Universe - Universe instance - """ - bounding_cell = openmc.Cell( - fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs)) - return openmc.Universe(cells=[bounding_cell]) + if cell_id not in self._cells: + self._cells[cell_id] = cell - @classmethod - def from_hdf5(cls, group): - """Create DAGMC universe from HDF5 group + def remove_cell(self, cell): + """Remove a cell from the universe. Parameters ---------- - group : h5py.Group - Group in HDF5 file - - Returns - ------- - openmc.DAGMCUniverse - DAGMCUniverse instance + cell : openmc.Cell + Cell to remove """ - id = int(group.name.split('/')[-1].lstrip('universe ')) - fname = group['filename'][()].decode() - name = group['name'][()].decode() if 'name' in group else None - out = cls(fname, universe_id=id, name=name) - - out.auto_geom_ids = bool(group.attrs['auto_geom_ids']) - out.auto_mat_ids = bool(group.attrs['auto_mat_ids']) - - return out - - @classmethod - def from_xml_element(cls, elem): - """Generate DAGMC universe from XML element + if not isinstance(cell, openmc.Cell): + msg = f'Unable to remove a Cell from Universe ID="{self._id}" ' \ + f'since "{cell}" is not a Cell' + raise TypeError(msg) - Parameters - ---------- - elem : lxml.etree._Element - `` element + # If the Cell is in the Universe's list of Cells, delete it + self._cells.pop(cell.id, None) - Returns - ------- - openmc.DAGMCUniverse - DAGMCUniverse instance + def create_xml_subelement(self, xml_element, memo=None): + if memo is None: + memo = set() - """ - id = int(get_text(elem, 'id')) - fname = get_text(elem, 'filename') + # Iterate over all Cells + for cell in self._cells.values(): - out = cls(fname, universe_id=id) + # If the cell was already written, move on + if cell in memo: + continue - name = get_text(elem, 'name') - if name is not None: - out.name = name + memo.add(cell) - out.auto_geom_ids = bool(elem.get('auto_geom_ids')) - out.auto_mat_ids = bool(elem.get('auto_mat_ids')) + # Create XML subelement for this Cell + cell_element = cell.create_xml_subelement(xml_element, memo) - return out + # Append the Universe ID to the subelement and add to Element + cell_element.set("universe", str(self._id)) + xml_element.append(cell_element) def _partial_deepcopy(self): - """Clone all of the openmc.DAGMCUniverse object's attributes except for - its cells, as they are copied within the clone function. This should - only to be used within the openmc.UniverseBase.clone() context. + """Clone all of the openmc.Universe object's attributes except for its cells, + as they are copied within the clone function. This should only to be + used within the openmc.UniverseBase.clone() context. """ - clone = openmc.DAGMCUniverse(name=self.name, filename=self.filename) + clone = openmc.Universe(name=self.name) clone.volume = self.volume - clone.auto_geom_ids = self.auto_geom_ids - clone.auto_mat_ids = self.auto_mat_ids return clone diff --git a/src/cell.cpp b/src/cell.cpp index 88876678706..d4d28fb70e6 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -252,12 +252,12 @@ void Cell::to_hdf5(hid_t cell_group) const // default constructor CSGCell::CSGCell() { - geom_type_ = GeometryType::CSG; + geom_type() = GeometryType::CSG; } CSGCell::CSGCell(pugi::xml_node cell_node) { - geom_type_ = GeometryType::CSG; + geom_type() = GeometryType::CSG; if (check_for_node(cell_node, "id")) { id_ = std::stoi(get_node_value(cell_node, "id")); diff --git a/src/dagmc.cpp b/src/dagmc.cpp index b79676c3626..1f73287ccce 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -72,6 +72,25 @@ DAGUniverse::DAGUniverse(pugi::xml_node node) adjust_material_ids_ = get_node_value_bool(node, "auto_mat_ids"); } + // get material assignment overloading + if (check_for_node(node, "material_overrides")) { + auto mat_node = node.child("material_overrides"); + // loop over all attributes (each attribute corresponds to a material) + for (pugi::xml_attribute attr = mat_node.first_attribute(); attr; + attr = attr.next_attribute()) { + // Store assignment reference name + std::string mat_ref_assignment = attr.name(); + + // Get mat name for each assignement instances + std::stringstream iss {attr.value()}; + vector instance_mats = split(iss.str(), ';'); + + // Store mat name for each instances + instance_material_overrides.insert( + std::make_pair(mat_ref_assignment, instance_mats)); + } + } + initialize(); } @@ -211,14 +230,17 @@ void DAGUniverse::init_geometry() if (mat_str == "graveyard") { graveyard = vol_handle; } - // material void checks if (mat_str == "void" || mat_str == "vacuum" || mat_str == "graveyard") { c->material_.push_back(MATERIAL_VOID); } else { - if (uses_uwuw()) { + if (instance_material_overrides.count("id_" + std::to_string(c->id_)) || + instance_material_overrides.count( + mat_str)) { // Check for material override + override_assign_material(mat_str, vol_handle, c); + } else if (uses_uwuw()) { // UWUW assignement uwuw_assign_material(vol_handle, c); - } else { + } else { // legacy assignement legacy_assign_material(mat_str, c); } } @@ -609,6 +631,30 @@ void DAGUniverse::uwuw_assign_material( fatal_error("DAGMC was not configured with UWUW."); #endif // OPENMC_UWUW } + +void DAGUniverse::override_assign_material(std::string key, + moab::EntityHandle vol_handle, std::unique_ptr& c) const +{ + + // if Cell ID matches an override key, use it to override the material + // assignment else if UWUW is used, get the material assignment from the DAGMC + // metadata + std::cout << "XML instance" + << instance_material_overrides.count("id_" + std::to_string(c->id_)) + << std::endl; + if (instance_material_overrides.count("id_" + std::to_string(c->id_))) { + key = "id_" + std::to_string(c->id_); + } else if (uses_uwuw()) { + key = dmd_ptr->volume_material_property_data_eh[vol_handle]; + } + + // Override the material assignment for each cell instance using the legacy + // assignement + for (auto mat_str_instance : instance_material_overrides.at(key)) { + legacy_assign_material(mat_str_instance, c); + } +} + //============================================================================== // DAGMC Cell implementation //============================================================================== @@ -616,7 +662,7 @@ void DAGUniverse::uwuw_assign_material( DAGCell::DAGCell(std::shared_ptr dag_ptr, int32_t dag_idx) : Cell {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx) { - geom_type_ = GeometryType::DAG; + geom_type() = GeometryType::DAG; }; std::pair DAGCell::distance( @@ -719,7 +765,7 @@ BoundingBox DAGCell::bounding_box() const DAGSurface::DAGSurface(std::shared_ptr dag_ptr, int32_t dag_idx) : Surface {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx) { - geom_type_ = GeometryType::DAG; + geom_type() = GeometryType::DAG; } // empty constructor moab::EntityHandle DAGSurface::mesh_handle() const @@ -818,12 +864,56 @@ int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ) return univp->cell_index(new_vol); } +extern "C" int openmc_dagmc_universe_get_cell_ids( + int32_t univ_id, int32_t* ids, size_t* n) +{ + // make sure the universe id is a DAGMC Universe + const auto& univ = model::universes[model::universe_map[univ_id]]; + if (univ->geom_type() != GeometryType::DAG) { + fatal_error( + "Universe " + std::to_string(univ_id) + " is not a DAGMC Universe!"); + } + + std::vector dag_cell_ids; + for (const auto& cell_index : univ->cells_) { + const auto& cell = model::cells[cell_index]; + if (cell->geom_type() == GeometryType::DAG) + dag_cell_ids.push_back(cell->id_); + } + std::copy(dag_cell_ids.begin(), dag_cell_ids.end(), ids); + *n = dag_cell_ids.size(); + return 0; +} + +extern "C" int openmc_dagmc_universe_get_num_cells(int32_t univ_id, size_t* n) +{ + // make sure the universe id is a DAGMC Universe + const auto& univ = model::universes[model::universe_map[univ_id]]; + if (univ->geom_type() != GeometryType::DAG) { + fatal_error( + "Universe " + std::to_string(univ_id) + " is not a DAGMC Universe"); + } + *n = univ->cells_.size(); + return 0; +} + } // namespace openmc #else namespace openmc { +extern "C" int openmc_dagmc_universe_get_cell_ids( + int32_t univ_id, int32_t* ids, size_t* n) +{ + fatal_error("OpenMC was not configured with DAGMC"); +}; + +extern "C" int openmc_dagmc_universe_get_num_cells(int32_t univ_id, size_t* n) +{ + fatal_error("OpenMC was not configured with DAGMC"); +}; + void read_dagmc_universes(pugi::xml_node node) { if (check_for_node(node, "dagmc_universe")) { diff --git a/src/particle.cpp b/src/particle.cpp index 64c50c9438f..0ea8650143d 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -533,7 +533,7 @@ void Particle::cross_surface(const Surface& surf) // if we're crossing a CSG surface, make sure the DAG history is reset #ifdef DAGMC - if (surf.geom_type_ == GeometryType::CSG) + if (surf.geom_type() == GeometryType::CSG) history().reset(); #endif @@ -548,7 +548,7 @@ void Particle::cross_surface(const Surface& surf) #ifdef DAGMC // in DAGMC, we know what the next cell should be - if (surf.geom_type_ == GeometryType::DAG) { + if (surf.geom_type() == GeometryType::DAG) { int32_t i_cell = next_cell(std::abs(surface()), cell_last(n_coord() - 1), lowest_coord().universe) - 1; @@ -668,7 +668,8 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u) // the lower universes. // (unless we're using a dagmc model, which has exactly one universe) n_coord() = 1; - if (surf.geom_type_ != GeometryType::DAG && !neighbor_list_find_cell(*this)) { + if (surf.geom_type() != GeometryType::DAG && + !neighbor_list_find_cell(*this)) { mark_as_lost("Couldn't find particle after reflecting from surface " + std::to_string(surf.id_) + "."); return; diff --git a/src/plot.cpp b/src/plot.cpp index 348138570c1..43f25a9a32f 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1301,7 +1301,7 @@ void ProjectionPlot::create_output() const int32_t i_surface = std::abs(p.surface()) - 1; if (i_surface > 0 && - model::surfaces[i_surface]->geom_type_ == GeometryType::DAG) { + model::surfaces[i_surface]->geom_type() == GeometryType::DAG) { #ifdef DAGMC int32_t i_cell = next_cell(i_surface, p.cell_last(p.n_coord() - 1), p.lowest_coord().universe); diff --git a/src/string_utils.cpp b/src/string_utils.cpp index 74f048e8d24..454f2dca0be 100644 --- a/src/string_utils.cpp +++ b/src/string_utils.cpp @@ -71,6 +71,31 @@ vector split(const std::string& in) return out; } +vector split(const std::string& in, char delim) +{ + vector out; + + for (int i = 0; i < in.size();) { + // Increment i until we find a non-delimiter character. + if (in[i] == delim) { + i++; + + } else { + // Find the next delimiter character at j. + int j = i + 1; + while (j < in.size() && in[j] != delim) { + j++; + } + + // Push-back everything between i and j. + out.push_back(in.substr(i, j - i)); + i = j + 1; // j is delimiter so leapfrog to j+1 + } + } + + return out; +} + bool ends_with(const std::string& value, const std::string& ending) { if (ending.size() > value.size()) diff --git a/src/surface.cpp b/src/surface.cpp index 50ef2a12830..dbcaf849848 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -165,9 +165,9 @@ void Surface::to_hdf5(hid_t group_id) const { hid_t surf_group = create_group(group_id, fmt::format("surface {}", id_)); - if (geom_type_ == GeometryType::DAG) { + if (geom_type() == GeometryType::DAG) { write_string(surf_group, "geom_type", "dagmc", false); - } else if (geom_type_ == GeometryType::CSG) { + } else if (geom_type() == GeometryType::CSG) { write_string(surf_group, "geom_type", "csg", false); if (bc_) { @@ -189,11 +189,11 @@ void Surface::to_hdf5(hid_t group_id) const CSGSurface::CSGSurface() : Surface {} { - geom_type_ = GeometryType::CSG; + geom_type() = GeometryType::CSG; }; CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface {surf_node} { - geom_type_ = GeometryType::CSG; + geom_type() = GeometryType::CSG; }; //============================================================================== diff --git a/tests/unit_tests/dagmc/test_model.py b/tests/unit_tests/dagmc/test_model.py new file mode 100644 index 00000000000..57c07699b9e --- /dev/null +++ b/tests/unit_tests/dagmc/test_model.py @@ -0,0 +1,139 @@ +import pkg_resources +from pathlib import Path + +import numpy as np +import pytest + +import openmc +from openmc import ZPlane, YPlane, XPlane, Cell + +pytestmark = pytest.mark.skipif( + not openmc.lib._dagmc_enabled(), + reason="DAGMC CAD geometry is not enabled.") + + +def set_dagmc_model(): + PITCH = 1.26 + + mats = {} + mats["no-void fuel"] = openmc.Material(1, "no-void fuel") + mats["no-void fuel"].add_nuclide("U235", 0.03) + mats["no-void fuel"].add_nuclide("U238", 0.97) + mats["no-void fuel"].add_nuclide("O16", 2.0) + mats["no-void fuel"].set_density("g/cm3", 10.0) + mats["no-void fuel"].name = "no-void fuel" + + mats["41"] = openmc.Material(name="41") + mats["41"].add_nuclide("H1", 2.0) + mats["41"].add_element("O", 1.0) + mats["41"].set_density("g/cm3", 1.0) + mats["41"].add_s_alpha_beta("c_H_in_H2O") + mats["41"].name = "41" + + p = pkg_resources.resource_filename(__name__, "dagmc.h5m") + + daguniv = openmc.DAGMCUniverse(p,auto_geom_ids=True,) + + def pattern(center, bc): + bc_ = { + "top": "transmission", + "bottom": "transmission", + "left": "transmission", + "right": "transmission", + } + bc_ |= bc + box = ( + -XPlane(center[0] + PITCH / 2, boundary_type=bc_["right"]) + & +XPlane(center[0] - PITCH / 2, boundary_type=bc_["left"]) + & -YPlane(center[1] + PITCH / 2, boundary_type=bc_["top"]) + & +YPlane(center[1] - PITCH / 2, boundary_type=bc_["bottom"]) + & -ZPlane(5, boundary_type="reflective") + & +ZPlane(-5, boundary_type="reflective") + ) + cell = Cell(region=box, fill=daguniv) + cell.translation = [*center, 0] + return [cell] + + root = openmc.Universe( + cells=[ + *pattern((-PITCH / 2, -PITCH / 2), + bc={"left": "reflective", "bottom": "reflective"}), + *pattern((-PITCH / 2, PITCH / 2), + bc={"left": "reflective", "top": "reflective"}), + *pattern((PITCH / 2, PITCH / 2), + bc={"right": "reflective", "top": "reflective"}), + *pattern((PITCH / 2, -PITCH / 2), + bc={"right": "reflective", "bottom": "reflective"}), + ] + ) + + point = openmc.stats.Point((0, 0, 0)) + source = openmc.IndependentSource(space=point) + + settings = openmc.Settings() + settings.source = source + settings.batches = 100 + settings.inactive = 10 + settings.particles = 1000 + + ll, ur = root.bounding_box + mat_vol = openmc.VolumeCalculation([mats["no-void fuel"]], 1000000, ll, ur) + cell_vol = openmc.VolumeCalculation( + list(root.cells.values()), 1000000, ll, ur) + settings.volume_calculations = [mat_vol, cell_vol] + + model = openmc.Model() + model.materials = openmc.Materials(mats.values()) + model.geometry = openmc.Geometry(root=root) + model.settings = settings + return model + + +def test_model_differentiate_depletable_with_DAGMC(): + model = set_dagmc_model() + + p = Path("differentiate_depletable_mats/divide_equally") + p.mkdir(parents=True, exist_ok=True) + model.init_lib() + model.sync_dagmc_universe() + model.calculate_volumes(cwd=p) + + # Get the volume of the no-void fuel material before differentiation + volume_before = np.sum([m.volume for m in model.materials if m.name == "no-void fuel"]) + + # Differentiate the depletable materials + model.differentiate_depletable_mats(diff_volume_method="divide equally") + # Get the volume of the no-void fuel material after differentiation + volume_after = np.sum([m.volume for m in model.materials if "fuel" in m.name]) + assert np.isclose(volume_before, volume_after) + + assert len(model.materials) == 4*2 +1 + model.finalize_lib() + + +def test_model_differentiate_with_DAGMC(): + + model = set_dagmc_model() + p = Path("differentiate_depletable_mats/divide_equally") + p.mkdir(parents=True, exist_ok=True) + model.init_lib() + model.sync_dagmc_universe() + model.calculate_volumes(cwd=p) + # Get the volume of the no-void fuel material before differentiation + volume_before = np.sum([m.volume for m in model.materials if m.name == "no-void fuel"]) + + # Differentiate all the materials + model.differentiate_mats(depletable_only=False) + + # Get the volume of the no-void fuel material after differentiation + mat_list = [m for m in model.materials] + mat_vol = openmc.VolumeCalculation(mat_list, 1000000, ll, ur) + cell_vol = openmc.VolumeCalculation(list(root.cells.values()), 1000000, ll, ur) + model.settings.volume_calculations = [mat_vol, cell_vol] + model.init_lib() # need to reinitialize the lib after differentiating the materials + model.calculate_volumes(cwd=p) + volume_after = np.sum([m.volume for m in model.materials if "fuel" in m.name]) + assert np.isclose(volume_before, volume_after) + + assert len(model.materials) == 4*2 + 4 + model.finalize_lib()