diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index dd4df21fb5d..60774871c01 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -68,8 +68,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: diff --git a/openmc/deplete/coupled_operator.py b/openmc/deplete/coupled_operator.py index 0d778c7450c..287a4e0d37b 100644 --- a/openmc/deplete/coupled_operator.py +++ b/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() diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index e145bce8763..0721535ca4f 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -5,8 +5,8 @@ """ 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 @@ -14,14 +14,30 @@ 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): diff --git a/openmc/lib/nuclide.py b/openmc/lib/nuclide.py index 399bb346520..8078882cf35 100644 --- a/openmc/lib/nuclide.py +++ b/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 ------- diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index 582c584654d..ad54026f014 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -51,6 +51,7 @@ 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') @@ -58,3 +59,53 @@ def test_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)