crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'refactor-deplete-operators' into transportoperator-subclass-notransport

Browse source
This commit is contained in:
yardasol 2022-07-19 10:08:20 -05:00
commit d69150db0c
7 changed files with 855 additions and 487 deletions
Download patch file Download diff file Expand all files Collapse all files

30
docs/source/pythonapi/deplete.rst
View file

@ -190,6 +190,35 @@ total system energy.
helpers.FissionYieldCutoffHelper
helpers.FluxCollapseHelper
Intermediate Classes
--------------------
Specific implementations of abstract base classes may utilize some of
the same methods and data structures. These methods and data are stored
in intermediate classes.
Methods common to tally-based implementation of :class:`FissionYieldHelper`
are stored in :class:`helpers.TalliedFissionYieldHelper`
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
helpers.TalliedFissionYieldHelper
Methods common to OpenMC-specific implementations of :class:`TransportOperator`
are stored in :class:`openmc_operator.OpenMCOperator`
.. autosummary::
:toctree: generated
:nosignatures:
:template: mycallable.rst
openmc_operator.OpenMCOperator
Abstract Base Classes
---------------------
@ -217,7 +246,6 @@ OpenMC simulations back on to the :class:`abc.TransportOperator`
abc.NormalizationHelper
abc.FissionYieldHelper
abc.ReactionRateHelper
abc.TalliedFissionYieldHelper
Custom integrators or depletion solvers can be developed by subclassing from
the following abstract base classes:

1
openmc/deplete/__init__.py
View file

@ -7,6 +7,7 @@ A depletion front-end tool.
from .nuclide import *
from .chain import *
from .openmc_operator import *
from .operator import *
from .flux_operator import *
from .reaction_rates import *

46
openmc/deplete/abc.py
View file

@ -10,6 +10,7 @@ from collections.abc import Iterable, Callable
from copy import deepcopy
from inspect import signature
from numbers import Real, Integral
from contextlib import contextmanager
import os
from pathlib import Path
import sys
@ -56,6 +57,24 @@ except AttributeError:
# Can't set __doc__ on properties on Python 3.4
pass
@contextmanager
def change_directory(output_dir):
"""
Helper function for managing the current directory.
Parameters
----------
output_dir : pathlib.Path
Directory to switch to.
"""
orig_dir = os.getcwd()
try:
output_dir.mkdir(exist_ok=True)
os.chdir(output_dir)
yield
finally:
os.chdir(orig_dir)
class TransportOperator(ABC):
"""Abstract class defining a transport operator
@ -87,9 +106,14 @@ class TransportOperator(ABC):
Initial atom density [atoms/cm^3] to add for nuclides that are zero
in initial condition to ensure they exist in the decay chain.
Only done for nuclides with reaction rates.
output_dir : pathlib.Path
Path to output directory to save results.
prev_res : Results or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
chain : openmc.deplete.Chain
The depletion chain information necessary to form matrices and tallies.
"""
def __init__(self, chain_file, fission_q=None, dilute_initial=1.0e3,
prev_results=None):
@ -133,18 +157,6 @@ class TransportOperator(ABC):
"""
def __enter__(self):
# Save current directory and move to specific output directory
self._orig_dir = os.getcwd()
self.output_dir.mkdir(exist_ok=True)
os.chdir(self.output_dir)
return self.initial_condition()
def __exit__(self, exc_type, exc_value, traceback):
self.finalize()
os.chdir(self._orig_dir)
@property
def output_dir(self):
return self._output_dir
@ -775,7 +787,8 @@ class Integrator(ABC):
.. versionadded:: 0.13.1
"""
with self.operator as conc:
with change_directory(self.operator.output_dir):
conc = self.operator.initial_condition()
t, self._i_res = self._get_start_data()
for i, (dt, source_rate) in enumerate(self):
@ -814,6 +827,8 @@ class Integrator(ABC):
source_rate, self._i_res + len(self), proc_time)
self.operator.write_bos_data(len(self) + self._i_res)
self.operator.finalize()
@add_params
class SIIntegrator(Integrator):
@ -931,7 +946,8 @@ class SIIntegrator(Integrator):
.. versionadded:: 0.13.1
"""
with self.operator as conc:
with change_directory(self.operator.output_dir):
conc = self.operator.initial_condition()
t, self._i_res = self._get_start_data()
for i, (dt, p) in enumerate(self):
@ -967,6 +983,8 @@ class SIIntegrator(Integrator):
p, self._i_res + len(self), proc_time)
self.operator.write_bos_data(self._i_res + len(self))
self.operator.finalize()
class DepSystemSolver(ABC):
r"""Abstract class for solving depletion equations

593
openmc/deplete/openmc_operator.py Normal file
View file

@ -0,0 +1,593 @@
"""OpenMC transport operator
This module implements functions used by both OpenMC transport operators as well as pure depletion operators.
"""
from abc import abstractmethod
from collections import OrderedDict
import numpy as np
import openmc
from openmc.mpi import comm
from .abc import TransportOperator, OperatorResult
from .atom_number import AtomNumber
from .reaction_rates import ReactionRates
__all__ = ["OpenMCOperator", "OperatorResult"]
def _distribute(items):
"""Distribute items across MPI communicator
Parameters
----------
items : list
List of items of distribute
Returns
-------
list
Items assigned to process that called
"""
min_size, extra = divmod(len(items), comm.size)
j = 0
for i in range(comm.size):
chunk_size = min_size + int(i < extra)
if comm.rank == i:
return items[j:j + chunk_size]
j += chunk_size
class OpenMCOperator(TransportOperator):
"""Abstract class holding OpenMC-specific functions for running
depletion calculations.
Specific classes for running coupled transport-depletion calculations or
depletion-only calculations are implemented as subclasses of OpenMCOperator
Parameters
----------
materials : openmc.Materials
List of all materials in the model
cross_sections : str or pandas.DataFrame
Path to continuous energy cross section library, or object containing
one-group cross-sections.
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the file
listed under ``depletion_chain`` in
:envvar:`OPENMC_CROSS_SECTIONS` environment variable.
prev_results : Results, optional
Results from a previous depletion calculation. If this argument is
specified, the depletion calculation will start from the latest state
in the previous results.
diff_burnable_mats : bool, optional
Whether to differentiate burnable materials with multiple instances.
Volumes are divided equally from the original material volume.
Default: False.
normalization_mode : str
Indicate how reaction rates should be normalized.
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV].
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 decay chain.
Only done for nuclides with reaction rates.
Defaults to 1.0e3.
fission_yield_mode : str
Key indicating what fission product yield scheme to use.
fission_yield_opts : dict of str to option, optional
Optional arguments to pass to the helper determined by
``fission_yield_mode``. Will be passed directly on to the
helper. Passing a value of None will use the defaults for
the associated helper.
reaction_rate_mode : str, optional
Indicate how one-group reaction rates should be calculated.
reaction_rate_opts : dict, optional
Keyword arguments that are passed to the reaction rate helper class.
reduce_chain : bool, optional
If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
depletion chain up to ``reduce_chain_level``. Default is False.
reduce_chain_level : int, optional
Depth of the search when reducing the depletion chain. Only used
if ``reduce_chain`` evaluates to true. The default value of
``None`` implies no limit on the depth.
Attributes
----------
materials : openmc.Materials
All materials present in the model
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 decay
chain. Only done for nuclides with reaction rates.
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
Whether or not to round output to OpenMC to 8 digits.
Useful in testing, as OpenMC is incredibly sensitive to exact values.
number : openmc.deplete.AtomNumber
Total number of atoms in simulation.
nuclides_with_data : set of str
A set listing all unique nuclides available from cross_sections.xml.
chain : openmc.deplete.Chain
The depletion chain information necessary to form matrices and tallies.
reaction_rates : openmc.deplete.ReactionRates
Reaction rates from the last operator step.
burnable_mats : list of str
All burnable material IDs
heavy_metal : float
Initial heavy metal inventory [g]
local_mats : list of str
All burnable material IDs being managed by a single process
prev_res : Results or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
"""
def __init__(
self,
materials=None,
cross_sections=None,
chain_file=None,
prev_results=None,
diff_burnable_mats=False,
normalization_mode=None,
fission_q=None,
dilute_initial=0.0,
fission_yield_mode=None,
fission_yield_opts=None,
reaction_rate_mode=None,
reaction_rate_opts=None,
reduce_chain=False,
reduce_chain_level=None):
super().__init__(chain_file, fission_q, dilute_initial, prev_results)
self.round_number = False
self.materials = materials
self.cross_sections = cross_sections
# Reduce the chain to only those nuclides present
if reduce_chain:
init_nuclides = set()
for material in self.materials:
if not material.depletable:
continue
for name, _dens_percent, _dens_type in material.nuclides:
init_nuclides.add(name)
self.chain = self.chain.reduce(init_nuclides, reduce_chain_level)
if diff_burnable_mats:
self._differentiate_burnable_mats()
# Determine which nuclides have cross section data
# This nuclides variables contains every nuclides
# for which there is an entry in the micro_xs parameter
openmc.reset_auto_ids()
self.burnable_mats, volumes, all_nuclides = self._get_burnable_mats()
self.local_mats = _distribute(self.burnable_mats)
self._mat_index_map = {
lm: self.burnable_mats.index(lm) for lm in self.local_mats}
if self.prev_res is not None:
self._load_previous_results()
self.nuclides_with_data = self._get_nuclides_with_data(
self.cross_sections)
# Select nuclides with data that are also in the chain
self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
if nuc.name in self.nuclides_with_data]
# Extract number densities from the geometry / previous depletion run
self._extract_number(self.local_mats,
volumes,
all_nuclides,
self.prev_res)
# Create reaction rates array
self.reaction_rates = ReactionRates(
self.local_mats, self._burnable_nucs, self.chain.reactions)
self._get_helper_classes(
reaction_rate_mode,
normalization_mode,
fission_yield_mode,
reaction_rate_opts,
fission_yield_opts)
@abstractmethod
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material"""
def _get_burnable_mats(self):
"""Determine depletable materials, volumes, and nuclides
Returns
-------
burnable_mats : list of str
List of burnable material IDs
volume : OrderedDict of str to float
Volume of each material in [cm^3]
nuclides : list of str
Nuclides in order of how they'll appear in the simulation.
"""
burnable_mats = set()
model_nuclides = set()
volume = OrderedDict()
self.heavy_metal = 0.0
# Iterate once through the geometry to get dictionaries
for mat in self.materials:
for nuclide in mat.get_nuclides():
model_nuclides.add(nuclide)
if mat.depletable:
burnable_mats.add(str(mat.id))
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
volume[str(mat.id)] = mat.volume
self.heavy_metal += mat.fissionable_mass
# Make sure there are burnable materials
if not burnable_mats:
raise RuntimeError(
"No depletable materials were found in the model.")
# Sort the sets
burnable_mats = sorted(burnable_mats, key=int)
model_nuclides = sorted(model_nuclides)
# Construct a global nuclide dictionary, burned first
nuclides = list(self.chain.nuclide_dict)
for nuc in model_nuclides:
if nuc not in nuclides:
nuclides.append(nuc)
return burnable_mats, volume, nuclides
@abstractmethod
def _load_previous_results(self):
"""Load results from a previous depletion simulation"""
@abstractmethod
def _get_nuclides_with_data(self, cross_sections):
"""Find nuclides with cross section data
Parameters
----------
cross_sections : str or pandas.DataFrame
Path to continuous energy cross section library, or object
containing one-group cross-sections.
Returns
-------
nuclides : set of str
Set of nuclide names that have cross secton data
"""
def _extract_number(self, local_mats, volume, all_nuclides, prev_res=None):
"""Construct AtomNumber using geometry
Parameters
----------
local_mats : list of str
Material IDs to be managed by this process
volume : OrderedDict of str to float
Volumes for the above materials in [cm^3]
all_nuclides : list of str
Nuclides to be used in the simulation.
prev_res : Results, optional
Results from a previous depletion calculation
"""
self.number = AtomNumber(
local_mats, all_nuclides, volume, len(
self.chain))
if self.dilute_initial != 0.0:
for nuc in self._burnable_nucs:
self.number.set_atom_density(
np.s_[:], nuc, self.dilute_initial)
# Now extract and store the number densities
# From the geometry if no previous depletion results
if prev_res is None:
for mat in self.materials:
if str(mat.id) in local_mats:
self._set_number_from_mat(mat)
# Else from previous depletion results
else:
for mat in self.materials:
if str(mat.id) in local_mats:
self._set_number_from_results(mat, prev_res)
def _set_number_from_mat(self, mat):
"""Extracts material and number densities from openmc.Material
Parameters
----------
mat : openmc.Material
The material to read from
"""
mat_id = str(mat.id)
for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items():
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
def _set_number_from_results(self, mat, prev_res):
"""Extracts material nuclides and number densities.
If the nuclide concentration's evolution is tracked, the densities come
from depletion results. Else, densities are extracted from the geometry
in the summary.
Parameters
----------
mat : openmc.Material
The material to read from
prev_res : Results
Results from a previous depletion calculation
"""
mat_id = str(mat.id)
# Get nuclide lists from geometry and depletion results
depl_nuc = prev_res[-1].nuc_to_ind
geom_nuc_densities = mat.get_nuclide_atom_densities()
# Merge lists of nuclides, with the same order for every calculation
geom_nuc_densities.update(depl_nuc)
for nuclide, atom_per_bcm in geom_nuc_densities.items():
if nuclide in depl_nuc:
concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1]
volume = prev_res[-1].volume[mat_id]
atom_per_cc = concentration / volume
else:
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
@abstractmethod
def _get_helper_classes(
self,
reaction_rate_mode,
normalization_mode,
fission_yield_mode,
reaction_rate_opts,
fission_yield_opts):
"""Create the ``_rate_helper``, ``_normalization_helper``, and
``_yield_helper`` objects.
Parameters
----------
reaction_rate_mode : str
Indicates the subclass of :class:`ReactionRateHelper` to
instantiate.
normalization_mode : str
Indicates the subclass of :class:`NormalizationHelper` to
instantiate.
fission_yield_mode : str
Indicates the subclass of :class:`FissionYieldHelper` to instatiate.
reaction_rate_opts : dict
Keyword arguments that are passed to the :class:`ReactionRateHelper`
subclass.
fission_yield_opts : dict
Keyword arguments that are passed to the :class:`FissionYieldHelper`
subclass.
"""
def initial_condition(self, materials):
"""Performs final setup and returns initial condition.
Parameters
----------
materials : list of str
list of material IDs
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
self._rate_helper.generate_tallies(materials, self.chain.reactions)
self._normalization_helper.prepare(
self.chain.nuclides, self.reaction_rates.index_nuc)
# Tell fission yield helper what materials this process is
# responsible for
self._yield_helper.generate_tallies(
materials, tuple(sorted(self._mat_index_map.values())))
# Return number density vector
return list(self.number.get_mat_slice(np.s_[:]))
def _update_materials_and_nuclides(self, vec):
"""Update the number density, material compositions, and nuclide
lists in helper objects
Parameters
----------
vec : list of numpy.ndarray
Total atoms.
"""
# Update the number densities regardless of the source rate
self.number.set_density(vec)
self._update_materials()
# Prevent OpenMC from complaining about re-creating tallies
openmc.reset_auto_ids()
# Update tally nuclides data in preparation for transport solve
nuclides = self._get_reaction_nuclides()
self._rate_helper.nuclides = nuclides
self._normalization_helper.nuclides = nuclides
self._yield_helper.update_tally_nuclides(nuclides)
@abstractmethod
def _update_materials(self):
"""Updates material compositions in OpenMC on all processes."""
def _get_reaction_nuclides(self):
"""Determine nuclides that should be tallied for reaction rates.
This method returns a list of all nuclides that have cross section data
and are listed in the depletion chain. Technically, we should count
nuclides that may not appear in the depletion chain because we still
need to get the fission reaction rate for these nuclides in order to
normalize power, but that is left as a future exercise.
Returns
-------
list of str
Nuclides with reaction rates
"""
nuc_set = set()
# Create the set of all nuclides in the decay chain in materials marked
# for burning in which the number density is greater than zero.
for nuc in self.number.nuclides:
if nuc in self.nuclides_with_data:
if np.sum(self.number[:, nuc]) > 0.0:
nuc_set.add(nuc)
# Communicate which nuclides have nonzeros to rank 0
if comm.rank == 0:
for i in range(1, comm.size):
nuc_newset = comm.recv(source=i, tag=i)
nuc_set |= nuc_newset
else:
comm.send(nuc_set, dest=0, tag=comm.rank)
if comm.rank == 0:
# Sort nuclides in the same order as self.number
nuc_list = [nuc for nuc in self.number.nuclides
if nuc in nuc_set]
else:
nuc_list = None
# Store list of nuclides on each process
nuc_list = comm.bcast(nuc_list)
return [nuc for nuc in nuc_list if nuc in self.chain]
def _calculate_reaction_rates(self, source_rate):
"""Unpack tallies from OpenMC and return an operator result
This method uses OpenMC's C API bindings to determine the k-effective
value and reaction rates from the simulation. The reaction rates are
normalized by a helper class depending on the method being used.
Parameters
----------
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
rates : openmc.deplete.ReactionRates
Reaction rates for nuclides
"""
rates = self.reaction_rates
rates.fill(0.0)
# Extract reaction nuclides
rxn_nuclides = self._rate_helper.nuclides
# Form fast map
nuc_ind = [rates.index_nuc[nuc] for nuc in rxn_nuclides]
react_ind = [rates.index_rx[react] for react in self.chain.reactions]
# Keep track of energy produced from all reactions in eV per source
# particle
self._normalization_helper.reset()
self._yield_helper.unpack()
# Store fission yield dictionaries
fission_yields = []
# Create arrays to store fission Q values, reaction rates, and nuclide
# numbers, zeroed out in material iteration
number = np.empty(rates.n_nuc)
fission_ind = rates.index_rx.get("fission")
# Extract results
for i, mat in enumerate(self.local_mats):
# Get tally index
mat_index = self._mat_index_map[mat]
# Zero out reaction rates and nuclide numbers
number.fill(0.0)
# Get new number densities
for nuc, i_nuc_results in zip(rxn_nuclides, nuc_ind):
number[i_nuc_results] = self.number[mat, nuc]
tally_rates = self._rate_helper.get_material_rates(
mat_index, nuc_ind, react_ind)
# Compute fission yields for this material
fission_yields.append(self._yield_helper.weighted_yields(i))
# Accumulate energy from fission
if fission_ind is not None:
self._normalization_helper.update(tally_rates[:, fission_ind])
# Divide by total number and store
rates[i] = self._rate_helper.divide_by_adens(number)
# Scale reaction rates to obtain units of reactions/sec
rates *= self._normalization_helper.factor(source_rate)
# Store new fission yields on the chain
self.chain.fission_yields = fission_yields
return rates
def get_results_info(self):
"""Returns volume list, material lists, and nuc lists.
Returns
-------
volume : dict of str float
Volumes corresponding to materials in full_burn_dict
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all material IDs to be burned. Used for sorting the simulation.
full_burn_list : list
List of all burnable material IDs
"""
nuc_list = self.number.burnable_nuclides
burn_list = self.local_mats
volume = {}
for i, mat in enumerate(burn_list):
volume[mat] = self.number.volume[i]
# Combine volume dictionaries across processes
volume_list = comm.allgather(volume)
volume = {k: v for d in volume_list for k, v in d.items()}
return volume, nuc_list, burn_list, self.burnable_mats

655
openmc/deplete/operator.py
View file

@ -8,7 +8,6 @@ densities is all done in-memory instead of through the filesystem.
"""
import copy
from collections import OrderedDict
import os
from warnings import warn
@ -21,10 +20,9 @@ from openmc.data import DataLibrary
from openmc.exceptions import DataError
import openmc.lib
from openmc.mpi import comm
from .abc import TransportOperator, OperatorResult
from .atom_number import AtomNumber
from .abc import OperatorResult
from .chain import _find_chain_file
from .reaction_rates import ReactionRates
from .openmc_operator import OpenMCOperator, _distribute
from .results import Results
from .helpers import (
DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper,
@ -35,31 +33,15 @@ from .helpers import (
__all__ = ["Operator", "OperatorResult"]
def _distribute(items):
"""Distribute items across MPI communicator
def _find_cross_sections(model):
"""Determine cross sections to use for depletion
Parameters
----------
items : list
List of items of distribute
Returns
-------
list
Items assigned to process that called
model : openmc.model.Model
Reactor model
"""
min_size, extra = divmod(len(items), comm.size)
j = 0
for i in range(comm.size):
chunk_size = min_size + int(i < extra)
if comm.rank == i:
return items[j:j + chunk_size]
j += chunk_size
def _find_cross_sections(model):
"""Determine cross sections to use for depletion"""
if model.materials and model.materials.cross_sections is not None:
# Prefer info from Model class if available
return model.materials.cross_sections
@ -74,7 +56,7 @@ def _find_cross_sections(model):
return cross_sections
class Operator(TransportOperator):
class Operator(OpenMCOperator):
"""OpenMC transport operator for depletion.
Instances of this class can be used to perform depletion using OpenMC as the
@ -233,8 +215,6 @@ class Operator(TransportOperator):
if fission_q is not None:
warn("Fission Q dictionary will not be used")
fission_q = None
super().__init__(chain_file, fission_q, dilute_initial, prev_results)
self.round_number = False
self.model = model
self.settings = model.settings
self.geometry = model.geometry
@ -244,66 +224,123 @@ class Operator(TransportOperator):
model.materials = openmc.Materials(
model.geometry.get_all_materials().values()
)
self.materials = model.materials
self.cleanup_when_done = True
# Reduce the chain before we create more materials
if reduce_chain:
all_isotopes = set()
for material in self.materials:
if not material.depletable:
continue
for name, _dens_percent, _dens_type in material.nuclides:
all_isotopes.add(name)
self.chain = self.chain.reduce(all_isotopes, reduce_chain_level)
super().__init__(
model.materials,
cross_sections,
chain_file,
prev_results,
diff_burnable_mats,
normalization_mode,
fission_q,
dilute_initial,
fission_yield_mode,
fission_yield_opts,
reaction_rate_mode,
reaction_rate_opts,
reduce_chain,
reduce_chain_level)
# Differentiate burnable materials with multiple instances
if diff_burnable_mats:
self._differentiate_burnable_mats()
self.materials = openmc.Materials(
model.geometry.get_all_materials().values()
)
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material"""
# Clear out OpenMC, create task lists, distribute
openmc.reset_auto_ids()
self.burnable_mats, volume, nuclides = self._get_burnable_mats()
self.local_mats = _distribute(self.burnable_mats)
# Count the number of instances for each cell and material
self.geometry.determine_paths(instances_only=True)
# Generate map from local materials => material index
self._mat_index_map = {
lm: self.burnable_mats.index(lm) for lm in self.local_mats}
# Extract all burnable materials which have multiple instances
distribmats = set(
[mat for mat in self.materials
if mat.depletable and mat.num_instances > 1])
if self.prev_res is not None:
# Reload volumes into geometry
prev_results[-1].transfer_volumes(self.model)
for mat in distribmats:
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
mat.volume /= mat.num_instances
# Store previous results in operator
# Distribute reaction rates according to those tracked
# on this process
if comm.size == 1:
self.prev_res = prev_results
else:
self.prev_res = Results()
mat_indexes = _distribute(range(len(self.burnable_mats)))
for res_obj in prev_results:
new_res = res_obj.distribute(self.local_mats, mat_indexes)
self.prev_res.append(new_res)
if distribmats:
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats:
mat = cell.fill
cell.fill = [mat.clone()
for i in range(cell.num_instances)]
# Determine which nuclides have incident neutron data
self.nuclides_with_data = self._get_nuclides_with_data(cross_sections)
self.materials = openmc.Materials(
self.model.geometry.get_all_materials().values()
)
# Select nuclides with data that are also in the chain
self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
if nuc.name in self.nuclides_with_data]
def _load_previous_results(self):
"""Load results from a previous depletion simulation"""
# Reload volumes into geometry
self.prev_res[-1].transfer_volumes(self.model)
# Extract number densities from the geometry / previous depletion run
self._extract_number(self.local_mats, volume, nuclides, self.prev_res)
# Store previous results in operator
# Distribute reaction rates according to those tracked
# on this process
if comm.size != 1:
prev_results = self.prev_res
self.prev_res = Results()
mat_indexes = _distribute(range(len(self.burnable_mats)))
for res_obj in prev_results:
new_res = res_obj.distribute(self.local_mats, mat_indexes)
self.prev_res.append(new_res)
# Create reaction rates array
self.reaction_rates = ReactionRates(
self.local_mats, self._burnable_nucs, self.chain.reactions)
def _get_nuclides_with_data(self, cross_sections):
"""Loads cross_sections.xml file to find nuclides with neutron data
Parameters
----------
cross_sections : str
Path to cross_sections.xml file
Returns
-------
nuclides : set of str
Set of nuclide names that have cross secton data
"""
nuclides = set()
data_lib = DataLibrary.from_xml(cross_sections)
for library in data_lib.libraries:
if library['type'] != 'neutron':
continue
for name in library['materials']:
if name not in nuclides:
nuclides.add(name)
return nuclides
def _get_helper_classes(
self,
reaction_rate_mode,
normalization_mode,
fission_yield_mode,
reaction_rate_opts,
fission_yield_opts):
"""Create the ``_rate_helper``, ``_normalization_helper``, and
``_yield_helper`` objects.
Parameters
----------
reaction_rate_mode : str
Indicates the subclass of :class:`ReactionRateHelper` to
instantiate.
normalization_mode : str
Indicates the subclass of :class:`NormalizationHelper` to
instatiate.
fission_yield_mode : str
Indicates the subclass of :class:`FissionYieldHelper` to instantiate.
reaction_rate_opts : dict
Keyword arguments that are passed to the :class:`ReactionRateHelper`
subclass.
fission_yield_opts : dict
Keyword arguments that are passed to the :class:`FissionYieldHelper`
subclass.
"""
# Get classes to assist working with tallies
if reaction_rate_mode == "direct":
self._rate_helper = DirectReactionRateHelper(
@ -342,6 +379,48 @@ class Operator(TransportOperator):
self._yield_helper = fission_helper.from_operator(
self, **fission_yield_opts)
def initial_condition(self):
"""Performs final setup and returns initial condition.
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
# Create XML files
if comm.rank == 0:
self.geometry.export_to_xml()
self.settings.export_to_xml()
self._generate_materials_xml()
# Initialize OpenMC library
comm.barrier()
if not openmc.lib.is_initialized:
openmc.lib.init(intracomm=comm)
# Generate tallies in memory
materials = [openmc.lib.materials[int(i)]
for i in self.burnable_mats]
return super().initial_condition(materials)
def _generate_materials_xml(self):
"""Creates materials.xml from self.number.
Due to uncertainty with how MPI interacts with OpenMC API, this
constructs the XML manually. The long term goal is to do this
through direct memory writing.
"""
# Sort nuclides according to order in AtomNumber object
nuclides = list(self.number.nuclides)
for mat in self.materials:
mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
self.materials.export_to_xml()
def __call__(self, vec, source_rate):
"""Runs a simulation.
@ -365,9 +444,7 @@ class Operator(TransportOperator):
# Reset results in OpenMC
openmc.lib.reset()
# Update the number densities regardless of the source rate
self.number.set_density(vec)
self._update_materials()
self._update_materials_and_nuclides(vec)
# If the source rate is zero, return zero reaction rates without running
# a transport solve
@ -376,235 +453,20 @@ class Operator(TransportOperator):
rates.fill(0.0)
return OperatorResult(ufloat(0.0, 0.0), rates)
# Prevent OpenMC from complaining about re-creating tallies
openmc.reset_auto_ids()
# Update tally nuclides data in preparation for transport solve
nuclides = self._get_tally_nuclides()
self._rate_helper.nuclides = nuclides
self._normalization_helper.nuclides = nuclides
self._yield_helper.update_tally_nuclides(nuclides)
# Run OpenMC
openmc.lib.run()
openmc.lib.reset_timers()
# Extract results
op_result = self._unpack_tallies_and_normalize(source_rate)
rates = self._calculate_reaction_rates(source_rate)
# Get k and uncertainty
keff = ufloat(*openmc.lib.keff())
op_result = OperatorResult(keff, rates)
return copy.deepcopy(op_result)
@staticmethod
def write_bos_data(step):
"""Write a state-point file with beginning of step data
Parameters
----------
step : int
Current depletion step including restarts
"""
openmc.lib.statepoint_write(
"openmc_simulation_n{}.h5".format(step),
write_source=False)
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material
"""
# Count the number of instances for each cell and material
self.geometry.determine_paths(instances_only=True)
# Extract all burnable materials which have multiple instances
distribmats = set(
[mat for mat in self.materials
if mat.depletable and mat.num_instances > 1])
for mat in distribmats:
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
mat.volume /= mat.num_instances
if distribmats:
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats:
mat = cell.fill
cell.fill = [mat.clone()
for i in range(cell.num_instances)]
def _get_burnable_mats(self):
"""Determine depletable materials, volumes, and nuclides
Returns
-------
burnable_mats : list of str
List of burnable material IDs
volume : OrderedDict of str to float
Volume of each material in [cm^3]
nuclides : list of str
Nuclides in order of how they'll appear in the simulation.
"""
burnable_mats = set()
model_nuclides = set()
volume = OrderedDict()
self.heavy_metal = 0.0
# Iterate once through the geometry to get dictionaries
for mat in self.materials:
for nuclide in mat.get_nuclides():
model_nuclides.add(nuclide)
if mat.depletable:
burnable_mats.add(str(mat.id))
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
volume[str(mat.id)] = mat.volume
self.heavy_metal += mat.fissionable_mass
# Make sure there are burnable materials
if not burnable_mats:
raise RuntimeError(
"No depletable materials were found in the model.")
# Sort the sets
burnable_mats = sorted(burnable_mats, key=int)
model_nuclides = sorted(model_nuclides)
# Construct a global nuclide dictionary, burned first
nuclides = list(self.chain.nuclide_dict)
for nuc in model_nuclides:
if nuc not in nuclides:
nuclides.append(nuc)
return burnable_mats, volume, nuclides
def _extract_number(self, local_mats, volume, nuclides, prev_res=None):
"""Construct AtomNumber using geometry
Parameters
----------
local_mats : list of str
Material IDs to be managed by this process
volume : OrderedDict of str to float
Volumes for the above materials in [cm^3]
nuclides : list of str
Nuclides to be used in the simulation.
prev_res : Results, optional
Results from a previous depletion calculation
"""
self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain))
if self.dilute_initial != 0.0:
for nuc in self._burnable_nucs:
self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial)
# Now extract and store the number densities
# From the geometry if no previous depletion results
if prev_res is None:
for mat in self.materials:
if str(mat.id) in local_mats:
self._set_number_from_mat(mat)
# Else from previous depletion results
else:
for mat in self.materials:
if str(mat.id) in local_mats:
self._set_number_from_results(mat, prev_res)
def _set_number_from_mat(self, mat):
"""Extracts material and number densities from openmc.Material
Parameters
----------
mat : openmc.Material
The material to read from
"""
mat_id = str(mat.id)
for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items():
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
def _set_number_from_results(self, mat, prev_res):
"""Extracts material nuclides and number densities.
If the nuclide concentration's evolution is tracked, the densities come
from depletion results. Else, densities are extracted from the geometry
in the summary.
Parameters
----------
mat : openmc.Material
The material to read from
prev_res : Results
Results from a previous depletion calculation
"""
mat_id = str(mat.id)
# Get nuclide lists from geometry and depletion results
depl_nuc = prev_res[-1].nuc_to_ind
geom_nuc_densities = mat.get_nuclide_atom_densities()
# Merge lists of nuclides, with the same order for every calculation
geom_nuc_densities.update(depl_nuc)
for nuclide, atom_per_bcm in geom_nuc_densities.items():
if nuclide in depl_nuc:
concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1]
volume = prev_res[-1].volume[mat_id]
atom_per_cc = concentration / volume
else:
atom_per_cc = atom_per_bcm * 1.0e24
self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
def initial_condition(self):
"""Performs final setup and returns initial condition.
Returns
-------
list of numpy.ndarray
Total density for initial conditions.
"""
# Create XML files
if comm.rank == 0:
self.geometry.export_to_xml()
self.settings.export_to_xml()
self._generate_materials_xml()
# Initialize OpenMC library
comm.barrier()
if not openmc.lib.is_initialized:
openmc.lib.init(intracomm=comm)
# Generate tallies in memory
materials = [openmc.lib.materials[int(i)]
for i in self.burnable_mats]
self._rate_helper.generate_tallies(materials, self.chain.reactions)
self._normalization_helper.prepare(
self.chain.nuclides, self.reaction_rates.index_nuc)
# Tell fission yield helper what materials this process is
# responsible for
self._yield_helper.generate_tallies(
materials, tuple(sorted(self._mat_index_map.values())))
# Return number density vector
return list(self.number.get_mat_slice(np.s_[:]))
def finalize(self):
"""Finalize a depletion simulation and release resources."""
if self.cleanup_when_done:
openmc.lib.finalize()
def _update_materials(self):
"""Updates material compositions in OpenMC on all processes."""
@ -631,192 +493,41 @@ class Operator(TransportOperator):
densities.append(val)
else:
# Only output warnings if values are significantly
# negative. CRAM does not guarantee positive values.
# negative. CRAM does not guarantee positive
# values.
if val < -1.0e-21:
print("WARNING: nuclide ", nuc, " in material ", mat,
" is negative (density = ", val, " at/barn-cm)")
print(
"WARNING: nuclide ",
nuc,
" in material ",
mat,
" is negative (density = ",
val,
" at/barn-cm)")
number_i[mat, nuc] = 0.0
# Update densities on C API side
mat_internal = openmc.lib.materials[int(mat)]
mat_internal.set_densities(nuclides, densities)
#TODO Update densities on the Python side, otherwise the
# TODO Update densities on the Python side, otherwise the
# summary.h5 file contains densities at the first time step
def _generate_materials_xml(self):
"""Creates materials.xml from self.number.
Due to uncertainty with how MPI interacts with OpenMC API, this
constructs the XML manually. The long term goal is to do this
through direct memory writing.
"""
# Sort nuclides according to order in AtomNumber object
nuclides = list(self.number.nuclides)
for mat in self.materials:
mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
self.materials.export_to_xml()
def _get_tally_nuclides(self):
"""Determine nuclides that should be tallied for reaction rates.
This method returns a list of all nuclides that have neutron data and
are listed in the depletion chain. Technically, we should tally nuclides
that may not appear in the depletion chain because we still need to get
the fission reaction rate for these nuclides in order to normalize
power, but that is left as a future exercise.
Returns
-------
list of str
Tally nuclides
"""
nuc_set = set()
# Create the set of all nuclides in the decay chain in materials marked
# for burning in which the number density is greater than zero.
for nuc in self.number.nuclides:
if nuc in self.nuclides_with_data:
if np.sum(self.number[:, nuc]) > 0.0:
nuc_set.add(nuc)
# Communicate which nuclides have nonzeros to rank 0
if comm.rank == 0:
for i in range(1, comm.size):
nuc_newset = comm.recv(source=i, tag=i)
nuc_set |= nuc_newset
else:
comm.send(nuc_set, dest=0, tag=comm.rank)
if comm.rank == 0:
# Sort nuclides in the same order as self.number
nuc_list = [nuc for nuc in self.number.nuclides
if nuc in nuc_set]
else:
nuc_list = None
# Store list of tally nuclides on each process
nuc_list = comm.bcast(nuc_list)
return [nuc for nuc in nuc_list if nuc in self.chain]
def _unpack_tallies_and_normalize(self, source_rate):
"""Unpack tallies from OpenMC and return an operator result
This method uses OpenMC's C API bindings to determine the k-effective
value and reaction rates from the simulation. The reaction rates are
normalized by a helper class depending on the method being used.
@staticmethod
def write_bos_data(step):
"""Write a state-point file with beginning of step data
Parameters
----------
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
openmc.deplete.OperatorResult
Eigenvalue and reaction rates resulting from transport operator
step : int
Current depletion step including restarts
"""
rates = self.reaction_rates
rates.fill(0.0)
openmc.lib.statepoint_write(
"openmc_simulation_n{}.h5".format(step),
write_source=False)
# Get k and uncertainty
keff = ufloat(*openmc.lib.keff())
# Extract tally bins
nuclides = self._rate_helper.nuclides
# Form fast map
nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
react_ind = [rates.index_rx[react] for react in self.chain.reactions]
# Keep track of energy produced from all reactions in eV per source
# particle
self._normalization_helper.reset()
self._yield_helper.unpack()
# Store fission yield dictionaries
fission_yields = []
# Create arrays to store fission Q values, reaction rates, and nuclide
# numbers, zeroed out in material iteration
number = np.empty(rates.n_nuc)
fission_ind = rates.index_rx.get("fission")
# Extract results
for i, mat in enumerate(self.local_mats):
# Get tally index
mat_index = self._mat_index_map[mat]
# Zero out reaction rates and nuclide numbers
number.fill(0.0)
# Get new number densities
for nuc, i_nuc_results in zip(nuclides, nuc_ind):
number[i_nuc_results] = self.number[mat, nuc]
tally_rates = self._rate_helper.get_material_rates(
mat_index, nuc_ind, react_ind)
# Compute fission yields for this material
fission_yields.append(self._yield_helper.weighted_yields(i))
# Accumulate energy from fission
if fission_ind is not None:
self._normalization_helper.update(tally_rates[:, fission_ind])
# Divide by total number and store
rates[i] = self._rate_helper.divide_by_adens(number)
# Scale reaction rates to obtain units of reactions/sec
rates *= self._normalization_helper.factor(source_rate)
# Store new fission yields on the chain
self.chain.fission_yields = fission_yields
return OperatorResult(keff, rates)
def _get_nuclides_with_data(self, cross_sections):
"""Loads cross_sections.xml file to find nuclides with neutron data"""
nuclides = set()
data_lib = DataLibrary.from_xml(cross_sections)
for library in data_lib.libraries:
if library['type'] != 'neutron':
continue
for name in library['materials']:
if name not in nuclides:
nuclides.add(name)
return nuclides
def get_results_info(self):
"""Returns volume list, material lists, and nuc lists.
Returns
-------
volume : dict of str float
Volumes corresponding to materials in full_burn_dict
nuc_list : list of str
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all material IDs to be burned. Used for sorting the simulation.
full_burn_list : list
List of all burnable material IDs
"""
nuc_list = self.number.burnable_nuclides
burn_list = self.local_mats
volume = {}
for i, mat in enumerate(burn_list):
volume[mat] = self.number.volume[i]
# Combine volume dictionaries across processes
volume_list = comm.allgather(volume)
volume = {k: v for d in volume_list for k, v in d.items()}
return volume, nuc_list, burn_list, self.burnable_mats
def finalize(self):
"""Finalize a depletion simulation and release resources."""
if self.cleanup_when_done:
openmc.lib.finalize()

12
openmc/universe.py
View file

@ -7,6 +7,7 @@ from pathlib import Path
from tempfile import TemporaryDirectory
from xml.etree import ElementTree as ET
import h5py
import numpy as np
import openmc
@ -630,6 +631,9 @@ class DAGMCUniverse(UniverseBase):
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps
in ID space between OpenMC and UWUW materials (False)
bounding_box : 2-tuple of numpy.array
Lower-left and upper-right coordinates of an axis-aligned bounding box
of the universe.
"""
def __init__(self,
@ -650,6 +654,14 @@ class DAGMCUniverse(UniverseBase):
string += '{: <16}=\t{}\n'.format('\tFile', self.filename)
return string
@property
def bounding_box(self):
with h5py.File(self.filename) as dagmc_file:
coords = dagmc_file['tstt']['nodes']['coordinates'][()]
lower_left_corner = coords.min(axis=0)
upper_right_corner = coords.max(axis=0)
return (lower_left_corner, upper_right_corner)
@property
def filename(self):
return self._filename

5
tests/regression_tests/dagmc/universes/test.py
View file

@ -46,6 +46,11 @@ class DAGMCUniverseTest(PyAPITestHarness):
# create the DAGMC universe
pincell_univ = openmc.DAGMCUniverse(filename='dagmc.h5m', auto_geom_ids=True)
# checks that the bounding box is calculated correctly
bounding_box = pincell_univ.bounding_box
assert bounding_box[0].tolist() == [-25., -25., -25.]
assert bounding_box[1].tolist() == [25., 25., 25.]
# create a 2 x 2 lattice using the DAGMC pincell
pitch = np.asarray((24.0, 24.0))
lattice = openmc.RectLattice()