Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

microxs from mg flux and chain file #2755

Merged
merged 19 commits into from
Dec 22, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion openmc/deplete/abc.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
14 changes: 8 additions & 6 deletions openmc/deplete/coupled_operator.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

import copy
from warnings import warn
from typing import Optional

import numpy as np
from uncertainties import ufloat
Expand All @@ -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")
Expand All @@ -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()
Expand Down
148 changes: 136 additions & 12 deletions openmc/deplete/microxs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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`

shimwell marked this conversation as resolved.
Show resolved Hide resolved
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(
paulromano marked this conversation as resolved.
Show resolved Hide resolved
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):
Expand Down
2 changes: 1 addition & 1 deletion openmc/lib/nuclide.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
-------
Expand Down
51 changes: 51 additions & 0 deletions tests/unit_tests/test_deplete_microxs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)