microxs from mg flux and chain file (openmc-dev#2755)

Co-authored-by: Jin Whan Bae <[email protected]>
Co-authored-by: shimwell <[email protected]>
Co-authored-by: Jonathan Shimwell <[email protected]>
Co-authored-by: Olek <[email protected]>
Co-authored-by: Paul Romano <[email protected]>

openmc/deplete/abc.py

@@ -71,8 +71,9 @@ def change_directory(output_dir):
         Directory to switch to.
     """
     orig_dir  = os.getcwd()
+    output_dir = Path(output_dir)
+    output_dir.mkdir(parents=True, exist_ok=True)
     try:
-        output_dir.mkdir(parents=True, exist_ok=True)
         os.chdir(output_dir)
         yield
     finally:

openmc/deplete/coupled_operator.py

@@ -10,6 +10,7 @@
 
 import copy
 from warnings import warn
+from typing import Optional
 
 import numpy as np
 from uncertainties import ufloat
@@ -33,18 +34,19 @@
 __all__ = ["CoupledOperator", "Operator", "OperatorResult"]
 
 
-def _find_cross_sections(model):
+def _find_cross_sections(model: Optional[str] = None):
     """Determine cross sections to use for depletion
 
     Parameters
     ----------
-    model : openmc.model.Model
+    model : openmc.model.Model, optional
         Reactor model
 
     """
-    if model.materials and model.materials.cross_sections is not None:
-        # Prefer info from Model class if available
-        return model.materials.cross_sections
+    if model:
+        if model.materials and model.materials.cross_sections is not None:
+            # Prefer info from Model class if available
+            return model.materials.cross_sections
 
     # otherwise fallback to environment variable
     cross_sections = openmc.config.get("cross_sections")
@@ -67,7 +69,7 @@ def _get_nuclides_with_data(cross_sections):
     Returns
     -------
     nuclides : set of str
-        Set of nuclide names that have cross secton data
+        Set of nuclide names that have cross section data
 
     """
     nuclides = set()

openmc/deplete/microxs.py

@@ -5,23 +5,39 @@
 """
 
 from __future__ import annotations
-import tempfile
-from typing import List, Tuple, Iterable, Optional, Union
+from tempfile import TemporaryDirectory
+from typing import List, Tuple, Iterable, Optional, Union, Sequence
 
 import pandas as pd
 import numpy as np
 
 from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike
 from openmc.exceptions import DataError
 from openmc import StatePoint
+from openmc.mgxs import GROUP_STRUCTURES
+from openmc.data import REACTION_MT
 import openmc
+from .abc import change_directory
 from .chain import Chain, REACTIONS
 from .coupled_operator import _find_cross_sections, _get_nuclides_with_data
+import openmc.lib
 
 _valid_rxns = list(REACTIONS)
 _valid_rxns.append('fission')
 
 
+def _resolve_chain_file_path(chain_file: str):
+    # Determine what reactions and nuclides are available in chain
+    if chain_file is None:
+        chain_file = openmc.config.get('chain_file')
+        if 'chain_file' in openmc.config:
+            raise DataError(
+                "No depletion chain specified and could not find depletion "
+                "chain in openmc.config['chain_file']"
+            )
+    return chain_file
+
+
 def get_microxs_and_flux(
         model: openmc.Model,
         domains,
@@ -69,13 +85,7 @@ def get_microxs_and_flux(
     original_tallies = model.tallies
 
     # 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_file = _resolve_chain_file_path(chain_file)
     chain = Chain.from_xml(chain_file)
     if reactions is None:
         reactions = chain.reactions
@@ -115,7 +125,7 @@ def get_microxs_and_flux(
     model.tallies = openmc.Tallies([rr_tally, flux_tally])
 
     # create temporary run
-    with tempfile.TemporaryDirectory() as temp_dir:
+    with TemporaryDirectory() as temp_dir:
         if run_kwargs is None:
             run_kwargs = {}
         else:
@@ -192,8 +202,122 @@ def __init__(self, data: np.ndarray, nuclides: List[str], reactions: List[str]):
         self._index_nuc = {nuc: i for i, nuc in enumerate(nuclides)}
         self._index_rx = {rx: i for i, rx in enumerate(reactions)}
 
-    # TODO: Add a classmethod for generating MicroXS directly from cross section
-    # data using a known flux spectrum
+    @classmethod
+    def from_multigroup_flux(
+        cls,
+        energies: Union[Sequence[float], str],
+        multigroup_flux: Sequence[float],
+        chain_file: Optional[PathLike] = None,
+        temperature: float = 293.6,
+        nuclides: Optional[Sequence[str]] = None,
+        reactions: Optional[Sequence[str]] = None,
+        **init_kwargs: dict,
+    ) -> MicroXS:
+        """Generated microscopic cross sections from a known flux.
+
+        The size of the MicroXS matrix depends on the chain file and cross
+        sections available. MicroXS entry will be 0 if the nuclide cross section
+        is not found.
+
+        .. versionadded:: 0.14.1
+
+        Parameters
+        ----------
+        energies : iterable of float or str
+            Energy group boundaries in [eV] or the name of the group structure
+        multi_group_flux : iterable of float
+            Energy-dependent multigroup flux values
+        chain_file : str, optional
+            Path to the depletion chain XML file that will be used in depletion
+            simulation.  Defaults to ``openmc.config['chain_file']``.
+        temperature : int, optional
+            Temperature for cross section evaluation in [K].
+        nuclides : list of str, optional
+            Nuclides to get cross sections for. If not specified, all burnable
+            nuclides from the depletion chain file are used.
+        reactions : list of str, optional
+            Reactions to get cross sections for. If not specified, all neutron
+            reactions listed in the depletion chain file are used.
+        **init_kwargs : dict
+            Keyword arguments passed to :func:`openmc.lib.init`
+
+        Returns
+        -------
+        MicroXS
+        """
+
+        check_type("temperature", temperature, (int, float))
+        # if energy is string then use group structure of that name
+        if isinstance(energies, str):
+            energies = GROUP_STRUCTURES[energies]
+        else:
+            # if user inputs energies check they are ascending (low to high) as
+            # some depletion codes use high energy to low energy.
+            if not np.all(np.diff(energies) > 0):
+                raise ValueError('Energy group boundaries must be in ascending order')
+
+        # check dimension consistency
+        if len(multigroup_flux) != len(energies) - 1:
+            raise ValueError('Length of flux array should be len(energies)-1')
+
+        chain_file_path = _resolve_chain_file_path(chain_file)
+        chain = Chain.from_xml(chain_file_path)
+
+        cross_sections = _find_cross_sections(model=None)
+        nuclides_with_data = _get_nuclides_with_data(cross_sections)
+
+        # If no nuclides were specified, default to all nuclides from the chain
+        if not nuclides:
+            nuclides = chain.nuclides
+            nuclides = [nuc.name for nuc in nuclides]
+
+        # Get reaction MT values. If no reactions specified, default to the
+        # reactions available in the chain file
+        if reactions is None:
+            reactions = chain.reactions
+        mts = [REACTION_MT[name] for name in reactions]
+
+        # Normalize multigroup flux
+        multigroup_flux = np.asarray(multigroup_flux)
+        multigroup_flux /= multigroup_flux.sum()
+
+        # Create 2D array for microscopic cross sections
+        microxs_arr = np.zeros((len(nuclides), len(mts)))
+
+        with TemporaryDirectory() as tmpdir:
+            # Create a material with all nuclides
+            mat_all_nucs = openmc.Material()
+            for nuc in nuclides:
+                if nuc in nuclides_with_data:
+                    mat_all_nucs.add_nuclide(nuc, 1.0)
+            mat_all_nucs.set_density("atom/b-cm", 1.0)
+
+            # Create simple model containing the above material
+            surf1 = openmc.Sphere(boundary_type="vacuum")
+            surf1_cell = openmc.Cell(fill=mat_all_nucs, region=-surf1)
+            model = openmc.Model()
+            model.geometry = openmc.Geometry([surf1_cell])
+            model.settings = openmc.Settings(
+                particles=1, batches=1, output={'summary': False})
+
+            with change_directory(tmpdir):
+                # Export model within temporary directory
+                model.export_to_model_xml()
+
+                with openmc.lib.run_in_memory(**init_kwargs):
+                    # For each nuclide and reaction, compute the flux-averaged
+                    # cross section
+                    for nuc_index, nuc in enumerate(nuclides):
+                        if nuc not in nuclides_with_data:
+                            continue
+                        lib_nuc = openmc.lib.nuclides[nuc]
+                        for mt_index, mt in enumerate(mts):
+                            xs = lib_nuc.collapse_rate(
+                                mt, temperature, energies, multigroup_flux
+                            )
+                            microxs_arr[nuc_index, mt_index] = xs
+
+        return cls(microxs_arr, nuclides, reactions)
 
     @classmethod
     def from_csv(cls, csv_file, **kwargs):

openmc/lib/nuclide.py

@@ -91,7 +91,7 @@ def collapse_rate(self, MT, temperature, energy, flux):
         energy : iterable of float
             Energy group boundaries in [eV]
         flux : iterable of float
-            Flux in each energt group (not normalized per eV)
+            Flux in each energy group (not normalized per eV)
 
         Returns
         -------

tests/unit_tests/test_deplete_microxs.py

@@ -51,10 +51,61 @@ def test_from_array():
                        r'match dimensions of data array of shape \(\d*\, \d*\)'):
         MicroXS(data[:, 0], nuclides, reactions)
 
+
 def test_csv():
     ref_xs = MicroXS.from_csv(ONE_GROUP_XS)
     ref_xs.to_csv('temp_xs.csv')
     temp_xs = MicroXS.from_csv('temp_xs.csv')
     assert np.all(ref_xs.data == temp_xs.data)
     remove('temp_xs.csv')
 
+
+def test_from_multigroup_flux():
+    energies = [0., 6.25e-1, 5.53e3, 8.21e5, 2.e7]
+    flux = [1.1e-7, 1.2e-6, 1.3e-5, 1.4e-4]
+    chain_file = Path(__file__).parents[1] / 'chain_simple.xml'
+    kwargs = {'multigroup_flux': flux, 'chain_file': chain_file}
+
+    # test with energy group structure from string
+    microxs = MicroXS.from_multigroup_flux(energies='CASMO-4', **kwargs)
+    assert isinstance(microxs, MicroXS)
+
+    # test with energy group structure as floats
+    microxs = MicroXS.from_multigroup_flux(energies=energies, **kwargs)
+    assert isinstance(microxs, MicroXS)
+
+    # test with nuclides provided
+    microxs = MicroXS.from_multigroup_flux(
+        energies=energies, nuclides=['Gd157', 'H1'], **kwargs
+    )
+    assert isinstance(microxs, MicroXS)
+    assert microxs.nuclides == ['Gd157', 'H1']
+
+    # test with reactions provided
+    microxs = MicroXS.from_multigroup_flux(
+        energies=energies, reactions=['fission', '(n,2n)'], **kwargs
+    )
+    assert isinstance(microxs, MicroXS)
+    assert microxs.reactions == ['fission', '(n,2n)']
+
+
+def test_multigroup_flux_same():
+    chain_file = Path(__file__).parents[1] / 'chain_simple.xml'
+
+    # Generate micro XS based on 4-group flux
+    energies = [0., 6.25e-1, 5.53e3, 8.21e5, 2.e7]
+    flux_per_ev = [0.3, 0.3, 1.0, 1.0]
+    flux = flux_per_ev * np.diff(energies)
+    microxs_4g = MicroXS.from_multigroup_flux(
+        energies=energies, multigroup_flux=flux, chain_file=chain_file)
+
+    # Generate micro XS based on 2-group flux, where the boundaries line up with
+    # the 4 group flux and have the same flux per eV across the full energy
+    # range
+    energies = [0., 5.53e3, 2.0e7]
+    flux_per_ev = [0.3, 1.0]
+    flux = flux_per_ev * np.diff(energies)
+    microxs_2g = MicroXS.from_multigroup_flux(
+        energies=energies, multigroup_flux=flux, chain_file=chain_file)
+
+    assert microxs_4g.data == pytest.approx(microxs_2g.data)