diff --git a/chains/chain_simple.xml b/chains/chain_simple.xml
new file mode 100644
index 000000000..345da2237
--- /dev/null
+++ b/chains/chain_simple.xml
@@ -0,0 +1,49 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 2.53000e-02
+
+ Gd157 Gd156 I135 Xe135 Xe136 Cs135
+ 1.093250e-04 2.087260e-04 2.780820e-02 6.759540e-03 2.392300e-02 4.356330e-05
+
+
+
+
+
+
+ 2.53000e-02
+
+ Gd157 Gd156 I135 Xe135 Xe136 Cs135
+ 6.142710e-5 1.483250e-04 0.0292737 0.002566345 0.0219242 4.9097e-6
+
+
+
+
+
+
+ 2.53000e-02
+
+ Gd157 Gd156 I135 Xe135 Xe136 Cs135
+ 4.141120e-04 7.605360e-04 0.0135457 0.00026864 0.0024432 3.7100E-07
+
+
+
+
diff --git a/chains/chain_test.xml b/chains/chain_test.xml
new file mode 100644
index 000000000..598570406
--- /dev/null
+++ b/chains/chain_test.xml
@@ -0,0 +1,23 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 0.0253
+
+ A B
+ 0.0292737 0.002566345
+
+
+
+
diff --git a/docs/source/pythonapi/deplete/index.rst b/docs/source/pythonapi/deplete/index.rst
new file mode 100644
index 000000000..55380c7a1
--- /dev/null
+++ b/docs/source/pythonapi/deplete/index.rst
@@ -0,0 +1,56 @@
+.. _api:
+
+=================
+API Documentation
+=================
+
+Integrators
+-----------
+
+.. toctree::
+ :maxdepth: 2
+
+ integrator.predictor
+ integrator.cecm
+
+Integrator Helper Functions
+---------------------------
+.. toctree::
+ :maxdepth: 2
+
+ integrator.CRAM16
+ integrator.CRAM48
+ integrator.save_results
+
+Metaclasses
+-----------
+
+.. autosummary::
+ :toctree: generated
+ :nosignatures:
+
+ opendeplete.Settings
+ opendeplete.Operator
+
+OpenMC Classes
+--------------
+
+.. autosummary::
+ :toctree: generated
+ :nosignatures:
+
+ opendeplete.OpenMCSettings
+ opendeplete.Materials
+ opendeplete.OpenMCOperator
+
+Data Classes
+------------
+.. autosummary::
+ :toctree: generated
+ :nosignatures:
+
+ opendeplete.AtomNumber
+ opendeplete.DepletionChain
+ opendeplete.Nuclide
+ opendeplete.ReactionRates
+ opendeplete.Results
diff --git a/docs/source/pythonapi/deplete/integrator.CRAM16.rst b/docs/source/pythonapi/deplete/integrator.CRAM16.rst
new file mode 100644
index 000000000..f9eba273e
--- /dev/null
+++ b/docs/source/pythonapi/deplete/integrator.CRAM16.rst
@@ -0,0 +1,6 @@
+integrator\.CRAM16
+==================
+
+.. currentmodule:: opendeplete.integrator
+
+.. autofunction:: CRAM16
diff --git a/docs/source/pythonapi/deplete/integrator.CRAM48.rst b/docs/source/pythonapi/deplete/integrator.CRAM48.rst
new file mode 100644
index 000000000..d7467a418
--- /dev/null
+++ b/docs/source/pythonapi/deplete/integrator.CRAM48.rst
@@ -0,0 +1,6 @@
+integrator\.CRAM48
+==================
+
+.. currentmodule:: opendeplete.integrator
+
+.. autofunction:: CRAM48
diff --git a/docs/source/pythonapi/deplete/integrator.cecm.rst b/docs/source/pythonapi/deplete/integrator.cecm.rst
new file mode 100644
index 000000000..507a638f6
--- /dev/null
+++ b/docs/source/pythonapi/deplete/integrator.cecm.rst
@@ -0,0 +1,6 @@
+integrator\.cecm
+=================
+
+.. currentmodule:: opendeplete.integrator
+
+.. autofunction:: cecm
diff --git a/docs/source/pythonapi/deplete/integrator.predictor.rst b/docs/source/pythonapi/deplete/integrator.predictor.rst
new file mode 100644
index 000000000..d6c0fd827
--- /dev/null
+++ b/docs/source/pythonapi/deplete/integrator.predictor.rst
@@ -0,0 +1,6 @@
+integrator\.predictor
+=====================
+
+.. currentmodule:: opendeplete.integrator
+
+.. autofunction:: predictor
diff --git a/docs/source/pythonapi/deplete/integrator.save_results.rst b/docs/source/pythonapi/deplete/integrator.save_results.rst
new file mode 100644
index 000000000..5c21dcb66
--- /dev/null
+++ b/docs/source/pythonapi/deplete/integrator.save_results.rst
@@ -0,0 +1,6 @@
+integrator\.save_results
+========================
+
+.. currentmodule:: opendeplete.integrator
+
+.. autofunction:: save_results
diff --git a/docs/source/pythonapi/deplete/opendeplete.Concentrations.rst b/docs/source/pythonapi/deplete/opendeplete.Concentrations.rst
new file mode 100644
index 000000000..6fa07a970
--- /dev/null
+++ b/docs/source/pythonapi/deplete/opendeplete.Concentrations.rst
@@ -0,0 +1,30 @@
+opendeplete.Concentrations
+==========================
+
+.. currentmodule:: opendeplete
+
+.. autoclass:: Concentrations
+
+
+ .. automethod:: __init__
+
+
+ .. rubric:: Methods
+
+ .. autosummary::
+
+ ~Concentrations.__init__
+ ~Concentrations.convert_nested_dict
+
+
+
+
+
+ .. rubric:: Attributes
+
+ .. autosummary::
+
+ ~Concentrations.n_cell
+ ~Concentrations.n_nuc
+
+
\ No newline at end of file
diff --git a/docs/source/pythonapi/deplete/opendeplete.ReactionRates.rst b/docs/source/pythonapi/deplete/opendeplete.ReactionRates.rst
new file mode 100644
index 000000000..99e048b56
--- /dev/null
+++ b/docs/source/pythonapi/deplete/opendeplete.ReactionRates.rst
@@ -0,0 +1,30 @@
+opendeplete.ReactionRates
+=========================
+
+.. currentmodule:: opendeplete
+
+.. autoclass:: ReactionRates
+
+
+ .. automethod:: __init__
+
+
+ .. rubric:: Methods
+
+ .. autosummary::
+
+ ~ReactionRates.__init__
+
+
+
+
+
+ .. rubric:: Attributes
+
+ .. autosummary::
+
+ ~ReactionRates.n_cell
+ ~ReactionRates.n_nuc
+ ~ReactionRates.n_react
+
+
\ No newline at end of file
diff --git a/docs/source/pythonapi/deplete/opendeplete.Results.rst b/docs/source/pythonapi/deplete/opendeplete.Results.rst
new file mode 100644
index 000000000..0ab8a1f71
--- /dev/null
+++ b/docs/source/pythonapi/deplete/opendeplete.Results.rst
@@ -0,0 +1,22 @@
+opendeplete.Results
+===================
+
+.. currentmodule:: opendeplete
+
+.. autoclass:: Results
+
+
+ .. automethod:: __init__
+
+
+ .. rubric:: Methods
+
+ .. autosummary::
+
+ ~Results.__init__
+
+
+
+
+
+
\ No newline at end of file
diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py
new file mode 100644
index 000000000..994a51e12
--- /dev/null
+++ b/openmc/deplete/__init__.py
@@ -0,0 +1,24 @@
+"""
+OpenDeplete
+===========
+
+A simple depletion front-end tool.
+"""
+
+from .dummy_comm import DummyCommunicator
+try:
+ from mpi4py import MPI
+ comm = MPI.COMM_WORLD
+ have_mpi = True
+except ImportError:
+ comm = DummyCommunicator()
+ have_mpi = False
+
+from .nuclide import *
+from .depletion_chain import *
+from .openmc_wrapper import *
+from .reaction_rates import *
+from .function import *
+from .results import *
+from .integrator import *
+from .utilities import *
diff --git a/openmc/deplete/atom_number.py b/openmc/deplete/atom_number.py
new file mode 100644
index 000000000..03bedbf53
--- /dev/null
+++ b/openmc/deplete/atom_number.py
@@ -0,0 +1,236 @@
+"""AtomNumber module.
+
+An ndarray to store atom densities with string, integer, or slice indexing.
+"""
+
+import numpy as np
+
+
+class AtomNumber(object):
+ """ AtomNumber module.
+
+ An ndarray to store atom densities with string, integer, or slice indexing.
+
+ Parameters
+ ----------
+ mat_to_ind : OrderedDict of str to int
+ A dictionary mapping material ID as string to index.
+ nuc_to_ind : OrderedDict of str to int
+ A dictionary mapping nuclide name as string to index.
+ volume : OrderedDict of int to float
+ Volume of geometry.
+ n_mat_burn : int
+ Number of materials to be burned.
+ n_nuc_burn : int
+ Number of nuclides to be burned.
+
+ Attributes
+ ----------
+ mat_to_ind : OrderedDict of str to int
+ A dictionary mapping cell ID as string to index.
+ nuc_to_ind : OrderedDict of str to int
+ A dictionary mapping nuclide name as string to index.
+ volume : numpy.array
+ Volume of geometry indexed by mat_to_ind. If a volume is not found,
+ it defaults to 1 so that reading density still works correctly.
+ n_mat_burn : int
+ Number of materials to be burned.
+ n_nuc_burn : int
+ Number of nuclides to be burned.
+ n_mat : int
+ Number of materials.
+ n_nuc : int
+ Number of nucs.
+ number : numpy.array
+ Array storing total atoms indexed by the above dictionaries.
+ burn_nuc_list : list of str
+ A list of all nuclide material names. Used for sorting the simulation.
+ burn_mat_list : list of str
+ A list of all burning material names. Used for sorting the simulation.
+ """
+
+ def __init__(self, mat_to_ind, nuc_to_ind, volume, n_mat_burn, n_nuc_burn):
+
+ self.mat_to_ind = mat_to_ind
+ self.nuc_to_ind = nuc_to_ind
+
+ self.volume = np.ones(self.n_mat)
+
+ for mat in volume:
+ if str(mat) in self.mat_to_ind:
+ ind = self.mat_to_ind[str(mat)]
+ self.volume[ind] = volume[mat]
+
+ self.n_mat_burn = n_mat_burn
+ self.n_nuc_burn = n_nuc_burn
+
+ self.number = np.zeros((self.n_mat, self.n_nuc))
+
+ # For performance, create storage for burn_nuc_list, burn_mat_list
+ self._burn_nuc_list = None
+ self._burn_mat_list = None
+
+ def __getitem__(self, pos):
+ """ Retrieves total atom number from AtomNumber.
+
+ Parameters
+ ----------
+ pos : tuple
+ A two-length tuple containing a material index and a nuc index.
+ These indexes can be strings (which get converted to integers via
+ the dictionaries), integers used directly, or slices.
+
+ Returns
+ -------
+ numpy.array
+ The value indexed from self.number.
+ """
+
+ mat, nuc = pos
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+ if isinstance(nuc, str):
+ nuc = self.nuc_to_ind[nuc]
+
+ return self.number[mat, nuc]
+
+ def __setitem__(self, pos, val):
+ """ Sets total atom number into AtomNumber.
+
+ Parameters
+ ----------
+ pos : tuple
+ A two-length tuple containing a material index and a nuc index.
+ These indexes can be strings (which get converted to integers via
+ the dictionaries), integers used directly, or slices.
+ val : float
+ The value to set the array to.
+ """
+
+ mat, nuc = pos
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+ if isinstance(nuc, str):
+ nuc = self.nuc_to_ind[nuc]
+
+ self.number[mat, nuc] = val
+
+ def get_atom_density(self, mat, nuc):
+ """ Accesses atom density instead of total number.
+
+ Parameters
+ ----------
+ mat : str, int or slice
+ Material index.
+ nuc : str, int or slice
+ Nuclide index.
+
+ Returns
+ -------
+ numpy.array
+ The density indexed.
+ """
+
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+ if isinstance(nuc, str):
+ nuc = self.nuc_to_ind[nuc]
+
+ return self[mat, nuc] / self.volume[mat]
+
+ def set_atom_density(self, mat, nuc, val):
+ """ Sets atom density instead of total number.
+
+ Parameters
+ ----------
+ mat : str, int or slice
+ Material index.
+ nuc : str, int or slice
+ Nuclide index.
+ val : numpy.array
+ Array of values to set.
+ """
+
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+ if isinstance(nuc, str):
+ nuc = self.nuc_to_ind[nuc]
+
+ self[mat, nuc] = val * self.volume[mat]
+
+ def get_mat_slice(self, mat):
+ """ Gets atom quantity indexed by mats for all burned nuclides
+
+ Parameters
+ ----------
+ mat : str, int or slice
+ Material index.
+
+ Returns
+ -------
+ numpy.array
+ The slice requested.
+ """
+
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+
+ return self[mat, 0:self.n_nuc_burn]
+
+ def set_mat_slice(self, mat, val):
+ """ Sets atom quantity indexed by mats for all burned nuclides
+
+ Parameters
+ ----------
+ mat : str, int or slice
+ Material index.
+ val : numpy.array
+ The slice to set.
+ """
+
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+
+ self[mat, 0:self.n_nuc_burn] = val
+
+ @property
+ def n_mat(self):
+ """Number of materials."""
+ return len(self.mat_to_ind)
+
+ @property
+ def n_nuc(self):
+ """Number of nuclides."""
+ return len(self.nuc_to_ind)
+
+ @property
+ def burn_nuc_list(self):
+ """ burn_nuc_list : list of str
+ A list of all nuclide material names. Used for sorting the simulation.
+ """
+
+ if self._burn_nuc_list is None:
+ self._burn_nuc_list = [None] * self.n_nuc_burn
+
+ for nuc in self.nuc_to_ind:
+ ind = self.nuc_to_ind[nuc]
+ if ind < self.n_nuc_burn:
+ self._burn_nuc_list[ind] = nuc
+
+ return self._burn_nuc_list
+
+ @property
+ def burn_mat_list(self):
+ """ burn_mat_list : list of str
+ A list of all burning material names. Used for sorting the simulation.
+ """
+
+ if self._burn_mat_list is None:
+ self._burn_mat_list = [None] * self.n_mat_burn
+
+ for mat in self.mat_to_ind:
+ ind = self.mat_to_ind[mat]
+ if ind < self.n_mat_burn:
+ self._burn_mat_list[ind] = mat
+
+ return self._burn_mat_list
diff --git a/openmc/deplete/depletion_chain.py b/openmc/deplete/depletion_chain.py
new file mode 100644
index 000000000..05cc9db43
--- /dev/null
+++ b/openmc/deplete/depletion_chain.py
@@ -0,0 +1,472 @@
+"""depletion_chain module.
+
+This module contains information about a depletion chain. A depletion chain is
+loaded from an .xml file and all the nuclides are linked together.
+"""
+
+from collections import OrderedDict, defaultdict
+from io import StringIO
+from itertools import chain
+import math
+import re
+import os
+
+from tqdm import tqdm
+import scipy.sparse as sp
+import openmc.data
+# Try to use lxml if it is available. It preserves the order of attributes and
+# provides a pretty-printer by default. If not available, use OpenMC function to
+# pretty print.
+try:
+ import lxml.etree as ET
+ _have_lxml = True
+except ImportError:
+ import xml.etree.ElementTree as ET
+ from openmc.clean_xml import clean_xml_indentation
+ _have_lxml = False
+
+from .nuclide import Nuclide, DecayTuple, ReactionTuple
+
+
+# tuple of (reaction name, possible MT values, (dA, dZ)) where dA is the change
+# in the mass number and dZ is the change in the atomic number
+_REACTIONS = [
+ ('(n,2n)', set(chain([16], range(875, 892))), (-1, 0)),
+ ('(n,3n)', {17}, (-2, 0)),
+ ('(n,4n)', {37}, (-3, 0)),
+ ('(n,gamma)', {102}, (1, 0)),
+ ('(n,p)', set(chain([103], range(600, 650))), (0, -1)),
+ ('(n,a)', set(chain([107], range(800, 850))), (-3, -2))
+]
+
+
+def _get_zai(s):
+ """Get ZAI value (10000*z + 10*A + metastable state) for sorting purposes"""
+ symbol, A, state = re.match(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)', s).groups()
+ Z = openmc.data.ATOMIC_NUMBER[symbol]
+ A = int(A)
+ state = int(state[2:]) if state else 0
+ return 10000*Z + 10*A + state
+
+
+def replace_missing(product, decay_data):
+ """Replace missing product with suitable decay daughter.
+
+ Parameters
+ ----------
+ product : str
+ Name of product in GND format, e.g. 'Y86_m1'.
+ decay_data : dict
+ Dictionary of decay data
+
+ Returns
+ -------
+ product : str
+ Replacement for missing product in GND format.
+
+ """
+
+ symbol, A, state = re.match(r'([A-Zn][a-z]*)(\d+)((?:_m\d+)?)',
+ product).groups()
+ Z = openmc.data.ATOMIC_NUMBER[symbol]
+ A = int(A)
+
+ # First check if ground state is available
+ if state:
+ metastable_state = int(state[2:])
+ product = '{}{}'.format(symbol, A)
+
+ # Find isotope with longest half-life
+ half_life = 0.0
+ for nuclide, data in decay_data.items():
+ m = re.match(r'{}(\d+)(?:_m\d+)?'.format(symbol), nuclide)
+ if m:
+ # If we find a stable nuclide, stop search
+ if data.nuclide['stable']:
+ mass_longest_lived = int(m.group(1))
+ break
+ if data.half_life.nominal_value > half_life:
+ mass_longest_lived = int(m.group(1))
+ half_life = data.half_life.nominal_value
+
+ # If mass number of longest-lived isotope is less than that of missing
+ # product, assume it undergoes beta-. Otherwise assume beta+.
+ beta_minus = (mass_longest_lived < A)
+
+ # Iterate until we find an existing nuclide
+ while product not in decay_data:
+ if Z > 98:
+ Z -= 2
+ A -= 4
+ else:
+ if beta_minus:
+ Z += 1
+ else:
+ Z -= 1
+ product = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
+
+ return product
+
+
+class DepletionChain(object):
+ """ The DepletionChain class.
+
+ This class contains a full representation of a depletion chain.
+
+ Attributes
+ ----------
+ n_nuclides : int
+ Number of nuclides in chain.
+ nuclides : list of Nuclide
+ List of nuclides in chain.
+ nuclide_dict : OrderedDict of str to int
+ Maps a nuclide name to an index in nuclides.
+ nuc_to_react_ind : OrderedDict of str to int
+ Dictionary mapping a nuclide name to an index in ReactionRates.
+ react_to_ind : OrderedDict of str to int
+ Dictionary mapping a reaction name to an index in ReactionRates.
+
+ """
+
+ def __init__(self):
+ self.nuclides = []
+ self.nuclide_dict = OrderedDict()
+ self.nuc_to_react_ind = OrderedDict()
+ self.react_to_ind = OrderedDict()
+
+ @property
+ def n_nuclides(self):
+ """Number of nuclides in chain."""
+ return len(self.nuclides)
+
+ @classmethod
+ def from_endf(cls, decay_files, fpy_files, neutron_files):
+ """Create a depletion chain from ENDF files.
+
+ Parameters
+ ----------
+ decay_files : list of str
+ List of ENDF decay sub-library files
+ fpy_files : list of str
+ List of ENDF neutron-induced fission product yield sub-library files
+ neutron_files : list of str
+ List of ENDF neutron reaction sub-library files
+
+ """
+ depl_chain = cls()
+
+ # Create dictionary mapping target to filename
+ reactions = {}
+ with tqdm(neutron_files) as pbar:
+ for f in pbar:
+ pbar.set_description('Processing {}'.format(os.path.basename(f)))
+ evaluation = openmc.data.endf.Evaluation(f)
+ name = evaluation.gnd_name
+ reactions[name] = {}
+ for mf, mt, nc, mod in evaluation.reaction_list:
+ if mf == 3:
+ file_obj = StringIO(evaluation.section[3, mt])
+ openmc.data.endf.get_head_record(file_obj)
+ q_value = openmc.data.endf.get_cont_record(file_obj)[1]
+ reactions[name][mt] = q_value
+
+ # Determine what decay and FPY nuclides are available
+ decay_data = {}
+ with tqdm(decay_files) as pbar:
+ for f in pbar:
+ pbar.set_description('Processing {}'.format(os.path.basename(f)))
+ data = openmc.data.Decay(f)
+ decay_data[data.nuclide['name']] = data
+
+ fpy_data = {}
+ with tqdm(fpy_files) as pbar:
+ for f in pbar:
+ pbar.set_description('Processing {}'.format(os.path.basename(f)))
+ data = openmc.data.FissionProductYields(f)
+ fpy_data[data.nuclide['name']] = data
+
+ print('Creating depletion_chain...')
+ missing_daughter = []
+ missing_rx_product = []
+ missing_fpy = []
+ missing_fp = []
+
+ reaction_index = 0
+ for idx, parent in enumerate(sorted(decay_data, key=_get_zai)):
+ data = decay_data[parent]
+
+ nuclide = Nuclide()
+ nuclide.name = parent
+
+ depl_chain.nuclides.append(nuclide)
+ depl_chain.nuclide_dict[parent] = idx
+
+ if not data.nuclide['stable'] and data.half_life.nominal_value != 0.0:
+ nuclide.half_life = data.half_life.nominal_value
+ nuclide.decay_energy = sum(E.nominal_value for E in
+ data.average_energies.values())
+ sum_br = 0.0
+ for i, mode in enumerate(data.modes):
+ type_ = ','.join(mode.modes)
+ if mode.daughter in decay_data:
+ target = mode.daughter
+ else:
+ print('missing {} {} {}'.format(parent, ','.join(mode.modes), mode.daughter))
+ target = replace_missing(mode.daughter, decay_data)
+
+ # Write branching ratio, taking care to ensure sum is unity
+ br = mode.branching_ratio.nominal_value
+ sum_br += br
+ if i == len(data.modes) - 1 and sum_br != 1.0:
+ br = 1.0 - sum(m.branching_ratio.nominal_value
+ for m in data.modes[:-1])
+
+ # Append decay mode
+ nuclide.decay_modes.append(DecayTuple(type_, target, br))
+
+ if parent in reactions:
+ reactions_available = set(reactions[parent].keys())
+ for name, mts, changes in _REACTIONS:
+ if mts & reactions_available:
+ delta_A, delta_Z = changes
+ A = data.nuclide['mass_number'] + delta_A
+ Z = data.nuclide['atomic_number'] + delta_Z
+ daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
+
+ if name not in depl_chain.react_to_ind:
+ depl_chain.react_to_ind[name] = reaction_index
+ reaction_index += 1
+
+ if daughter not in decay_data:
+ missing_rx_product.append((parent, name, daughter))
+
+ # Store Q value
+ for mt in sorted(mts):
+ if mt in reactions[parent]:
+ q_value = reactions[parent][mt]
+ break
+ else:
+ q_value = 0.0
+
+ nuclide.reactions.append(ReactionTuple(
+ name, daughter, q_value, 1.0))
+
+ if any(mt in reactions_available for mt in [18, 19, 20, 21, 38]):
+ if parent in fpy_data:
+ q_value = reactions[parent][18]
+ nuclide.reactions.append(
+ ReactionTuple('fission', 0, q_value, 1.0))
+
+ if 'fission' not in depl_chain.react_to_ind:
+ depl_chain.react_to_ind['fission'] = reaction_index
+ reaction_index += 1
+ else:
+ missing_fpy.append(parent)
+
+ if parent in fpy_data:
+ fpy = fpy_data[parent]
+
+ if fpy.energies is not None:
+ nuclide.yield_energies = fpy.energies
+ else:
+ nuclide.yield_energies = [0.0]
+
+ for E, table in zip(nuclide.yield_energies, fpy.independent):
+ yield_replace = 0.0
+ yields = defaultdict(float)
+ for product, y in table.items():
+ # Handle fission products that have no decay data available
+ if product not in decay_data:
+ daughter = replace_missing(product, decay_data)
+ product = daughter
+ yield_replace += y.nominal_value
+
+ yields[product] += y.nominal_value
+
+ if yield_replace > 0.0:
+ missing_fp.append((parent, E, yield_replace))
+
+ nuclide.yield_data[E] = []
+ for k in sorted(yields, key=_get_zai):
+ nuclide.yield_data[E].append((k, yields[k]))
+
+ # Display warnings
+ if missing_daughter:
+ print('The following decay modes have daughters with no decay data:')
+ for mode in missing_daughter:
+ print(' {}'.format(mode))
+ print('')
+
+ if missing_rx_product:
+ print('The following reaction products have no decay data:')
+ for vals in missing_rx_product:
+ print('{} {} -> {}'.format(*vals))
+ print('')
+
+ if missing_fpy:
+ print('The following fissionable nuclides have no fission product yields:')
+ for parent in missing_fpy:
+ print(' ' + parent)
+ print('')
+
+ if missing_fp:
+ print('The following nuclides have fission products with no decay data:')
+ for vals in missing_fp:
+ print(' {}, E={} eV (total yield={})'.format(*vals))
+
+ return depl_chain
+
+ @classmethod
+ def xml_read(cls, filename):
+ """Reads a depletion chain XML file.
+
+ Parameters
+ ----------
+ filename : str
+ The path to the depletion chain XML file.
+
+ Todo
+ ----
+ Allow for branching on capture, etc.
+ """
+ depl_chain = cls()
+
+ # Load XML tree
+ try:
+ root = ET.parse(filename)
+ except:
+ if filename is None:
+ print("No chain specified, either manually or in environment variable OPENDEPLETE_CHAIN.")
+ else:
+ print('Decay chain "', filename, '" is invalid.')
+ raise
+
+ reaction_index = 0
+ for i, nuclide_elem in enumerate(root.findall('nuclide_table')):
+ nuc = Nuclide.xml_read(nuclide_elem)
+ depl_chain.nuclide_dict[nuc.name] = i
+
+ # Check for reaction paths
+ for rx in nuc.reactions:
+ if rx.type not in depl_chain.react_to_ind:
+ depl_chain.react_to_ind[rx.type] = reaction_index
+ reaction_index += 1
+
+ depl_chain.nuclides.append(nuc)
+
+ return depl_chain
+
+ def xml_write(self, filename):
+ """Writes a depletion chain XML file.
+
+ Parameters
+ ----------
+ filename : str
+ The path to the depletion chain XML file.
+
+ """
+
+ root_elem = ET.Element('depletion')
+ for nuclide in self.nuclides:
+ root_elem.append(nuclide.xml_write())
+
+ tree = ET.ElementTree(root_elem)
+ if _have_lxml:
+ tree.write(filename, encoding='utf-8', pretty_print=True)
+ else:
+ clean_xml_indentation(root_elem, spaces_per_level=2)
+ tree.write(filename, encoding='utf-8')
+
+ def form_matrix(self, rates):
+ """ Forms depletion matrix.
+
+ Parameters
+ ----------
+ rates : numpy.ndarray
+ 2D array indexed by nuclide then by cell.
+
+ Returns
+ -------
+ scipy.sparse.csr_matrix
+ Sparse matrix representing depletion.
+ """
+
+ matrix = defaultdict(float)
+ reactions = set()
+
+ for i, nuc in enumerate(self.nuclides):
+
+ if nuc.n_decay_modes != 0:
+ # Decay paths
+ # Loss
+ decay_constant = math.log(2) / nuc.half_life
+
+ if decay_constant != 0.0:
+ matrix[i, i] -= decay_constant
+
+ # Gain
+ for _, target, branching_ratio in nuc.decay_modes:
+ # Allow for total annihilation for debug purposes
+ if target != 'Nothing':
+ branch_val = branching_ratio * decay_constant
+
+ if branch_val != 0.0:
+ k = self.nuclide_dict[target]
+ matrix[k, i] += branch_val
+
+ if nuc.name in self.nuc_to_react_ind:
+ # Extract all reactions for this nuclide in this cell
+ nuc_ind = self.nuc_to_react_ind[nuc.name]
+ nuc_rates = rates[nuc_ind, :]
+
+ for r_type, target, _, br in nuc.reactions:
+ # Extract reaction index, and then final reaction rate
+ r_id = self.react_to_ind[r_type]
+ path_rate = nuc_rates[r_id]
+
+ # Loss term -- make sure we only count loss once for
+ # reactions with branching ratios
+ if r_type not in reactions:
+ reactions.add(r_type)
+ if path_rate != 0.0:
+ matrix[i, i] -= path_rate
+
+ # Gain term; allow for total annihilation for debug purposes
+ if target != 'Nothing':
+ if r_type != 'fission':
+ if path_rate != 0.0:
+ k = self.nuclide_dict[target]
+ matrix[k, i] += path_rate * br
+ else:
+ # Assume that we should always use thermal fission
+ # yields. At some point it would be nice to account
+ # for the energy-dependence..
+ energy, data = sorted(nuc.yield_data.items())[0]
+ for product, y in data:
+ yield_val = y * path_rate
+ if yield_val != 0.0:
+ k = self.nuclide_dict[product]
+ matrix[k, i] += yield_val
+
+ # Clear set of reactions
+ reactions.clear()
+
+ # Use DOK matrix as intermediate representation, then convert to CSR and return
+ matrix_dok = sp.dok_matrix((self.n_nuclides, self.n_nuclides))
+ dict.update(matrix_dok, matrix)
+ return matrix_dok.tocsr()
+
+ def nuc_by_ind(self, ind):
+ """ Extracts nuclides from the list by dictionary key.
+
+ Parameters
+ ----------
+ ind : str
+ Name of nuclide.
+
+ Returns
+ -------
+ Nuclide
+ Nuclide object that corresponds to ind.
+ """
+ return self.nuclides[self.nuclide_dict[ind]]
diff --git a/openmc/deplete/dummy_comm.py b/openmc/deplete/dummy_comm.py
new file mode 100644
index 000000000..b3fa27264
--- /dev/null
+++ b/openmc/deplete/dummy_comm.py
@@ -0,0 +1,27 @@
+class DummyCommunicator(object):
+ rank = 0
+ size = 1
+
+ def allgather(self, sendobj):
+ return [sendobj]
+
+ def allreduce(self, sendobj, op=None):
+ return sendobj
+
+ def barrier(self):
+ pass
+
+ def bcast(self, obj, root=0):
+ return obj
+
+ def gather(self, sendobj, root=0):
+ return [sendobj]
+
+ def py2f(self):
+ return 0
+
+ def reduce(self, sendobj, op=None, root=0):
+ return sendobj
+
+ def scatter(self, sendobj, root=0):
+ return sendobj[0]
diff --git a/openmc/deplete/function.py b/openmc/deplete/function.py
new file mode 100644
index 000000000..74eb92422
--- /dev/null
+++ b/openmc/deplete/function.py
@@ -0,0 +1,114 @@
+"""function module.
+
+This module contains the Operator class, which is then passed to an integrator
+to run a full depletion simulation.
+"""
+
+from abc import ABCMeta, abstractmethod
+
+class Settings(object):
+ """ The Settings class.
+
+ Contains all parameters necessary for the integrator.
+
+ Attributes
+ ----------
+ dt_vec : numpy.array
+ Array of time steps to take.
+ output_dir : str
+ Path to output directory to save results.
+ """
+
+ def __init__(self):
+ # Integrator specific
+ self.dt_vec = None
+ self.output_dir = None
+
+class Operator(metaclass=ABCMeta):
+ """ The Operator metaclass.
+
+ This defines all functions that the integrator needs to operate.
+
+ Attributes
+ ----------
+ settings : Settings
+ Settings object.
+ """
+
+ def __init__(self, settings):
+ self.settings = settings
+
+ @abstractmethod
+ def initial_condition(self):
+ """ Performs final setup and returns initial condition.
+
+ Returns
+ -------
+ list of numpy.array
+ Total density for initial conditions.
+ """
+
+ pass
+
+ @abstractmethod
+ def eval(self, vec, print_out=True):
+ """ Runs a simulation.
+
+ Parameters
+ ----------
+ vec : list of numpy.array
+ Total atoms to be used in function.
+ print_out : bool, optional
+ Whether or not to print out time.
+
+ Returns
+ -------
+ k : float
+ Eigenvalue of the problem.
+ rates : ReactionRates
+ Reaction rates from this simulation.
+ seed : int
+ Seed for this simulation.
+ """
+
+ pass
+
+ @abstractmethod
+ def get_results_info(self):
+ """ Returns volume list, cell lists, and nuc lists.
+
+ Returns
+ -------
+ volume : list of float
+ Volumes corresponding to materials in burn_list
+ 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 cell IDs to be burned. Used for sorting the simulation.
+ full_burn_list : list of int
+ All burnable materials in the geometry.
+ """
+
+ pass
+
+ @abstractmethod
+ def form_matrix(self, y, mat):
+ """ Forms the f(y) matrix in y' = f(y)y.
+
+ Nominally a depletion matrix, this is abstracted on the off chance
+ that the function f has nothing to do with depletion at all.
+
+ Parameters
+ ----------
+ y : numpy.ndarray
+ An array representing y.
+ mat : int
+ Material id.
+
+ Returns
+ -------
+ scipy.sparse.csr_matrix
+ Sparse matrix representing f(y).
+ """
+
+ pass
diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py
new file mode 100644
index 000000000..607650dc6
--- /dev/null
+++ b/openmc/deplete/integrator/__init__.py
@@ -0,0 +1,11 @@
+"""
+Integrator
+===========
+
+The integrator subcomponents.
+"""
+
+from .cecm import *
+from .cram import *
+from .predictor import *
+from .save_results import *
diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py
new file mode 100644
index 000000000..4d9baebb2
--- /dev/null
+++ b/openmc/deplete/integrator/cecm.py
@@ -0,0 +1,133 @@
+""" The CE/CM integrator."""
+
+import copy
+from itertools import repeat
+import os
+from multiprocessing import Pool
+import time
+
+from .. import comm
+from .cram import CRAM48, cram_wrapper
+from .save_results import save_results
+
+
+def cecm(operator, print_out=True):
+ """The CE/CM integrator.
+
+ Implements the second order CE/CM Predictor-Corrector algorithm [ref]_.
+ This algorithm is mathematically defined as:
+
+ .. math::
+ y' &= A(y, t) y(t)
+
+ A_p &= A(y_n, t_n)
+
+ y_m &= \\text{expm}(A_p h/2) y_n
+
+ A_c &= A(y_m, t_n + h/2)
+
+ y_{n+1} &= \\text{expm}(A_c h) y_n
+
+ .. [ref]
+ Isotalo, Aarno. "Comparison of Neutronics-Depletion Coupling Schemes
+ for Burnup Calculations—Continued Study." Nuclear Science and
+ Engineering 180.3 (2015): 286-300.
+
+ Parameters
+ ----------
+ operator : Operator
+ The operator object to simulate on.
+ print_out : bool, optional
+ Whether or not to print out time.
+ """
+
+ # Save current directory
+ dir_home = os.getcwd()
+
+ # Move to folder
+ os.makedirs(operator.settings.output_dir, exist_ok=True)
+ os.chdir(operator.settings.output_dir)
+
+ # Generate initial conditions
+ vec = operator.initial_condition()
+
+ n_mats = len(vec)
+
+ t = 0.0
+
+ for i, dt in enumerate(operator.settings.dt_vec):
+ # Create vectors
+ x = [copy.deepcopy(vec)]
+ seeds = []
+ eigvls = []
+ rates_array = []
+
+ eigvl, rates, seed = operator.eval(x[0])
+
+ eigvls.append(eigvl)
+ seeds.append(seed)
+ rates_array.append(rates)
+
+ t_start = time.time()
+
+ chains = repeat(operator.chain, n_mats)
+ vecs = (x[0][i] for i in range(n_mats))
+ rates = (rates_array[0][i, :, :] for i in range(n_mats))
+ dts = repeat(dt/2, n_mats)
+
+ with Pool() as pool:
+ iters = zip(chains, vecs, rates, dts)
+ x_result = list(pool.starmap(cram_wrapper, iters))
+
+ t_end = time.time()
+ if comm.rank == 0:
+ if print_out:
+ print("Time to matexp: ", t_end - t_start)
+
+ x.append(x_result)
+
+ eigvl, rates, seed = operator.eval(x[1])
+
+ eigvls.append(eigvl)
+ seeds.append(seed)
+ rates_array.append(rates)
+
+ t_start = time.time()
+
+ chains = repeat(operator.chain, n_mats)
+ vecs = (x[0][i] for i in range(n_mats))
+ rates = (rates_array[1][i, :, :] for i in range(n_mats))
+ dts = repeat(dt, n_mats)
+
+ with Pool() as pool:
+ iters = zip(chains, vecs, rates, dts)
+ x_result = list(pool.starmap(cram_wrapper, iters))
+
+ t_end = time.time()
+ if comm.rank == 0:
+ if print_out:
+ print("Time to matexp: ", t_end - t_start)
+
+ # Create results, write to disk
+ save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i)
+
+ t += dt
+ vec = copy.deepcopy(x_result)
+
+ # Perform one last simulation
+ x = [copy.deepcopy(vec)]
+ seeds = []
+ eigvls = []
+ rates_array = []
+ eigvl, rates, seed = operator.eval(x[0])
+
+ eigvls.append(eigvl)
+ seeds.append(seed)
+ rates_array.append(rates)
+
+ # Create results, write to disk
+ save_results(operator, x, rates_array, eigvls, seeds, [t, t],
+ len(operator.settings.dt_vec))
+
+ # Return to origin
+ os.chdir(dir_home)
diff --git a/openmc/deplete/integrator/cram.py b/openmc/deplete/integrator/cram.py
new file mode 100644
index 000000000..a18d8450c
--- /dev/null
+++ b/openmc/deplete/integrator/cram.py
@@ -0,0 +1,185 @@
+""" Chebyshev Rational Approximation Method module
+
+Implements two different forms of CRAM for use in opendeplete.
+"""
+
+import numpy as np
+import scipy.sparse as sp
+import scipy.sparse.linalg as sla
+
+
+def cram_wrapper(chain, n0, rates, dt):
+ """Wraps depletion matrix creation / CRAM solve for multiprocess execution
+
+ Parameters
+ ----------
+ chain : DepletionChain
+ Depletion chain used to construct the burnup matrix
+ n0 : numpy.array
+ Vector to operate a matrix exponent on.
+ rates : numpy.ndarray
+ 2D array indexed by nuclide then by cell.
+ dt : float
+ Time to integrate to.
+
+ Returns
+ -------
+ numpy.array
+ Results of the matrix exponent.
+ """
+ A = chain.form_matrix(rates)
+ return CRAM48(A, n0, dt)
+
+
+def CRAM16(A, n0, dt):
+ """ Chebyshev Rational Approximation Method, order 16
+
+ Algorithm is the 16th order Chebyshev Rational Approximation Method,
+ implemented in the more stable incomplete partial fraction (IPF) form
+ [cram16]_.
+
+ .. [cram16]
+ Pusa, Maria. "Higher-Order Chebyshev Rational Approximation Method and
+ Application to Burnup Equations." Nuclear Science and Engineering 182.3
+ (2016).
+
+ Parameters
+ ----------
+ A : scipy.linalg.csr_matrix
+ Matrix to take exponent of.
+ n0 : numpy.array
+ Vector to operate a matrix exponent on.
+ dt : float
+ Time to integrate to.
+
+ Returns
+ -------
+ numpy.array
+ Results of the matrix exponent.
+ """
+
+ alpha = np.array([+2.124853710495224e-16,
+ +5.464930576870210e+3 - 3.797983575308356e+4j,
+ +9.045112476907548e+1 - 1.115537522430261e+3j,
+ +2.344818070467641e+2 - 4.228020157070496e+2j,
+ +9.453304067358312e+1 - 2.951294291446048e+2j,
+ +7.283792954673409e+2 - 1.205646080220011e+5j,
+ +3.648229059594851e+1 - 1.155509621409682e+2j,
+ +2.547321630156819e+1 - 2.639500283021502e+1j,
+ +2.394538338734709e+1 - 5.650522971778156e+0j],
+ dtype=np.complex128)
+ theta = np.array([+0.0,
+ +3.509103608414918 + 8.436198985884374j,
+ +5.948152268951177 + 3.587457362018322j,
+ -5.264971343442647 + 16.22022147316793j,
+ +1.419375897185666 + 10.92536348449672j,
+ +6.416177699099435 + 1.194122393370139j,
+ +4.993174737717997 + 5.996881713603942j,
+ -1.413928462488886 + 13.49772569889275j,
+ -10.84391707869699 + 19.27744616718165j],
+ dtype=np.complex128)
+
+ n = A.shape[0]
+
+ alpha0 = 2.124853710495224e-16
+
+ k = 8
+
+ y = np.array(n0, dtype=np.float64)
+ for l in range(1, k+1):
+ y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
+
+ y *= alpha0
+ return y
+
+
+def CRAM48(A, n0, dt):
+ """ Chebyshev Rational Approximation Method, order 48
+
+ Algorithm is the 48th order Chebyshev Rational Approximation Method,
+ implemented in the more stable incomplete partial fraction (IPF) form
+ [cram48]_.
+
+ .. [cram48]
+ Pusa, Maria. "Higher-Order Chebyshev Rational Approximation Method and
+ Application to Burnup Equations." Nuclear Science and Engineering 182.3
+ (2016).
+
+ Parameters
+ ----------
+ A : scipy.linalg.csr_matrix
+ Matrix to take exponent of.
+ n0 : numpy.array
+ Vector to operate a matrix exponent on.
+ dt : float
+ Time to integrate to.
+
+ Returns
+ -------
+ numpy.array
+ Results of the matrix exponent.
+ """
+
+ theta_r = np.array([-4.465731934165702e+1, -5.284616241568964e+0,
+ -8.867715667624458e+0, +3.493013124279215e+0,
+ +1.564102508858634e+1, +1.742097597385893e+1,
+ -2.834466755180654e+1, +1.661569367939544e+1,
+ +8.011836167974721e+0, -2.056267541998229e+0,
+ +1.449208170441839e+1, +1.853807176907916e+1,
+ +9.932562704505182e+0, -2.244223871767187e+1,
+ +8.590014121680897e-1, -1.286192925744479e+1,
+ +1.164596909542055e+1, +1.806076684783089e+1,
+ +5.870672154659249e+0, -3.542938819659747e+1,
+ +1.901323489060250e+1, +1.885508331552577e+1,
+ -1.734689708174982e+1, +1.316284237125190e+1])
+ theta_i = np.array([+6.233225190695437e+1, +4.057499381311059e+1,
+ +4.325515754166724e+1, +3.281615453173585e+1,
+ +1.558061616372237e+1, +1.076629305714420e+1,
+ +5.492841024648724e+1, +1.316994930024688e+1,
+ +2.780232111309410e+1, +3.794824788914354e+1,
+ +1.799988210051809e+1, +5.974332563100539e+0,
+ +2.532823409972962e+1, +5.179633600312162e+1,
+ +3.536456194294350e+1, +4.600304902833652e+1,
+ +2.287153304140217e+1, +8.368200580099821e+0,
+ +3.029700159040121e+1, +5.834381701800013e+1,
+ +1.194282058271408e+0, +3.583428564427879e+0,
+ +4.883941101108207e+1, +2.042951874827759e+1])
+ theta = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
+
+ alpha_r = np.array([+6.387380733878774e+2, +1.909896179065730e+2,
+ +4.236195226571914e+2, +4.645770595258726e+2,
+ +7.765163276752433e+2, +1.907115136768522e+3,
+ +2.909892685603256e+3, +1.944772206620450e+2,
+ +1.382799786972332e+5, +5.628442079602433e+3,
+ +2.151681283794220e+2, +1.324720240514420e+3,
+ +1.617548476343347e+4, +1.112729040439685e+2,
+ +1.074624783191125e+2, +8.835727765158191e+1,
+ +9.354078136054179e+1, +9.418142823531573e+1,
+ +1.040012390717851e+2, +6.861882624343235e+1,
+ +8.766654491283722e+1, +1.056007619389650e+2,
+ +7.738987569039419e+1, +1.041366366475571e+2])
+ alpha_i = np.array([-6.743912502859256e+2, -3.973203432721332e+2,
+ -2.041233768918671e+3, -1.652917287299683e+3,
+ -1.783617639907328e+4, -5.887068595142284e+4,
+ -9.953255345514560e+3, -1.427131226068449e+3,
+ -3.256885197214938e+6, -2.924284515884309e+4,
+ -1.121774011188224e+3, -6.370088443140973e+4,
+ -1.008798413156542e+6, -8.837109731680418e+1,
+ -1.457246116408180e+2, -6.388286188419360e+1,
+ -2.195424319460237e+2, -6.719055740098035e+2,
+ -1.693747595553868e+2, -1.177598523430493e+1,
+ -4.596464999363902e+3, -1.738294585524067e+3,
+ -4.311715386228984e+1, -2.777743732451969e+2])
+ alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
+ n = A.shape[0]
+
+ alpha0 = 2.258038182743983e-47
+
+ k = 24
+
+ y = np.array(n0, dtype=np.float64)
+ for l in range(k):
+ y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
+
+ y *= alpha0
+ return y
diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py
new file mode 100644
index 000000000..6c9d538fd
--- /dev/null
+++ b/openmc/deplete/integrator/predictor.py
@@ -0,0 +1,100 @@
+""" The Predictor algorithm."""
+
+import copy
+from itertools import repeat
+import os
+from multiprocessing import Pool
+import time
+
+from .. import comm
+from .cram import CRAM48, cram_wrapper
+from .save_results import save_results
+
+
+def predictor(operator, print_out=True):
+ """The basic predictor integrator.
+
+ Implements the first order predictor algorithm. This algorithm is
+ mathematically defined as:
+
+ .. math::
+ y' &= A(y, t) y(t)
+
+ A_p &= A(y_n, t_n)
+
+ y_{n+1} &= \\text{expm}(A_p h) y_n
+
+ Parameters
+ ----------
+ operator : Operator
+ The operator object to simulate on.
+ print_out : bool, optional
+ Whether or not to print out time.
+ """
+
+ # Save current directory
+ dir_home = os.getcwd()
+
+ # Move to folder
+ os.makedirs(operator.settings.output_dir, exist_ok=True)
+ os.chdir(operator.settings.output_dir)
+
+ # Generate initial conditions
+ vec = operator.initial_condition()
+
+ n_mats = len(vec)
+
+ t = 0.0
+
+ for i, dt in enumerate(operator.settings.dt_vec):
+ # Create vectors
+ x = [copy.deepcopy(vec)]
+ seeds = []
+ eigvls = []
+ rates_array = []
+
+ eigvl, rates, seed = operator.eval(x[0])
+
+ eigvls.append(eigvl)
+ seeds.append(seed)
+ rates_array.append(rates)
+
+ # Create results, write to disk
+ save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i)
+
+ t_start = time.time()
+
+ chains = repeat(operator.chain, n_mats)
+ vecs = (x[0][i] for i in range(n_mats))
+ rates = (rates_array[0][i, :, :] for i in range(n_mats))
+ dts = repeat(dt, n_mats)
+
+ with Pool() as pool:
+ iters = zip(chains, vecs, rates, dts)
+ x_result = list(pool.starmap(cram_wrapper, iters))
+
+ t_end = time.time()
+ if comm.rank == 0:
+ if print_out:
+ print("Time to matexp: ", t_end - t_start)
+
+ t += dt
+ vec = copy.deepcopy(x_result)
+
+ # Perform one last simulation
+ x = [copy.deepcopy(vec)]
+ seeds = []
+ eigvls = []
+ rates_array = []
+ eigvl, rates, seed = operator.eval(x[0])
+
+ eigvls.append(eigvl)
+ seeds.append(seed)
+ rates_array.append(rates)
+
+ # Create results, write to disk
+ save_results(operator, x, rates_array, eigvls, seeds, [t, t],
+ len(operator.settings.dt_vec))
+
+ # Return to origin
+ os.chdir(dir_home)
diff --git a/openmc/deplete/integrator/save_results.py b/openmc/deplete/integrator/save_results.py
new file mode 100644
index 000000000..35cbc7f3f
--- /dev/null
+++ b/openmc/deplete/integrator/save_results.py
@@ -0,0 +1,46 @@
+""" Generic result saving code for integrators.
+
+"""
+from opendeplete.results import Results, write_results
+
+def save_results(op, x, rates, eigvls, seeds, t, step_ind):
+ """ Creates and writes results to disk
+
+ Parameters
+ ----------
+ op : Function
+ The operator used to generate these results.
+ x : list of list of numpy.array
+ The prior x vectors. Indexed [i][cell] using the above equation.
+ rates : list of ReactionRates
+ The reaction rates for each substep.
+ eigvls : list of float
+ Eigenvalue for each substep
+ seeds : list of int
+ Seeds for each substep.
+ t : list of float
+ Time indices.
+ step_ind : int
+ Step index.
+ """
+
+ # Get indexing terms
+ vol_list, nuc_list, burn_list, full_burn_list = op.get_results_info()
+
+ # Create results
+ stages = len(x)
+ results = Results()
+ results.allocate(vol_list, nuc_list, burn_list, full_burn_list, stages)
+
+ n_mat = len(burn_list)
+
+ for i in range(stages):
+ for mat_i in range(n_mat):
+ results[i, mat_i, :] = x[i][mat_i][:]
+
+ results.k = eigvls
+ results.seeds = seeds
+ results.time = t
+ results.rates = rates
+
+ write_results(results, "results.h5", step_ind)
diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
new file mode 100644
index 000000000..1208a9b3c
--- /dev/null
+++ b/openmc/deplete/nuclide.py
@@ -0,0 +1,178 @@
+"""Nuclide module.
+
+Contains the per-nuclide components of a depletion chain.
+"""
+
+from collections import namedtuple
+try:
+ import lxml.etree as ET
+except ImportError:
+ import xml.etree.ElementTree as ET
+
+DecayTuple = namedtuple('DecayTuple', 'type target branching_ratio')
+ReactionTuple = namedtuple('ReactionTuple', 'type target Q branching_ratio')
+
+
+class Nuclide(object):
+ """The Nuclide class.
+
+ Contains everything in a depletion chain relating to a single nuclide.
+
+ Attributes
+ ----------
+ name : str
+ Name of nuclide.
+ half_life : float
+ Half life of nuclide in s^-1.
+ decay_energy : float
+ Energy deposited from decay in eV.
+ n_decay_modes : int
+ Number of decay pathways.
+ decay_modes : list of DecayTuple
+ Decay mode information. Each element of the list is a named tuple with
+ attributes 'type', 'target', and 'branching_ratio'.
+ n_reaction_paths : int
+ Number of possible reaction pathways.
+ reactions : list of ReactionTuple
+ Reaction information. Each element of the list is a named tuple with
+ attribute 'type', 'target', 'Q', and 'branching_ratio'.
+ yield_data : dict of float to list
+ Maps tabulated energy to list of (product, yield) for all
+ neutron-induced fission products.
+ yield_energies : list of float
+ Energies at which fission product yiels exist
+
+ """
+
+ def __init__(self):
+ # Information about the nuclide
+ self.name = None
+ self.half_life = None
+ self.decay_energy = 0.0
+
+ # Decay paths
+ self.decay_modes = []
+
+ # Reaction paths
+ self.reactions = []
+
+ # Neutron fission yields, if present
+ self.yield_data = {}
+ self.yield_energies = []
+
+ @property
+ def n_decay_modes(self):
+ """Number of decay modes."""
+ return len(self.decay_modes)
+
+ @property
+ def n_reaction_paths(self):
+ """Number of possible reaction pathways."""
+ return len(self.reactions)
+
+ @classmethod
+ def xml_read(cls, element):
+ """Read nuclide from an XML element.
+
+ Parameters
+ ----------
+ element : xml.etree.ElementTree.Element
+ XML element to write nuclide data to
+
+ Returns
+ -------
+ nuc : Nuclide
+ Instance of a nuclide
+
+ """
+ nuc = cls()
+ nuc.name = element.get('name')
+
+ # Check for half-life
+ if 'half_life' in element.attrib:
+ nuc.half_life = float(element.get('half_life'))
+ nuc.decay_energy = float(element.get('decay_energy', '0'))
+
+ # Check for decay paths
+ for decay_elem in element.iter('decay_type'):
+ d_type = decay_elem.get('type')
+ target = decay_elem.get('target')
+ branching_ratio = float(decay_elem.get('branching_ratio'))
+ nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))
+
+ # Check for reaction paths
+ for reaction_elem in element.iter('reaction_type'):
+ r_type = reaction_elem.get('type')
+ Q = float(reaction_elem.get('Q', '0'))
+ branching_ratio = float(reaction_elem.get('branching_ratio', '1'))
+
+ # If the type is not fission, get target and Q value, otherwise
+ # just set null values
+ if r_type != 'fission':
+ target = reaction_elem.get('target')
+ else:
+ target = None
+
+ # Append reaction
+ nuc.reactions.append(ReactionTuple(
+ r_type, target, Q, branching_ratio))
+
+ fpy_elem = element.find('neutron_fission_yields')
+ if fpy_elem is not None:
+ for yields_elem in fpy_elem.iter('fission_yields'):
+ E = float(yields_elem.get('energy'))
+ products = yields_elem.find('products').text.split()
+ yields = [float(y) for y in
+ yields_elem.find('data').text.split()]
+ nuc.yield_data[E] = list(zip(products, yields))
+ nuc.yield_energies = list(sorted(nuc.yield_data.keys()))
+
+ return nuc
+
+ def xml_write(self):
+ """Write nuclide to XML element.
+
+ Returns
+ -------
+ elem : xml.etree.ElementTree.Element
+ XML element to write nuclide data to
+
+ """
+ elem = ET.Element('nuclide_table')
+ elem.set('name', self.name)
+
+ if self.half_life is not None:
+ elem.set('half_life', str(self.half_life))
+ elem.set('decay_modes', str(len(self.decay_modes)))
+ elem.set('decay_energy', str(self.decay_energy))
+ for mode, daughter, br in self.decay_modes:
+ mode_elem = ET.SubElement(elem, 'decay_type')
+ mode_elem.set('type', mode)
+ mode_elem.set('target', daughter)
+ mode_elem.set('branching_ratio', str(br))
+
+ elem.set('reactions', str(len(self.reactions)))
+ for rx, daughter, Q, br in self.reactions:
+ rx_elem = ET.SubElement(elem, 'reaction_type')
+ rx_elem.set('type', rx)
+ rx_elem.set('Q', str(Q))
+ if rx != 'fission':
+ rx_elem.set('target', daughter)
+ if br != 1.0:
+ rx_elem.set('branching_ratio', str(br))
+
+ if self.yield_data:
+ fpy_elem = ET.SubElement(elem, 'neutron_fission_yields')
+ energy_elem = ET.SubElement(fpy_elem, 'energies')
+ energy_elem.text = ' '.join(str(E) for E in self.yield_energies)
+
+ for E in self.yield_energies:
+ yields_elem = ET.SubElement(fpy_elem, 'fission_yields')
+ yields_elem.set('energy', str(E))
+
+ products_elem = ET.SubElement(yields_elem, 'products')
+ products_elem.text = ' '.join(x[0] for x in self.yield_data[E])
+ data_elem = ET.SubElement(yields_elem, 'data')
+ data_elem.text = ' '.join(str(x[1]) for x in self.yield_data[E])
+
+ return elem
diff --git a/openmc/deplete/openmc_wrapper.py b/openmc/deplete/openmc_wrapper.py
new file mode 100644
index 000000000..347dc7185
--- /dev/null
+++ b/openmc/deplete/openmc_wrapper.py
@@ -0,0 +1,853 @@
+""" The OpenMC wrapper module.
+
+This module implements the OpenDeplete -> OpenMC linkage.
+"""
+
+import copy
+from collections import OrderedDict
+import os
+import random
+import sys
+import time
+try:
+ import lxml.etree as ET
+ _have_lxml = True
+except ImportError:
+ import xml.etree.ElementTree as ET
+ from openmc.clean_xml import clean_xml_indentation
+ _have_lxml = False
+
+import h5py
+import numpy as np
+import openmc
+import openmc.capi
+
+from . import comm
+from .atom_number import AtomNumber
+from .depletion_chain import DepletionChain
+from .reaction_rates import ReactionRates
+from .function import Settings, Operator
+
+
+_JOULE_PER_EV = 1.6021766208e-19
+
+
+def chunks(items, n):
+ min_size, extra = divmod(len(items), n)
+ j = 0
+ chunk_list = []
+ for i in range(n):
+ chunk_size = min_size + int(i < extra)
+ chunk_list.append(items[j:j + chunk_size])
+ j += chunk_size
+ return chunk_list
+
+
+class OpenMCSettings(Settings):
+ """The OpenMCSettings class.
+
+ Extends Settings to provide information OpenMC needs to run.
+
+ Attributes
+ ----------
+ dt_vec : numpy.array
+ Array of time steps to take. (From Settings)
+ tol : float
+ Tolerance for adaptive time stepping. (From Settings)
+ output_dir : str
+ Path to output directory to save results. (From Settings)
+ chain_file : str
+ Path to the depletion chain xml file. Defaults to the environment
+ variable "OPENDEPLETE_CHAIN" if it exists.
+ openmc_call : str
+ OpenMC executable path. Defaults to "openmc".
+ particles : int
+ Number of particles to simulate per batch.
+ batches : int
+ Number of batches.
+ inactive : int
+ Number of inactive batches.
+ lower_left : list of float
+ Coordinate of lower left of bounding box of geometry.
+ upper_right : list of float
+ Coordinate of upper right of bounding box of geometry.
+ entropy_dimension : list of int
+ Grid size of entropy.
+ dilute_initial : float, default 1.0e3
+ Initial atom density 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.
+ 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.
+ constant_seed : int
+ If present, all runs will be performed with this seed.
+ power : float
+ Power of the reactor in W. For a 2D problem, the power can be given in
+ W/cm as long as the "volume" assigned to a depletion material is
+ actually an area in cm^2.
+ """
+
+ def __init__(self):
+ super().__init__()
+ # OpenMC specific
+ try:
+ self.chain_file = os.environ["OPENDEPLETE_CHAIN"]
+ except KeyError:
+ self.chain_file = None
+ self.openmc_call = "openmc"
+ self.particles = None
+ self.batches = None
+ self.inactive = None
+ self.lower_left = None
+ self.upper_right = None
+ self.entropy_dimension = None
+ self.dilute_initial = 1.0e3
+
+ # OpenMC testing specific
+ self.round_number = False
+ self.constant_seed = None
+
+ # Depletion problem specific
+ self.power = None
+
+
+class Materials(object):
+ """The Materials class.
+
+ Contains information about cross sections for a cell.
+
+ Attributes
+ ----------
+ temperature : float
+ Temperature in Kelvin for each region.
+ sab : str or list of str
+ ENDF S(a,b) name for a region that needs S(a,b) data. Not set if no
+ S(a,b) needed for region.
+ """
+
+ def __init__(self):
+ self.temperature = None
+ self.sab = None
+
+
+class OpenMCOperator(Operator):
+ """The OpenMC Operator class.
+
+ Provides Operator functions for OpenMC.
+
+ Parameters
+ ----------
+ geometry : openmc.Geometry
+ The OpenMC geometry object.
+ settings : OpenMCSettings
+ Settings object.
+
+ Attributes
+ ----------
+ settings : OpenMCSettings
+ Settings object. (From Operator)
+ geometry : openmc.Geometry
+ The OpenMC geometry object.
+ materials : list of Materials
+ Materials to be used for this simulation.
+ seed : int
+ The RNG seed used in last OpenMC run.
+ number : AtomNumber
+ Total number of atoms in simulation.
+ participating_nuclides : set of str
+ A set listing all unique nuclides available from cross_sections.xml.
+ chain : DepletionChain
+ The depletion chain information necessary to form matrices and tallies.
+ reaction_rates : ReactionRates
+ Reaction rates from the last operator step.
+ power : OrderedDict of str to float
+ Material-by-Material power. Indexed by material ID.
+ mat_name : OrderedDict of str to int
+ The name of region each material is set to. Indexed by material ID.
+ burn_mat_to_id : OrderedDict of str to int
+ Dictionary mapping material ID (as a string) to an index in reaction_rates.
+ burn_nuc_to_id : OrderedDict of str to int
+ Dictionary mapping nuclide name (as a string) to an index in
+ reaction_rates.
+ n_nuc : int
+ Number of nuclides considered in the decay chain.
+ mat_tally_ind : OrderedDict of str to int
+ Dictionary mapping material ID to index in tally.
+ """
+
+ def __init__(self, geometry, settings):
+ super().__init__(settings)
+
+ self.geometry = geometry
+ self.seed = 0
+ self.number = None
+ self.participating_nuclides = None
+ self.reaction_rates = None
+ self.power = None
+ self.mat_name = OrderedDict()
+ self.burn_mat_to_ind = OrderedDict()
+ self.burn_nuc_to_ind = None
+
+ # Read depletion chain
+ self.chain = DepletionChain.xml_read(settings.chain_file)
+
+ # Clear out OpenMC, create task lists, distribute
+ if comm.rank == 0:
+ clean_up_openmc()
+ mat_burn_list, mat_not_burn_list, volume, self.mat_tally_ind, \
+ nuc_dict = self.extract_mat_ids()
+ else:
+ # Dummy variables
+ mat_burn_list = None
+ mat_not_burn_list = None
+ volume = None
+ nuc_dict = None
+ self.mat_tally_ind = None
+
+ mat_burn = comm.scatter(mat_burn_list)
+ mat_not_burn = comm.scatter(mat_not_burn_list)
+ nuc_dict = comm.bcast(nuc_dict)
+ volume = comm.bcast(volume)
+ self.mat_tally_ind = comm.bcast(self.mat_tally_ind)
+
+ # Load participating nuclides
+ self.load_participating()
+
+ # Extract number densities from the geometry
+ self.extract_number(mat_burn, mat_not_burn, volume, nuc_dict)
+
+ # Create reaction rate tables
+ self.initialize_reaction_rates()
+
+ def __del__(self):
+ openmc.capi.finalize()
+
+ def extract_mat_ids(self):
+ """ Extracts materials and assigns them to processes.
+
+ Returns
+ -------
+ mat_burn_lists : list of list of int
+ List of burnable materials indexed by rank.
+ mat_not_burn_lists : list of list of int
+ List of non-burnable materials indexed by rank.
+ volume : OrderedDict of str to float
+ Volume of each cell
+ mat_tally_ind : OrderedDict of str to int
+ Dictionary mapping material ID to index in tally.
+ nuc_dict : OrderedDict of str to int
+ Nuclides in order of how they'll appear in the simulation.
+ """
+
+ mat_burn = set()
+ mat_not_burn = set()
+ nuc_set = set()
+
+ volume = OrderedDict()
+
+ # Iterate once through the geometry to get dictionaries
+ cells = self.geometry.get_all_material_cells()
+ for cell in cells.values():
+ name = cell.name
+
+ if isinstance(cell.fill, openmc.Material):
+ mat = cell.fill
+ for nuclide in mat.get_nuclide_densities():
+ nuc_set.add(nuclide)
+ if mat.depletable:
+ mat_burn.add(str(mat.id))
+ volume[str(mat.id)] = mat.volume
+ else:
+ mat_not_burn.add(str(mat.id))
+ self.mat_name[mat.id] = name
+ else:
+ for mat in cell.fill:
+ for nuclide in mat.get_nuclide_densities():
+ nuc_set.add(nuclide)
+ if mat.depletable:
+ mat_burn.add(str(mat.id))
+ volume[str(mat.id)] = mat.volume
+ else:
+ mat_not_burn.add(str(mat.id))
+ self.mat_name[mat.id] = name
+
+ need_vol = []
+
+ for mat_id in volume:
+ if volume[mat_id] is None:
+ need_vol.append(mat_id)
+
+ if need_vol:
+ exit("Need volumes for materials: " + str(need_vol))
+
+ # Sort the sets
+ mat_burn = sorted(mat_burn, key=int)
+ mat_not_burn = sorted(mat_not_burn, key=int)
+ nuc_set = sorted(nuc_set)
+
+ # Construct a global nuclide dictionary, burned first
+ nuc_dict = copy.deepcopy(self.chain.nuclide_dict)
+
+ i = len(nuc_dict)
+
+ for nuc in nuc_set:
+ if nuc not in nuc_dict:
+ nuc_dict[nuc] = i
+ i += 1
+
+ # Decompose geometry
+ mat_burn_lists = chunks(mat_burn, comm.size)
+ mat_not_burn_lists = chunks(mat_not_burn, comm.size)
+
+ mat_tally_ind = OrderedDict()
+
+ for i, mat in enumerate(mat_burn):
+ mat_tally_ind[mat] = i
+
+ return mat_burn_lists, mat_not_burn_lists, volume, mat_tally_ind, nuc_dict
+
+ def extract_number(self, mat_burn, mat_not_burn, volume, nuc_dict):
+ """ Construct self.number read from geometry
+
+ Parameters
+ ----------
+ mat_burn : list of int
+ Materials to be burned managed by this thread.
+ mat_not_burn
+ Materials not to be burned managed by this thread.
+ volume : OrderedDict of str to float
+ Volumes for the above materials.
+ nuc_dict : OrderedDict of str to int
+ Nuclides to be used in the simulation.
+ """
+
+ # Same with materials
+ mat_dict = OrderedDict()
+ self.burn_mat_to_ind = OrderedDict()
+ i = 0
+ for mat in mat_burn:
+ mat_dict[mat] = i
+ self.burn_mat_to_ind[mat] = i
+ i += 1
+
+ for mat in mat_not_burn:
+ mat_dict[mat] = i
+ i += 1
+
+ n_mat_burn = len(mat_burn)
+ n_nuc_burn = len(self.chain.nuclide_dict)
+
+ self.number = AtomNumber(mat_dict, nuc_dict, volume, n_mat_burn, n_nuc_burn)
+
+ if self.settings.dilute_initial != 0.0:
+ for nuc in self.burn_nuc_to_ind:
+ self.number.set_atom_density(np.s_[:], nuc, self.settings.dilute_initial)
+
+ # Now extract the number densities and store
+ cells = self.geometry.get_all_material_cells()
+ for cell in cells.values():
+ if isinstance(cell.fill, openmc.Material):
+ if str(cell.fill.id) in mat_dict:
+ self.set_number_from_mat(cell.fill)
+ else:
+ for mat in cell.fill:
+ if str(mat.id) in mat_dict:
+ self.set_number_from_mat(mat)
+
+ def set_number_from_mat(self, mat):
+ """ Extracts material and number densities from openmc.Material
+
+ Parameters
+ ----------
+ mat : openmc.Materials
+ The material to read from
+ """
+
+ mat_id = str(mat.id)
+ mat_ind = self.number.mat_to_ind[mat_id]
+
+ nuc_dens = mat.get_nuclide_atom_densities()
+ for nuclide in nuc_dens:
+ name = nuclide.name
+ number = nuc_dens[nuclide][1] * 1.0e24
+ self.number.set_atom_density(mat_id, name, number)
+
+ def initialize_reaction_rates(self):
+ """ Create reaction rates object. """
+ self.reaction_rates = ReactionRates(
+ self.burn_mat_to_ind,
+ self.burn_nuc_to_ind,
+ self.chain.react_to_ind)
+
+ self.chain.nuc_to_react_ind = self.burn_nuc_to_ind
+
+ def eval(self, vec, print_out=True):
+ """ Runs a simulation.
+
+ Parameters
+ ----------
+ vec : list of numpy.array
+ Total atoms to be used in function.
+ print_out : bool, optional
+ Whether or not to print out time.
+
+ Returns
+ -------
+ mat : list of scipy.sparse.csr_matrix
+ Matrices for the next step.
+ k : float
+ Eigenvalue of the problem.
+ rates : ReactionRates
+ Reaction rates from this simulation.
+ seed : int
+ Seed for this simulation.
+ """
+
+ # Prevent OpenMC from complaining about re-creating tallies
+ clean_up_openmc()
+
+ # Update status
+ self.set_density(vec)
+
+ time_start = time.time()
+
+ # Update material compositions and tally nuclides
+ self._update_materials()
+ openmc.capi.tallies[1].nuclides = self._get_tally_nuclides()
+
+ # Run OpenMC
+ openmc.capi.reset()
+ openmc.capi.run()
+
+ time_openmc = time.time()
+
+ # Extract results
+ k = self.unpack_tallies_and_normalize()
+
+ if comm.rank == 0:
+ time_unpack = time.time()
+
+ if print_out:
+ print("Time to openmc: ", time_openmc - time_start)
+ print("Time to unpack: ", time_unpack - time_openmc)
+
+ return k, copy.deepcopy(self.reaction_rates), self.seed
+
+ def form_matrix(self, y, mat):
+ """ Forms the depletion matrix.
+
+ Parameters
+ ----------
+ y : numpy.ndarray
+ An array representing reaction rates for this cell.
+ mat : int
+ Material id.
+
+ Returns
+ -------
+ scipy.sparse.csr_matrix
+ Sparse matrix representing the depletion matrix.
+ """
+
+ return copy.deepcopy(self.chain.form_matrix(y[mat, :, :]))
+
+ def initial_condition(self):
+ """ Performs final setup and returns initial condition.
+
+ Returns
+ -------
+ list of numpy.array
+ Total density for initial conditions.
+ """
+
+ # Create XML files
+ if comm.rank == 0:
+ self.geometry.export_to_xml()
+ self.generate_settings_xml()
+ self.generate_materials_xml()
+
+ # Initialize OpenMC library
+ comm.barrier()
+ openmc.capi.init(comm)
+
+ # Generate tallies in memory
+ self.generate_tallies()
+
+ # Return number density vector
+ return self.total_density_list()
+
+ def _update_materials(self):
+ """Updates material compositions in OpenMC on all processes."""
+
+ for rank in range(comm.size):
+ number_i = comm.bcast(self.number, root=rank)
+
+ for mat in number_i.mat_to_ind:
+ nuclides = []
+ densities = []
+ for nuc in number_i.nuc_to_ind:
+ if nuc in self.participating_nuclides:
+ val = 1.0e-24 * number_i.get_atom_density(mat, nuc)
+
+ # If nuclide is zero, do not add to the problem.
+ if val > 0.0:
+ if self.settings.round_number:
+ val_magnitude = np.floor(np.log10(val))
+ val_scaled = val / 10**val_magnitude
+ val_round = round(val_scaled, 8)
+
+ val = val_round * 10**val_magnitude
+
+ nuclides.append(nuc)
+ densities.append(val)
+ else:
+ # Only output warnings if values are significantly
+ # 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)")
+ number_i[mat, nuc] = 0.0
+
+ mat_internal = openmc.capi.materials[int(mat)]
+ mat_internal.set_densities(nuclides, densities)
+
+ 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.
+ """
+
+ materials = openmc.Materials(self.geometry.get_all_materials()
+ .values())
+
+ # Sort nuclides according to order in AtomNumber object
+ nuclides = list(self.number.nuc_to_ind.keys())
+ for mat in materials:
+ mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
+
+ materials.export_to_xml()
+
+ def generate_settings_xml(self):
+ """ Generates settings.xml.
+
+ This function creates settings.xml using the value of the settings
+ variable.
+
+ Todo
+ ----
+ Rewrite to generalize source box.
+ """
+
+ batches = self.settings.batches
+ inactive = self.settings.inactive
+ particles = self.settings.particles
+
+ # Just a generic settings file to get it running.
+ settings_file = openmc.Settings()
+ settings_file.batches = batches
+ settings_file.inactive = inactive
+ settings_file.particles = particles
+ settings_file.source = openmc.Source(space=openmc.stats.Box(
+ self.settings.lower_left, self.settings.upper_right))
+
+ if self.settings.entropy_dimension is not None:
+ entropy_mesh = openmc.Mesh()
+ entropy_mesh.lower_left = self.settings.lower_left
+ entropy_mesh.upper_right = self.settings.upper_right
+ entropy_mesh.dimension = self.settings.entropy_dimension
+ settings_file.entropy_mesh = entropy_mesh
+
+ # Set seed
+ if self.settings.constant_seed is not None:
+ seed = self.settings.constant_seed
+ else:
+ seed = random.randint(1, sys.maxsize-1)
+
+ settings_file.seed = self.seed = seed
+
+ settings_file.export_to_xml()
+
+ def _get_tally_nuclides(self):
+ nuc_set = set()
+
+ # Create the set of all nuclides in the decay chain in cells marked for
+ # burning in which the number density is greater than zero.
+ for nuc in self.number.nuc_to_ind:
+ if nuc in self.participating_nuclides:
+ 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.nuc_to_ind
+ if nuc in nuc_set]
+ else:
+ nuc_list = None
+
+ # Store list of tally nuclides on each process
+ nuc_list = comm.bcast(nuc_list, root=0)
+ tally_nuclides = [nuc for nuc in nuc_list
+ if nuc in self.chain.nuclide_dict]
+
+ return tally_nuclides
+
+ def generate_tallies(self):
+ """Generates depletion tallies.
+
+ Using information from self.depletion_chain as well as the nuclides
+ currently in the problem, this function automatically generates a
+ tally.xml for the simulation.
+ """
+
+ # Create tallies for depleting regions
+ materials = [openmc.capi.materials[int(i)]
+ for i in self.mat_tally_ind]
+ mat_filter = openmc.capi.MaterialFilter(materials, 1)
+
+ # Set up a tally that has a material filter covering each depletable
+ # material and scores corresponding to all reactions that cause
+ # transmutation. The nuclides for the tally are set later when eval() is
+ # called.
+ tally_dep = openmc.capi.Tally(1)
+ tally_dep.scores = self.chain.react_to_ind.keys()
+ tally_dep.filters = [mat_filter]
+
+ def total_density_list(self):
+ """ Returns a list of total density lists.
+
+ This list is in the exact same order as depletion_matrix_list, so that
+ matrix exponentiation can be done easily.
+
+ Returns
+ -------
+ list of numpy.array
+ A list of np.arrays containing total atoms of each cell.
+ """
+
+ total_density = [self.number.get_mat_slice(i) for i in range(self.number.n_mat_burn)]
+
+ return total_density
+
+ def set_density(self, total_density):
+ """ Sets density.
+
+ Sets the density in the exact same order as total_density_list outputs,
+ allowing for internal consistency
+
+ Parameters
+ ----------
+ total_density : list of numpy.array
+ Total atoms.
+ """
+
+ # Fill in values
+ for i in range(self.number.n_mat_burn):
+ self.number.set_mat_slice(i, total_density[i])
+
+ def unpack_tallies_and_normalize(self):
+ """ Unpack tallies from OpenMC
+
+ This function reads the tallies generated by OpenMC (from the tally.xml
+ file generated in generate_tally_xml) normalizes them so that the total
+ power generated is new_power, and then stores them in the reaction rate
+ database.
+
+ Returns
+ -------
+ k : float
+ Eigenvalue of the last simulation.
+
+ Todo
+ ----
+ Provide units for power
+ """
+
+ rates = self.reaction_rates
+ rates[:, :, :] = 0.0
+
+ k_combined = openmc.capi.keff()[0]
+
+ # Extract tally bins
+ materials = list(self.mat_tally_ind.keys())
+ nuclides = openmc.capi.tallies[1].nuclides
+ reactions = list(self.chain.react_to_ind.keys())
+
+ # Form fast map
+ nuc_ind = [rates.nuc_to_ind[nuc] for nuc in nuclides]
+ react_ind = [rates.react_to_ind[react] for react in reactions]
+
+ # Compute fission power
+ # TODO : improve this calculation
+
+ # Keep track of energy produced from all reactions in eV per source
+ # particle
+ energy = 0.0
+
+ # Create arrays to store fission Q values, reaction rates, and nuclide
+ # numbers
+ fission_Q = np.zeros(rates.n_nuc)
+ rates_expanded = np.zeros((rates.n_nuc, rates.n_react))
+ number = np.zeros(rates.n_nuc)
+
+ fission_ind = rates.react_to_ind["fission"]
+
+ for nuclide in self.chain.nuclides:
+ if nuclide.name in rates.nuc_to_ind:
+ for rx in nuclide.reactions:
+ if rx.type == 'fission':
+ ind = rates.nuc_to_ind[nuclide.name]
+ fission_Q[ind] = rx.Q
+ break
+
+ # Extract results
+ for i, mat in enumerate(self.number.burn_mat_list):
+ # Get tally index
+ slab = materials.index(mat)
+
+ # Get material results hyperslab
+ results = openmc.capi.tallies[1].results[slab, :, 1]
+
+ # Zero out reaction rates and nuclide numbers
+ rates_expanded[:] = 0.0
+ number[:] = 0.0
+
+ # Expand into our memory layout
+ j = 0
+ for nuc, i_nuc_results in zip(nuclides, nuc_ind):
+ number[i_nuc_results] = self.number[mat, nuc]
+ for react in react_ind:
+ rates_expanded[i_nuc_results, react] = results[j]
+ j += 1
+
+ # Accumulate energy from fission
+ energy += np.dot(rates_expanded[:, fission_ind], fission_Q)
+
+ # Divide by total number and store
+ for i_nuc_results in nuc_ind:
+ if number[i_nuc_results] != 0.0:
+ for react in react_ind:
+ rates_expanded[i_nuc_results, react] /= number[i_nuc_results]
+
+ rates.rates[i, :, :] = rates_expanded
+
+ # Reduce energy produced from all processes
+ energy = comm.allreduce(energy)
+
+ # Determine power in eV/s
+ power = self.settings.power / _JOULE_PER_EV
+
+ # Scale reaction rates to obtain units of reactions/sec
+ rates[:, :, :] *= power / energy
+
+ return k_combined
+
+ def load_participating(self):
+ """ Loads a cross_sections.xml file to find participating nuclides.
+
+ This allows for nuclides that are important in the decay chain but not
+ important neutronically, or have no cross section data.
+ """
+
+ # Reads cross_sections.xml to create a dictionary containing
+ # participating (burning and not just decaying) nuclides.
+
+ try:
+ filename = os.environ["OPENMC_CROSS_SECTIONS"]
+ except KeyError:
+ filename = None
+
+ self.participating_nuclides = set()
+
+ try:
+ tree = ET.parse(filename)
+ except:
+ if filename is None:
+ msg = "No cross_sections.xml specified in materials."
+ else:
+ msg = 'Cross section file "{}" is invalid.'.format(filename)
+ raise IOError(msg)
+
+ root = tree.getroot()
+ self.burn_nuc_to_ind = OrderedDict()
+ nuc_ind = 0
+
+ for nuclide_node in root.findall('library'):
+ mats = nuclide_node.get('materials')
+ if not mats:
+ continue
+ for name in mats.split():
+ # Make a burn list of the union of nuclides in cross_sections.xml
+ # and nuclides in depletion chain.
+ if name not in self.participating_nuclides:
+ self.participating_nuclides.add(name)
+ if name in self.chain.nuclide_dict:
+ self.burn_nuc_to_ind[name] = nuc_ind
+ nuc_ind += 1
+
+ @property
+ def n_nuc(self):
+ """Number of nuclides considered in the decay chain."""
+ return len(self.chain.nuclides)
+
+ def get_results_info(self):
+ """ Returns volume list, cell 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 cell IDs to be burned. Used for sorting the simulation.
+ full_burn_dict : OrderedDict of str to int
+ Maps cell name to index in global geometry.
+ """
+
+ nuc_list = self.number.burn_nuc_list
+ burn_list = self.number.burn_mat_list
+
+ 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.mat_tally_ind
+
+def density_to_mat(dens_dict):
+ """ Generates an OpenMC material from a cell ID and self.number_density.
+ Parameters
+ ----------
+ m_id : int
+ Cell ID.
+ Returns
+ -------
+ openmc.Material
+ The OpenMC material filled with nuclides.
+ """
+
+ mat = openmc.Material()
+ for key in dens_dict:
+ mat.add_nuclide(key, 1.0e-24*dens_dict[key])
+ mat.set_density('sum')
+
+ return mat
+
+def clean_up_openmc():
+ """ Resets all automatic indexing in OpenMC, as these get in the way. """
+ openmc.reset_auto_ids()
diff --git a/openmc/deplete/reaction_rates.py b/openmc/deplete/reaction_rates.py
new file mode 100644
index 000000000..7b934027a
--- /dev/null
+++ b/openmc/deplete/reaction_rates.py
@@ -0,0 +1,113 @@
+"""ReactionRates module.
+
+An ndarray to store reaction rates with string, integer, or slice indexing.
+"""
+
+import numpy as np
+
+
+class ReactionRates(object):
+ """ ReactionRates class.
+
+ An ndarray to store reaction rates with string, integer, or slice indexing.
+
+ Parameters
+ ----------
+ mat_to_ind : OrderedDict of str to int
+ A dictionary mapping material ID as string to index.
+ nuc_to_ind : OrderedDict of str to int
+ A dictionary mapping nuclide name as string to index.
+ react_to_ind : OrderedDict of str to int
+ A dictionary mapping reaction name as string to index.
+
+ Attributes
+ ----------
+ mat_to_ind : OrderedDict of str to int
+ A dictionary mapping cell ID as string to index.
+ nuc_to_ind : OrderedDict of str to int
+ A dictionary mapping nuclide name as string to index.
+ react_to_ind : OrderedDict of str to int
+ A dictionary mapping reaction name as string to index.
+ n_mat : int
+ Number of materials.
+ n_nuc : int
+ Number of nucs.
+ n_react : int
+ Number of reactions.
+ rates : numpy.array
+ Array storing rates indexed by the above dictionaries.
+ """
+
+ def __init__(self, mat_to_ind, nuc_to_ind, react_to_ind):
+
+ self.mat_to_ind = mat_to_ind
+ self.nuc_to_ind = nuc_to_ind
+ self.react_to_ind = react_to_ind
+
+ self.rates = np.zeros((self.n_mat, self.n_nuc, self.n_react))
+
+ def __getitem__(self, pos):
+ """ Retrieves an item from reaction_rates.
+
+ Parameters
+ ----------
+ pos : tuple
+ A three-length tuple containing a material index, a nuc index, and a
+ reaction index. These indexes can be strings (which get converted
+ to integers via the dictionaries), integers used directly, or
+ slices.
+
+ Returns
+ -------
+ numpy.array
+ The value indexed from self.rates.
+ """
+
+ mat, nuc, react = pos
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+ if isinstance(nuc, str):
+ nuc = self.nuc_to_ind[nuc]
+ if isinstance(react, str):
+ react = self.react_to_ind[react]
+
+ return self.rates[mat, nuc, react]
+
+ def __setitem__(self, pos, val):
+ """ Sets an item from reaction_rates.
+
+ Parameters
+ ----------
+ pos : tuple
+ A three-length tuple containing a material index, a nuc index, and a
+ reaction index. These indexes can be strings (which get converted
+ to integers via the dictionaries), integers used directly, or
+ slices.
+ val : float
+ The value to set the array to.
+ """
+
+ mat, nuc, react = pos
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+ if isinstance(nuc, str):
+ nuc = self.nuc_to_ind[nuc]
+ if isinstance(react, str):
+ react = self.react_to_ind[react]
+
+ self.rates[mat, nuc, react] = val
+
+ @property
+ def n_mat(self):
+ """Number of cells."""
+ return len(self.mat_to_ind)
+
+ @property
+ def n_nuc(self):
+ """Number of nucs."""
+ return len(self.nuc_to_ind)
+
+ @property
+ def n_react(self):
+ """Number of reactions."""
+ return len(self.react_to_ind)
diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py
new file mode 100644
index 000000000..c0ec1627e
--- /dev/null
+++ b/openmc/deplete/results.py
@@ -0,0 +1,454 @@
+""" The results module.
+
+Contains results generation and saving capabilities.
+"""
+
+from collections import OrderedDict
+import copy
+
+import numpy as np
+import h5py
+
+from . import comm, have_mpi
+from .reaction_rates import ReactionRates
+
+RESULTS_VERSION = 2
+
+class Results(object):
+ """ Contains output of opendeplete.
+
+ Attributes
+ ----------
+ k : list of float
+ Eigenvalue for each substep.
+ seeds : list of int
+ Seeds for each substep.
+ time : list of float
+ Time at beginning, end of step, in seconds.
+ n_mat : int
+ Number of mats.
+ n_nuc : int
+ Number of nuclides.
+ rates : list of ReactionRates
+ The reaction rates for each substep.
+ volume : OrderedDict of int to float
+ Dictionary mapping mat id to volume.
+ mat_to_ind : OrderedDict of str to int
+ A dictionary mapping mat ID as string to index.
+ nuc_to_ind : OrderedDict of str to int
+ A dictionary mapping nuclide name as string to index.
+ mat_to_hdf5_ind : OrderedDict of str to int
+ A dictionary mapping mat ID as string to global index.
+ n_hdf5_mats : int
+ Number of materials in entire geometry.
+ n_stages : int
+ Number of stages in simulation.
+ data : numpy.array
+ Atom quantity, stored by stage, mat, then by nuclide.
+ """
+
+ def __init__(self):
+ self.k = None
+ self.seeds = None
+ self.time = None
+ self.p_terms = None
+ self.rates = None
+ self.volume = None
+
+ self.mat_to_ind = None
+ self.nuc_to_ind = None
+ self.mat_to_hdf5_ind = None
+
+ self.data = None
+
+ def allocate(self, volume, nuc_list, burn_list, full_burn_dict, stages):
+ """ Allocates memory of Results.
+
+ Parameters
+ ----------
+ 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 mat IDs to be burned. Used for sorting the simulation.
+ full_burn_dict : dict of str to int
+ Map of material name to id in global geometry.
+ stages : int
+ Number of stages in simulation.
+ """
+
+ self.volume = copy.deepcopy(volume)
+ self.nuc_to_ind = OrderedDict()
+ self.mat_to_ind = OrderedDict()
+ self.mat_to_hdf5_ind = copy.deepcopy(full_burn_dict)
+
+ for i, mat in enumerate(burn_list):
+ self.mat_to_ind[mat] = i
+
+ for i, nuc in enumerate(nuc_list):
+ self.nuc_to_ind[nuc] = i
+
+ # Create storage array
+ self.data = np.zeros((stages, self.n_mat, self.n_nuc))
+
+ @property
+ def n_mat(self):
+ """Number of mats."""
+ return len(self.mat_to_ind)
+
+ @property
+ def n_nuc(self):
+ """Number of nuclides."""
+ return len(self.nuc_to_ind)
+
+ @property
+ def n_hdf5_mats(self):
+ """Number of materials in entire geometry."""
+ return len(self.mat_to_hdf5_ind)
+
+ @property
+ def n_stages(self):
+ """Number of stages in simulation."""
+ return self.data.shape[0]
+
+ def __getitem__(self, pos):
+ """ Retrieves an item from results.
+
+ Parameters
+ ----------
+ pos : tuple
+ A three-length tuple containing a stage index, mat index and a nuc
+ index. All can be integers or slices. The second two can be
+ strings corresponding to their respective dictionary.
+
+ Returns
+ -------
+ float
+ The atoms for stage, mat, nuc
+ """
+
+ stage, mat, nuc = pos
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+ if isinstance(nuc, str):
+ nuc = self.nuc_to_ind[nuc]
+
+ return self.data[stage, mat, nuc]
+
+ def __setitem__(self, pos, val):
+ """ Sets an item from results.
+
+ Parameters
+ ----------
+ pos : tuple
+ A three-length tuple containing a stage index, mat index and a nuc
+ index. All can be integers or slices. The second two can be
+ strings corresponding to their respective dictionary.
+
+ val : float
+ The value to set data to.
+ """
+
+ stage, mat, nuc = pos
+ if isinstance(mat, str):
+ mat = self.mat_to_ind[mat]
+ if isinstance(nuc, str):
+ nuc = self.nuc_to_ind[nuc]
+
+ self.data[stage, mat, nuc] = val
+
+ def create_hdf5(self, handle):
+ """ Creates file structure for a blank HDF5 file.
+
+ Parameters
+ ----------
+ handle : h5py.File or h5py.Group
+ An hdf5 file or group type to store this in.
+ """
+
+ # Create and save the 5 dictionaries:
+ # quantities
+ # self.mat_to_ind -> self.volume (TODO: support for changing volumes)
+ # self.nuc_to_ind
+ # reactions
+ # self.rates[0].nuc_to_ind (can be different from above, above is superset)
+ # self.rates[0].react_to_ind
+ # these are shared by every step of the simulation, and should be deduplicated.
+
+ # Store concentration mat and nuclide dictionaries (along with volumes)
+
+ handle.create_dataset("version", data=RESULTS_VERSION)
+
+ mat_int = sorted([int(mat) for mat in self.mat_to_hdf5_ind])
+ mat_list = [str(mat) for mat in mat_int]
+ nuc_list = sorted(self.nuc_to_ind.keys())
+ rxn_list = sorted(self.rates[0].react_to_ind.keys())
+
+ n_mats = self.n_hdf5_mats
+ n_nuc_number = len(nuc_list)
+ n_nuc_rxn = len(self.rates[0].nuc_to_ind)
+ n_rxn = len(rxn_list)
+ n_stages = self.n_stages
+
+ mat_group = handle.create_group("cells")
+
+ for mat in mat_list:
+ mat_single_group = mat_group.create_group(mat)
+ mat_single_group.attrs["index"] = self.mat_to_hdf5_ind[mat]
+ mat_single_group.attrs["volume"] = self.volume[mat]
+
+ nuc_group = handle.create_group("nuclides")
+
+ for nuc in nuc_list:
+ nuc_single_group = nuc_group.create_group(nuc)
+ nuc_single_group.attrs["atom number index"] = self.nuc_to_ind[nuc]
+ if nuc in self.rates[0].nuc_to_ind:
+ nuc_single_group.attrs["reaction rate index"] = self.rates[0].nuc_to_ind[nuc]
+
+ rxn_group = handle.create_group("reactions")
+
+ for rxn in rxn_list:
+ rxn_single_group = rxn_group.create_group(rxn)
+ rxn_single_group.attrs["index"] = self.rates[0].react_to_ind[rxn]
+
+ # Construct array storage
+
+ handle.create_dataset("number", (1, n_stages, n_mats, n_nuc_number),
+ maxshape=(None, n_stages, n_mats, n_nuc_number),
+ chunks=(1, 1, n_mats, n_nuc_number),
+ dtype='float64')
+
+ handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
+ maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
+ chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
+ dtype='float64')
+
+ handle.create_dataset("eigenvalues", (1, n_stages),
+ maxshape=(None, n_stages), dtype='float64')
+
+ handle.create_dataset("seeds", (1, n_stages), maxshape=(None, n_stages), dtype='int64')
+
+ handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
+
+ def to_hdf5(self, handle, index):
+ """ Converts results object into an hdf5 object.
+
+ Parameters
+ ----------
+ handle : h5py.File or h5py.Group
+ An hdf5 file or group type to store this in.
+ index : int
+ What step is this?
+ """
+
+ if "/number" not in handle:
+ comm.barrier()
+ self.create_hdf5(handle)
+
+ comm.barrier()
+
+ # Grab handles
+ number_dset = handle["/number"]
+ rxn_dset = handle["/reaction rates"]
+ eigenvalues_dset = handle["/eigenvalues"]
+ seeds_dset = handle["/seeds"]
+ time_dset = handle["/time"]
+
+ # Get number of results stored
+ number_shape = list(number_dset.shape)
+ number_results = number_shape[0]
+
+ new_shape = index + 1
+
+ if number_results < new_shape:
+ # Extend first dimension by 1
+ number_shape[0] = new_shape
+ number_dset.resize(number_shape)
+
+ rxn_shape = list(rxn_dset.shape)
+ rxn_shape[0] = new_shape
+ rxn_dset.resize(rxn_shape)
+
+ eigenvalues_shape = list(eigenvalues_dset.shape)
+ eigenvalues_shape[0] = new_shape
+ eigenvalues_dset.resize(eigenvalues_shape)
+
+ seeds_shape = list(seeds_dset.shape)
+ seeds_shape[0] = new_shape
+ seeds_dset.resize(seeds_shape)
+
+ time_shape = list(time_dset.shape)
+ time_shape[0] = new_shape
+ time_dset.resize(time_shape)
+
+ # If nothing to write, just return
+ if len(self.mat_to_ind) == 0:
+ return
+
+ # Add data
+ # Note, for the last step, self.n_stages = 1, even if n_stages != 1.
+ n_stages = self.n_stages
+ inds = [self.mat_to_hdf5_ind[mat] for mat in self.mat_to_ind]
+ low = min(inds)
+ high = max(inds)
+ for i in range(n_stages):
+ number_dset[index, i, low:high+1, :] = self.data[i, :, :]
+ rxn_dset[index, i, low:high+1, :, :] = self.rates[i][:, :, :]
+ if comm.rank == 0:
+ eigenvalues_dset[index, i] = self.k[i]
+ seeds_dset[index, i] = self.seeds[i]
+ if comm.rank == 0:
+ time_dset[index, :] = self.time
+
+ def from_hdf5(self, handle, index):
+ """ Loads results object from HDF5.
+
+ Parameters
+ ----------
+ handle : h5py.File or h5py.Group
+ An hdf5 file or group type to load from.
+ index : int
+ What step is this?
+ """
+
+ # Grab handles
+ number_dset = handle["/number"]
+ eigenvalues_dset = handle["/eigenvalues"]
+ seeds_dset = handle["/seeds"]
+ time_dset = handle["/time"]
+
+ self.data = number_dset[index, :, :, :]
+ self.k = eigenvalues_dset[index, :]
+ self.seeds = seeds_dset[index, :]
+ self.time = time_dset[index, :]
+
+ # Reconstruct dictionaries
+ self.volume = OrderedDict()
+ self.mat_to_ind = OrderedDict()
+ self.nuc_to_ind = OrderedDict()
+ rxn_nuc_to_ind = OrderedDict()
+ rxn_to_ind = OrderedDict()
+
+ for mat in handle["/cells"]:
+ mat_handle = handle["/cells/" + mat]
+ vol = mat_handle.attrs["volume"]
+ ind = mat_handle.attrs["index"]
+
+ self.volume[mat] = vol
+ self.mat_to_ind[mat] = ind
+
+ for nuc in handle["/nuclides"]:
+ nuc_handle = handle["/nuclides/" + nuc]
+ ind_atom = nuc_handle.attrs["atom number index"]
+ self.nuc_to_ind[nuc] = ind_atom
+
+ if "reaction rate index" in nuc_handle.attrs:
+ rxn_nuc_to_ind[nuc] = nuc_handle.attrs["reaction rate index"]
+
+ for rxn in handle["/reactions"]:
+ rxn_handle = handle["/reactions/" + rxn]
+ rxn_to_ind[rxn] = rxn_handle.attrs["index"]
+
+ self.rates = []
+ # Reconstruct reactions
+ for i in range(self.n_stages):
+ rate = ReactionRates(self.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind)
+
+ rate.rates = handle["/reaction rates"][index, i, :, :, :]
+ self.rates.append(rate)
+
+
+def get_dict(number):
+ """ Given an operator nested dictionary, output indexing dictionaries.
+
+ These indexing dictionaries map mat IDs and nuclide names to indices
+ inside of Results.data.
+
+ Parameters
+ ----------
+ number : AtomNumber
+ The object to extract dictionaries from
+
+ Returns
+ -------
+ mat_to_ind : OrderedDict of str to int
+ Maps mat strings to index in array.
+ nuc_to_ind : OrderedDict of str to int
+ Maps nuclide strings to index in array.
+ """
+ mat_to_ind = OrderedDict()
+ nuc_to_ind = OrderedDict()
+
+ for nuc in number.nuc_to_ind:
+ nuc_ind = number.nuc_to_ind[nuc]
+ if nuc_ind < number.n_nuc_burn:
+ nuc_to_ind[nuc] = nuc_ind
+
+ for mat in number.mat_to_ind:
+ mat_ind = number.mat_to_ind[mat]
+ if mat_ind < number.n_mat_burn:
+ mat_to_ind[mat] = mat_ind
+
+ return mat_to_ind, nuc_to_ind
+
+
+def write_results(result, filename, index):
+ """ Outputs result to an .hdf5 file.
+
+ Parameters
+ ----------
+ result : Results
+ Object to be stored in a file.
+ filename : String
+ Target filename.
+ index : int
+ What step is this?
+ """
+
+ if have_mpi and h5py.get_config().mpi:
+ kwargs = {'driver': 'mpio', 'comm': comm}
+ else:
+ kwargs = {}
+
+ kwargs['mode'] = "w" if index == 0 else "a"
+
+ with h5py.File(filename, **kwargs) as handle:
+ result.to_hdf5(handle, index)
+
+
+def read_results(filename):
+ """ Reads out a list of results objects from an hdf5 file.
+
+ Parameters
+ ----------
+ filename : str
+ The filename to read from.
+
+ Returns
+ -------
+ results : list of Results
+ The result objects.
+ """
+
+ file = h5py.File(filename, "r")
+
+ assert file["/version"].value == RESULTS_VERSION
+
+ # Grab handles
+ number_dset = file["/number"]
+
+ # Get number of results stored
+ number_shape = list(number_dset.shape)
+ number_results = number_shape[0]
+
+ results = []
+
+ for i in range(number_results):
+ result = Results()
+ result.from_hdf5(file, i)
+ results.append(result)
+
+ file.close()
+
+ return results
diff --git a/openmc/deplete/utilities.py b/openmc/deplete/utilities.py
new file mode 100644
index 000000000..54632ed9c
--- /dev/null
+++ b/openmc/deplete/utilities.py
@@ -0,0 +1,98 @@
+""" The utilities module.
+
+Contains functions that can be used to post-process objects that come out of
+the results module.
+"""
+
+import numpy as np
+
+def evaluate_single_nuclide(results, cell, nuc):
+ """ Evaluates a single nuclide in a single cell from a results list.
+
+ Parameters
+ ----------
+ results : list of results
+ The results to extract data from. Must be sorted and continuous.
+ cell : str
+ Cell name to evaluate
+ nuc : str
+ Nuclide name to evaluate
+
+ Returns
+ -------
+ time : numpy.array
+ Time vector.
+ concentration : numpy.array
+ Total number of atoms in the cell.
+ """
+
+ n_points = len(results)
+ time = np.zeros(n_points)
+ concentration = np.zeros(n_points)
+
+ # Evaluate value in each region
+ for i, result in enumerate(results):
+ time[i] = result.time[0]
+ concentration[i] = result[0, cell, nuc]
+
+ return time, concentration
+
+def evaluate_reaction_rate(results, cell, nuc, rxn):
+ """ Evaluates a single nuclide reaction rate in a single cell from a results list.
+
+ Parameters
+ ----------
+ results : list of Results
+ The results to extract data from. Must be sorted and continuous.
+ cell : str
+ Cell name to evaluate
+ nuc : str
+ Nuclide name to evaluate
+ rxn : str
+ Reaction rate to evaluate
+
+ Returns
+ -------
+ time : numpy.array
+ Time vector.
+ rate : numpy.array
+ Reaction rate.
+ """
+
+ n_points = len(results)
+ time = np.zeros(n_points)
+ rate = np.zeros(n_points)
+ # Evaluate value in each region
+ for i, result in enumerate(results):
+ time[i] = result.time[0]
+ rate[i] = result.rates[0][cell, nuc, rxn] * result[0, cell, nuc]
+
+ return time, rate
+
+def evaluate_eigenvalue(results):
+ """ Evaluates the eigenvalue from a results list.
+
+ Parameters
+ ----------
+ results : list of Results
+ The results to extract data from. Must be sorted and continuous.
+
+ Returns
+ -------
+ time : numpy.array
+ Time vector.
+ eigenvalue : numpy.array
+ Eigenvalue.
+ """
+
+ n_points = len(results)
+ time = np.zeros(n_points)
+ eigenvalue = np.zeros(n_points)
+
+ # Evaluate value in each region
+ for i, result in enumerate(results):
+
+ time[i] = result.time[0]
+ eigenvalue[i] = result.k[0]
+
+ return time, eigenvalue
diff --git a/scripts/example_geometry.py b/scripts/example_geometry.py
new file mode 100644
index 000000000..9afcc0d46
--- /dev/null
+++ b/scripts/example_geometry.py
@@ -0,0 +1,358 @@
+"""An example file showing how to make a geometry.
+
+This particular example creates a 3x3 geometry, with 8 regular pins and one
+Gd-157 2 wt-percent enriched. All pins are segmented.
+"""
+
+from collections import OrderedDict
+import math
+
+import numpy as np
+import openmc
+
+from opendeplete import density_to_mat
+
+
+def generate_initial_number_density():
+ """ Generates initial number density.
+
+ These results were from a CASMO5 run in which the gadolinium pin was
+ loaded with 2 wt percent of Gd-157.
+ """
+
+ # Concentration to be used for all fuel pins
+ fuel_dict = OrderedDict()
+ fuel_dict['U235'] = 1.05692e21
+ fuel_dict['U234'] = 1.00506e19
+ fuel_dict['U238'] = 2.21371e22
+ fuel_dict['O16'] = 4.62954e22
+ fuel_dict['O17'] = 1.127684e20
+ fuel_dict['I135'] = 1.0e10
+ fuel_dict['Xe135'] = 1.0e10
+ fuel_dict['Xe136'] = 1.0e10
+ fuel_dict['Cs135'] = 1.0e10
+ fuel_dict['Gd156'] = 1.0e10
+ fuel_dict['Gd157'] = 1.0e10
+ # fuel_dict['O18'] = 9.51352e19 # Does not exist in ENDF71, merged into 17
+
+ # Concentration to be used for the gadolinium fuel pin
+ fuel_gd_dict = OrderedDict()
+ fuel_gd_dict['U235'] = 1.03579e21
+ fuel_gd_dict['U238'] = 2.16943e22
+ fuel_gd_dict['Gd156'] = 3.95517E+10
+ fuel_gd_dict['Gd157'] = 1.08156e20
+ fuel_gd_dict['O16'] = 4.64035e22
+ fuel_dict['I135'] = 1.0e10
+ fuel_dict['Xe136'] = 1.0e10
+ fuel_dict['Xe135'] = 1.0e10
+ fuel_dict['Cs135'] = 1.0e10
+ # There are a whole bunch of 1e-10 stuff here.
+
+ # Concentration to be used for cladding
+ clad_dict = OrderedDict()
+ clad_dict['O16'] = 3.07427e20
+ clad_dict['O17'] = 7.48868e17
+ clad_dict['Cr50'] = 3.29620e18
+ clad_dict['Cr52'] = 6.35639e19
+ clad_dict['Cr53'] = 7.20763e18
+ clad_dict['Cr54'] = 1.79413e18
+ clad_dict['Fe54'] = 5.57350e18
+ clad_dict['Fe56'] = 8.74921e19
+ clad_dict['Fe57'] = 2.02057e18
+ clad_dict['Fe58'] = 2.68901e17
+ clad_dict['Cr50'] = 3.29620e18
+ clad_dict['Cr52'] = 6.35639e19
+ clad_dict['Cr53'] = 7.20763e18
+ clad_dict['Cr54'] = 1.79413e18
+ clad_dict['Ni58'] = 2.51631e19
+ clad_dict['Ni60'] = 9.69278e18
+ clad_dict['Ni61'] = 4.21338e17
+ clad_dict['Ni62'] = 1.34341e18
+ clad_dict['Ni64'] = 3.43127e17
+ clad_dict['Zr90'] = 2.18320e22
+ clad_dict['Zr91'] = 4.76104e21
+ clad_dict['Zr92'] = 7.27734e21
+ clad_dict['Zr94'] = 7.37494e21
+ clad_dict['Zr96'] = 1.18814e21
+ clad_dict['Sn112'] = 4.67352e18
+ clad_dict['Sn114'] = 3.17992e18
+ clad_dict['Sn115'] = 1.63814e18
+ clad_dict['Sn116'] = 7.00546e19
+ clad_dict['Sn117'] = 3.70027e19
+ clad_dict['Sn118'] = 1.16694e20
+ clad_dict['Sn119'] = 4.13872e19
+ clad_dict['Sn120'] = 1.56973e20
+ clad_dict['Sn122'] = 2.23076e19
+ clad_dict['Sn124'] = 2.78966e19
+
+ # Gap concentration
+ # Funny enough, the example problem uses air.
+ gap_dict = OrderedDict()
+ gap_dict['O16'] = 7.86548e18
+ gap_dict['O17'] = 2.99548e15
+ gap_dict['N14'] = 3.38646e19
+ gap_dict['N15'] = 1.23717e17
+
+ # Concentration to be used for coolant
+ # No boron
+ cool_dict = OrderedDict()
+ cool_dict['H1'] = 4.68063e22
+ cool_dict['O16'] = 2.33427e22
+ cool_dict['O17'] = 8.89086e18
+
+ # Store these dictionaries in the initial conditions dictionary
+ initial_density = OrderedDict()
+ initial_density['fuel_gd'] = fuel_gd_dict
+ initial_density['fuel'] = fuel_dict
+ initial_density['gap'] = gap_dict
+ initial_density['clad'] = clad_dict
+ initial_density['cool'] = cool_dict
+
+ # Set up libraries to use
+ temperature = OrderedDict()
+ sab = OrderedDict()
+
+ # Toggle betweeen MCNP and NNDC data
+ MCNP = False
+
+ if MCNP:
+ temperature['fuel_gd'] = 900.0
+ temperature['fuel'] = 900.0
+ # We approximate temperature of everything as 600K, even though it was
+ # actually 580K.
+ temperature['gap'] = 600.0
+ temperature['clad'] = 600.0
+ temperature['cool'] = 600.0
+ else:
+ temperature['fuel_gd'] = 293.6
+ temperature['fuel'] = 293.6
+ temperature['gap'] = 293.6
+ temperature['clad'] = 293.6
+ temperature['cool'] = 293.6
+
+ sab['cool'] = 'c_H_in_H2O'
+
+ # Set up burnable materials
+ burn = OrderedDict()
+ burn['fuel_gd'] = True
+ burn['fuel'] = True
+ burn['gap'] = False
+ burn['clad'] = False
+ burn['cool'] = False
+
+ return temperature, sab, initial_density, burn
+
+def segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad):
+ """ Calculates a segmented pin.
+
+ Separates a pin with n_rings and n_wedges. All cells have equal volume.
+ Pin is centered at origin.
+ """
+
+ # Calculate all the volumes of interest
+ v_fuel = math.pi * r_fuel**2
+ v_gap = math.pi * r_gap**2 - v_fuel
+ v_clad = math.pi * r_clad**2 - v_fuel - v_gap
+ v_ring = v_fuel / n_rings
+ v_segment = v_ring / n_wedges
+
+ # Compute ring radiuses
+ r_rings = np.zeros(n_rings)
+
+ for i in range(n_rings):
+ r_rings[i] = math.sqrt(1.0/(math.pi) * v_ring * (i+1))
+
+ # Compute thetas
+ theta = np.linspace(0, 2*math.pi, n_wedges + 1)
+
+ # Compute surfaces
+ fuel_rings = [openmc.ZCylinder(x0=0, y0=0, R=r_rings[i])
+ for i in range(n_rings)]
+
+ fuel_wedges = [openmc.Plane(A=math.cos(theta[i]), B=math.sin(theta[i]))
+ for i in range(n_wedges)]
+
+ gap_ring = openmc.ZCylinder(x0=0, y0=0, R=r_gap)
+ clad_ring = openmc.ZCylinder(x0=0, y0=0, R=r_clad)
+
+ # Create cells
+ fuel_cells = []
+ if n_wedges == 1:
+ for i in range(n_rings):
+ cell = openmc.Cell(name='fuel')
+ if i == 0:
+ cell.region = -fuel_rings[0]
+ else:
+ cell.region = +fuel_rings[i-1] & -fuel_rings[i]
+ fuel_cells.append(cell)
+ else:
+ for i in range(n_rings):
+ for j in range(n_wedges):
+ cell = openmc.Cell(name='fuel')
+ if i == 0:
+ if j != n_wedges-1:
+ cell.region = (-fuel_rings[0]
+ & +fuel_wedges[j]
+ & -fuel_wedges[j+1])
+ else:
+ cell.region = (-fuel_rings[0]
+ & +fuel_wedges[j]
+ & -fuel_wedges[0])
+ else:
+ if j != n_wedges-1:
+ cell.region = (+fuel_rings[i-1]
+ & -fuel_rings[i]
+ & +fuel_wedges[j]
+ & -fuel_wedges[j+1])
+ else:
+ cell.region = (+fuel_rings[i-1]
+ & -fuel_rings[i]
+ & +fuel_wedges[j]
+ & -fuel_wedges[0])
+ fuel_cells.append(cell)
+
+ # Gap ring
+ gap_cell = openmc.Cell(name='gap')
+ gap_cell.region = +fuel_rings[-1] & -gap_ring
+ fuel_cells.append(gap_cell)
+
+ # Clad ring
+ clad_cell = openmc.Cell(name='clad')
+ clad_cell.region = +gap_ring & -clad_ring
+ fuel_cells.append(clad_cell)
+
+ # Moderator
+ mod_cell = openmc.Cell(name='cool')
+ mod_cell.region = +clad_ring
+ fuel_cells.append(mod_cell)
+
+ # Form universe
+ fuel_u = openmc.Universe()
+ fuel_u.add_cells(fuel_cells)
+
+ return fuel_u, v_segment, v_gap, v_clad
+
+def generate_geometry(n_rings, n_wedges):
+ """ Generates example geometry.
+
+ This function creates the initial geometry, a 9 pin reflective problem.
+ One pin, containing gadolinium, is discretized into sectors.
+
+ In addition to what one would do with the general OpenMC geometry code, it
+ is necessary to create a dictionary, volume, that maps a cell ID to a
+ volume. Further, by naming cells the same as the above materials, the code
+ can automatically handle the mapping.
+
+ Parameters
+ ----------
+ n_rings : int
+ Number of rings to generate for the geometry
+ n_wedges : int
+ Number of wedges to generate for the geometry
+ """
+
+ pitch = 1.26197
+ r_fuel = 0.412275
+ r_gap = 0.418987
+ r_clad = 0.476121
+
+ n_pin = 3
+
+ # This table describes the 'fuel' to actual type mapping
+ # It's not necessary to do it this way. Just adjust the initial conditions
+ # below.
+ mapping = ['fuel', 'fuel', 'fuel',
+ 'fuel', 'fuel_gd', 'fuel',
+ 'fuel', 'fuel', 'fuel']
+
+ # Form pin cell
+ fuel_u, v_segment, v_gap, v_clad = segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad)
+
+ # Form lattice
+ all_water_c = openmc.Cell(name='cool')
+ all_water_u = openmc.Universe(cells=(all_water_c, ))
+
+ lattice = openmc.RectLattice()
+ lattice.pitch = [pitch]*2
+ lattice.lower_left = [-pitch*n_pin/2, -pitch*n_pin/2]
+ lattice_array = [[fuel_u for i in range(n_pin)] for j in range(n_pin)]
+ lattice.universes = lattice_array
+ lattice.outer = all_water_u
+
+ # Bound universe
+ x_low = openmc.XPlane(x0=-pitch*n_pin/2, boundary_type='reflective')
+ x_high = openmc.XPlane(x0=pitch*n_pin/2, boundary_type='reflective')
+ y_low = openmc.YPlane(y0=-pitch*n_pin/2, boundary_type='reflective')
+ y_high = openmc.YPlane(y0=pitch*n_pin/2, boundary_type='reflective')
+ z_low = openmc.ZPlane(z0=-10, boundary_type='reflective')
+ z_high = openmc.ZPlane(z0=10, boundary_type='reflective')
+
+ # Compute bounding box
+ lower_left = [-pitch*n_pin/2, -pitch*n_pin/2, -10]
+ upper_right = [pitch*n_pin/2, pitch*n_pin/2, 10]
+
+ root_c = openmc.Cell(fill=lattice)
+ root_c.region = (+x_low & -x_high
+ & +y_low & -y_high
+ & +z_low & -z_high)
+ root_u = openmc.Universe(universe_id=0, cells=(root_c, ))
+ geometry = openmc.Geometry(root_u)
+
+ v_cool = pitch**2 - (v_gap + v_clad + n_rings * n_wedges * v_segment)
+
+ # Store volumes for later usage
+ volume = {'fuel': v_segment, 'gap':v_gap, 'clad':v_clad, 'cool':v_cool}
+
+ return geometry, volume, mapping, lower_left, upper_right
+
+def generate_problem(n_rings=5, n_wedges=8):
+ """ Merges geometry and materials.
+
+ This function initializes the materials for each cell using the dictionaries
+ provided by generate_initial_number_density. It is assumed a cell named
+ 'fuel' will have further region differentiation (see mapping).
+
+ Parameters
+ ----------
+ n_rings : int, optional
+ Number of rings to generate for the geometry
+ n_wedges : int, optional
+ Number of wedges to generate for the geometry
+ """
+
+ # Get materials dictionary, geometry, and volumes
+ temperature, sab, initial_density, burn = generate_initial_number_density()
+ geometry, volume, mapping, lower_left, upper_right = generate_geometry(n_rings, n_wedges)
+
+ # Apply distribmats, fill geometry
+ cells = geometry.root_universe.get_all_cells()
+ for cell_id in cells:
+ cell = cells[cell_id]
+ if cell.name == 'fuel':
+
+ omc_mats = []
+
+ for cell_type in mapping:
+ omc_mat = density_to_mat(initial_density[cell_type])
+
+ if cell_type in sab:
+ omc_mat.add_s_alpha_beta(sab[cell_type])
+ omc_mat.temperature = temperature[cell_type]
+ omc_mat.depletable = burn[cell_type]
+ omc_mat.volume = volume['fuel']
+
+ omc_mats.append(omc_mat)
+
+ cell.fill = omc_mats
+ elif cell.name != '':
+ omc_mat = density_to_mat(initial_density[cell.name])
+
+ if cell.name in sab:
+ omc_mat.add_s_alpha_beta(sab[cell.name])
+ omc_mat.temperature = temperature[cell.name]
+ omc_mat.depletable = burn[cell.name]
+ omc_mat.volume = volume[cell.name]
+
+ cell.fill = omc_mat
+
+ return geometry, lower_left, upper_right
diff --git a/scripts/example_plot.py b/scripts/example_plot.py
new file mode 100644
index 000000000..d2c6ee9d6
--- /dev/null
+++ b/scripts/example_plot.py
@@ -0,0 +1,46 @@
+"""An example file showing how to plot data from a simulation."""
+
+import matplotlib.pyplot as plt
+
+from opendeplete import read_results, \
+ evaluate_single_nuclide, \
+ evaluate_reaction_rate, \
+ evaluate_eigenvalue
+
+# Set variables for where the data is, and what we want to read out.
+result_folder = "test"
+
+# Load data
+results = read_results(result_folder + "/results.h5")
+
+cell = "5"
+nuc = "Gd157"
+rxn = "(n,gamma)"
+
+# Total number of nuclides
+plt.figure()
+# Pointwise data
+x, y = evaluate_single_nuclide(results, cell, nuc)
+plt.semilogy(x, y)
+
+plt.xlabel("Time, s")
+plt.ylabel("Total Number")
+plt.savefig("number.pdf")
+
+# Reaction rate
+plt.figure()
+x, y = evaluate_reaction_rate(results, cell, nuc, rxn)
+plt.plot(x, y)
+plt.xlabel("Time, s")
+plt.ylabel("Reaction Rate, 1/s")
+
+plt.savefig("rate.pdf")
+
+# Eigenvalue
+plt.figure()
+x, y = evaluate_eigenvalue(results)
+plt.plot(x, y)
+plt.xlabel("Time, s")
+plt.ylabel("Eigenvalue")
+
+plt.savefig("eigvl.pdf")
diff --git a/scripts/example_run.py b/scripts/example_run.py
new file mode 100644
index 000000000..bb80f6582
--- /dev/null
+++ b/scripts/example_run.py
@@ -0,0 +1,39 @@
+"""An example file showing how to run a simulation."""
+
+import numpy as np
+import opendeplete
+
+import example_geometry
+
+# Load geometry from example
+geometry, lower_left, upper_right = example_geometry.generate_problem()
+
+# Create dt vector for 5.5 months with 15 day timesteps
+dt1 = 15*24*60*60 # 15 days
+dt2 = 5.5*30*24*60*60 # 5.5 months
+N = np.floor(dt2/dt1)
+
+dt = np.repeat([dt1], N)
+
+# Create settings variable
+settings = opendeplete.OpenMCSettings()
+
+settings.openmc_call = "openmc"
+# An example for mpiexec:
+# settings.openmc_call = ["mpiexec", "openmc"]
+settings.particles = 1000
+settings.batches = 100
+settings.inactive = 40
+settings.lower_left = lower_left
+settings.upper_right = upper_right
+settings.entropy_dimension = [10, 10, 1]
+
+joule_per_mev = 1.6021766208e-13
+settings.power = 2.337e15*4*joule_per_mev # MeV/second cm from CASMO
+settings.dt_vec = dt
+settings.output_dir = 'test'
+
+op = opendeplete.OpenMCOperator(geometry, settings)
+
+# Perform simulation using the MCNPX/MCNP6 algorithm
+opendeplete.integrator.cecm(op)
diff --git a/scripts/make_chain.py b/scripts/make_chain.py
new file mode 100644
index 000000000..2e0d9d3bc
--- /dev/null
+++ b/scripts/make_chain.py
@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+
+import glob
+import os
+from zipfile import ZipFile
+
+import requests
+from tqdm import tqdm
+import opendeplete
+
+
+urls = [
+ 'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-neutrons.zip',
+ 'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-decay.zip',
+ 'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-nfy.zip'
+]
+
+
+def download_file(url):
+ response = requests.get(url, stream=True)
+ filesize = int(response.headers.get('content-length'))
+
+ # Check if file already downloaded
+ basename = url.split('/')[-1]
+ if os.path.exists(basename):
+ if os.path.getsize(basename) == filesize:
+ return basename
+ else:
+ overwrite = input('Overwrite {}? ([y]/n) '.format(basename))
+ if overwrite.lower().startswith('n'):
+ return basename
+
+ with open(basename, 'wb') as f:
+ with tqdm(desc='Downloading {}'.format(basename),
+ total=filesize, unit='B', unit_scale=True) as pbar:
+ for i, chunk in enumerate(response.iter_content(chunk_size=4096)):
+ pbar.update(4096)
+ if chunk:
+ f.write(chunk)
+
+ return basename
+
+
+def main():
+ for url in urls:
+ basename = download_file(url)
+ with ZipFile(basename, 'r') as zf:
+ print('Extracting {}...'.format(basename))
+ zf.extractall()
+
+ decay_files = glob.glob(os.path.join('decay', '*.endf'))
+ nfy_files = glob.glob(os.path.join('nfy', '*.endf'))
+ neutron_files = glob.glob(os.path.join('neutrons', '*.endf'))
+
+ chain = opendeplete.DepletionChain.from_endf(decay_files, nfy_files, neutron_files)
+ chain.xml_write('chain_endfb71.xml')
+
+
+if __name__ == '__main__':
+ main()
diff --git a/tests/deplete_tests/__init__.py b/tests/deplete_tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/deplete_tests/dummy_geometry.py b/tests/deplete_tests/dummy_geometry.py
new file mode 100644
index 000000000..610151941
--- /dev/null
+++ b/tests/deplete_tests/dummy_geometry.py
@@ -0,0 +1,165 @@
+""" The OpenMC wrapper module.
+
+This module implements the OpenDeplete -> OpenMC linkage.
+"""
+
+import numpy as np
+import scipy.sparse as sp
+
+from opendeplete.reaction_rates import ReactionRates
+from opendeplete.function import Operator
+
+class DummyGeometry(Operator):
+ """ This is a dummy geometry class with no statistical uncertainty.
+
+ y_1' = sin(y_2) y_1 + cos(y_1) y_2
+ y_2' = -cos(y_2) y_1 + sin(y_1) y_2
+
+ y_1(0) = 1
+ y_2(0) = 1
+
+ y_1(1.5) ~ 2.3197067076743316
+ y_2(1.5) ~ 3.1726475740397628
+
+ """
+
+ def __init__(self, settings):
+ Operator.__init__(self, settings)
+
+ @property
+ def chain(self):
+ return self
+
+ def eval(self, vec, print_out=False):
+ """ Evaluates F(y)
+
+ Parameters
+ ----------
+ vec : list of numpy.array
+ Total atoms to be used in function.
+ print_out : bool, optional, ignored
+ Whether or not to print out time.
+
+ Returns
+ -------
+ k : float
+ Zero.
+ rates : ReactionRates
+ Reaction rates from this simulation.
+ seed : int
+ Zero.
+ """
+
+ cell_to_ind = {"1" : 0}
+ nuc_to_ind = {"1" : 0, "2" : 1}
+ react_to_ind = {"1" : 0}
+
+ reaction_rates = ReactionRates(cell_to_ind, nuc_to_ind, react_to_ind)
+
+ reaction_rates[0, 0, 0] = vec[0][0]
+ reaction_rates[0, 1, 0] = vec[0][1]
+
+ # Create a fake rates object
+
+ return 0.0, reaction_rates, 0
+
+ def form_matrix(self, rates):
+ """ Forms the f(y) matrix in y' = f(y)y.
+
+ Nominally a depletion matrix, this is abstracted on the off chance
+ that the function f has nothing to do with depletion at all.
+
+ Parameters
+ ----------
+ rates : numpy.ndarray
+ Slice of reaction rates for a single material
+
+ Returns
+ -------
+ scipy.sparse.csr_matrix
+ Sparse matrix representing f(y).
+ """
+
+ y_1 = rates[0, 0]
+ y_2 = rates[1, 0]
+
+ mat = np.zeros((2, 2))
+ a11 = np.sin(y_2)
+ a12 = np.cos(y_1)
+ a21 = -np.cos(y_2)
+ a22 = np.sin(y_1)
+
+ return sp.csr_matrix(np.array([[a11, a12], [a21, a22]]))
+
+ @property
+ def volume(self):
+ """
+ volume : dict of str float
+ Volumes of material
+ """
+
+ return {"1": 0.0}
+
+ @property
+ def nuc_list(self):
+ """
+ nuc_list : list of str
+ A list of all nuclide names. Used for sorting the simulation.
+ """
+
+ return ["1", "2"]
+
+ @property
+ def burn_list(self):
+ """
+ burn_list : list of str
+ A list of all cell IDs to be burned. Used for sorting the simulation.
+ """
+
+ return ["1"]
+
+ @property
+ def mat_tally_ind(self):
+ """Maps cell name to index in global geometry."""
+ return {"1": 0}
+
+
+ @property
+ def reaction_rates(self):
+ """
+ reaction_rates : ReactionRates
+ Reaction rates from the last operator step.
+ """
+ cell_to_ind = {"1" : 0}
+ nuc_to_ind = {"1" : 0, "2" : 1}
+ react_to_ind = {"1" : 0}
+
+ return ReactionRates(cell_to_ind, nuc_to_ind, react_to_ind)
+
+ def initial_condition(self):
+ """ Returns initial vector.
+
+ Returns
+ -------
+ list of numpy.array
+ Total density for initial conditions.
+ """
+
+ return [np.array((1.0, 1.0))]
+
+ def get_results_info(self):
+ """ Returns volume list, cell 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 cell IDs to be burned. Used for sorting the simulation.
+ full_burn_dict : OrderedDict of str to int
+ Maps cell name to index in global geometry.
+ """
+
+ return self.volume, self.nuc_list, self.burn_list, self.mat_tally_ind
diff --git a/tests/deplete_tests/example_geometry.py b/tests/deplete_tests/example_geometry.py
new file mode 120000
index 000000000..1071aabc0
--- /dev/null
+++ b/tests/deplete_tests/example_geometry.py
@@ -0,0 +1 @@
+../../scripts/example_geometry.py
\ No newline at end of file
diff --git a/tests/deplete_tests/test_atom_number.py b/tests/deplete_tests/test_atom_number.py
new file mode 100644
index 000000000..9a17230f8
--- /dev/null
+++ b/tests/deplete_tests/test_atom_number.py
@@ -0,0 +1,180 @@
+""" Tests for atom_number.py. """
+
+import unittest
+
+import numpy as np
+
+from opendeplete import atom_number
+
+class TestAtomNumber(unittest.TestCase):
+ """ Tests for the AtomNumber class. """
+
+ def test_indexing(self):
+ """Tests the __getitem__ and __setitem__ routines simultaneously."""
+
+ mat_to_ind = {"10000" : 0, "10001" : 1, "10002" : 2}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "U234" : 2}
+ volume = {"10000" : 0.38, "10001" : 0.21}
+
+ number = atom_number.AtomNumber(mat_to_ind, nuc_to_ind, volume, 2, 2)
+
+ number["10000", "U238"] = 1.0
+ number["10001", "U238"] = 2.0
+ number["10000", "U235"] = 3.0
+ number["10001", "U235"] = 4.0
+
+ # String indexing
+ self.assertEqual(number["10000", "U238"], 1.0)
+ self.assertEqual(number["10001", "U238"], 2.0)
+ self.assertEqual(number["10000", "U235"], 3.0)
+ self.assertEqual(number["10001", "U235"], 4.0)
+
+ # Int indexing
+ self.assertEqual(number[0, 0], 1.0)
+ self.assertEqual(number[1, 0], 2.0)
+ self.assertEqual(number[0, 1], 3.0)
+ self.assertEqual(number[1, 1], 4.0)
+
+ number[0, 0] = 5.0
+
+ self.assertEqual(number[0, 0], 5.0)
+ self.assertEqual(number["10000", "U238"], 5.0)
+
+ def test_n_mat(self):
+ """ Test number of materials property. """
+ mat_to_ind = {"10000" : 0, "10001" : 1}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
+ volume = {"10000" : 0.38, "10001" : 0.21}
+
+ number = atom_number.AtomNumber(mat_to_ind, nuc_to_ind, volume, 2, 2)
+
+ self.assertEqual(number.n_mat, 2)
+
+ def test_n_nuc(self):
+ """ Test number of nuclides property. """
+ mat_to_ind = {"10000" : 0, "10001" : 1}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
+ volume = {"10000" : 0.38, "10001" : 0.21}
+
+ number = atom_number.AtomNumber(mat_to_ind, nuc_to_ind, volume, 2, 2)
+
+ self.assertEqual(number.n_nuc, 3)
+
+ def test_burn_nuc_list(self):
+ """ Test the list of burned nuclides property """
+ mat_to_ind = {"10000" : 0, "10001" : 1}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
+ volume = {"10000" : 0.38, "10001" : 0.21}
+
+ number = atom_number.AtomNumber(mat_to_ind, nuc_to_ind, volume, 2, 2)
+
+ self.assertEqual(number.burn_nuc_list, ["U238", "U235"])
+
+ def test_burn_mat_list(self):
+ """ Test the list of burned nuclides property """
+ mat_to_ind = {"10000" : 0, "10001" : 1, "10002" : 2}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
+ volume = {"10000" : 0.38, "10001" : 0.21}
+
+ number = atom_number.AtomNumber(mat_to_ind, nuc_to_ind, volume, 2, 2)
+
+ self.assertEqual(number.burn_mat_list, ["10000", "10001"])
+
+ def test_density_indexing(self):
+ """Tests the get and set_atom_density routines simultaneously."""
+
+ mat_to_ind = {"10000" : 0, "10001" : 1, "10002" : 2}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "U234" : 2}
+ volume = {"10000" : 0.38, "10001" : 0.21}
+
+ number = atom_number.AtomNumber(mat_to_ind, nuc_to_ind, volume, 2, 2)
+
+ number.set_atom_density("10000", "U238", 1.0)
+ number.set_atom_density("10001", "U238", 2.0)
+ number.set_atom_density("10002", "U238", 3.0)
+ number.set_atom_density("10000", "U235", 4.0)
+ number.set_atom_density("10001", "U235", 5.0)
+ number.set_atom_density("10002", "U235", 6.0)
+ number.set_atom_density("10000", "U234", 7.0)
+ number.set_atom_density("10001", "U234", 8.0)
+ number.set_atom_density("10002", "U234", 9.0)
+
+ # String indexing
+ self.assertEqual(number.get_atom_density("10000", "U238"), 1.0)
+ self.assertEqual(number.get_atom_density("10001", "U238"), 2.0)
+ self.assertEqual(number.get_atom_density("10002", "U238"), 3.0)
+ self.assertEqual(number.get_atom_density("10000", "U235"), 4.0)
+ self.assertEqual(number.get_atom_density("10001", "U235"), 5.0)
+ self.assertEqual(number.get_atom_density("10002", "U235"), 6.0)
+ self.assertEqual(number.get_atom_density("10000", "U234"), 7.0)
+ self.assertEqual(number.get_atom_density("10001", "U234"), 8.0)
+ self.assertEqual(number.get_atom_density("10002", "U234"), 9.0)
+
+ # Int indexing
+ self.assertEqual(number.get_atom_density(0, 0), 1.0)
+ self.assertEqual(number.get_atom_density(1, 0), 2.0)
+ self.assertEqual(number.get_atom_density(2, 0), 3.0)
+ self.assertEqual(number.get_atom_density(0, 1), 4.0)
+ self.assertEqual(number.get_atom_density(1, 1), 5.0)
+ self.assertEqual(number.get_atom_density(2, 1), 6.0)
+ self.assertEqual(number.get_atom_density(0, 2), 7.0)
+ self.assertEqual(number.get_atom_density(1, 2), 8.0)
+ self.assertEqual(number.get_atom_density(2, 2), 9.0)
+
+
+ number.set_atom_density(0, 0, 5.0)
+
+ self.assertEqual(number.get_atom_density(0, 0), 5.0)
+
+ # Verify volume is used correctly
+ self.assertEqual(number[0, 0], 5.0 * 0.38)
+ self.assertEqual(number[1, 0], 2.0 * 0.21)
+ self.assertEqual(number[2, 0], 3.0 * 1.0)
+ self.assertEqual(number[0, 1], 4.0 * 0.38)
+ self.assertEqual(number[1, 1], 5.0 * 0.21)
+ self.assertEqual(number[2, 1], 6.0 * 1.0)
+ self.assertEqual(number[0, 2], 7.0 * 0.38)
+ self.assertEqual(number[1, 2], 8.0 * 0.21)
+ self.assertEqual(number[2, 2], 9.0 * 1.0)
+
+ def test_get_mat_slice(self):
+ """Tests getting slices."""
+
+ mat_to_ind = {"10000" : 0, "10001" : 1, "10002" : 2}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "U234" : 2}
+ volume = {"10000" : 0.38, "10001" : 0.21}
+
+ number = atom_number.AtomNumber(mat_to_ind, nuc_to_ind, volume, 2, 2)
+
+ number.number = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])
+
+ sl = number.get_mat_slice(0)
+
+ np.testing.assert_array_equal(sl, np.array([1.0, 2.0]))
+
+ sl = number.get_mat_slice("10000")
+
+ np.testing.assert_array_equal(sl, np.array([1.0, 2.0]))
+
+ def test_set_mat_slice(self):
+ """Tests getting slices."""
+
+ mat_to_ind = {"10000" : 0, "10001" : 1, "10002" : 2}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "U234" : 2}
+ volume = {"10000" : 0.38, "10001" : 0.21}
+
+ number = atom_number.AtomNumber(mat_to_ind, nuc_to_ind, volume, 2, 2)
+
+ number.set_mat_slice(0, [1.0, 2.0])
+
+ self.assertEqual(number[0, 0], 1.0)
+ self.assertEqual(number[0, 1], 2.0)
+
+ number.set_mat_slice("10000", [3.0, 4.0])
+
+ self.assertEqual(number[0, 0], 3.0)
+ self.assertEqual(number[0, 1], 4.0)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_cecm_regression.py b/tests/deplete_tests/test_cecm_regression.py
new file mode 100644
index 000000000..23a634200
--- /dev/null
+++ b/tests/deplete_tests/test_cecm_regression.py
@@ -0,0 +1,69 @@
+""" Regression tests for cecm.py"""
+
+import os
+import unittest
+
+import numpy as np
+
+import opendeplete
+from opendeplete import results
+from opendeplete import utilities
+import test.dummy_geometry as dummy_geometry
+
+
+class TestCECMRegression(unittest.TestCase):
+ """ Regression tests for opendeplete.integrator.cecm algorithm.
+
+ These tests integrate a simple test problem described in dummy_geometry.py.
+ """
+
+ @classmethod
+ def setUpClass(cls):
+ """ Save current directory in case integrator crashes."""
+ cls.cwd = os.getcwd()
+ cls.results = "test_integrator_regression"
+
+ def test_cecm(self):
+ """ Integral regression test of integrator algorithm using CE/CM. """
+
+ settings = opendeplete.Settings()
+ settings.dt_vec = [0.75, 0.75]
+ settings.output_dir = self.results
+
+ op = dummy_geometry.DummyGeometry(settings)
+
+ # Perform simulation using the MCNPX/MCNP6 algorithm
+ opendeplete.cecm(op, print_out=False)
+
+ # Load the files
+ res = results.read_results(settings.output_dir + "/results.h5")
+
+ _, y1 = utilities.evaluate_single_nuclide(res, "1", "1")
+ _, y2 = utilities.evaluate_single_nuclide(res, "1", "2")
+
+ # Mathematica solution
+ s1 = [1.86872629872102, 1.395525772416039]
+ s2 = [2.18097439443550, 2.69429754646747]
+
+ tol = 1.0e-13
+
+ self.assertLess(np.absolute(y1[1] - s1[0]), tol)
+ self.assertLess(np.absolute(y2[1] - s1[1]), tol)
+
+ self.assertLess(np.absolute(y1[2] - s2[0]), tol)
+ self.assertLess(np.absolute(y2[2] - s2[1]), tol)
+
+ @classmethod
+ def tearDownClass(cls):
+ """ Clean up files"""
+
+ os.chdir(cls.cwd)
+
+ opendeplete.comm.barrier()
+ if opendeplete.comm.rank == 0:
+ os.remove(os.path.join(cls.results, "results.h5"))
+ os.rmdir(cls.results)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_cram.py b/tests/deplete_tests/test_cram.py
new file mode 100644
index 000000000..2744adbf4
--- /dev/null
+++ b/tests/deplete_tests/test_cram.py
@@ -0,0 +1,48 @@
+""" Tests for cram.py """
+
+import unittest
+
+import numpy as np
+import scipy.sparse as sp
+
+from opendeplete.integrator import CRAM16, CRAM48
+
+class TestCram(unittest.TestCase):
+ """ Tests for cram.py
+
+ Compares a few Mathematica matrix exponentials to CRAM16/CRAM48.
+ """
+
+ def test_CRAM16(self):
+ """ Test 16-term CRAM. """
+ x = np.array([1.0, 1.0])
+ mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
+ dt = 0.1
+
+ z = CRAM16(mat, x, dt)
+
+ # Solution from mathematica
+ z0 = np.array((0.904837418035960, 0.576799023327476))
+
+ tol = 1.0e-15
+
+ self.assertLess(np.linalg.norm(z - z0), tol)
+
+ def test_CRAM48(self):
+ """ Test 48-term CRAM. """
+ x = np.array([1.0, 1.0])
+ mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
+ dt = 0.1
+
+ z = CRAM48(mat, x, dt)
+
+ # Solution from mathematica
+ z0 = np.array((0.904837418035960, 0.576799023327476))
+
+ tol = 1.0e-15
+
+ self.assertLess(np.linalg.norm(z - z0), tol)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_depletion_chain.py b/tests/deplete_tests/test_depletion_chain.py
new file mode 100644
index 000000000..216d1e68f
--- /dev/null
+++ b/tests/deplete_tests/test_depletion_chain.py
@@ -0,0 +1,197 @@
+""" Tests for depletion_chain.py"""
+
+from collections import OrderedDict
+import os
+import unittest
+
+import numpy as np
+
+from opendeplete import comm, depletion_chain, reaction_rates, nuclide
+
+
+class TestDepletionChain(unittest.TestCase):
+ """ Tests for DepletionChain class."""
+
+ def test__init__(self):
+ """ Test depletion chain initialization."""
+ dep = depletion_chain.DepletionChain()
+
+ self.assertIsInstance(dep.nuclides, list)
+ self.assertIsInstance(dep.nuclide_dict, OrderedDict)
+ self.assertIsInstance(dep.react_to_ind, OrderedDict)
+
+ def test_n_nuclides(self):
+ """ Test depletion chain n_nuclides parameter. """
+ dep = depletion_chain.DepletionChain()
+
+ dep.nuclides = ["NucA", "NucB", "NucC"]
+
+ self.assertEqual(dep.n_nuclides, 3)
+
+ def test_from_endf(self):
+ """Test depletion chain building from ENDF. Empty at the moment until we figure
+ out a good way to unit-test this."""
+ pass
+
+ def test_xml_read(self):
+ """ Read chain_test.xml and ensure all values are correct. """
+ # Unfortunately, this routine touches a lot of the code, but most of
+ # the components external to depletion_chain.py are simple storage
+ # types.
+
+ dep = depletion_chain.DepletionChain.xml_read("chains/chain_test.xml")
+
+ # Basic checks
+ self.assertEqual(dep.n_nuclides, 3)
+
+ # A tests
+ nuc = dep.nuclides[dep.nuclide_dict["A"]]
+
+ self.assertEqual(nuc.name, "A")
+ self.assertEqual(nuc.half_life, 2.36520E+04)
+ self.assertEqual(nuc.n_decay_modes, 2)
+ modes = nuc.decay_modes
+ self.assertEqual([m.target for m in modes], ["B", "C"])
+ self.assertEqual([m.type for m in modes], ["beta1", "beta2"])
+ self.assertEqual([m.branching_ratio for m in modes], [0.6, 0.4])
+ self.assertEqual(nuc.n_reaction_paths, 1)
+ self.assertEqual([r.target for r in nuc.reactions], ["C"])
+ self.assertEqual([r.type for r in nuc.reactions], ["(n,gamma)"])
+ self.assertEqual([r.branching_ratio for r in nuc.reactions], [1.0])
+
+ # B tests
+ nuc = dep.nuclides[dep.nuclide_dict["B"]]
+
+ self.assertEqual(nuc.name, "B")
+ self.assertEqual(nuc.half_life, 3.29040E+04)
+ self.assertEqual(nuc.n_decay_modes, 1)
+ modes = nuc.decay_modes
+ self.assertEqual([m.target for m in modes], ["A"])
+ self.assertEqual([m.type for m in modes], ["beta"])
+ self.assertEqual([m.branching_ratio for m in modes], [1.0])
+ self.assertEqual(nuc.n_reaction_paths, 1)
+ self.assertEqual([r.target for r in nuc.reactions], ["C"])
+ self.assertEqual([r.type for r in nuc.reactions], ["(n,gamma)"])
+ self.assertEqual([r.branching_ratio for r in nuc.reactions], [1.0])
+
+ # C tests
+ nuc = dep.nuclides[dep.nuclide_dict["C"]]
+
+ self.assertEqual(nuc.name, "C")
+ self.assertEqual(nuc.n_decay_modes, 0)
+ self.assertEqual(nuc.n_reaction_paths, 3)
+ self.assertEqual([r.target for r in nuc.reactions], [None, "A", "B"])
+ self.assertEqual([r.type for r in nuc.reactions], ["fission", "(n,gamma)", "(n,gamma)"])
+ self.assertEqual([r.branching_ratio for r in nuc.reactions], [1.0, 0.7, 0.3])
+
+ # Yield tests
+ self.assertEqual(nuc.yield_energies, [0.0253])
+ self.assertEqual(list(nuc.yield_data.keys()), [0.0253])
+ self.assertEqual(nuc.yield_data[0.0253],
+ [("A", 0.0292737), ("B", 0.002566345)])
+
+ def test_xml_write(self):
+ """Test writing a depletion chain to XML."""
+
+ # Prevent different MPI ranks from conflicting
+ filename = 'test%u.xml' % comm.rank
+
+ A = nuclide.Nuclide()
+ A.name = "A"
+ A.half_life = 2.36520e4
+ A.decay_modes = [
+ nuclide.DecayTuple("beta1", "B", 0.6),
+ nuclide.DecayTuple("beta2", "C", 0.4)
+ ]
+ A.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)]
+
+ B = nuclide.Nuclide()
+ B.name = "B"
+ B.half_life = 3.29040e4
+ B.decay_modes = [nuclide.DecayTuple("beta", "A", 1.0)]
+ B.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)]
+
+ C = nuclide.Nuclide()
+ C.name = "C"
+ C.reactions = [
+ nuclide.ReactionTuple("fission", None, 2.0e8, 1.0),
+ nuclide.ReactionTuple("(n,gamma)", "A", 0.0, 0.7),
+ nuclide.ReactionTuple("(n,gamma)", "B", 0.0, 0.3)
+ ]
+ C.yield_energies = [0.0253]
+ C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]}
+
+ chain = depletion_chain.DepletionChain()
+ chain.nuclides = [A, B, C]
+ chain.xml_write(filename)
+
+ original = open('chains/chain_test.xml', 'r').read()
+ chain_xml = open(filename, 'r').read()
+ self.assertEqual(original, chain_xml)
+
+ os.remove(filename)
+
+ def test_form_matrix(self):
+ """ Using chain_test, and a dummy reaction rate, compute the matrix. """
+ # Relies on test_xml_read passing.
+
+ dep = depletion_chain.DepletionChain.xml_read("chains/chain_test.xml")
+
+ cell_ind = {"10000": 0, "10001": 1}
+ nuc_ind = {"A": 0, "B": 1, "C": 2}
+ react_ind = dep.react_to_ind
+
+ react = reaction_rates.ReactionRates(cell_ind, nuc_ind, react_ind)
+
+ dep.nuc_to_react_ind = nuc_ind
+
+ react["10000", "C", "fission"] = 1.0
+ react["10000", "A", "(n,gamma)"] = 2.0
+ react["10000", "B", "(n,gamma)"] = 3.0
+ react["10000", "C", "(n,gamma)"] = 4.0
+
+ mat = dep.form_matrix(react[0, :, :])
+ # Loss A, decay, (n, gamma)
+ mat00 = -np.log(2) / 2.36520E+04 - 2
+ # A -> B, decay, 0.6 branching ratio
+ mat10 = np.log(2) / 2.36520E+04 * 0.6
+ # A -> C, decay, 0.4 branching ratio + (n,gamma)
+ mat20 = np.log(2) / 2.36520E+04 * 0.4 + 2
+
+ # B -> A, decay, 1.0 branching ratio
+ mat01 = np.log(2)/3.29040E+04
+ # Loss B, decay, (n, gamma)
+ mat11 = -np.log(2)/3.29040E+04 - 3
+ # B -> C, (n, gamma)
+ mat21 = 3
+
+ # C -> A fission, (n, gamma)
+ mat02 = 0.0292737 * 1.0 + 4.0 * 0.7
+ # C -> B fission, (n, gamma)
+ mat12 = 0.002566345 * 1.0 + 4.0 * 0.3
+ # Loss C, fission, (n, gamma)
+ mat22 = -1.0 - 4.0
+
+ self.assertEqual(mat[0, 0], mat00)
+ self.assertEqual(mat[1, 0], mat10)
+ self.assertEqual(mat[2, 0], mat20)
+ self.assertEqual(mat[0, 1], mat01)
+ self.assertEqual(mat[1, 1], mat11)
+ self.assertEqual(mat[2, 1], mat21)
+ self.assertEqual(mat[0, 2], mat02)
+ self.assertEqual(mat[1, 2], mat12)
+ self.assertEqual(mat[2, 2], mat22)
+
+ def test_nuc_by_ind(self):
+ """ Test nuc_by_ind converter function. """
+ dep = depletion_chain.DepletionChain()
+
+ dep.nuclides = ["NucA", "NucB", "NucC"]
+ dep.nuclide_dict = {"NucA" : 0, "NucB" : 1, "NucC" : 2}
+
+ self.assertEqual("NucA", dep.nuc_by_ind("NucA"))
+ self.assertEqual("NucB", dep.nuc_by_ind("NucB"))
+ self.assertEqual("NucC", dep.nuc_by_ind("NucC"))
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_full.py b/tests/deplete_tests/test_full.py
new file mode 100644
index 000000000..f9a6c7493
--- /dev/null
+++ b/tests/deplete_tests/test_full.py
@@ -0,0 +1,119 @@
+""" Full system test suite. """
+
+import shutil
+import unittest
+
+import numpy as np
+
+import opendeplete
+from opendeplete import results
+from opendeplete import utilities
+import test.example_geometry as example_geometry
+
+
+class TestFull(unittest.TestCase):
+ """ Full system test suite.
+
+ Runs an entire OpenMC simulation with depletion coupling and verifies
+ that the outputs match a reference file. Sensitive to changes in
+ OpenMC.
+ """
+
+ def test_full(self):
+ """
+ This test runs a complete OpenMC simulation and tests the outputs.
+ It will take a while.
+ """
+
+ n_rings = 2
+ n_wedges = 4
+
+ # Load geometry from example
+ geometry, lower_left, upper_right = \
+ example_geometry.generate_problem(n_rings=n_rings, n_wedges=n_wedges)
+
+ # Create dt vector for 3 steps with 15 day timesteps
+ dt1 = 15*24*60*60 # 15 days
+ dt2 = 1.5*30*24*60*60 # 1.5 months
+ N = np.floor(dt2/dt1)
+
+ dt = np.repeat([dt1], N)
+
+ # Create settings variable
+ settings = opendeplete.OpenMCSettings()
+
+ settings.chain_file = "chains/chain_simple.xml"
+ settings.openmc_call = "openmc"
+ settings.openmc_npernode = 2
+ settings.particles = 100
+ settings.batches = 100
+ settings.inactive = 40
+ settings.lower_left = lower_left
+ settings.upper_right = upper_right
+ settings.entropy_dimension = [10, 10, 1]
+
+ settings.round_number = True
+ settings.constant_seed = 1
+
+ joule_per_mev = 1.6021766208e-13
+ settings.power = 2.337e15*4*joule_per_mev # MeV/second cm from CASMO
+ settings.dt_vec = dt
+ settings.output_dir = "test_full"
+
+ op = opendeplete.OpenMCOperator(geometry, settings)
+
+ # Perform simulation using the predictor algorithm
+ opendeplete.integrator.predictor(op)
+
+ # Load the files
+ res_test = results.read_results(settings.output_dir + "/results.h5")
+
+ # Load the reference
+ res_old = results.read_results("test/test_reference.h5")
+
+ # Assert same mats
+ for mat in res_old[0].mat_to_ind:
+ self.assertIn(mat, res_test[0].mat_to_ind,
+ msg="Cell " + mat + " not in new results.")
+ for nuc in res_old[0].nuc_to_ind:
+ self.assertIn(nuc, res_test[0].nuc_to_ind,
+ msg="Nuclide " + nuc + " not in new results.")
+
+ for mat in res_test[0].mat_to_ind:
+ self.assertIn(mat, res_old[0].mat_to_ind,
+ msg="Cell " + mat + " not in old results.")
+ for nuc in res_test[0].nuc_to_ind:
+ self.assertIn(nuc, res_old[0].nuc_to_ind,
+ msg="Nuclide " + nuc + " not in old results.")
+
+ for mat in res_test[0].mat_to_ind:
+ for nuc in res_test[0].nuc_to_ind:
+ _, y_test = utilities.evaluate_single_nuclide(res_test, mat, nuc)
+ _, y_old = utilities.evaluate_single_nuclide(res_old, mat, nuc)
+
+ # Test each point
+
+ tol = 1.0e-6
+
+ correct = True
+ for i, ref in enumerate(y_old):
+ if ref != y_test[i]:
+ if ref != 0.0:
+ if np.abs(y_test[i] - ref) / ref > tol:
+ correct = False
+ else:
+ correct = False
+
+ self.assertTrue(correct,
+ msg="Discrepancy in mat " + mat + " and nuc " + nuc
+ + "\n" + str(y_old) + "\n" + str(y_test))
+
+ def tearDown(self):
+ """ Clean up files"""
+ opendeplete.comm.barrier()
+ if opendeplete.comm.rank == 0:
+ shutil.rmtree("test_full", ignore_errors=True)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_integrator.py b/tests/deplete_tests/test_integrator.py
new file mode 100644
index 000000000..7e121ce16
--- /dev/null
+++ b/tests/deplete_tests/test_integrator.py
@@ -0,0 +1,116 @@
+""" Tests for integrator.py """
+
+import copy
+import os
+import unittest
+from unittest.mock import MagicMock
+
+import numpy as np
+
+from opendeplete import integrator, ReactionRates, results, comm
+
+
+class TestIntegrator(unittest.TestCase):
+ """ Tests for integrator.py
+
+ It is worth noting that opendeplete.integrate is extremely complex, to
+ the point I am unsure if it can be reasonably unit-tested. For the time
+ being, it will be left unimplemented and testing will be done via
+ regression (in test_integrator_regression.py)
+ """
+
+ def test_save_results(self):
+ """ Test data save module """
+
+ stages = 3
+
+ np.random.seed(comm.rank)
+
+ # Mock geometry
+ op = MagicMock()
+
+ vol_dict = {}
+ full_burn_dict = {}
+
+ j = 0
+ for i in range(comm.size):
+ vol_dict[str(2*i)] = 1.2
+ vol_dict[str(2*i + 1)] = 1.2
+ full_burn_dict[str(2*i)] = j
+ full_burn_dict[str(2*i + 1)] = j + 1
+ j += 2
+
+ burn_list = [str(i) for i in range(2*comm.rank, 2*comm.rank + 2)]
+ nuc_list = ["na", "nb"]
+
+ op.get_results_info.return_value = vol_dict, nuc_list, burn_list, full_burn_dict
+
+ # Construct x
+ x1 = []
+ x2 = []
+
+ for i in range(stages):
+ x1.append([np.random.rand(2), np.random.rand(2)])
+ x2.append([np.random.rand(2), np.random.rand(2)])
+
+ # Construct r
+ cell_dict = {s:i for i, s in enumerate(burn_list)}
+ r1 = ReactionRates(cell_dict, {"na":0, "nb":1}, {"ra":0, "rb":1})
+ r1.rates = np.random.rand(2, 2, 2)
+
+ rate1 = []
+ rate2 = []
+
+ for i in range(stages):
+ rate1.append(copy.deepcopy(r1))
+ r1.rates = np.random.rand(2, 2, 2)
+ rate2.append(copy.deepcopy(r1))
+ r1.rates = np.random.rand(2, 2, 2)
+
+ # Create global terms
+ eigvl1 = np.random.rand(stages)
+ eigvl2 = np.random.rand(stages)
+ seed1 = [np.random.randint(100) for i in range(stages)]
+ seed2 = [np.random.randint(100) for i in range(stages)]
+
+ eigvl1 = comm.bcast(eigvl1, root=0)
+ eigvl2 = comm.bcast(eigvl2, root=0)
+ seed1 = comm.bcast(seed1, root=0)
+ seed2 = comm.bcast(seed2, root=0)
+
+ t1 = [0.0, 1.0]
+ t2 = [1.0, 2.0]
+
+ integrator.save_results(op, x1, rate1, eigvl1, seed1, t1, 0)
+ integrator.save_results(op, x2, rate2, eigvl2, seed2, t2, 1)
+
+ # Load the files
+ res = results.read_results("results.h5")
+
+ for i in range(stages):
+ for mat_i, mat in enumerate(burn_list):
+
+ for nuc_i, nuc in enumerate(nuc_list):
+ self.assertEqual(res[0][i, mat, nuc], x1[i][mat_i][nuc_i])
+ self.assertEqual(res[1][i, mat, nuc], x2[i][mat_i][nuc_i])
+ np.testing.assert_array_equal(res[0].rates[i][mat, nuc, :],
+ rate1[i][mat, nuc, :])
+ np.testing.assert_array_equal(res[1].rates[i][mat, nuc, :],
+ rate2[i][mat, nuc, :])
+
+ np.testing.assert_array_equal(res[0].k, eigvl1)
+ np.testing.assert_array_equal(res[0].seeds, seed1)
+ np.testing.assert_array_equal(res[0].time, t1)
+
+ np.testing.assert_array_equal(res[1].k, eigvl2)
+ np.testing.assert_array_equal(res[1].seeds, seed2)
+ np.testing.assert_array_equal(res[1].time, t2)
+
+ # Delete files
+ comm.barrier()
+ if comm.rank == 0:
+ os.remove("results.h5")
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_nuclide.py b/tests/deplete_tests/test_nuclide.py
new file mode 100644
index 000000000..c5439b2aa
--- /dev/null
+++ b/tests/deplete_tests/test_nuclide.py
@@ -0,0 +1,121 @@
+""" Tests for nuclide.py. """
+
+import unittest
+import xml.etree.ElementTree as ET
+
+from opendeplete import nuclide
+
+
+class TestNuclide(unittest.TestCase):
+ """ Tests for the nuclide class. """
+
+ def test_n_decay_modes(self):
+ """ Test the decay mode count parameter. """
+
+ nuc = nuclide.Nuclide()
+
+ nuc.decay_modes = [
+ nuclide.DecayTuple("beta1", "a", 0.5),
+ nuclide.DecayTuple("beta2", "b", 0.3),
+ nuclide.DecayTuple("beta3", "c", 0.2)
+ ]
+
+ self.assertEqual(nuc.n_decay_modes, 3)
+
+ def test_n_reaction_paths(self):
+ """ Test the reaction path count parameter. """
+
+ nuc = nuclide.Nuclide()
+
+ nuc.reactions = [
+ nuclide.ReactionTuple("(n,2n)", "a", 0.0, 1.0),
+ nuclide.ReactionTuple("(n,3n)", "b", 0.0, 1.0),
+ nuclide.ReactionTuple("(n,4n)", "c", 0.0, 1.0)
+ ]
+
+ self.assertEqual(nuc.n_reaction_paths, 3)
+
+ def test_xml_read(self):
+ """Test reading nuclide data from an XML element."""
+
+ data = """
+
+
+
+
+
+
+
+
+
+ 0.0253
+
+ Te134 Zr100 Xe138
+ 0.062155 0.0497641 0.0481413
+
+
+
+ """
+
+ element = ET.fromstring(data)
+ u235 = nuclide.Nuclide.xml_read(element)
+
+ self.assertEqual(u235.decay_modes, [
+ nuclide.DecayTuple('sf', 'U235', 7.2e-11),
+ nuclide.DecayTuple('alpha', 'Th231', 1 - 7.2e-11)
+ ])
+ self.assertEqual(u235.reactions, [
+ nuclide.ReactionTuple('(n,2n)', 'U234', -5297781.0, 1.0),
+ nuclide.ReactionTuple('(n,3n)', 'U233', -12142300.0, 1.0),
+ nuclide.ReactionTuple('(n,4n)', 'U232', -17885600.0, 1.0),
+ nuclide.ReactionTuple('(n,gamma)', 'U236', 6545200.0, 1.0),
+ nuclide.ReactionTuple('fission', None, 193405400.0, 1.0),
+ ])
+ self.assertEqual(u235.yield_energies, [0.0253])
+ self.assertEqual(u235.yield_data, {
+ 0.0253: [('Te134', 0.062155), ('Zr100', 0.0497641),
+ ('Xe138', 0.0481413)]
+ })
+
+ def test_xml_write(self):
+ """Test writing nuclide data to an XML element."""
+
+ C = nuclide.Nuclide()
+ C.name = "C"
+ C.half_life = 0.123
+ C.decay_modes = [
+ nuclide.DecayTuple('beta-', 'B', 0.99),
+ nuclide.DecayTuple('alpha', 'D', 0.01)
+ ]
+ C.reactions = [
+ nuclide.ReactionTuple('fission', None, 2.0e8, 1.0),
+ nuclide.ReactionTuple('(n,gamma)', 'A', 0.0, 1.0)
+ ]
+ C.yield_energies = [0.0253]
+ C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]}
+ element = C.xml_write()
+
+ self.assertEqual(element.get("half_life"), "0.123")
+
+ decay_elems = element.findall("decay_type")
+ self.assertEqual(len(decay_elems), 2)
+ self.assertEqual(decay_elems[0].get("type"), "beta-")
+ self.assertEqual(decay_elems[0].get("target"), "B")
+ self.assertEqual(decay_elems[0].get("branching_ratio"), "0.99")
+ self.assertEqual(decay_elems[1].get("type"), "alpha")
+ self.assertEqual(decay_elems[1].get("target"), "D")
+ self.assertEqual(decay_elems[1].get("branching_ratio"), "0.01")
+
+ rx_elems = element.findall("reaction_type")
+ self.assertEqual(len(rx_elems), 2)
+ self.assertEqual(rx_elems[0].get("type"), "fission")
+ self.assertEqual(float(rx_elems[0].get("Q")), 2.0e8)
+ self.assertEqual(rx_elems[1].get("type"), "(n,gamma)")
+ self.assertEqual(rx_elems[1].get("target"), "A")
+ self.assertEqual(float(rx_elems[1].get("Q")), 0.0)
+
+ self.assertIsNotNone(element.find('neutron_fission_yields'))
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_predictor_regression.py b/tests/deplete_tests/test_predictor_regression.py
new file mode 100644
index 000000000..c72ae8a47
--- /dev/null
+++ b/tests/deplete_tests/test_predictor_regression.py
@@ -0,0 +1,68 @@
+""" Regression tests for predictor.py"""
+
+import os
+import unittest
+
+import numpy as np
+
+import opendeplete
+from opendeplete import results
+from opendeplete import utilities
+import test.dummy_geometry as dummy_geometry
+
+class TestPredictorRegression(unittest.TestCase):
+ """ Regression tests for opendeplete.integrator.predictor algorithm.
+
+ These tests integrate a simple test problem described in dummy_geometry.py.
+ """
+
+ @classmethod
+ def setUpClass(cls):
+ """ Save current directory in case integrator crashes."""
+ cls.cwd = os.getcwd()
+ cls.results = "test_integrator_regression"
+
+ def test_predictor(self):
+ """ Integral regression test of integrator algorithm using CE/CM. """
+
+ settings = opendeplete.Settings()
+ settings.dt_vec = [0.75, 0.75]
+ settings.output_dir = self.results
+
+ op = dummy_geometry.DummyGeometry(settings)
+
+ # Perform simulation using the predictor algorithm
+ opendeplete.predictor(op, print_out=False)
+
+ # Load the files
+ res = results.read_results(settings.output_dir + "/results.h5")
+
+ _, y1 = utilities.evaluate_single_nuclide(res, "1", "1")
+ _, y2 = utilities.evaluate_single_nuclide(res, "1", "2")
+
+ # Mathematica solution
+ s1 = [2.46847546272295, 0.986431226850467]
+ s2 = [4.11525874568034, -0.0581692232513460]
+
+ tol = 1.0e-13
+
+ self.assertLess(np.absolute(y1[1] - s1[0]), tol)
+ self.assertLess(np.absolute(y2[1] - s1[1]), tol)
+
+ self.assertLess(np.absolute(y1[2] - s2[0]), tol)
+ self.assertLess(np.absolute(y2[2] - s2[1]), tol)
+
+ @classmethod
+ def tearDownClass(cls):
+ """ Clean up files"""
+
+ os.chdir(cls.cwd)
+
+ opendeplete.comm.barrier()
+ if opendeplete.comm.rank == 0:
+ os.remove(os.path.join(cls.results, "results.h5"))
+ os.rmdir(cls.results)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_reaction_rates.py b/tests/deplete_tests/test_reaction_rates.py
new file mode 100644
index 000000000..4821ec18c
--- /dev/null
+++ b/tests/deplete_tests/test_reaction_rates.py
@@ -0,0 +1,86 @@
+""" Tests for reaction_rates.py. """
+
+import unittest
+
+from opendeplete import reaction_rates
+
+
+class TestReactionRates(unittest.TestCase):
+ """ Tests for the ReactionRates class. """
+
+ def test_indexing(self):
+ """Tests the __getitem__ and __setitem__ routines simultaneously."""
+
+ mat_to_ind = {"10000" : 0, "10001" : 1}
+ nuc_to_ind = {"U238" : 0, "U235" : 1}
+ react_to_ind = {"fission" : 0, "(n,gamma)" : 1}
+
+ rates = reaction_rates.ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
+
+ rates["10000", "U238", "fission"] = 1.0
+ rates["10001", "U238", "fission"] = 2.0
+ rates["10000", "U235", "fission"] = 3.0
+ rates["10001", "U235", "fission"] = 4.0
+ rates["10000", "U238", "(n,gamma)"] = 5.0
+ rates["10001", "U238", "(n,gamma)"] = 6.0
+ rates["10000", "U235", "(n,gamma)"] = 7.0
+ rates["10001", "U235", "(n,gamma)"] = 8.0
+
+ # String indexing
+ self.assertEqual(rates["10000", "U238", "fission"], 1.0)
+ self.assertEqual(rates["10001", "U238", "fission"], 2.0)
+ self.assertEqual(rates["10000", "U235", "fission"], 3.0)
+ self.assertEqual(rates["10001", "U235", "fission"], 4.0)
+ self.assertEqual(rates["10000", "U238", "(n,gamma)"], 5.0)
+ self.assertEqual(rates["10001", "U238", "(n,gamma)"], 6.0)
+ self.assertEqual(rates["10000", "U235", "(n,gamma)"], 7.0)
+ self.assertEqual(rates["10001", "U235", "(n,gamma)"], 8.0)
+
+ # Int indexing
+ self.assertEqual(rates[0, 0, 0], 1.0)
+ self.assertEqual(rates[1, 0, 0], 2.0)
+ self.assertEqual(rates[0, 1, 0], 3.0)
+ self.assertEqual(rates[1, 1, 0], 4.0)
+ self.assertEqual(rates[0, 0, 1], 5.0)
+ self.assertEqual(rates[1, 0, 1], 6.0)
+ self.assertEqual(rates[0, 1, 1], 7.0)
+ self.assertEqual(rates[1, 1, 1], 8.0)
+
+ rates[0, 0, 0] = 5.0
+
+ self.assertEqual(rates[0, 0, 0], 5.0)
+ self.assertEqual(rates["10000", "U238", "fission"], 5.0)
+
+ def test_n_mat(self):
+ """ Test number of materials property. """
+ mat_to_ind = {"10000" : 0, "10001" : 1}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
+ react_to_ind = {"fission" : 0, "(n,gamma)" : 1, "(n,2n)" : 2, "(n,3n)" : 3}
+
+ rates = reaction_rates.ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
+
+ self.assertEqual(rates.n_mat, 2)
+
+ def test_n_nuc(self):
+ """ Test number of nuclides property. """
+ mat_to_ind = {"10000" : 0, "10001" : 1}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
+ react_to_ind = {"fission" : 0, "(n,gamma)" : 1, "(n,2n)" : 2, "(n,3n)" : 3}
+
+ rates = reaction_rates.ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
+
+ self.assertEqual(rates.n_nuc, 3)
+
+ def test_n_react(self):
+ """ Test number of reactions property. """
+ mat_to_ind = {"10000" : 0, "10001" : 1}
+ nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
+ react_to_ind = {"fission" : 0, "(n,gamma)" : 1, "(n,2n)" : 2, "(n,3n)" : 3}
+
+ rates = reaction_rates.ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
+
+ self.assertEqual(rates.n_react, 4)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/deplete_tests/test_reference.h5 b/tests/deplete_tests/test_reference.h5
new file mode 100644
index 000000000..ef3ae0090
Binary files /dev/null and b/tests/deplete_tests/test_reference.h5 differ
diff --git a/tests/deplete_tests/test_utilities.py b/tests/deplete_tests/test_utilities.py
new file mode 100644
index 000000000..faa1a500b
--- /dev/null
+++ b/tests/deplete_tests/test_utilities.py
@@ -0,0 +1,68 @@
+""" Full system test suite. """
+
+import unittest
+
+import numpy as np
+
+from opendeplete import results
+from opendeplete import utilities
+
+
+class TestUtilities(unittest.TestCase):
+ """ Tests the utilities classes.
+
+ This also tests the results read/write code.
+ """
+
+ def test_evaluate_single_nuclide(self):
+ """ Tests evaluating single nuclide utility code.
+ """
+
+ # Load the reference
+ res = results.read_results("test/test_reference.h5")
+
+ x, y = utilities.evaluate_single_nuclide(res, "1", "Xe135")
+
+ x_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
+ y_ref = [6.6747328233649218e+08, 3.5519299354458244e+14,
+ 3.4599104054580338e+14, 3.3821165110278112e+14]
+
+ np.testing.assert_array_equal(x, x_ref)
+ np.testing.assert_array_equal(y, y_ref)
+
+ def test_evaluate_reaction_rate(self):
+ """ Tests evaluating reaction rate utility code.
+ """
+
+ # Load the reference
+ res = results.read_results("test/test_reference.h5")
+
+ x, y = utilities.evaluate_reaction_rate(res, "1", "Xe135", "(n,gamma)")
+
+ x_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
+ xe_ref = np.array([6.6747328233649218e+08, 3.5519299354458244e+14,
+ 3.4599104054580338e+14, 3.3821165110278112e+14])
+ r_ref = np.array([4.0643598574337784e-05, 4.1457730544386974e-05,
+ 3.4121248544056681e-05, 3.9204686657643301e-05])
+
+ np.testing.assert_array_equal(x, x_ref)
+ np.testing.assert_array_equal(y, xe_ref * r_ref)
+
+ def test_evaluate_eigenvalue(self):
+ """ Tests evaluating eigenvalue
+ """
+
+ # Load the reference
+ res = results.read_results("test/test_reference.h5")
+
+ x, y = utilities.evaluate_eigenvalue(res)
+
+ x_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
+ y_ref = [1.1921986054449838, 1.1712785643938586, 1.1927099024502694, 1.2269183590698847]
+
+ np.testing.assert_array_equal(x, x_ref)
+ np.testing.assert_array_equal(y, y_ref)
+
+
+if __name__ == '__main__':
+ unittest.main()