Remove initial dilute nuclides in MicroXS (#2579)

openmc/deplete/microxs.py

@@ -5,14 +5,13 @@
 """
 
 import tempfile
-from copy import deepcopy
 
 from pandas import DataFrame, read_csv, Series
 import numpy as np
 
 from openmc.checkvalue import check_type, check_value, check_iterable_type
 from openmc.exceptions import DataError
-from openmc import StatePoint, Materials
+from openmc import StatePoint
 import openmc
 from .chain import Chain, REACTIONS
 from .coupled_operator import _find_cross_sections, _get_nuclides_with_data
@@ -31,9 +30,10 @@ class MicroXS(DataFrame):
     @classmethod
     def from_model(cls,
                    model,
-                   reaction_domain,
+                   domain,
+                   nuclides=None,
+                   reactions=None,
                    chain_file=None,
-                   dilute_initial=1.0e3,
                    energy_bounds=(0, 20e6),
                    run_kwargs=None):
         """Generate a one-group cross-section dataframe using OpenMC.
@@ -44,17 +44,19 @@ def from_model(cls,
         ----------
         model : openmc.Model
             OpenMC model object. Must contain geometry, materials, and settings.
-        reaction_domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh
+        domain : openmc.Material or openmc.Cell or openmc.Universe
             Domain in which to tally reaction rates.
+        nuclides : list of str
+            Nuclides to get cross sections for. If not specified, all burnable
+            nuclides from the depletion chain file are used.
+        reactions : list of str
+            Reactions to get cross sections for. If not specified, all neutron
+            reactions listed in the depletion chain file are used.
         chain_file : str, optional
             Path to the depletion chain XML file that will be used in depletion
             simulation. Used to determine cross sections for materials not
             present in the inital composition. Defaults to
             ``openmc.config['chain_file']``.
-        dilute_initial : float, optional
-            Initial atom density [atoms/cm^3] to add for nuclides that
-            are zero in initial condition to ensure they exist in the cross
-            section data. Only done for nuclides with reaction rates.
         energy_bound : 2-tuple of float, optional
             Bounds for the energy group.
         run_kwargs : dict, optional
@@ -66,29 +68,40 @@ def from_model(cls,
             Cross section data in [b]
 
         """
-        # Set up the reaction tallies
+        # Save any original tallies on the model
         original_tallies = model.tallies
-        original_materials = deepcopy(model.materials)
-        xs = {}
-        reactions, burnable_nucs, diluted_materials = cls._add_dilute_nuclides(
-            chain_file, model, dilute_initial)
-        model.materials = diluted_materials
 
+        # Determine what reactions and nuclides are available in chain
+        if chain_file is None:
+            chain_file = openmc.config.get('chain_file')
+            if chain_file is None:
+                raise DataError(
+                    "No depletion chain specified and could not find depletion "
+                    "chain in openmc.config['chain_file']"
+                )
+        chain = Chain.from_xml(chain_file)
+        if reactions is None:
+            reactions = chain.reactions
+        if not nuclides:
+            cross_sections = _find_cross_sections(model)
+            nuclides_with_data = _get_nuclides_with_data(cross_sections)
+            nuclides = [nuc.name for nuc in chain.nuclides
+                        if nuc.name in nuclides_with_data]
+
+        # Set up the reaction rate and flux tallies
         energy_filter = openmc.EnergyFilter(energy_bounds)
-        if isinstance(reaction_domain, openmc.Material):
-            domain_filter = openmc.MaterialFilter([reaction_domain])
-        elif isinstance(reaction_domain, openmc.Cell):
-            domain_filter = openmc.CellFilter([reaction_domain])
-        elif isinstance(reaction_domain, openmc.Universe):
-            domain_filter = openmc.UniverseFilter([reaction_domain])
+        if isinstance(domain, openmc.Material):
+            domain_filter = openmc.MaterialFilter([domain])
+        elif isinstance(domain, openmc.Cell):
+            domain_filter = openmc.CellFilter([domain])
+        elif isinstance(domain, openmc.Universe):
+            domain_filter = openmc.UniverseFilter([domain])
         else:
-            raise ValueError(f"Unsupported domain type: {type(reaction_domain)}")
+            raise ValueError(f"Unsupported domain type: {type(domain)}")
 
-        # TODO: Right now, we use all nuclides from the material but it probably
-        # should be based on the burnable nuclides
         rr_tally = openmc.Tally(name='MicroXS RR')
         rr_tally.filters = [domain_filter, energy_filter]
-        rr_tally.nuclides = reaction_domain.get_nuclides()
+        rr_tally.nuclides = nuclides
         rr_tally.multiply_density = False
         rr_tally.scores = reactions
 
@@ -126,68 +139,8 @@ def from_model(cls,
 
         # Revert to the original tallies and materials
         model.tallies = original_tallies
-        model.materials = original_materials
-
-        return cls(series)
-
-    @classmethod
-    def _add_dilute_nuclides(cls, chain_file, model, dilute_initial):
-        """
-        Add nuclides not present in burnable materials that have neutron data
-        and are present in the depletion chain to those materials. This allows
-        us to tally those specific nuclides for reactions to create one-group
-        cross sections.
-
-        Parameters
-        ----------
-        chain_file : str
-            Path to the depletion chain XML file that will be used in depletion
-            simulation. Used to determine cross sections for materials not
-            present in the inital composition.
-        model : openmc.Model
-            Model object
-        dilute_initial : float
-            Initial atom density [atoms/cm^3] to add for nuclides that
-            are zero in initial condition to ensure they exist in the cross
-            section data. Only done for nuclides with reaction rates.
 
-        Returns
-        -------
-        reactions : list of str
-            List of reaction names
-        diluted_materials : openmc.Materials
-            :class:`openmc.Materials` object with nuclides added to burnable
-            materials.
-        """
-        if chain_file is None:
-            chain_file = openmc.config.get('chain_file')
-            if chain_file is None:
-                raise DataError(
-                    "No depletion chain specified and could not find depletion "
-                    "chain in openmc.config['chain_file']"
-                )
-        chain = Chain.from_xml(chain_file)
-        reactions = chain.reactions
-        cross_sections = _find_cross_sections(model)
-        nuclides_with_data = _get_nuclides_with_data(cross_sections)
-        burnable_nucs = [nuc.name for nuc in chain.nuclides
-                         if nuc.name in nuclides_with_data]
-        diluted_materials = Materials()
-        for material in model.materials:
-            if material.depletable:
-                nuc_densities = material.get_nuclide_atom_densities()
-                dilute_density = 1.0e-24 * dilute_initial
-                material.set_density('sum')
-                for nuc, density in nuc_densities.items():
-                    material.remove_nuclide(nuc)
-                    material.add_nuclide(nuc, density)
-                for burn_nuc in burnable_nucs:
-                    if burn_nuc not in nuc_densities:
-                        material.add_nuclide(burn_nuc,
-                                                     dilute_density)
-            diluted_materials.append(material)
-
-        return reactions, burnable_nucs, diluted_materials
+        return cls(series).rename_axis('nuclide')
 
     @classmethod
     def from_array(cls, nuclides, reactions, data):

tests/regression_tests/microxs/test.py

@@ -40,14 +40,17 @@ def model():
 
     settings = openmc.Settings()
     settings.particles = 1000
-    settings.inactive = 10
-    settings.batches = 50
+    settings.inactive = 5
+    settings.batches = 10
 
     return openmc.Model(geometry, materials, settings)
 
 
 def test_from_model(model):
-    test_xs = MicroXS.from_model(model, model.materials[0], CHAIN_FILE)
+    fuel = model.materials[0]
+    nuclides = ['U234', 'U235', 'U238', 'U236', 'O16', 'O17', 'I135', 'Xe135',
+                'Xe136', 'Cs135', 'Gd157', 'Gd156']
+    test_xs = MicroXS.from_model(model, fuel, nuclides, chain_file=CHAIN_FILE)
     if config['update']:
         test_xs.to_csv('test_reference.csv')
 

tests/regression_tests/microxs/test_reference.csv

@@ -1,13 +1,13 @@
 nuclide,"(n,gamma)",fission
-U234,20.548033586079335,0.4951725071956495
-U235,10.593745111766133,48.86980740247932
-U238,0.8607296097035912,0.10623994948321437
-U236,8.697176401063281,0.32148140073986475
-O16,7.503456435273737e-05,0.0
-O17,0.0004107265933745623,0.0
-I135,6.896228129273278,0.0
-Xe135,229100.9245987756,0.0
-Xe136,0.02336047367105298,0.0
-Cs135,2.055822714073886,0.0
-Gd157,12927.465334134899,0.0
-Gd156,3.500756543915523,0.0
+U234,21.418670317831197,0.5014588470882195
+U235,10.343944102483244,47.46718472611891
+U238,0.8741166723597251,0.10829568455139126
+U236,9.083486784689326,0.3325287927011428
+O16,7.548646353912453e-05,0.0
+O17,0.0004018486221310307,0.0
+I135,6.6912565089429235,0.0
+Xe135,223998.64185667288,0.0
+Xe136,0.022934362666193576,0.0
+Cs135,2.28453952223533,0.0
+Gd157,12582.079620036275,0.0
+Gd156,2.9421127515332417,0.0