microxs from mg flux and chain file

openmc/deplete/microxs.py

@@ -6,7 +6,7 @@
 
 from __future__ import annotations
 import tempfile
-from typing import List, Tuple, Iterable, Optional, Union
+from typing import List, Tuple, Iterable, Optional, Union, Sequence
 
 import pandas as pd
 import numpy as np
@@ -22,6 +22,29 @@
 _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:
+        if 'chain_file' in openmc.config:
+            raise DataError(
+                "No depletion chain specified and could not find depletion "
+                "chain in openmc.config['chain_file']"
+            )
+        else:
+            chain_file = openmc.config['chain_file']
+    return chain_file
+
+
+def _resolve_openmc_data_path(openmc_data_path:str=None):
+    if not openmc_data_path:
+        if 'cross_sections' not in openmc.config:
+            raise ValueError("`openmc_data_path` is not defined nor `openmc.config['cross_sections']` is defined.")
+        else:
+            openmc_data_path = openmc.config['cross_sections']
+    return openmc_data_path
+
+
+
 def get_microxs_and_flux(
         model: openmc.Model,
         domains,
@@ -69,13 +92,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
@@ -192,8 +209,93 @@ 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_multi_group_flux(
+        cls,
+        energies: Union[Iterable[float], str],
+        multi_group_flux: Sequence[float],
+        chain_file: Optional[PathLike] = None,
+        openmc_data_path: Optional[PathLike] = None,
+        temperature:int=294
+    ):
+        """Generated MicroXS object from a known flux and a chain file. The size of the MicroXs matrix depends
+        on the chain file.
+
+        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.
+        openmc_dat_path: str, optional
+            Path to the cross section XML file that contains data library paths.
+        temperature: int, optional
+            Temperature for cross section evaluation in [K].
+
+        Returns
+        -------
+        MicroXS
+        """
+        # check energies
+        if isinstance(energies, str):
+            energies = openmc.EnergyFilter.from_group_structure(energies).values
+        else:
+            # if user inputs energies check they are ascending (low to high) as some
+            # depletion codes use high energy to low energy.
+            if not all(energies[i] <= energies[i+1] for i in range(len(energies) - 1)):
+                raise ValueError('Energy bin must be in ascending order')
+
+        # check dimension consistency
+        if not len(multi_group_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 = openmc.deplete.Chain.from_xml(chain_file_path)
+
+        openmc_data_path = _resolve_openmc_data_path(openmc_data_path)
+        data_lib = openmc.data.DataLibrary.from_xml(openmc_data_path)
+
+        # get reactions and nuclides from chain file
+        nuclides, reactions = chain.nuclides, chain.reactions
+        nuclides = [nuc.name for nuc in nuclides]
+        if 'fission' in reactions:
+            # send to back
+            reactions.append(reactions.pop(reactions.index('fission')))
+        # convert reactions to mt file
+        #! I feel like this zero indexing will backfire
+        #! There has to be a better way for this
+        mts = [list(REACTIONS[reaction].mts)[0] for reaction in reactions if reaction != 'fission']
+        if 'fission' in reactions:
+            mts.append(18) # fission is not in the REACTIONS
+
+        # normalize flux and get energy midpoint
+        norm_multi_group_flux = np.array(multi_group_flux) / sum(multi_group_flux)
+        energies_midpoint = [(energies[i]+energies[i+1])/2 for i in range(len(energies)-1)]
+        temperature_key = f'{temperature}K'
+        microxs_arr = np.zeros((len(nuclides), len(mts)))
+
+        for nuc_indx, nuc in enumerate(nuclides):
+            mat_lib = data_lib.get_by_material(nuc, data_type='neutron')
+            if not mat_lib: # file does not exist
+                microxs_arr[nuc_indx, :] = 0
+                continue
+            else:
+                hdf5_path = mat_lib['path']
+            nuc_data = openmc.data.IncidentNeutron.from_hdf5(hdf5_path)
+            for mt_indx, mt in enumerate(mts):
+                try: # sometimes it fails cause there's no entry
+                    # perhaps a multigroup mgxs should be made as chosing the xs at the energy bin center is not always fair
+                    total_xs = np.array(nuc_data[mt].xs[temperature_key](energies_midpoint))
+                    # collapse
+                    total_norm = sum(total_xs * norm_multi_group_flux)
+                    microxs_arr[nuc_indx, mt_indx] = total_norm
+                except KeyError:
+                    microxs_arr[nuc_indx, mt_indx] = 0.0
+
+        return cls(nuclides=nuclides, reactions=reactions, data=microxs_arr)
 
     @classmethod
     def from_csv(cls, csv_file, **kwargs):

tests/unit_tests/test_deplete_microxs.py

@@ -58,3 +58,20 @@ def test_csv():
     assert np.all(ref_xs.data == temp_xs.data)
     remove('temp_xs.csv')
 
+def test_from_multi_group_flux():
+
+    microxs = MicroXS.from_multi_group_flux(
+        energies='CASMO-4',
+        multi_group_flux=[1.1e-7, 1.2e-6, 1.3e-5, 1.4e-4],
+        temperature='293',
+        chain_file=Path(__file__).parents[1] / 'chain_simple.xml'
+    )
+    assert isinstance(microxs, MicroXS)
+
+    microxs = MicroXS.from_multi_group_flux(
+        energies=[0., 6.25e-1, 5.53e3, 8.21e5, 2.e7],
+        multi_group_flux=[1.1e-7, 1.2e-6, 1.3e-5, 1.4e-4],
+        temperature='293',
+        chain_file=Path(__file__).parents[1] / 'chain_simple.xml'
+    )
+    assert isinstance(microxs, MicroXS)