Merge pull request #11 from openmc-dev/develop

Update with latest develop before open pr

Author: church89

include/openmc/dagmc.h

@@ -37,6 +37,8 @@ class DAGSurface : public Surface {
 public:
   DAGSurface(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx);
 
+  moab::EntityHandle mesh_handle() const;
+
   double evaluate(Position r) const override;
   double distance(Position r, Direction u, bool coincident) const override;
   Direction normal(Position r) const override;
@@ -45,7 +47,7 @@ class DAGSurface : public Surface {
   inline void to_hdf5_inner(hid_t group_id) const override {};
 
   // Accessor methods
-  const std::shared_ptr<moab::DagMC>& dagmc_ptr() const { return dagmc_ptr_; }
+  moab::DagMC* dagmc_ptr() const { return dagmc_ptr_.get(); }
   int32_t dag_index() const { return dag_index_; }
 
 private:
@@ -57,6 +59,8 @@ class DAGCell : public Cell {
 public:
   DAGCell(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx);
 
+  moab::EntityHandle mesh_handle() const;
+
   bool contains(Position r, Direction u, int32_t on_surface) const override;
 
   std::pair<double, int32_t> distance(
@@ -67,7 +71,7 @@ class DAGCell : public Cell {
   void to_hdf5_inner(hid_t group_id) const override;
 
   // Accessor methods
-  const std::shared_ptr<moab::DagMC>& dagmc_ptr() const { return dagmc_ptr_; }
+  moab::DagMC* dagmc_ptr() const { return dagmc_ptr_.get(); }
   int32_t dag_index() const { return dag_index_; }
 
 private:
@@ -122,6 +126,16 @@ class DAGUniverse : public Universe {
   void legacy_assign_material(
     std::string mat_string, std::unique_ptr<DAGCell>& 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
+  int32_t cell_index(moab::EntityHandle vol) const;
+
+  //! Return the index into the model surfaces vector for a given DAGMC surface
+  //! handle in the universe
+  //! \param[in] surf MOAB handle to the DAGMC surface set
+  int32_t surface_index(moab::EntityHandle surf) const;
+
   //! Generate a string representing the ranges of IDs present in the DAGMC
   //! model. Contiguous chunks of IDs are represented as a range (i.e. 1-10). If
   //! there is a single ID a chunk, it will be represented as a single number
@@ -142,6 +156,7 @@ class DAGUniverse : public Universe {
                             //!< universe in OpenMC's surface vector
 
   // Accessors
+  moab::DagMC* dagmc_ptr() const { return dagmc_instance_.get(); }
   bool has_graveyard() const { return has_graveyard_; }
 
 private:

openmc/deplete/chain.py

@@ -686,16 +686,16 @@ def form_matrix(self, rates, fission_yields=None):
         dict.update(matrix_dok, matrix)
         return matrix_dok.tocsr()
 
-    def form_rr_term(self, transfer_rates, materials):
+    def form_rr_term(self, tr_rates, mats):
         """Function to form the transfer rate term matrices.
 
         .. versionadded:: 0.13.4
 
         Parameters
         ----------
-        transfer_rates : openmc.deplete.TransferRates
+        tr_rates : openmc.deplete.TransferRates
             Instance of openmc.deplete.TransferRates
-        materials : string or two-tuple of strings
+        mats : string or two-tuple of strings
             Two cases are possible:
 
             1) Material ID as string:
@@ -719,25 +719,27 @@ def form_rr_term(self, transfer_rates, materials):
         """
         matrix = defaultdict(float)
 
-        for i, nuclide in enumerate(self.nuclides):
-            element = re.split(r'\d+', nuclide.name)[0]
+        for i, nuc in enumerate(self.nuclides):
+            elm = re.split(r'\d+', nuc.name)[0]
             # Build transfer terms matrices
-            if isinstance(materials, str):
-                material = materials
-                components = transfer_rates.get_components(material)
-                if element in components:
-                    matrix[i, i] = transfer_rates.get_transfer_rate(material, element)
-                elif nuclide.name in components:
-                    matrix[i, i] = transfer_rates.get_transfer_rate(material, nuclide.name)
+            if isinstance(mats, str):
+                mat = mats
+                components = tr_rates.get_components(mat)
+                if elm in components:
+                    matrix[i, i] = sum(tr_rates.get_transfer_rate(mat, elm))
+                elif nuc.name in components:
+                    matrix[i, i] = sum(tr_rates.get_transfer_rate(mat, nuc.name))
                 else:
                     matrix[i, i] = 0.0
             #Build transfer terms matrices
-            elif isinstance(materials, tuple):
-                destination_material, material = materials
-                if transfer_rates.get_destination_material(material, element) == destination_material:
-                    matrix[i, i] = transfer_rates.get_transfer_rate(material, element)
-                elif transfer_rates.get_destination_material(material, nuclide.name) == destination_material:
-                    matrix[i, i] = transfer_rates.get_transfer_rate(material, nuclide.name)
+            elif isinstance(mats, tuple):
+                dest_mat, mat = mats
+                if dest_mat in tr_rates.get_destination_material(mat, elm):
+                    dest_mat_idx = tr_rates.get_destination_material(mat, elm).index(dest_mat)
+                    matrix[i, i] = tr_rates.get_transfer_rate(mat, elm)[dest_mat_idx]
+                elif dest_mat in tr_rates.get_destination_material(mat, nuc.name):
+                    dest_mat_idx = tr_rates.get_destination_material(mat, nuc.name).index(dest_mat)
+                    matrix[i, i] = tr_rates.get_transfer_rate(mat, nuc.name)[dest_mat_idx]
                 else:
                     matrix[i, i] = 0.0
             #Nothing else is allowed

openmc/deplete/coupled_operator.py

@@ -416,7 +416,7 @@ def _generate_materials_xml(self):
         for mat in self.materials:
             mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
 
-        self.materials.export_to_xml()
+        self.materials.export_to_xml(nuclides_to_ignore=self._decay_nucs)
 
     def __call__(self, vec, source_rate):
         """Runs a simulation.

openmc/deplete/openmc_operator.py

@@ -6,6 +6,7 @@
 """
 
 from abc import abstractmethod
+from warnings import warn
 from typing import List, Tuple, Dict
 
 import numpy as np
@@ -133,6 +134,18 @@ def __init__(
         # This nuclides variables contains every nuclides
         # for which there is an entry in the micro_xs parameter
         openmc.reset_auto_ids()
+
+        self.nuclides_with_data = self._get_nuclides_with_data(
+            self.cross_sections)
+
+        # Select nuclides with data that are also in the chain
+        self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
+                               if nuc.name in self.nuclides_with_data]
+
+        # Select nuclides without data that are also in the chain
+        self._decay_nucs = [nuc.name for nuc in self.chain.nuclides
+                            if nuc.name not in self.nuclides_with_data]
+
         self.burnable_mats, volumes, all_nuclides = self._get_burnable_mats()
         self.local_mats = _distribute(self.burnable_mats)
 
@@ -142,13 +155,6 @@ def __init__(
         if self.prev_res is not None:
             self._load_previous_results()
 
-        self.nuclides_with_data = self._get_nuclides_with_data(
-            self.cross_sections)
-
-        # Select nuclides with data that are also in the chain
-        self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
-                               if nuc.name in self.nuclides_with_data]
-
         # Extract number densities from the geometry / previous depletion run
         self._extract_number(self.local_mats,
                              volumes,
@@ -171,7 +177,7 @@ def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]:
         Returns
         -------
         burnable_mats : list of str
-            List of burnable material IDs
+            list of burnable material IDs
         volume : dict of str to float
             Volume of each material in [cm^3]
         nuclides : list of str
@@ -188,7 +194,13 @@ def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]:
         # Iterate once through the geometry to get dictionaries
         for mat in self.materials:
             for nuclide in mat.get_nuclides():
-                model_nuclides.add(nuclide)
+                if nuclide in self.nuclides_with_data or self._decay_nucs:
+                    model_nuclides.add(nuclide)
+                else:
+                    msg = (f"Nuclilde {nuclide} in material {mat.id} is not "
+                           "present in the depletion chain and has no cross "
+                           "section data.")
+                    raise warn(msg)
             if mat.depletable:
                 burnable_mats.add(str(mat.id))
                 if mat.volume is None:

openmc/deplete/transfer_rates.py

@@ -91,13 +91,13 @@ def get_transfer_rate(self, material, component):
 
         Returns
         -------
-        transfer_rate : float
-            Transfer rate value
+        transfer_rate : list of floats
+            Transfer rate values
 
         """
         material_id = self._get_material_id(material)
         check_type('component', component, str)
-        return self.transfer_rates[material_id][component][0]
+        return [i[0] for i in self.transfer_rates[material_id][component]]
 
     def get_destination_material(self, material, component):
         """Return destination material for given material and
@@ -112,15 +112,17 @@ def get_destination_material(self, material, component):
 
         Returns
         -------
-        destination_material_id : str
+        destination_material_id : list of str
             Depletable material ID to where the element or nuclide gets
             transferred
 
         """
         material_id = self._get_material_id(material)
         check_type('component', component, str)
         if component in self.transfer_rates[material_id]:
-            return self.transfer_rates[material_id][component][1]
+            return [i[1] for i in self.transfer_rates[material_id][component]]
+        else:
+            return []
 
     def get_components(self, material):
         """Extract removing elements and/or nuclides for a given material
@@ -217,7 +219,11 @@ def set_transfer_rate(self, material, components, transfer_rate,
                                          f'where element {element} already has '
                                          'a transfer rate.')
 
-            self.transfer_rates[material_id][component] = \
-                transfer_rate / unit_conv, destination_material_id
+            if component in self.transfer_rates[material_id]:
+                self.transfer_rates[material_id][component].append(
+                    (transfer_rate / unit_conv, destination_material_id))
+            else:
+                self.transfer_rates[material_id][component] = [
+                    (transfer_rate / unit_conv, destination_material_id)]
             if destination_material_id is not None:
                 self.index_transfer.add((destination_material_id, material_id))

openmc/material.py

@@ -1250,7 +1250,7 @@ def clone(self, memo: Optional[dict] = None) -> Material:
 
         return memo[self]
 
-    def _get_nuclide_xml(self, nuclide: str) -> ET.Element:
+    def _get_nuclide_xml(self, nuclide: NuclideTuple) -> ET.Element:
         xml_element = ET.Element("nuclide")
         xml_element.set("name", nuclide.name)
 
@@ -1267,15 +1267,28 @@ def _get_macroscopic_xml(self, macroscopic: str) -> ET.Element:
 
         return xml_element
 
-    def _get_nuclides_xml(self, nuclides: typing.Iterable[str]) -> List[ET.Element]:
+    def _get_nuclides_xml(
+            self, nuclides: typing.Iterable[NuclideTuple],
+            nuclides_to_ignore: Optional[typing.Iterable[str]] = None)-> List[ET.Element]:
         xml_elements = []
-        for nuclide in nuclides:
-            xml_elements.append(self._get_nuclide_xml(nuclide))
+
+        # Remove any nuclides to ignore from the XML export
+        if nuclides_to_ignore:
+            nuclides = [nuclide for nuclide in nuclides if nuclide.name not in nuclides_to_ignore]
+
+        xml_elements = [self._get_nuclide_xml(nuclide) for nuclide in nuclides]
+
         return xml_elements
 
-    def to_xml_element(self) -> ET.Element:
+    def to_xml_element(
+            self, nuclides_to_ignore: Optional[typing.Iterable[str]] = None) -> ET.Element:
         """Return XML representation of the material
 
+        Parameters
+        ----------
+        nuclides_to_ignore : list of str
+            Nuclides to ignore when exporting to XML.
+
         Returns
         -------
         element : lxml.etree._Element
@@ -1320,7 +1333,8 @@ def to_xml_element(self) -> ET.Element:
 
         if self._macroscopic is None:
             # Create nuclide XML subelements
-            subelements = self._get_nuclides_xml(self._nuclides)
+            subelements = self._get_nuclides_xml(self._nuclides,
+                                                 nuclides_to_ignore=nuclides_to_ignore)
             for subelement in subelements:
                 element.append(subelement)
         else:
@@ -1576,7 +1590,8 @@ def make_isotropic_in_lab(self):
         for material in self:
             material.make_isotropic_in_lab()
 
-    def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_indent=True):
+    def _write_xml(self, file, header=True, level=0, spaces_per_level=2,
+                   trailing_indent=True, nuclides_to_ignore=None):
         """Writes XML content of the materials to an open file handle.
 
         Parameters
@@ -1591,6 +1606,8 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in
             Number of spaces per indentation
         trailing_indentation : bool
             Whether or not to write a trailing indentation for the materials element
+        nuclides_to_ignore : list of str
+            Nuclides to ignore when exporting to XML.
 
         """
         indentation = level*spaces_per_level*' '
@@ -1611,7 +1628,7 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in
 
         # Write the <material> elements.
         for material in sorted(self, key=lambda x: x.id):
-            element = material.to_xml_element()
+            element = material.to_xml_element(nuclides_to_ignore=nuclides_to_ignore)
             clean_indentation(element, level=level+1)
             element.tail = element.tail.strip(' ')
             file.write((level+1)*spaces_per_level*' ')
@@ -1626,13 +1643,16 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in
         if trailing_indent:
             file.write(indentation)
 
-    def export_to_xml(self, path: PathLike = 'materials.xml'):
+    def export_to_xml(self, path: PathLike = 'materials.xml',
+                      nuclides_to_ignore: Optional[typing.Iterable[str]] = None):
         """Export material collection to an XML file.
 
         Parameters
         ----------
         path : str
             Path to file to write. Defaults to 'materials.xml'.
+        nuclides_to_ignore : list of str
+            Nuclides to ignore when exporting to XML.
 
         """
         # Check if path is a directory
@@ -1645,7 +1665,7 @@ def export_to_xml(self, path: PathLike = 'materials.xml'):
         # one go.
         with open(str(p), 'w', encoding='utf-8',
                   errors='xmlcharrefreplace') as fh:
-            self._write_xml(fh)
+            self._write_xml(fh, nuclides_to_ignore=nuclides_to_ignore)
 
     @classmethod
     def from_xml_element(cls, elem) -> Material:

openmc/region.py

@@ -4,7 +4,6 @@
 
 import numpy as np
 
-from .checkvalue import check_type
 from .bounding_box import BoundingBox
 
 
@@ -362,6 +361,9 @@ class Intersection(Region, MutableSequence):
 
     def __init__(self, nodes):
         self._nodes = list(nodes)
+        for node in nodes:
+            if not isinstance(node, Region):
+                raise ValueError('Intersection operands must be of type Region')
 
     def __and__(self, other):
         new = Intersection(self)
@@ -450,6 +452,9 @@ class Union(Region, MutableSequence):
 
     def __init__(self, nodes):
         self._nodes = list(nodes)
+        for node in nodes:
+            if not isinstance(node, Region):
+                raise ValueError('Union operands must be of type Region')
 
     def __or__(self, other):
         new = Union(self)
@@ -566,7 +571,8 @@ def node(self):
 
     @node.setter
     def node(self, node):
-        check_type('node', node, Region)
+        if not isinstance(node, Region):
+            raise ValueError('Complement operand must be of type Region')
         self._node = node
 
     @property