From f8764416d2c727c5b1693c297c3f7994c8314d1b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 9 Feb 2018 13:29:04 -0600 Subject: [PATCH] Copy OpenDeplete files from commit 2d804c227a --- chains/chain_simple.xml | 49 + chains/chain_test.xml | 23 + docs/source/pythonapi/deplete/index.rst | 56 ++ .../pythonapi/deplete/integrator.CRAM16.rst | 6 + .../pythonapi/deplete/integrator.CRAM48.rst | 6 + .../pythonapi/deplete/integrator.cecm.rst | 6 + .../deplete/integrator.predictor.rst | 6 + .../deplete/integrator.save_results.rst | 6 + .../deplete/opendeplete.Concentrations.rst | 30 + .../deplete/opendeplete.ReactionRates.rst | 30 + .../pythonapi/deplete/opendeplete.Results.rst | 22 + openmc/deplete/__init__.py | 24 + openmc/deplete/atom_number.py | 236 +++++ openmc/deplete/depletion_chain.py | 472 ++++++++++ openmc/deplete/dummy_comm.py | 27 + openmc/deplete/function.py | 114 +++ openmc/deplete/integrator/__init__.py | 11 + openmc/deplete/integrator/cecm.py | 133 +++ openmc/deplete/integrator/cram.py | 185 ++++ openmc/deplete/integrator/predictor.py | 100 ++ openmc/deplete/integrator/save_results.py | 46 + openmc/deplete/nuclide.py | 178 ++++ openmc/deplete/openmc_wrapper.py | 853 ++++++++++++++++++ openmc/deplete/reaction_rates.py | 113 +++ openmc/deplete/results.py | 454 ++++++++++ openmc/deplete/utilities.py | 98 ++ scripts/example_geometry.py | 358 ++++++++ scripts/example_plot.py | 46 + scripts/example_run.py | 39 + scripts/make_chain.py | 60 ++ tests/deplete_tests/__init__.py | 0 tests/deplete_tests/dummy_geometry.py | 165 ++++ tests/deplete_tests/example_geometry.py | 1 + tests/deplete_tests/test_atom_number.py | 180 ++++ tests/deplete_tests/test_cecm_regression.py | 69 ++ tests/deplete_tests/test_cram.py | 48 + tests/deplete_tests/test_depletion_chain.py | 197 ++++ tests/deplete_tests/test_full.py | 119 +++ tests/deplete_tests/test_integrator.py | 116 +++ tests/deplete_tests/test_nuclide.py | 121 +++ .../test_predictor_regression.py | 68 ++ tests/deplete_tests/test_reaction_rates.py | 86 ++ tests/deplete_tests/test_reference.h5 | Bin 0 -> 165384 bytes tests/deplete_tests/test_utilities.py | 68 ++ 44 files changed, 5025 insertions(+) create mode 100644 chains/chain_simple.xml create mode 100644 chains/chain_test.xml create mode 100644 docs/source/pythonapi/deplete/index.rst create mode 100644 docs/source/pythonapi/deplete/integrator.CRAM16.rst create mode 100644 docs/source/pythonapi/deplete/integrator.CRAM48.rst create mode 100644 docs/source/pythonapi/deplete/integrator.cecm.rst create mode 100644 docs/source/pythonapi/deplete/integrator.predictor.rst create mode 100644 docs/source/pythonapi/deplete/integrator.save_results.rst create mode 100644 docs/source/pythonapi/deplete/opendeplete.Concentrations.rst create mode 100644 docs/source/pythonapi/deplete/opendeplete.ReactionRates.rst create mode 100644 docs/source/pythonapi/deplete/opendeplete.Results.rst create mode 100644 openmc/deplete/__init__.py create mode 100644 openmc/deplete/atom_number.py create mode 100644 openmc/deplete/depletion_chain.py create mode 100644 openmc/deplete/dummy_comm.py create mode 100644 openmc/deplete/function.py create mode 100644 openmc/deplete/integrator/__init__.py create mode 100644 openmc/deplete/integrator/cecm.py create mode 100644 openmc/deplete/integrator/cram.py create mode 100644 openmc/deplete/integrator/predictor.py create mode 100644 openmc/deplete/integrator/save_results.py create mode 100644 openmc/deplete/nuclide.py create mode 100644 openmc/deplete/openmc_wrapper.py create mode 100644 openmc/deplete/reaction_rates.py create mode 100644 openmc/deplete/results.py create mode 100644 openmc/deplete/utilities.py create mode 100644 scripts/example_geometry.py create mode 100644 scripts/example_plot.py create mode 100644 scripts/example_run.py create mode 100644 scripts/make_chain.py create mode 100644 tests/deplete_tests/__init__.py create mode 100644 tests/deplete_tests/dummy_geometry.py create mode 120000 tests/deplete_tests/example_geometry.py create mode 100644 tests/deplete_tests/test_atom_number.py create mode 100644 tests/deplete_tests/test_cecm_regression.py create mode 100644 tests/deplete_tests/test_cram.py create mode 100644 tests/deplete_tests/test_depletion_chain.py create mode 100644 tests/deplete_tests/test_full.py create mode 100644 tests/deplete_tests/test_integrator.py create mode 100644 tests/deplete_tests/test_nuclide.py create mode 100644 tests/deplete_tests/test_predictor_regression.py create mode 100644 tests/deplete_tests/test_reaction_rates.py create mode 100644 tests/deplete_tests/test_reference.h5 create mode 100644 tests/deplete_tests/test_utilities.py diff --git a/chains/chain_simple.xml b/chains/chain_simple.xml new file mode 100644 index 0000000000..345da2237d --- /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 0000000000..5985704063 --- /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 0000000000..55380c7a1c --- /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 0000000000..f9eba273ed --- /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 0000000000..d7467a418a --- /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 0000000000..507a638f69 --- /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 0000000000..d6c0fd827c --- /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 0000000000..5c21dcb664 --- /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 0000000000..6fa07a970b --- /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 0000000000..99e048b565 --- /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 0000000000..0ab8a1f711 --- /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 0000000000..994a51e12c --- /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 0000000000..03bedbf531 --- /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 0000000000..05cc9db435 --- /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 0000000000..b3fa272648 --- /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 0000000000..74eb92422b --- /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 0000000000..607650dc69 --- /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 0000000000..4d9baebb2e --- /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 0000000000..a18d8450c3 --- /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 0000000000..6c9d538fd6 --- /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 0000000000..35cbc7f3f1 --- /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 0000000000..1208a9b3c3 --- /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 0000000000..347dc71857 --- /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 0000000000..7b934027a5 --- /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 0000000000..c0ec1627ea --- /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 0000000000..54632ed9c2 --- /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 0000000000..9afcc0d464 --- /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 0000000000..d2c6ee9d6a --- /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 0000000000..bb80f65820 --- /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 0000000000..2e0d9d3bc4 --- /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 0000000000..e69de29bb2 diff --git a/tests/deplete_tests/dummy_geometry.py b/tests/deplete_tests/dummy_geometry.py new file mode 100644 index 0000000000..6101519411 --- /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 0000000000..1071aabc05 --- /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 0000000000..9a17230f86 --- /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 0000000000..23a6342000 --- /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 0000000000..2744adbf48 --- /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 0000000000..216d1e68f7 --- /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 0000000000..f9a6c7493d --- /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 0000000000..7e121ce167 --- /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 0000000000..c5439b2aa3 --- /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 0000000000..c72ae8a475 --- /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 0000000000..4821ec18cd --- /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 0000000000000000000000000000000000000000..ef3ae0090943bc7ecdc01a1afaa70c0216e029f0 GIT binary patch literal 165384 zcmeEP2|QHY`yXo(5<+PakxKS8?ltZmOC%~qJC&lclonD_+O=p=+9Z{fii*%WC}|TF zX%{0zWi2i8ztcT;&TD?}^uBuY)9?S@>60_hxzByhbH3;KKIhz-bMI`hXW2?i_LgLD zUlI}wQHJdIx743j@RC?7{Jn~jspB5tjSwi;gEE0sX9#`&$6zRf`X0bOzn<$D8yf~g zm_ga6N^lxOPn3LT1}aqZ$QC6i1-kryjexz4wF|d{$)J}3po~tWA`evj;zEcaPDC*A z0?i09$cUp_6(Qo8`(Bo)CXX<=+6*y5;?@fb34d3WAQ-@XBSIMf+FV`kOYRSLM;eDt zU@#)d1Hakdv7?+>LVw5-L1LtX-#baXjRi85M!pXkQEz{e()jZZQTZXiAE`%aHv?ID z09nK1|0qi1&+Yy0wSw{@KYnChlogm4W#!49dVmfmW8|AKAg7Ne-^c|EF4&)#} zq+%e4tsvjj0(t!w@{KHjJMsW2t5V1H$T#lbco>_MwZQS9Gvu3OaJ=|2DHj9TrG}Ji zft>q|lx4yEkuOaipvJwp4=LM|YXADbLPVzFBBBhCkOjnZS0F2r{{kTZGh^{K$6{2z z!3qpRo3p}l46WN#&;pGBjR1`RjR1`Rjllmj0@QUC8!QCuK^c9L?3SSQg|abuoLXN9 z((?cdH%gBJi(J<-1g>i-J-YmNB0$BZJLo=nZOve$%eBRWTCcjat^R3=y1u2_V|O7< zQSAiDl5#ZggKAe5)Khw42BcoB5m{3Gu;_GjLuPh5~WkV*h!Ie7zWq z9;#O0lM^Pxs2jtmRfIf)I{$6SQ@4qVXRP2A5@oS%DJXyCDv^50K%LH8ZE>Tewg?Ey zTLB=BC~Yb}sr*6J|7=O|O|K_u!3DI}UcOO&br#>uR7n46zWr<@{JYK*#J7^(ZE>LD z`d9NUR+(&v4ccoj->A>uI*V@!{Yd|5zO{`PTKaPY1o6!U#1VB?NAXRz57~kqXs^9| z1It5V?b$i&No6(Cf0}Ppq|oJ$5D>&S_JFo$brj#i`jhR%g7(_WH?VD^)A+VtgY=)~ zTibY{r9VeN5a0Ab98qU=6yIvW{2&Y3YcJo#fPATVnG~H z=XMm|JT%D`!a#fN<(oKgu(SB)txY;W^NorWy8ICWg7_v2;)pu8qxhCTm~5eTC|S0b zZ{X{gPFqjr>yv)de4`?TE`Nl8Aijm^wmqw(_-5}$w&P9gDCWM}%QvvSw8QyEeecQy zKs3jZ{!+N89wuexQL?lL2fX9k>O}%MIf1Nq0eZ#9+v+h+kg_V+=RoaWjhxt4&pnBh z(?LBc{l0&x?}3=ZNFO~woV9n{B|sNDdR;?*pZfd8((mAzp!~@O`zk0OsXXCARs6$} z;t~Hoj`s4XEAWHn(LWpz`t*N%y?pBiy46{HV-h>|xG$P-Bo6f7A0r@$Z%+)`p={q)j|pTC2cYnq7u5AFKMw89 zThgGPot3w|?MMgcy!D$0+;8HZpu8mu;)rsS;=c=3@efOiZ~XbBy?m1aesmV!ia~y( z`SuS3;$L#QAijn1&66F)H~x9Ky?pE5$$YaPMLIz9?O#Hef7neye5-Y8Yq+EMmh4Ej zlMmW!FW+Q=o1HbkF~*V((0pqfFSPXM2ngbv2YBulb=I%u8+9GL9yG8X@Q!MSpNEk^ zuNllIl%A?Hd4LK0jRtzjfCrRbF_4?bkqQ(~CBXfrlpg>6r@-8ATTttAejM7Hw|amc zc2?dB7*9Gt=dIs1)P5_+1m!Jd5J!}=6t7*Viho#AeB;k2?d2Qw7?955TQS&IMf2?+ z2E@PQbU}Pu5B6nJFmx2(`1d=umv3O(cBk#XV4o!sqWC5V&tkrPn!g%PYMiK2&~sxL zQ``DW?{lOOqsu=-Ku|mgNAZm@n`}W9wAWs~ zQIDPMEWSzjlK#_tqauYae}sS_z9r9VdsavBEnp7WP9$iry?j#!Zgv*m*z-vTXuh?L z7h3vr1O)L-HK48jujU(d9n1#J6ocy^svZ7$FBj<91AbC^{Cv#^dg)-^q4Y|C%mDo= z26_Pt$TyT;G?3F50)^kapw{L5{cCUD>IM4QS$Ruu5$OP(w|?`0`%T;vl(&lccuVox zc+?k^a+sYa1`L z^ydf&;+xirw)(%CZ`5^g323GSTnAC@@XvddK+gs6lhWhot1{3F0&gfi{@*L91$t~S zKT&#>K$hU^KT!Ye3Tj=>?{|CiR$tJw&RQp!ttK6y^VV+{V!x$Bg7Q{Lcw2WVUc0od z{%MKgTkJ~G1RJ#1UcU7M9&{Gp64sFp(0u#pQ21A!D2Q(^AdV<7I*Ml98x<*Z`6C1b@r}K)?O7egw=i&nWh`i~y?j#zZgv*m^0$%>(0pqf zFSPXM2ngbv-j=rdznX8m^E$e}0qz2iyUVDLsB( z%Yp;BV1A$V zHN|U}w$(o^QG82?Bu(Uk_S(xgHQ+&K@hy25=>W~QpALn8)ro@m<__YB0;8k&ro5AE z!3?z5UcRXV2Rn;zs(VQXXueUALYF^6KoH+@_q08$qxiOdH`z`CXs^9|8wlL&EWQ=T zkPgs%Ya1`L^ydf&;+xrlw)(%CZ{&5bI71m24g>2-svZ7yrYbm43$B+aJ^uMo3motO zJf`&cd94QyZ4CMSnWW6eYs}gL$qV#0J1`$2#3m~BP zs|spe&hK}7^Ogp1sI&5xeH`fkowt5fWc>fOBPegx0zOi`r+DZ>Rs6$};u}Bz+siji z;74ciE&V9z0L{057!d!G(*^O(1H=&pLr3wA|94T_%eO(m!Or5F)^XASnr~F3(B+R1 z5X86qgtljO6yNxNU%0(|(*kaG7T-!HkPgs%Ya1`L^ydf&;v4g1Tm4_nH_Ey`cw+!c zLr@xl(ioH`pfm-g87R#`sjo-A(FX?gfkAy>P#+l72L|r0rpIQJriKh1lThH_Dq016JXC2*fRz8Oo2TrE`vZE{I?4FIbQ!Qr2jYD zAwTC6V?=@gb-zjWgNoNMAjjS&>!ZPe+&gXcYJn_!kF1XcdS>};^*n$acAuC*kJ{=nf$Z^^tnaM9)aMTFNu&qq6Unl@>#`xB1D&-l(>qT(kPQ5w z^9KNeF8&k&LF=;O)VAt?z6_VP)jus!e2YFunn(xjwU=*0fd`$%x7ZBQ0h(_=9SZ-d z69w_j9>fs^Mn~~YB8_Z83$)i>zG(voJBx3!mq`a`zEP1vmp?*45Z}@-wLPn&_!e}5 zY$qDD*IvF218#N}-*Pia2WY;vjTc(_a|8tOO^e-D|5x*kybkWd@BoGbz`BxZhku>v z4Gxrm>m^E$e|`)A2V4M;DLsB(2Y~}@Fh5axl|YsN{SN~Ng0e_FD7{!9v#*i$>%oD_ z>uvR9b4Xb)mn^C80o=h2nUuXSAjb|M_5SW(YF*CncYE`e4rs8m@|GFEjm}$tcjW)w zdV=y+3E(4z4#h7Qs^TA(6yNyy-(J3fpONdd=OiV7d_?o@9|pv~Qnr~F3(B+R15X3h&_+E%Qx1;#R|2xAJ@{UX z<{K3$bonC$1o2G|#1VCFNAZpS``q^O4QzYpG`=;1@0V%5QISHIKSDqd-(o==QRj9P z-}v`Mw3ly&z`@R1PkI-T4$yp~B84u0gn%Hv$$~hd&h03^@$WBcFW|4dnH&$ogDx zAi1=yUNMlHIb?l4&{KWgR?i;D-j!s1G0=;wYO9wH!YvpR}z{OjfR@(tYn z)@ghz0qb{~Z&akv<&O{$#JBYzj;M1xif{bao$ckDIdHJE)|1h1Ne5`YQISHIKSDqd z-i_(^t!fsB|v60 zlJ&JfPphe|o(qrzK9cpn`(NsFC;rdD+MBnAgN}Aq-jV?0L+7pEJ>q^3{{-c&WRM3a zFR6UuLRI|3lHwbG-f1u2EPx-K#kU}k2Wh_j!+`jgoGyrOsvwRi7&?k?{PTBv`DWS4 ze5(ZO37T(Iq|oJ$5D>(-NWOJxNAZn+9o$~NA)U-O_b;RaG~a#{W%Mx$0YQ9YfH(Ng42%n%w|=(>`#t>;l()1%98q3VEOwzP{$WY+jX$5Xmv2nqM`!V^ z7UVaYZ~rhL{w1di;#)M|JlRovp;~>#wo<9W)w2+CV}AdaXWP`v(Cd5hv3e?Dn1-&mc@x8|NCGBn?) zNTJIgAs~ovu^^79b32M}9`dB+FwkCm`DOzg?5z3CyBEzjDn{t?uMrT$H(3xzl!cDs zTfQP`xmKAh+sijw;9zI*O|>7*H!4Qx@~;sP#J8|M6vIdxl=Q3ClhpGs^uTqi8JOp& zcKGK#GoY6M_(|#U^OXtolzWpJlpYhv9-v=L08PT+w*AWoaxM4)9DAUr41O?xXAj68 z)B}Qn-haEG*5&;DYj57N1CDom-V*x$&Ac*bkMh46oQ42Fv)sX^%l0XbTWtf6?8JGiZ0Es$lk$$C}bICEHAJ#Qed*Cp$J zwtxSv^Qd`?A1^5elOaktA;gd)?){fwcrNt!RrOosJKbx6s*m4HuOQ%{EHMD&7%;!r zP`wrb$DtjC?{*e=dCpkq47++ z&ZXk;SM!d4K8JkyIo@&SO708V`8nQ!%YMSr-<$6kdL%+L@4mxN|M?391o19`Z@t=4 zyyIUVL%#eR@2KuVJ3q%eYQ6e*c^70z^NzxvF8>+jBy5r-+$7I zD8K(i3{jeK&vo}A6}Tn$Jjy>0?!gXtQTwH|fbDc}T}rjff1Roa^z^_wn$qK6SAXve zp_>5aeM*mi9l!+kR6(Ai^!W3ZJp5fZo`%4`IQtMY#tPrA9^0|EX#stVk)2m72TFSY!%XM6hkdV+0^ zj9ET`fyA}~UjCx;0j2*}CB4t-uMqloALoxtD!-N3x19$l=~vC~*nZ=X{v z-~}7_LFLU_u=41O(+BcRn6dto~JbhvEx=-e@mhD8C>R_!E7J+plQUjbYRZzW5a*YPrXS$?evJr21VJ z5f;;PCXI6ca{qIi<~I9%3XlWVxCh1?hfZJ+`zC+Ii<=j?IfUvLRWcRYXFX^N={$8F zDgCt#?|W^}qpl17>P&knQj=;Aey{%HGqpyvZGk$E@|S=)f04^4FiY{jc!-KlkOhfB*bcNQ~T>@mKqDs9hb@{uM#{IVipM@6*u& z`*g%}7m&uef4TpueL7szyzL_YxnD=?dw>4!A3^(ae)ji9l*37Ppu+ZL(wseW(=f%~HNX;Ip|0U;05e`y401ZV_k1ZV_k z1ZV_k1ZV_k1ZV_k1ZV{QUn7um!>D9Myf8|2a`B6x;HIs;c(VHu-EXlTSv=Y5+y*VB z21TBn*j;Fy#vpyvcjxr2UXmPS)^jyO6Hj$?f89 zAN>L!C1)`!s=WY}RVU8Y>_8~l9WiCto$bN1zcOyvi-KG!o;={u(^(r2UggQAV^$x0 z|4M-;e{qaIr#O8Gnm;|W>kz$iL{5LwhoUKJXm(MN*02W{YIg9w#0}SKM0UsQmm@#1 z5siph5f_#T^Zb#LXC{pqxt-_FNTtkUo4&U4+VC;eUKm;=?|aH|>^o%c%TvvZo!N-K zc=G!MQRvTm1BSYN-8NqPfr|pqelGdMlWR6iy7O&WU!Fffx4jkY%h)_QY@Xsc19LcD z%EO*#Sp++w0eMSzmU)&VTehuncwnoB-pRU`v;7E;o=hyU`JDF}so-3*8|9gagdDvg zJuy#&htDwf!qQ?_I9_I39KW=RHuL=1p?Jw-%6iy;jU^fy*#V1q{^VYEpAk0<;(2Mx z*}k5b5&G3|-_)~OWyrRx9-~7VmqS zjJu;CUs%mx@Dt6iF?8v16wjYoHN9F6bp%y9^3>-;H)d{3HjAO ze8#-J_7XgQWOVv@EJWpb^0d00)>jlEK1Wq{x#cOr_(2bw+XR2nL=}$e?mr|{jugfp z8f_P>j#|!`x%*fe3vKfAoEo{i8tL21aN>tIY(#EA*3SJd(4XM3hGI_#L3}=>xt8~} zhWTgdMUV8Nvk)(dGrf*Ah+N|NH+jNn`H!n0zb3IX2RtxE&{Ubp9-|g<5M0&wVb~-! z^ovlx>(j?$sQsD(c^QQoWJvCFwNw6VBm>{HMqd*0YmCO{8PW|fAKKVvPEWrh!fWsG z>%IEDXF8+^>c?vAF-)Ud%>F+tnN25iltHK7L7o`dY zHC)8dGjEN(iWI64L(N8C{0kd7rytVGwi)`<%c$4KgW@orv#s3jEbwdK@on3M)Y&Ry z;d~K(vLV1W2j)Y&1H;`oJ}~|-gamE;QZ^7RVkmm|IA4n7AAjXJ-$)TXuOohpxlgm}~Uin4Fvfp)s@$ zx4BO~LVO}u-qy(IZ;bB#cvy7r`!a-uJ!H)@QAaN{3d`N_wL$%{THon2YLI6~S}(+! zvXPz4lWz(UaUS2VjO{-${Vv4kef(f`p^Gpc@(awzt;vIYU;QrCB4`!#Cw$#*e{ols z4;vcK)rIb1qMy36x4TU)M{3tJ<@NY75G`K+pvS3f8`N}`hKf<;dt`6x*1i30G7-P( z5jXdohkRf6lCd$YZxku!an7Q?+m5u^SxA2%d@WRWhxk2cX>{oQ0yGy z{1x&`;#%ED`PXo~W?(rlBi2BFR*VfPX_1lR`7>s<(!I%6knfQbB*O}YVSYZhO?tk} z9t(77Rp0~VFb+}>we_Q%hdOE&5*Q`c9Y?jh9Ng^HuLk+jHL|iioQawOQsW199THPp4mOFdl?M+3$xkHOzpB9mJJ!{Ukc#@E$1x2$Ht@lsu!-FF}d z`r|cJuR-q|%zskNm!j(~LVOG>3^+OUFh6hDEu=b62#(jFZ~L+=Rtd zj9zfQc$szT)7;x|J#SKXA=q_JAP=uJ-PV{I>o2@`bN5QS7bByKR^CF~A22wGOwQ>d zuQl>$V(i7niK#=;wQh5ksRUIbvi6s^Z>nJ<3+7!taLf4y&)xu!ae={agL(GUXAF|i z=D>WTXYZDJ{ttX!Wi(~>>>+S|zj7ugPFM%#gUc1;GN7D|Bj8V1aAbajJibpM1 zM+f@pdpW$q(fdO$%`0fAK^nVmmos%^BYV9AWro^UqJ&hPdvxpNYIVSaXt z+dpRHcp0Amhn|Fwe&YcBNi-9>lJ*tG&zKd~Vb%LAP@Q?r$9%_dkXVzXp|%^<&|@d< zTN_{E=*@(e$KtwIAqg&T2cX}sBG2BGI;y{c{zQeG=+0Dw{yZAm!{?>}TyN|Q9MU~o z4f6dP=YY~t3z&c8*B@W%mIm{og3|lTYi}5!xi=K(0DG_TRd3$u)tEte*3!pQm{jUVi@4QaMoG0_MZ2^DiuGbKv}Sv#x2&=kBoozSa(* zyBEWJc=uSLLs%pnugaC;JCbk<^jd@HqshZKh@xWN+!411ptYrb51M=6=z^K6Gt;n2 zWMD$o#n*aQk?nWI8f)jn@tW!|H2n1q$oF`b@g~X1Fka-Yo}H0x(v645oO44TS>A+v zx7g$&;_x2MC$WoP>|XN92tCl6Shki}ALfK#ynNnV4b9TLb}&jBLk-Ux9az4!3Naej z!Z|Z96Pc!QO5W`;oKGDb_wDbx9P<6S^pX>y`X6|F6PhghxK}OYd;gd#$d@$`AIqp6 zD|@eo_@tPY`C3e}Ky|xpiA;6mAl(eQj%H0(L-z-vuN7=?)N17&R&SwNPkBjHaK&YaPyd7C9?s0_!o#EL^nTwLM1>P% z^K+HO=aN_^A+lA%I<>Uj7(nHbeo<;ZIA=)O)SCrgdHPmmxh5n4Hw}Qb2pF+u<7fwbAXfUTPcURv=c9NgqELvypou&mv4dzUrpfV+haQ$#2mE?Of7Fn-3wk7XVaf#YR< z?Bw>?H!wdt_D)&qq6+&jQnsLKwj{*ok_JmoE)@DR>8ZKS)YJNCy*FOk9Lqrl%|AI? z`;{6xGWwHL&uKV%@m$223jHc%d(J)g;%YYX8huzeO&9u8;c`*jG8f{r^h4x`=h<-m zvvZfmg*+jcf4V1xF)nw3`T2r;b;RjYFrF3amgb%6qmAzIQC$6E5C`!!i?ti3tB%I3 zJ{maJoP|p5DRX4*c#B-IzK|HEoryRopN+rd1mpR}_+6=nS0O(AKdm%=TLba&LvC9O z<-_#``#`v^M-BAHQhS2#j40^urLOqNC4@&`txki z@>w5N9mX3F1F;(G4+U(bB_hP^VmG*ce#BI9QdtJ)i_bfD7=O)x^H*r{u9QBL;duEp zJlVWB4X($IYE+$GS^J58)u_iQB?AZVW(YC76aC-h!cO&uvHS zwQE3%ZR(|tyR(rkK6T4p2Ez4Eb<_Cy7e~SIdUy6IqVy4d-m_5mnNrVgknj2vto`fm z!}&|u#bURtEu6oQ6WVS`s}VF*)8=V_6$hy~+pj)UbRg<@^SKL5A1`WIMp&Jh4BVh8tzA7W(dgHL6 zR8Z9jcs!zyQ)yAGIhuFj_SMClGGt#Lp;s8O{^@F0v(ph}p(g(KopUvs5L4%v0bhHw z5jz!MH%$W=KcW-H$2EJy_=&mOWAQsDh>y=E>&+Iba6R2pK7TH2Tof-K*u6qrv;TXr7`*UxL*uk`X=2iNzzh6snfItj;Xz&4e|B?3x$;zWlsC*jhCAd64>-g8Ghjf3n-zJnWBd_`~(Zr4%b{Pqryq-B)Y)*ZUkqsPb-AqSrul z+=;znLC!Yl^};cy7~>m|b*(3mPp+BBu+>wBci#l@(en>Wh&yXZF|{lDD~pSGL0-Z{O^-#Z8LeeY=Rw`#ZGd~tbVmj@maO}uzYF*&pI zi5K+G+|zhUi}@a2{O0-eUCVg(GM405HS31++FvZpwj8Mm<5xv_Q`lN< zDDVDyswwvgjOXfSLc>P;!2H7;IOp6qgQ4iq$q^GL5TDO49QifgYPbe^aLw&7p_eRl zf=9*9z%O-3-K>w<`;#(}`m!+j%g^9^A>vvvGE^G!%VE~MCgvZmmo@B%YQLz3*E@P8 z>za!n!SSlE&Fgiz7LHepifCGVg&BIrg*oEf?J~r@QSr-Bdv%nfI`d)pI~;Xcyyd8o z@H^znRiXZe>#iUZTe^#HcnQbLP(>nK%n|ZyVVbO@@g+E)Hzrgjo*D}C&mNV#S^HnY z{4+5?S*db5JpNj5)+WP6hUlL2ehbvcbC98%mk+|28fbTkWZxxW_UNM>l9%pRHXyjC z$2xnpOk`ey?mp|uU3u51H+BWxsXYes;WCq6?umYoUl$n*@0vJ5zK?v7x^R%vHQsnL zMVup7Plfq8d*ZUByNXQoU6pf0@Tqb{N%Ns8+eQ=hKltiO;VWBoO}xb>iNZRh)@A8e zvE)p|aM{koM{OX#mQ1|tx~UcN>)h-N%|4-!Pp0#}wJ=p6-uBlLoITxr$7+ zo}oItFN_~e@0&;V*~9UY6xsN2#R@oox%ebKj~xfEhh3K@RUM6n__S;)H=p|y=0iLC zdzR1LEYVwjYc^X==OCv|cW&vcrHSe(J+R-r(H2!sn|n#m?E?~fF{odv1{;~JvuaMV z37jupzKv*++6?hY8eDMx=^s8ntVWLv)3kus%jw@fdM$ee@tM?YpI=!8`895v{NQb- z*61iTb9pVi9Chuz(1ASya^=5&7(8}x0;M%9ydYmq5Mr8O2;GLg{-CNCep6Xs{l z#o-U+@4)$8UF*v8%U|L8%g$K$Sll1JAG|wU^hECiaQzcfuTWbb0^`f)icahnClmBg z#hD!+kC!1K9}a#@U=Bd-&9cuYFT&8xU&?N@nARdu-8T>TDslxmHd|}aQ+>F8E_kJW z{PrvOd4gZXIKB_Nn39mi?^K)sg^{)HNpg&uy-FmE?0R6FI@AOvasfXSu+%ROI1_w!4 z{A^~IS;TrQ{-sd!3LDhAW%U8SGj+(;#c3f0k!&P9f16KVHuPsiDK3@thwls2iwcX5 z4T1bVk*yoFXfoupXm5__{aEPFgR72)9AD^9$e1NglFRkb%gFPdJx6m8&E1#2y7_CM zx!3ya-D+=(HkRf`UYOW~NH~vA&{>{|2zLqMNZf?}ywviWh#EqEP0l}Rc`yv}JtA*e zt@j2PucLdnEL_Zi{>TkR-=PxFpQ@r-7t38{sE)f(UD!-we@ktte)F4w=yY@THSbL} z=$xIeQ&Wljyrz$W{;(uA@^t+9Z=nW|U!R2^7Y-_cypwf4t zwsIGV{Stnw=M1<5{fV)@vZd7l`qLVi8xs)*_lJDZ*yp@)CA^+b*Y192!f}XCQt{ee zrZ1pBtM52{JF^c*Q${o(!~=@hBAflNML|zSj%rk7tR)ipc@OynNusI@NS%5ad^R<(0ba!Vn*(-sNS5e$b!! zN#iHqPq9Gfw_2@N8&QUcrzY1dHd069XH?uTEN7t(TZI*+Pt+mLPV1N39=d{rZ`!7# zlMClBzn-oSzKnHCbX+AirwoZ5<1}^OU1d~7Vc3CrpA1mcPV&x~$L|qMxwZPnEH?5mO)}i|)hnL8 z45h^7=Rb$>@?q=L0MqCuXurm3I?`h^?Ej1&(;mpR!1bZ(%^L>>Pld-F9~)dvd2fJ* z-rOpWsViR;7GEPqWj%JM*;(Ka?{{R6E3p@;^=Fshg8utGMn`)i@2oI1ouIz&4} zX*9%V&*>1=1sfqgX8mt8IQYZ-d{Xi0!d^LWywc5*TUK_1`0ZHNH)zrYGqeg%JE?oP z43XSb_ANa{9X+!o!Eu?J4La-Iu8}KG)*?${RR*6PaurGX`a0G7ER3H=Jy>rmpF)3H z<#%&rkHYqsR_;93*bA;#&wY!CpPvH#i3s-DF6{&TsTiU$t2xmSz440WWBrDMI9WF~ zX&xSk?%0`C6&7a{1+Qn~1&EoRUkYTw^tuLa+KKV^p3 zym^5{rQf0%^L%D{ao>J_E%5T3v5@$I6{^B-k?(Y`1q=pNue!i=UNgDpEcEg9;u@mv3uJPIpLqS}{ix*qn?l5c+qp0P`%(GtJ8l2{rQ~ztMTAA% zPmxV<|8oC-zlHHX?#-ugKLJTnm6PQV#FK7iwDrf;cvjK1s%&n3j ztZ-jiwq)t?U1BgI(|&A^S|t}}{aBA$52^}UP=fY4s|Z+>vh(~75qx&Ev&dz=Mr_ED z+xsH+dEf?OTeFs(W#Z4ye(TC{O2eM&Cz_=mD8sV*-n($)Kpg&wGc#s}sSJzt*;m|x z_`S#HQfb>1J__Nj@yzIfj*Xb-eVyqBC#K+CVs5SrV>5BfC6OiXKc!=PpUxhs8c~XU z={t1Y1;WnY+>K$^tYlbui^s^C6Lt#BuR2cKEQFtz)_Po^(TH_@Vrts-crrd3-MdmZ zmWfMF+xh7BrgY3T{_fWd$x@6h|IzU}VaGGSad5V^3`^ekQeUn=#uhtf9m^2LHD5l< zIAqm`wTKV6W;c2o-hH_4({mS?_}CGlwH}AlvG<#-9W^Y=u#ky;9n%RrtZJ#j^%gR$ z&<&A!ri7jPUhn#JwXMN+GM5L5#?)g@;dA;MbhpGEmpwlqFNxyjRdcU-Mit@y3M1sx zv>sp&t$k0*5Wf#OELD1b&=YA^xXNU;6yo6dmaw+%zV&*Nm^QssJlEF zUl}ju@N6;@|6CBv)~w9HrdF!F zMDdGiTYOIKZp1<+R97AC>WQx|TO_kU!W!@Wwq`;6+6*k#*|$h-AqR^Wk>R6m z<*05f!?IoayqgJOM}O$PuFi78_{&qfEnIY)um-Usx%O+O;QOa-?XhmDHU3IjvviZh z1?&GM!e3wt~=IQl~m~GZ>CH1Wycvwo2#6uTr{P8F=+mxC# zY}bYyNss4c*uJ!|m`jA6v)OUS-lAJmzG5{~nu) zaL+ycDq_Aoz?c{89ps5Pd~~EaZrCGfmfbap9;rkeo_zM!b7+tVUUf*`XX5lG?CceM zl!2)SzVCLXOyvn{JZMmvocZz#*o~Y$GMm4aVPf6BZs+DK-^`Z!S&i83x3l_M>P*2U%>!1e7BTVOS3(cU=VV~R#~zYq*py;JVzw>G zChXMgl8?zG`qgb%>{Kp3b-p)o8)s2mw6JXI#L10Vqi@UMPg6YcGZRmpF?-C!FJ1C# z(b$}ZEpTwXdUFs5lN6SJpG(+zB%5$%9}zzf)0WQW^4)IZ&CNCOqIl1zOV7PJ+k~B0 znDYEWfG2L$&G5KJy)_=wb^TXnKsx5OJg>>}JO@)re7lj$_oa8@-$@$CuY? zUL)*Oe5$$qjIh(=xPzi4VQ2Nd^P@^MMDfcGo6nqG+K6pAcd>iyU{8FdoBAlbIwtNv zYh;+(w=}Hg`IYr3;(6`S7AaSE*1V8CD5%O)z)7g!gRFS3E3=r;IPXG+Vh5 zGnQ;Rn|j?7e>|CWLg+9PkC(h`5v!YqJyHLZTDX>jITbVuar@=FEl5n5h(oE7`I1~a zDeTd6pQsDtC-O&3d-Aven<=#(-^!Ybd*2Qy%W*8lN}p`h z$|39|F}rp@NQ}E|M)g>Pu(Ruw-ReiS#JE?_7CKwkfDPO;Kk~HmWc*gqohX-uOgv-s zB-V#f7qGB|gMD)clwwc&ryS(w!y>Je0Uks=&x&>$|vju%)QOUdFwrY>tsZRbz;rlb=>#~ zt?xI^w^bP5Jh*;EIWylxI*n7{S#Iw-nXJ_ z`Io{X{Pgwsgx7l>V3ugeLb_Ym{)y*=|!ogv~-iCMQi>4gy9Kjh6+tBOWU z<=G3haavRGrWc1RFQC@=+Nh6}$%ijs2A6L4(|lcuNsK}~x%1AJ8wT|eL>_w`>!fT% z_%kHkW1mN{Fm9byCH6k05t}jdrQE*R({T4(Tm4tPiFujPZ+-7s8CdVeh`OQu%doxv z`&isOrq&wdr$?+ej{2IJaK~NSFlWY~#lpD4Xh)H2(;BgA_oFv&rcA|O7#ABPJz(Nr zHJ`k=vG4*`xu8V^tt!QQ(eMx4`SkusyId0jXS;-Yb1ps=hMM`~tVHm#y`yf3FKfh_ z!V`*|hI!yG%N*W}dB?voGN z|0ZEaVVD$iAc2p}I3+gMPEDP=QSQt-Ec-~{VY@H&*y7!uiRtkye9g*Mo45B-e9gG< zamBmx@z762GCEfaF}b}@Jtc`ajJK-QHM}p)nir&X{Tva8ujYS=8-G$5kKg`PKA)IR z-4;uHOVF5x&sa6)i`Nq-p0_en-!U@-Tb9;W$XLA$Gahs!>nh>T!dAJ^A(k>M^|kM# zx#MM}5&Nx_A&OV7x;>%%TqCxj`Sc{YdJlZa%}I;R`dj1qYj=vydy$THbDt(7XT!m^ z4exn}I}ge|ofz&-tPcY}_gKZ{`_`#__MeUt!yjMKPQSFh3ESUb&VC&_9gm+9YBu

+|}q z5>43s-LK1bewdD94>o_4oo9`|UK#mty~KH}->b+C3!}<0RmBslx&2CN>bGOzP#Ko9 z{CY2L{74*&o9|jDjBk+IeLKgi0o$-?d3;b$58QoV>Cr)wOgwj_uIjo~X_)?Vzbo_2 z%CHTWKThV3yZs^cFXxE7v8w3fXRe)Xa||Bphl%2x*M`m!$&J|YtGm}5ukpm)p82R$ z_Or%SbZz>6TXh~g%|7X}wWSQp3w`$a4q?ac+-!xpM1F2Oe5@HI>?nWf<(rR);&{DJ z!mIQ~Oy`v8tc`}Ac=@eynQjBE@sal`o(+&q$42#PP;m3%U`LLtr(GfJtT>q5Z!j@m z2#-%|xeQu4GiMi{THpLl+>RwGs*oICf5<}`fHp*zB3jxlk=eqR@C)k?<{ddrup zG?Zf9=M@TZc`^4|%dYXnxG&gYYsJkQ?wdkpehRL~j`Wdiz2@0~SveXmD@?G(C0F)y zx*2YT7vC8^I5;N{x4sh-anz**yhF{`q+C;%8%G z{k(gx^E2-CMNtFxDx+E$pKDlPc6x3j7Vnv@;h8=ScRr|8aqcA(zu$ah@Iqof?Xt<= zU~)+*mg(4g(=Eb|e&UL}Xkxq~5^tn(=PxVK{m*v~=z^Dxy;-o|h`3I9dZyu+?RGvcEc@xF* zf2b|Pa#ECh%I(+ri*Cu%9Ae!fGe+~$aAF<4jTLX6F%7R1iwcfQW8$n(DIuYZ^B8ma zq_U*XrC7!r-}@Paox`{qgS&3Y3u3Xjy!hNdEZ%pt82+$MvtstOMy%E(Bwo>DI(~P} zwE}MgYrM(+Yu`%W^O*JHb^CYJaIopS`sH%h&tca&o^y4Gb+~nP3^)HIdly>g3=zez zo5w_#2{&TwVw=FOEKhu!CAQgXD-)kHu64qoyfkdfb~D*qOE?${J9(1p&y|C|wFihe z9JEzFiQBKLF-FzRTZQnKehrU8Bpb2C{$;z28Yko3u@Hm78=1Jqkm-F?Q_?W+6LM~u z$*-{Nxq7l(UMx{|n;%T%=QW3CXjl<;5cXPy%{%L`)5jmy+*sOx84kUC`qEw&-p4*> z(S{f+{8s#%iid^;c-^4_&MvG7yB&YPU5bdqWy*8T6cvQorvc)Os!*o$Rt2S(_^nGW zt_ugd(qDBQcOJZb{!zXyVMpiD=CNFU70r3%88bi(f6)C-hW*ke%w^&5F%69A_>C$h zQP~{Aj?%qU+j|7xvtsfN3@pdI*?Tv0<8ZO6V&8+hGOSw?V}rPSKXOJbW0Jopo}_o> z;i&LN?BXc*`e%zhaXtOz@hgdW=kx9@)d2_6u#?KyGwyhDun}U~r@7-U`(^iaDZ-!7 z+i^|YJa*z)zv1=2EG)^PX6YBFAIo2 zxG5c7?0v$L(BWWO(yJ|?%!AzL(CV4$308u#!r8d@U?CWL~&%Q zt;f`ljhI)goYTW;p19$5>p4t!YkdEGw(fADGz_bn@2InmgB?F@(~Wz5kN|BtAzf{AB*@k_Y}m!972g7_aUjLt9Sxr(w^7Q+3|2 zD#K*v8qLWk?ATA7-dmB#KdS3`8gTQo?ocrs#cB1}XVGBiF=-9h{n+!7j>@*UlYjlS z!=l7GDkVEq*f}2`sVOHH&MCs=LI%2U*N54KHM)<8e0Zk|zKR=%wOVI}ch@7AfHwi_Q%5a~jUUVmI#|$sMl*Vm;wb@GI=8 zd>D7U6wX@q#f*t{OL$Ak%jiZdv}}*64iV2wGZqWkN)mRo`p?VhlZNHSlq8L+D8s5y z!>Qlrofw_-B7_~z(1QuBJw9j~sT4fo60qWHcN-NovSjaaZ% zx3v4|o_NI42Ztrytnn^QyWf-+pT}4eM+cl)#ldE%-_PdqOEbuI*LZyyR`=)Ma4x@Y z%$T;K;T?l_ziHq1E?2yE!28Lsosnh*je_?}E}n9F-ZMjZKWWZ%=JoyE^iU^9##kBR z`7o}JE?9-`Q%A>zu`~8Wu+YtZ*S$D8?~#c^T>UcAGm-blC*@4O1n=KG-uzi+krBLq zv*7aHDnn~{e{0D3lJ^%T!25r-`xjar=>wmKkRLN%U10>g-!rb_@y^@JnW(w{Q_-ul zImp9ejll9IHFVsk=6B0R;poKiGkZFRzd?p&p0<58B@@}3^K8u+A^1Fm!gE6ph;4<> zhnc=BYHwu%e4d8F8?~e}?y&!#J}z1jl?v~-&7UDQ+w~9kbGvs@817lCk7hVojO|vz zK^)W~(5*@8=w)?e;hOm@RJ*~{(C+CwWXN7j&S_O9BA@2#=qe8Hx3zm#yLQhB_&g2a z+Gn_Q?~lCufiortPg^Ji@o`gJZgI{X;^Y6RaP4O^h>v~?J8bDzBlK?Xxbrr|^MTil zvt7JkusWJ}IZw*@3x^T!qvXA(xJtW+B6N;Kq9m;PY6z^?N(x|ZI+)h^ygJ`mCVzb(4V4{FIS&Z{R7_* zm-TTT1n<}0*1Noq=UMoC@?}e8#UB6RdAwb^t=v=?YKn%{DTROZ;~;L8T3$oV2BNDj z3wQIojHA}CJ{+F3s0I-!xZzj%gpE8OUclHt2l})4jrql8!=XQ~GqEX?`ayi=+{B|R^HZu-OX|3y!F>Ftq%`S@n~y9b$G zh^R&;`3%}AS(=GhCaWZl@rC}(f0fhk+Z~8cAJ&HZTQ))d7U;~;mhabx$43Rvc{5#u z;e3&mdSy+48^otcG2VKmoIaZ5n<200P>zVb?>2Yk2z8XXF|OLKgoUbawwE2J@*cUe zLFW16j7&u5?Sh;Qjc~q*DAYCXY7gg&8;|Z<^oab!^Q-5ds*lly@l@~d@@*F;jGv;b z%Er0+a6H$3%yQYDG#uUM`8BeOQ#o>8{j99pR`%~^O zA)X%?p_jNY7V`aKe&D^E=ON!u9sQhE*atp8V}G~Y16F@{exC)0^=i@fFdklGj1Kfj zx((;wtg3Al%k~T^AJc~3ZAElU}qUt9+8^e=_42v2&OdXf<@_%I<@D5YPXQy4KsGStoA^IG@A3WY6RxY0A-@hybCar7 z8-h+t-0Jw+iCAxl`M3(tP)8&Ew$?FYS!l?9{X9d_*T{l}9QCO)GLamExVaNjp+CE` z-#nd{2KmLvD1Kai1;%gL!;?y;MvyZzk4%I2nsvMKb=ylq+cZLid_b|HL2*_h`(|9Wo?9$$mi{Jrly%e7uV?C%v$}3e}0j{_OM6 z@zq_0tpGeUn)y`hc{|X@&U_)}<#Q14gW^1IYY`AH%DQ2A=m z($vm)#|bO@C^^M^d}*hoRWfRYyr!s6v^Lj$jv~H#>@og4Y7GXs!j_u?0^p! zZ8hw7oPj;6Fxuy))*c|gj^`o9`L1Qp$tifzvl4 zzHbwtMWNYBjsrq4?I3dVy8!}!rJZE$z~1vxA-!kU zRd4$VkpcCpVAd?^VF>g~^S)*{P2soCM_-7G0xX*VZ$e6s=zI?CW|k7O$C9i34T>S(1XO~3V z$$t*aD?D)_~;g?E5HA3s%8dNcw1WHH-sO7G487cl#in?iu-p#T1sgbTyj4W8Gou4%!IUiZ5z zM<$^2s!~Yd0W97VLg~{3fqjVO$ZO6vz&{d~V+nqyb_VEk zFMs;YYfpLDp-3U#po##^xRu%zyJ71=-ScSgRTM_k6~*;vW})jYV#26A9uj}dC9!)R z=p*~>CzK37;KLmK5VP$eV<^=eMS)?HB%Canc zK0JFf{4xOwD&MB5t`~wUdA}mOrCRXdUwy82l}X4(>e9hMUpz#4(MFBW8}K2A#-gW8 zG~gedC8J-X8NgpWmtVP9yMpy|&tWdfo?w8VTlKoXt9pTbmYW+&59!Inaw4h;vxmo^ zH|KgUaFhwbK^JN{UY4TpHL-r-;msN7+272_3s3OS%3@pR_5om@$e~GJDh_}j)uU;p{ATPu-X_@xb~;-^)Dfq zi@iYZ3pQ`oL1*J<`P?jY$4g*&c?b__d)Xg-s0i#cD*MIbk3P^R;dJ)>MMeoChS%g8*A$cG6KFHI+Y(gB;leUIYH zw+$Ki7tibJ&h!Q7;IoyECG7nXtJllxddi?a%d?M^PdHsb zds_KQtjOdwx0hjB*as4Bq`x2Zh|-?98Yyxn@{j{tnEc|^zO5CV6E znutVB|Ak0OZIk8#cxXE7l4R{qP(KU#t+$?N2KI4wS8wJp1pM<$#UWUC2dpm?Y-6ru zxqy1LB6{ZRUo+sZyyijPuL>$KbF>F#ukjePSvDr z(2{@Vb9yJ%JvWRkP`)kKU*WY{ z$jzjnWxo^Nl?v$;=A)NFn)H33e~s)>ahgAsgR@~ zU>~pD12LE006Y&gy51E%3Fe^-v#@7HO@R4tU#8-^-##aZn59TD^Q-gmn?}tloWp^-LWkF)VTE7&)*{;_cm(jtmg*=lD;frXN6k zQ9rzbu%`ihPbRRrjdiB)<0q$Hpn&oju;)&ay0t8=AdK2i2G<4>pcI3u9J)SkxK+;Y ziRiQ<+?Ks^;%D44bpF7rOa*H^G?(N4o<(GOU*5ux!W$1&BlhtV`N96ZL>G|%DfVws z-VNZ74YMnE&YlJF_39CtJNXpE*QHpwT4sF+PRvx(Snwf0J%%iE1K9qL)NS(fen(Mw zXrnm#_=OFqDu20YbG{av&$(5s_8G*xgzo}7X%obIOv}+4W-q`uRmVO(%bEcB`>bU2 z5Kj%LH!9eNY;oFv57kvNpJ?Od;k4c#UW8`^C_Fc1eFl5q=kqB~QMQX%eaMnG<#=Kd zYVD5P;BCZ18xOwtJk|z$2z|DFbN>V2!*`M>-zQFxKh+Y+=j_(NdhkwJ$+^OKfS;zG zu`k`LK%XZ+L{_ruPQqvTCuLbrj6u2$(~Tzgh2h)p%&ok5wBS{Op*F+fBGfx$>rwft z4mxsBU#Rvh$nSrevr|=#yVyKlpId__QGQzRfl7bdgMXHx>qk?x#pG(C=yxGX zOKCu#B>kovT6I953G#}MbHX5BydR=BoAL&D+GR;`j~M~;uEotS+?uok_&V-N195*- zf}#7@&OXNW6YT!7;0h{4VKmH9x{gg3eq`43>sP`8q&v^{pPwKeTIdZtergEtAuj_T z+IRrWZ-2r`j6B~8)_<)e!S4UA0scL`W>J@F`E@@YM(GD|R#V{LvkR=1_|B`q51BYD z`W*<6Y3xde%?A7v zB2yF6JPFo!H*Lr{kvWK%OCXL5#HU;N}tm;*be8@D)jujVpx2g+T>H>~l+9JVQ7a z_)GS`z@CUJAfMKmG708df_!ndr0n9o+aUfPGk=`DHl+mH#a!SKg~y;nSNf*Iu=jN4 z53^XvsB6RM%M_ckTPqO7+|@^`UuvP;xcUKGE5L_jxtU$+)*xS;l{V=ddfzSl&^+RL~j3 z^Uv4o!lu{$$zR!gS515Yes)dDnOm{_jcbXO5eiayoxHFNUEPr*0 z_&<2)$qg~GnV-O4W4gIn<^@0>KC@Xvdutqy;pD&zrvwyQ%r-2R!TpF?AFUpf-_>%#G6=@D!m@&DCNt;A49 zT`4KJ{hr*~CK~~=rX86+<%QYjyx`><20CzrQ^z|Z#~L(2&Pc-#VduG%$_`I5fWORq zQ(Zg%nYZ8eax@^zmvUeJ-=E`DYBvEt3sZfTemM^G3Gm~6enIq~cn=SFr|v3RcYe1jtJOvi+`X|WDGTy^Yq$KOTBVBNPXDYq61q~v7#!vgf_IC7~drWg3@ zaeDBjIz{lF7n{b$r;EdYpQ*-Yl8*0!eBtN(`QH5=z)z736SIPuFudeb2?zTVAWsD^ z<2Y=6nqAE8FN zgj)w@keId)jb4RVWwJ#4-0_fw9i-(m0qnznO6Jc0_!0$eH9XtQ( z7X8*$x(FRr>#AgP#zSm}o)??{1ncDvjlT>3C4qc;)R^Kc_j6DmCQf-XS`32ymjBq! zit$~4vH$xKKC~SHVL*RE*wOi8`q+E!AIoaXouJf7klJcF{ z=2^(V`9YDOdM(t^vp)a8WPD%V@g3_rK^5>m%k_}%!z^oHy)#gpnpK$x@S0xm;pJl# z`0J8=e)GXJ;II6AbD{8HX*hS5J4>4ud*2Ypu4+;!3NOBFjw^M*_$TBUSAGlj9#c%E zpH?OwvaY(7!(|WhY5h%eoc*n*`}_f-Iqz1sI8&^A4(#_<<#`my&V&8=D(1iC&eNdYs2q%#Su+9tJY**0Sb0eSPO0Rh zpbH*_vh+W%QLhTa+qWOf3r_36<1)g&5m@=~b7Kv3bVL!aAqnaYs=RjB7b&3L zaK73r#wic*n@k@3fVU0kL%(z6@031>uL!$a zO+H@)R?5q}wv5f&e|)5Kd$V>Cq7=9@HpYsFu9Gv?d`|`OrD^#~|4Rn&*T=IGCm#Rv zUh{AUO>6l#upXRj-}?0c2Kr1%YicU40Q@}gscp-?rUI)JSmxcu&M!`vF z2kQ+>51d?|62K37TTb#_0Km_Y=+m!p8Ykha?`YM(OOHWBKV4UIULiQr5MHp4!S;Vh z*(9?d5+vYzQK5er4>cYAR6=|V);mm(kK_K9gZ0iGqp$A_dBA&cbsTqTs#d{%=C37N z*6e@I2Onn+&?}n)e--#ma&B#_!y)I^Zf>Pv>o30a$c|cJSU=w6r)9SmtbTxuT(x}x z>gnfmSjFC3{Bzjk_SIj&UmJvc&7l!spK0!667A7T`|*E#`TCHd6NtC<)5oU#t<(4U zAY*z4cK+vmVNZ<{st)(G;h|q$iSINB(1uVLr70IL+&E-p+AFIF2a{gNCNwQU8a-S| z16OOIwH10?tt+SxAK$>^?~p)!D8Tocw`LyXPf5JhMn@^Y&(h@^V^t|&zpkp>H9EH) z@Kdm<%B179Q*h1tN1WjY0_2Qio@KKyNNh|ACh^^X#w=V=^|-}e`ZO~Xf^5$L0klBqj^9D{9+D{A6` z3D9`Zr z6C&`A_X)bc_jTYqqcn#1w2$LjS3D)f&8TxXLmY|4EReOib{^FT2lYF z3p=PmUURRW2x>qzSIDxhWG_&npEOvRd)^U|PfV?E2Xfue(BxdfYu}B~D9vNv3WUp$ zOjD({rXnJ;JVmD5`vU6=^^B>ya^>H=fn@^ z;8PluC}btxHIPd}0t-g9-G8~DvTr~9b`v*7x#{{c+?7vG&~FAl=w<;TKuncT7s~9vd7Zd zx9at=zR5sEj=V$?B1(3Nfg#@sU3dJ}DobvRo>27|?&_~XUNBkYd&?{#i4*o)d*3xF zM4-F#7pG;HEB~Dw`z{)`KNDQ8Sl{-Hmc7$53F)T|oES&k(AU=Wl zv%DpjjM2@m7p0TsWr&^2w8KUv5mC;&x%C~BbN6`5EPWqN>#sP=^q!rovp;?~ltqdD z3M^tV_(MblyATJ&&K0d|lKkq_WQ6L!`74^|U53cLadZixAtIMrKG3vca@0%L1DS_$ zTESIGxdEvJaQHjv1r2B_GZ(R#qa(3WE`E zFDxN^RPXJ+VREiZ_HMDy;IuN6)X4VaFi{z-SMn|(CW73)2|6Ss4W+1MQ8Ypo$oLBb zx{c5QNYm(%Vi(#j@g@0p76Cc@h?7GL!(l;Fkl``xj&aPD@P-f!hX%_1L06N=Q5}cw z=VI7B>)x`$1hak@l!56}S81U!x=WHeYFSo=gsS(aewQF3Qf)=HjhH@u;;qCk?7LnC z2>&lGe!1~YndniV@mpmVV+%=GKgXNqf<n3rr(&{=Tsx@(eO{07n;Wd;`ruxqf_1koev#d2n zQ-k{iBg zB9fV^0kuEJc+6_snC^lxiYLN>A+8mOyNT+3l@CPZ;mMDEd;Ie*IR0A@hC{vfJc&Iy z@cpE)6mcrFmwIqCigy_i*xuyQ&2&Xg?m3j4AR3|9%d>GopUMzZDs@FkULsN=GH15u z_X5?kkGuzQTB1ki+xGA*KTu##R!)hE7MEB35LiYQuD){Z`{9am7M@WLO)*B%v$D@B z*~^gujr)X-3q)kMeN1S}9d-ZUw6;40m-qbs7hR;H$NKdXM?RPJxRQ|4_lnM= zyl!Z}-*2uDUyM*j7REXvhG)0QALhI_iOB0KJUw-o9A-^NIRi}2aB6_(o}E?X4ag50 zlB35T4{a&l#ril`uS~Okazan4$#CCfHAaWSTV|6*D-aV$)^$&|CFHeppWqKnPIE>r z^+Ak(uxFxr_r5EQGJ-TFGmAWU6BwMWv4jk+5gsJt&!RRq6Q&ogpF>>=w6_HnyU|lI zRBCMbqeyHxSCJBiL$!J-hKh2W7Ec5D;9U%dW}WreGhvjdlkQl4=mQe6Wj>rFMCFE3 z^Av0lZH>|2afkZlu|Cc}hX|Q6=ZHwhi=^B=`?zc@@fZJx(>jplIk~5gfUtihc_IaB zUn@W@e`Xo+>&+K+vUWvZj3$+LMH{0s6=ih7ofU|cx(*eNpNK3f1bp4|doEoFzcc19 zqG+)Bo*ab?rJU6~EZ%<>bXnb7MzWMn9xPgNLw^>~Sql~!qbk$91{deckt|+hwVn*? z1Gz3U*^kM2J4vgdf#I;Z>InXS_^GVuV4NUF@8u1oY|{Bv_D=GLB^gjMZ>CI|%ze`2Q1K7{2R z>I*#HFI-Ug?Bv+jA4X`s%YfvM#R?>(w7HKuZVB=K6rS`8lT$Qu2F)Yjv>pYLccYT>DR?mC+boFJLO5zFdw($H+EM?JOY@WqoO# zn4I4ySa94}9%Q3A@^X(4+qPpuR+A}E)3hGsQ7Q?M%B~4lU~omVPC{&o)qY1f0?h!g+!z(BlGaH za#xfgSIz8sml3L?lcv(lU4eWT6uRCD6Or~1FY`T~hft-mEX-eW+?V_I^f|6ELJvvL zBb=7c<|W@^egEGRzA~LXhw7y>JS`b8M6IEv$fV^iG+t~?m4$MP* zdYI&7j1te>KI>XkjpV8SidElULS$L~Qtjb6CTi|uBW9lfYPFwx-xX%wy&*kGj`DPc zgoHROBbPhx{Zu*Tf;Qp3%95WMV{CRCi zl6v`0{hl25L%Q4Y7|*kMAN#eZ&&@gZ2lsfWQ03=!zvW;OQuvX@EY8jy{lKiLAuef* z{yBEk-%h>^i6pP~j&~*@MjR^3d%SV#LGAFhcAVCSht?8%a@>Z@Wn@oKp*a>b4ozP% z{^7I>+@*9wll9Ab8bplI=wzmG=7BNAZm=h9}0+h9bbw#KaZO{g3W(Px|K?| zC#L}pU_SF5r}fMyFnce)Msm7UxyC6_lHTEG&vJ>#Ic~YS*pqIkh?&OKtMf)^O-r(I zsZ#|qgTJ-j<4#249lTWc;!6S-WyOQ>kGyZC#2%ihPDaO;6PA(woU4E2vAQMU?_fa2 z*cG&R%wQ=v+z7@0{hWrofJcvOgvU84x=d#4sZgNZp1+-MS7CE3GG$O#em9iPxH3A)-xw9PlG~C0QGsM! zxS2frl!&-Q(01(AKh9idqsp-P6Src&hV0c*l=Dm}C$WCw%_<=m4GceT6LjgU^xRRV z&gv|OabwhRLon@NSOsExeGq41kM*+^y#L&X$vLIY;M|S%htZ~%c0pa^*eHyJ#4yXgXO`CC+$jkjZwZNj{de&Wk|T%NgZaaU-Xj| zMdqFy6H4Sv|0%!Q|BFf3TtB!EjjK4L?3eLKJrf%w%txqBr1oyunUJs8IYG z^l+jY_59t+d+6l^64wy)UI)YB;=_b15ihX%r?bN;48tK`SqtZ-dMfn$g5Kx?HZNrK zRs4UZzV7II8ZyW8#>VL7gd@X}^`*#*5ssE?PFR0t?<(1zK6AvIM@=z1FUi~V?fG4( zu+i9qpBlAV6u0egT}BAY5A{(tchpBDgFu~Xj6TY>Uu*qRj^v0;opi;K5Te|zrT;KF z1fI0&Ls;Lx9JKAW_g#KD3l;f)DA9TEw;b;ayuFBDyT}doND?C^o%#>Ty~+)a(cL>};}i`$N`mDm0;$$gKR1giuhwyLr*r z1D$d>XK+dEJUX#^^)99;+JsxZ3BE51{ zrA8MYM_bO`Tt)&PJqc7Ua%Li6IDV*b{L<5;xUYla;T=C_?&RSDN#&~tGY2ttOiA@zM_*z1Rmj5chJZMHjR8o#{`oDKaWy)W1 zNl+T2US*Zc_Vi`QpwDPT%S#f{`1h|^J0>UeomrhwD^APGt-N$E-YYqHBh6$f(bsl4 zC#(;XkRp>yb1z9Q=%~`@sPljkYLFm!Lk=!OYPET2zGf~VNhdy2?D>nYF>dMymcQod zBGu?IIqITKf6~2{keZJ#9Y6RlA$7u|a+Sg+=rhNSmY$b}sA|V3<@jhldcH!|=yb^s{a#o3+U+@sHU0LQ{I@t#$>{OfG%jkBx*pC21D# z@sH4?&y9Rcj*9R-%RLB^KNyyTMfbl3M=TO1p^1*L7t?$E5IMJDo2Hqbf$%F=q)1^IUA3@5R?@U?AJ82bAc{xKeeyCkc^kg<9@5V*3eP zWbJGg#^^m^nbfAnGUPx|*UwKchzJwosMwx;x<&hNemz({ufZw4$Ipk@IK+D-mXQn6 zE{k6f67u@5x;TNw6jiVei#E%6-PuzZ?$Wa7h~-`U1FpMA@wKv$C_xgX^ak-Dy^M%kBcsE1XKV#tOO z%I1>wotvuy!7)4zm<=N$x47hO_VVC+{i?eMG5a8kXIJ<9B_pX~b>Sf;dKTL!c;mZ_ zTrn#+_VkE5>ML2sJN^ROPYAAzA791hfrQ1!)ZftzF)q3 z`mFJ9bk1QoER0Rf*{h@MjP-({(_laNDf9Zh=s|EEH!?rMBEM9`7tQ?nU>GwJpuR)k=e$uM-%3I@4fq-WeyS)p38ieioK70N<}jN zpZ)1i4oOjUR^UEY=C^$M<_NG~FHEatOj!!f3)GGey}RdyU)~FU_UtEhLn{1 zqAFZn$5rs;!6;4$M^?0wp}i;tYjOMyNV1*}$T zv*7%hYwd{_#Yu3VqaDeXWH?2jq{D27o?a??}ey76bb@(3W4$ zs94-TZ#kwZQm%a(=p)Eua#!9dYTsTd-h{pVzi9v%;5KUYHuQTzml!Tx=D;g_brd9H%{+CoJp_ma5* zo@MOt52!`K{c-l`5*I@^U>~~A_YFM~!1-9d7T>}lSq0e8u&cZC5CPhKp4DY`8tebI zdHs{g6@|x*{xgx`nTAqT?WVI$Ya!8agZ`mra6Xp7Q%Y#>1?OV|!wkv(M_2ak^G$eM zMeQE&$DMn@dN%3cJi(v5lqnzv#NRWfx}uC%Ct;=?ccVqD|62*AJNI!|@&< z-eX>>rg`ZCKB%%voXf@fAqZKtO!L_J>Tz{V-~#_N9U$0;546*(K=kxDpZvi#V;A zd{hz5<@ufh3Y z3x}5yyA<%3p8Hz0B0U?j;f)rCn|^m@cT8x(^S+BFN`8}&QIuU^JFXUj7aF84lL0(a)((unGX&>j zVVtc_UjN)Di@YeXc}EfW>y}HMojn^kUp+B?FQMQMz?Uy28ynF;9gd=N^Hjy|n`{mK zd}&}M2*b=J+=ago*mmYWqCq9rzl#fLy3dY>6gCfTw><~;QD?fV^5QbE&x`seKbEWj zp67-Ol6frwo~RXmFeRaY4?8W3hH1tDAC|TH#YpI4{c{A@Sz*IbNa5b>U&1F2_yRfp zQoz^=SmdYML4m7d(19?%2`$YUXzdpL+OeR%eLURM=y(v<9=y*t=RDf~^PUFxxgyUU zyqom`@QZeShzIjA;4gxbT2f~q;KMgA(m$nM>cBOc1hNCr7*uNN<01B55Wcr6nR!VS zg*~4iUHk$~Lju2~g3htvq4Bd@QSa~o&&2;+=9n6ReN-~&3m(1ze5mp+AcpFn`x0>+ zO1Fo+06!-r3$`;P0)3!$%_G$@~}SBdYQ4{@shq`}R`vVB`r;1Lx1t-}E1~CINkD zquPH7KL+vMkpJL}?RyRQ(-X}ZGrBRTvCPTvL836cR)0z3UZ57t`Pw@tBX1fK3sq{~ ze2m=}5KF=bwg7z^s%{>7o(JMx&^>T6w+-9}swrOMm%If0HAy-2GCmmGXV;8=A$&;` z=o99&qaj)Z!Od`+=@`^=lDhV&x-e|BF>~i4mQOt%ef6wB=Aq@4hc0m+@erC=DY^Oq*eCz@ zl8rD8&*m)J<*&8C{o~hCA^~d? zAb%Xhzw}8t60vVTE=pF_Uj;z_@W8%j`g?2_S~<78ka zmlO-3qEU!0RI=@hek~MAmqOru(Yr5ix3pwpB?;7rgpIp?7EgivuQrS9<-;LS+@)7h*S}_Xs}+^w|vZ zUKn~0;$8pk)-`76+0^V1nRaGxvhx30~OI;b}m4x0XEy#w-T z)0Mu%B4aA>Y}bX=63H=W-KoO<wV#FlM6q3we2j-wq;9dz zhk$%~R?py1#57nBieH=h$S6d%@6RW{LwlGn@$CQYua?%3TLJPHdEVx7?pF|B9ycW> z4HPwCd9}>#k1Yg9_Pbw9q@@V#`o;?L=HBj4^^l9hEoOaeCX(Qd>#)I3(6Fan#uZ4}JT#79jG^eQitkQy&Q50Dj!k*v`3F zgZtWR40vm9C(8Zbe;cE5$bbQ?cLI|O2(379f6Bb$vz_Q~B{*t5;Y|L=Q7Gm~X+s;< z&nI14>K%VU2j&&e)A)O43W{j(8sTli;(bPkZHofjpSs~^evXG3)Q7_C(K_x=!Fo`l z+&bZ{F5qV^?#HwYQy`yGr-LC6fra@Jbf$2$tDR*%^g=4wJ@mD1>0Mm?WiO z?0DeIa4&6@#0->9r6b0%S_>^p1V>+z1MyDC4RQ0<1^&80*+3H30Q}SUpy8FJKgh4= z&vNgYDMjwv>+9G`ReKMp5Bp40)do9cVVa)5>u!WmC_7*oO~?98x?c%YD{RZcQ5t16 z6Rwlcm&Q<)Ywv2HRGj1=2~JQSjbjDe9sZP;@+we(}_9Ars*qKrWX582c+ z-Md-_?)$6or_^ND#5n?_kxvuZdjh+kcz6W<6|M!Rl>W3YyfzQDj#yO*wA4bKe;Rh)Gk|#i zrvIiBSqJtJZi>)tQ3bz0^i^|>_XV(z_2^N1KW5;s!*9JSm`Z>?AqM!xq9tvZ$DqH7 zA!Z!AmtDyHNm&@a61@0*t5OTjVx&3Kt~>?F;}qsnvHn6?od{=!UEr^sl;<_+|IE9R zoM5}?`W4`Lt)gE?y@6#PzvriqdfzSs`UJ5AQu>_$`jDy2>@;1IhIa=h1qKzzpsQ;P z6ztSOuuT+vbYTaDGqaYiDmyGd;Y-XX7zlXi*|F!pQA%K+hZ^E@wd;VN4Zr5}O38!u z)}y{+-Owt4Z}(5%ZpIKn{c|9YpGnIC#QVKZFq@r^3>-u|rf)4r!1{@jA5yRg!6xx| zyJl+yeo4A6gm^7N1#)6N;xV<*w?3<(cz=+;aKDD%{N@1hKK%BF39$yOztVIW9jbc) zzK3kCu?uwpKDf*0ZLV8%jjJbft$u3Fgge5!2A`d z0_eaJWRY-Y+!xzF^hs{)tTG08&V8O@%5I#zzrN{uv-$QUH}IFt6ysZVZot=v-v+rD z4g-Ct8uV!|+qit(V|S#6m0Nl6slGZN%Ne>3y}M=d0l@Y>c? z73^>LUol?LO9XhP73|hLt^n%8c6owqQn%2){&RWc5*y9HUsdY94a`3PAC6R-TAt;T zfxo|rjJC)6Q=S(_9@`2byvLx;;X&o-S)SC_E)Kt7vpvWfqi_X42V+FYVeuRPfL_# z1c<@z#I!(y5NsLb_95_=7VPFQM4r2_0Bvl14{Coz zU4Kv>tS3&z&uZJg1p2Tqk=Z^O0sgwlvi7p;H_+#2MGNf?wE`SIWb(m6e+;_M7Z+Li zR0J-n@{+6#!S+`V@oeCQ7NGg6Uk7vE;Gx7@pUA`gfj)~*KmWS>Pkk8LrzPB&1L}=A zZn^_eCqO(;jm5CuUIzLU{X8NeeG2%?7Cq;^*dqlu$*$@7cwzV1i_WT4VE5S%Nhs-V zmuti6i{nOvGYgPX7Q@N3m0HNg>v@e?7tm)o{V;a78u*KT>Dcc31P zJ3SVmF9GTw)b(3Fxgq%XC&STMeXnHT4}}&ohyAhs9)Iejwx2@q@rR>SvDiF0>i{!{ z_gMdK;NvZR}exOPMe9b`ak#^k{ z;Ms;Iilxj6*k|q|1K>sZUBc>i+4{BjX@{AUJbeCEDZavjvn!-(T1s~shi(RFJZ+scgp5rJT&|@ zLEP~Q@K^g%j_!pmfajaDZ-?F5KtD~)*DI#qdjVcgDzO%K9(lVT|3Ru;?~mUDeAp>< z#AHlZ0e-dP3FXR+LY^Okv|}WB;SUL|vgaQv!QMlEo+s%qLORcG`erKDK&IUsnmVC? z56RJE|E(PXd`M>*s=77|=EZIg#S#hnAYW#)>Fn^=fcl5bDWrYw6v*%HHrl)A)zn~F zvUM}FWCGMHZq#j;65C z1>i^OjMH|^9k4(4yXruBcPQ}Jcs%8!*?;CMh{umg>i(lo7t3*!tnCDxc1q~x`EYFB zXB1mIzkndT&iQA+;0p?KEOxchon3?i$$GmC&GC@!U7Qi(1nLcWlNI8W9ni;o=t0eB zE@?l1{CetRR(uN7^WKW*Z@u9I=cjk)r+#Pmf%tL|c}zR$j=<|J&92ME1js4-j5ICQ z-{Z~r{cnpC3fFt}J3r~1g4$f1bRM0>_J2M%QH?r-_;RSZ@mT=74ETTjrFBNLTZg`b zc*mEXOIa!b{;E#%h)_=bC%z60JrSS={FD%U`Fs+G65PnaZb-&H2KncnOkO%F3KO@7 zpRVTUzzK|zmA7Mw(0cXx)FF#n$hTp$jBX9YmvzXw_e!-OzR;<9D(-*I}#o;d;Q+j1FRIX4DXOg=cUgZ1yCCH__qWp!XX zlGDek&;`iRyDc&jo3FWL?ya!A0rVjstfD=67wB_HnZl{L5!mO6t@25oA)rq^`+a-T z5{UPmGE`{12gKh%$EC=PBt=-tK6gpu00D9jN>wS{6o$)i9S5XcwPDM6zZdO)mY||v zLWJ6HHPB3Tsi^sD5bu#}oND)|fj;&UG3(zuK>d!3^$mQ*4)ilrWiM&;0{RH(!!;$0 zz&@|Sx`&Vy3HVV)xQ-%9fL>)>V!pvB1eXlW8LnN^hM^@rvxtMskjBcVJD%CK(Aaim zm4E`Uk5Qq>TwgoTM^1cIBkiAgc0S>|CFB=H_WkpD@|jq_J_|6hyc$kmv>;jmFqmcU-yuTGIJaAfc6r!OUHDV72L|9=81US|eQ!g_z! zsW|_RLf1s2DieyiV5iG84(ilOF#q{dxGZ)FQZ>5X@RzX`+8lUFP}K)~sN0e!ylMpa zuv;+M(o+h=_kV9r55Eck^-%s)VuN)nSnsTM%|$%>r=OUG#BrISSqUyYo+=@$MS#*5 zbF`A&vH3UQZVq=mwczB>s}DRZ7oc{c>W0oMtY6mNqQyoT*yllSNDS*P;GcvFdaf%K zfDaYAJ}{VV0zO>kc-hWr0OG68=|1~qdVu#SFGFWitP)&*{%+}AY=1s+=&U&Mp9uW^ zq4P_RM>;T9U;F3Lx_O92Y+{+I5)b)#6s&iB1p6%`Ly140v;zD%8Xf(uYXS7pXkNDT z9PBmM(ScQU`ycR4 z5+V8Laz2JOwNU#2;&j3g#C!TqbQjHWV4s&KxY0ulp#HHqLQHx@FS>7^-H+iWCyqAk z|9!Kp`ir6vSbup&*&Y05rwXSB7f`8V^RAMjVmsx?MB#_q{?8A4>B8CbUs7AT7NGY% z50u(|;30#RUzeErK)-bby}85jOtAja?|CS^4T1gUI0l}NS^xAij>_k9hyOFrn|1l| ztnf9!hm`Y#m9;&lV?eQ_#7Xa{t_2~{g(&clVj)agj5bXt^RrFmlmD*8|*&G zXH&+V040EDkGNZMhp&KmU+nfoUGf2bD&lEyYKMWnf}@|s+`R_yTrg$Pp1=s~^Iuar z{@|gLaKz*T@?v!Y^l;Fu%4tFbwjADC+PtX?cPEMRI+!g%t6JLGePnn@mG6h$bQ!SE zeN7J=3U#nvE)=P~dRh{k54@K}d-b!yelH{O?#gf%;KNuoO|66f%=hJX$hsMMNFFwR z{MzDf6gDrPKeIxzR~Y`_Jx-j)>eZ1o-TH;j6-eYKj+_;n4|E}s_U=7fV4qt3& z@{EFA|8MdAcz9Ku4j6X_%T1DhI>bHg9%%e}D%x`GH zD6=nObYK~h6vHi8V)uVWMP$PM;|KiXak+Qt(Fowf8!yVDsSbht9>wb?uWoLGeEQJV zcxFu&*eA7cOPp5=*e4#Rxq0@eF3fl-a@K)t9I8KVdad@580_byNgw=C7bZEMjF*mB zfjV-ep2TD4x#Tvxr~P4&FE;08Xc9+(eQv2lZ*ZbK^~MKP{JHFXAG;4KP9X_rZ)1I34(aKsVOz+@@I*^O3ij^cbk8Zh z2%MIooqkCRcE2@^HMsA73k52vRJbbNNkTkh-$i{S^Ff1j2ad&Qo1t+s61!jCmLlca zMHfYJTgXAQpL`yZ^I>|vfF=s-D^;VF>B8j5HAt>!VDA_9<)g1Fs7Z)>68hhPGd^g% z;PtTSdQ&unEAX7aa5=IQa**M-!WJUhSC~M=r+1A@Q&#UTA_D)^4m%)diIo~ zb^60O#JD10>2|)62G7UT{I)opmbbYad>*?edq8EvGn|_Og>IX{O4$3mtBC>?M>lSx zq0GlR<_b(vEoR|usn~Mli;C`ar`;yPEk5Ya7+Ip z75Y?wW=t`jgxpgcm9{g!L$98QUuwQriX_pd*54;>BMA={%zp2E*AE$D z`)!<7s=Cus6DFr|Fq-#m6*+p&VKA$Hfr$8s989;J^F~p6-+Rno&Cs>E$8jEIm59QH z2XT(ZTZr@(3U~~Y6TJOQqUSzN>*#{)u~tkDGpd>sb)OPde%X#Md5HBxkak{)kNTot zj|(1Tv^B%t{jy(O#fB>HnAKc7mbH!KtXts6F*%}RHQuH^*n3?pDr9?dmL-1Sgp(*x zb$jtM-?&JKkcsHwlo1Sv@_Y#zuBK>n-z5%P?7iD#BGpDHaSJI(-c0_VzxqvXiHG8} zNcO%Xd*4+Q!kLjIM~+&m-4A!c`Vbi@qdhy*e9-pqjjZkSrf98!OmXr5tLnPLxqPF( zEwdsdNf{9#JL7rqJYSWq$SNx%-w44qP&5Nld$^k@JLBkKm25TjPY9!GaNgc*sZ%%g2hjmoZR@hiB+{d%NPT4 z{8q#A#BQ<>L=TE9k#We+^Nn>Ri4rGu$`EaW!CbpEvt-ebH{N3T-j8O{3}^YIIo|44 ziS5b0r^@89foaE`i0K4!qKfv?1Ova9du1Ia+sD=|zVn;ZJZ3k0hf9W+gpIDR318bZ z!=;;#OH8Do#(9jNB|ZXu*2*dQ=jm7c@ZPw6 zH@Fha@#Xu+I~T#+q?(8(;qSdItjN8>X9md8T*GQOe8Ak~GbSUloy(Z>pRB*4!U+MS zdagDS#uWXSRyN2F_t82_ErXil>TE$2CW7Ud?-rZc7rre_eqZ?j5y;8HB8<)kf;sX= zH8rxF|H_?C{@|m+IsI&fgTdTem~cz=wY$DJcOgZ*1F&DgXeu~m3?;&9yWPI8-&-&s_N1nVY1^Q?- zB&B7Clwv2^Cr!AFHnAo9z7#TkF1s=g>i{{H3 zdp#=bhg)bfk5%iL;Xbd>Ks=)qJ2Uq6tH$OgW^}DPfb6eUp`*9kZ?F(Nq^7o;fgIaW zF`Esb&mP&+Sp#6-TQ;ZY1Y?~y9@@$p(VS(5*I#ZIO7^M1qOFf-J4tO|i;?v&dw?AG zRZbDfM=S&$xhL^$K+bq5&FzBXKUlImb<V|)-<`D{C3jEfeZ8@yoD zhUfp*4oW4>V2y(D=2rj?NjWxu;(;7=qr*!F;ING6UQ*j-N<58{GbrK0B4#tU*luL& zi{HFiEUZ#)hCg%*-kuIE#<+fV2|VZB#JB>2tk!@&UmIWaR|ESTJy7=O2hfLYJacYp z7nma_9nw<;a|}lwuL<)00REb-xRsz{j`Q1k^JprRW7ah#Unl1`F}LV|uVjC5yNugr zfjHe~I2uH@&sz%0w-13n{x+={z14KKy$o06q6g0?UUBd!)Q|8ToAcC!FJ9L!PCB>O#ioYumTg9`6&5VmcR zTf}HP+d}s<`Qj!5qSWflW;mnq;;Rk%QjC*UczIvX26ooCr>dl@HtqT@ zaz34nwFu1zdo-t=dF(PgNZ9bB+>OTTzW9sigJz6iVaZztd$tne5{x4>@Jd7aCU!R3 zEMN%8N&54%uLb16lPAnF$$5~zTt27o6%9@&^OK1a>|J*LramSVDVdhi`&&%x;(Nm!y`^wHc>Z}V(m2f8{yiUn8+Jxhe06c=dSA~$^rjePvNyB^G2Lm zphLam99BXX#!0Wfh;8m>k)?Ar#bYn_uB99^#qB=549Wfd6<4&8x#IHj54K=m_tY8S z(9lER*)-svtz292vjB(ZB^=wWGr(ThyyXy%kDyL?Snjm{qz|qzr0;tZ2lT-r9 z*l+zZTS|9s11pGngj_!+x{{g;>B7Bxc^ z{qUttHXUXvb3Al&_(;>$63lp4ql=u_7Pi~w77aO{YDBRWTnlC)1outmk#Sf!X2Vg( zMvZScUu}}}1$&^UXjV2)`r_$(xwvVp&2fg(*Vadv$}p3I#AWotCdO=}jgtBB;TOlp z-#}ah6MSmOaVjX%Kf!R08vo7l?u@(_csHVnu3zM%FW&L3pyXnbIewLNuCee#8P;*x z+Q$9x7H0Zo)Pu|$d1)P+dqPLPaWrZrQtsXu<4 zg@qxj)(kh;tV!v3P=Xajo$flrxQ#i!e7Q=FuVhnik#An$T?BgZ4zhhzxRaprMDmC zukzMaLvmbbG}SF|m(t)}zcLC`%ShOlkRS2LY5=a%6dqa5VvYv{w(|4-DZ);fKl8cu zeH*I|IqXiZqZ;{nf^OKc5DMIeOvrH&bJl&aeU=ul*KZ#``*#T=DeK%&sP@NcYlvw zB|#tHaN-Na{cZ3b22TRR<5_zW#=9@ntGdz`KO3anGwfrIiw=id$W<3(6?oN^f$U9e z>{p-{x$gb%g!BNJ$C$K?FO&1vb%8GmB9>r2RW-Y}Je-87FA)3XR{e0V%J*T-1arLn zHARh+YcYm(J-O>#v4xEpe5NMHyUuH6$4QV+Llu|)kmLR3Pd%!V2`YTbf>pfDbrDN> z?B{gE#~)9B|7YeJ$X|CTGu4XEm0|gAc8ZTEx3O)>uQp^J8;Pqt?%@Hx3vRhbkn^B> zwz#?T6gAGee_KczydU6l$7k%KMgXowpK>`J{9WG~oa57ZN-*nZ0cmf;x3NF>9`=y! z{J=0=N5L7?8`r)RlJWdKpzbu_!xVRtTWqX9 z{>N@|Y>dksFF}V1EZ}=5{;H|L-uqkFD^o{qDlpeIIaqA&aE*n~z;O5PPavmHO{KAV zj2frf7v7%GyomM0azr}44#0;$t35NkYlinA_f(yNi?Gw3o;MS&1APkn!^nSEB%4%! z_+>C}JrPAD7%!2((a z>eL@^VZLX#X*z%$jr^6_eZXI%lWn)jaq-@mc{^wS3O1+^A8C4Z5t}t;{60yt!+XzC zD&_$kdPfKBDYwqXD+4+sju5kmepHm?v`MYglUlM`MB06)!t??jUG#m_4Rxn~DyaNlH)JljALR;dy< zJ?kBa(5${J*-z+zOEWcI*VV7HBoF3y8*_(CXZ-OWeU$VOqULz1s6<$`X)#tX zGV)c4VjB~96jm}1nh;6Jfn7~V}6C{rXx|@QyaMB-^A=^i^C?d0!k_Nx{ z;PP7S0EkmbqnkdQet2MDI8FjNd-M*6{t96s z9RF@|@jH+sBK^*x^VKr;nB=4XOMDTdlQqtM$Zw0wjXEvRtefC=rJE_Gx^MB9Jbr%% zzRqJo-!|I301oALuW4?(u@E8^!%E0}*nm{o(M?j}Iz$QP4GHjGnqD!%2G}Q7muEX> zuQ@K4n!8z2TZ|p>ZW;~WyNzWPt>}^Kjf2~b={i78!EkmyIllgSB*;G;qQPmJ>O1d) zyraq@B7dJ9#A_pc-su3Ty_e;a&vTQ@tYvc>x08oS+xLhhL3=>}S#oU1Fq`$G5SYkjw{frpL9 z$#}MB)28AM0QJW0a!N9O`VzxUEAP^PJ@hH+e*+Efvkb}Ic2kFd8f_slV=2q zE`2PJ(3JPZUrAxXE?+?GmlbO@QB6`3ae&d@N3xLC^-7dBrz|P12 z{+iGRIHVS%tGvxmg1Nc*r~lt0QIwV zs#tG)W)b$8n;^)iw~gIarhQ1p55<4A5xs!Ntmo)Ik$L0TrN{#npq@`>h*hslUc?qk zrnaM0{qZ{DAwmwbIsS7Gle#dN*XMLG|6JL%h5e}LF(=2V#-Hg(BH-uxypOzO96s^> zuSTW^%-LTZPFEl-VX?Pbu+scMJm$l4sco`3E;D-Wiu+^{_SNFuP2cptSWaZ(6>^<0 za$h&|#YGkZKBqQ9mg65$(;NKHKF+S0_of9B@I2E^;caYzJv^U&|M|$1bFJ|F)m8V- z=Qy*W{hFl$uRFW-mf z_y2t`XhPfI`RhY-t#33TR6FO-Hv6@nH~rJ^CgxvzeyReV2WKtRnYMJm^YUi*;^4Ll zY4mK9d9h3%n3uh)z)R2(K&jrkmk*@js1p5mN|%7&h={Jui9{~&-gk||?TCN&^Of6; zR@u`+eMBcl?@o2W`(`^YnIyM7sNX@fNH@hB-q$TFP8Mf`!hVAEg#xCVVkb~jX}1xt z`y6yH~2!t3heh^ zx45bHm;>%By>%`Bv&KK?>n!{O*6p9|_Svn5&`Jk8j@wy_sPv|($a;C}$ z-tSVj{3mfS3*HA+&0mg2pT+I?vzR8&InM>&Kcoa$?xEilK=(;iOIunEBcaTpnTeqP zbqTq>r|+#4I(u+~LM81ta^Ramy0Kd&V)Wt3{qK8+cjVEv=DeY-fc@J1iQoQa+d%(t zg$XqgzQg-Cu4^<>3Cj>aMTZxVj(_GwZfzYFtcym_=|8t}sho&Nvoyns8t7Mwz?A;Z z1Q1Xb^|nx5nMq`^GE;x>R29OxZ~h!<55%)jP3o!WQiz|BmMZ1W`|v)2S7%cwNDbmO zT%B^nr~}4V>iX+g$p=v1a@kaYVXLEPufQWIH(4UmG+d@u)GC1XxU6{p=_a7IM>H>H z2}~jhIjf9m)0N1ofBNsaP6=pVxoi zp+1IHQ9&=dV7w>95ZNXTa5To%el`uv9~`xIx%bVTA04$Hz1vM7paq>C9;>Dkh!wHh z9|8NrW;5O1+2z3dmRRvVDv42OpH1n%HZ6+KJ_TnrM%^5tKA+g+(v6LwKHF7Ze)T7z zeToK?gc9_`QGvO-51Q>nv1JdNydd63o&;+7%r`rqejn>@X%k$tdVl3J|zZpUtjpS^#7TJb$Tfp)f9j1aDX_w=Rj z;c;_E&?`57+d zKd5lTrOJKhA%4bWnfir~L4VD^SWo!$UK%~pH5Hbf0_Lqf8IAgs1<_k9TTQ>E2xyze z%>CKcIpp2*OToHKHHbV(TvLJ#`m1xVlEcXh&|gmVc^tv(@V-~aNWw*v4feNw&Y`Nz;7y^r%wM}>-!lcG znS%;YABTLa6|M<<p?hdb)*(h=dyM`>e%})-=a+-v;l6)GC^>DA7(Jwpk&E`tAXJ z?#S5nvp_!7@@MxJm4y7GDekm&BzAGfe;d1=bIa*NJYR@BuBMRy^(lDXtUuTe^~w0s zx3;fS4i)d;T^{j*h2kAA$Ix57;Bf z-U{Pul6zIn#UJY9VyM-%wgU5+iEA;Jf)>1=XBrAK7xTH#hu`d${(X{?V1L#*OdOR4>(FyHyaNcUv{Riw*|f23Rl`=!|-p2TJ@;C92%KHIy4Jat&1ztpZv=_LfHqVeo8r!&WiNEm(D<=k>X zH20ADzOrNjx{Gvd?CJCvB4?4`%hX+kh-S!0N8WAPk(tS)mMcC-HrX~W&roktsUUwD6ED4qTb``&0hlf?ba*8lcN#P zK2mAXD=%7NcH(RPwbwUsJ;+yjkE|l~4F%EH$M40lDh(r)Lm^X~fS*icKO`M>kwULt zTTPxNf%@={hVKi`O2p^ecLolX!5w)QN%{tE1+YFWYsl_ynMyKkUq;6{h2yxj;gdRUV7<6$9`;d3)~jpAl%EiHDhlMP7)X zfGv9_^MCfsriF6k0CNe}bO*#IkaN9b4CTS)$c5qe zk?XFo|HWdd@pPpyjIZ?r>U_*ukpBXIsB_$PfOu-046J`$2>m5FmbANv6Z*@T&fp{# zD}kPpwzw|*a~P42kw++}j-pq%L>OJE!2E${MD3Fru&?yjr!`StkWW1>sfT$&en$5V zb;lfr{z@D86YnSs^H;v?66e={=Bb571QxiLpug1F80@?E!+hFU=UKygOBmIOZzx=S zKtvb|LoxIbKbnYDB1h73l(T?N+w#B|a)8-yY%Z<}Irhr*SQ{h6kISL4Z&x*8KBc-9 zDL1$X=h1yAZSPoxK|ajCRGPRZ3-kM?{v1`DG1NyCycN>aiK4soRcsoPN05NWAC$Ot zgwOD1a$b#;JUuoI6~_>KFsM}gQQi~hsQ@jerEqQsNwGp^~tJ!jP(V>di6S2 zdGneB#Lw^SzYOE!&^`<@ZGJq7P@la*Jbot(&!Ep+PUeP>f`0b1ucb#?`BBN-&iWcI z0=nNmDs;c{B*JfZNwWdW$JaQ!5>~C$c z-sL&y0{L*PNu3m?2Ki8q_OjQ-4K>u%@zB`UdLq)fZ_SD{a~M@KHgd`)sH34Cr`8d6iQq3JNC@Iaiq(~ z|6)LDIb!rXHrUHaa{4H=kEUDb554tK$P7T>u z@z_7_yL^#ge|JXw6xvcP-em?4lvK9QZ&-K_{Yx_Ju*x!1GhGXgM0?glePy6LAe}wupC>^KnI|bvt zPc5Cmh}-EsU-(g*%F$Pgkp#4+gyzLL;TfcA>E;=cV>O6&5nbA^ za;Q(bsYHQ;3XFH0n`1b(9pc5KaoGEZKD5`(Wu5vzUtxU3JbY5{oCND1fvf79}L=A5R>;8s8XYJ3fK@y4pg6)l?vV^Q^6j0T9n2Y*PcS zd=Sr1V|n^^Pr>iyNPpgx;MBc4{w>s$)w&`5Y3FAP)=tKSWoRF!h*rtb4mtE`N>*E7 z7!mp26X^9WTL7JStowp^Ndt{KB*lMuX$ILhYd0~*3*Im3Z{K2VfqdAkt>Jll1lot; zeO;G;1H_Bf?Xm77rcfU{iOearf96F_=o1_cmcx8euB&MMn^y^yvR&oPEFdDH?anj) z`-RY#;gZ(8V85eY&DviTjd|qi-DX7&u%AyuwR^4QpZ(3(u4jCFP7C`biIqb`J%&&p zA6xsaw^UG{f?NG}xL2UR99LNB^{znvX}!@LwRRQklg@Zc$~Xu1xi;#AohuhW52gJQ zU6ml9vUdiU6W+`tQFttpDeIf(DVM`;*rp zKWZ2hzmVjC^}LNqWn=6j^w-1PQgOq)Vki=lFCSh_MC#=GZ>1#(qUy?#hTqR>pud<8 zoPHv?i0r?V@$P_dHS(W_tu@^g^jG}t{`?PLpue(Z7@cJnVE#H}Z6ZvNfqZB9)!&{< z2-@fN@Aj_M3lPtk=zl{$Z%LqHf+c)E+K5QC8Aq)Y6F>Txb2*Y#2S_jQk;jA;+UMQ}m#4eqVZCA2Jk>4p9-e0`irhFzKL+a$A?kR9pY{C?o<-)J z|E@7Zz7$zDea&BvqSqh43!->5g!q`#)HdYvp%lJ(4u`i@(D<97g{eYwh*42K-@f)T z#LLuiO?rG~M_xg4?W0=bs2zSLDY(!JU4`d|uZ`o2xXdAbTDyN4Hn~851qi9R-p_;n zivP0O_f1v<_3QJ`5?UQbiW6A+62Uy*-49aYJ)quD`e>bh@%IeE%l1*ZceD~I-uSCP z(*xuEBBOC^;Xm`u{faLwD=1;SSL?lU-rEi9Ly5&+>AlfVpU#$k*I^=zFM;6t;SZfs zsCMt>y{}?KQk3a z$k2&~coy`+*`;s5^OmBhy}jm9@OvQbD$%TA-r^t`t|$1foF;wfh@8C4_QaYXs@#|%65{Ij&NQ$nvH|2(ay3U!Nz z`C*kIx`{Rk`a?ETEi2=n^IraI`&0Q-A$}@$oqWpLDUS9YwlG7!5)sl?>GNO@@V*O8 zKjkUF&owy~lMgrNky6Va#b>~M%8^$x%iCNq-aDv|OCJq}@t!d+{UY=kte2B0o>q)C zL4T>=nbtFIL|y`CeZT5f`HEG z>(4k~i-;$UYQeq4N@QWrC)=z4=;PQMYjDdH@?oyw^O#?MpnXz&=yldb;Q3))JeT9o zOP_b*GjzgR_^kze&nPRT1@(Z^>PNe*MvJ@uO(JtMKvBSZl2p@!)m+x|X>Q1{tP zXV_LEF7ctlU8(T<%WGxD>wFmG!#F%iyeI?u%Y2(!&S(bK&&vI3hr&<3-LcQP7UjpP z!SMSiocc|y2cn2dWaaod-ytG}Z)A=5t{p<17d?+(^;JgcC@5#6Nb^YjudI^1v(-q_ z{_~I6vmqb8?hVto-g9rq9-|89jbHS@d|EJjy7Z+1tcTJc+laf>!+5u?sBgt=?EC4m805nTx93g>tB+1+gAqiVN9;;pO>tV57p`t=*Gc( z$^V@<$TY_OzC8$kUs6DmxVd=U4t_Mg7Bfp8g8v_RonoQ7ObShEFm<#JAtK?N3L?K< z1W|tN%3%{94fLha#;|GG0`kLu#-b3dM$C_3S>ETNeMCrlF%#R6e@q0`74lLce)l;1 zmrfxG`3JQJZLtR+{{$po7gah4E2iq=`n4ENyu>7}Okgfk^DzejwX zKAJTT^?4g?aKlwe5jENU?!OXn9{yR9C+Np#us`+15Bo11r_f-1@dH+3O9)q#OwP@Z zRS5N1>-^Uu=r2)&=rx@QsLvPLnz!GDR(J4o<7&=&&jVN=`g!RumaD=0WjvPEZtpJ0 zhf>$Fk9;_-f>tnOh8wRC5q6#92A0c5P>)rfD?yzY`cjiSwBgMRa=)!CwaKj-`E*B~ zwx9&^AumVv_TICQ4{t6ySwFc*y@OYubuEYQ523y4)Gl51*#qmtg1!T-D_7w6_2HEA z@N#ouw1vTdS|ec?ab^$tubYw&9U499mTszoc3Y~rM=>oTWn0Us3muinSL^egkDft3 zj0#$sYDZvw=J34Ah5Vd*TRX?q)fp#gBu{ww^AkT=z;O6JmNJ53*4ng-cAr{K;21K6>w;_aUVOU0s{BylLeZjF^Jpt8A|eU4*>>F%L`N_Y ztu@wD=n5TH^B-Xv`K5GTf7gR51f%v|9i)cyh2{oh${%VWo*x((JG5Cr{-I<}I;61; z@idV7QBIo%#(Ve|UJ9fL=8I(c?!c%pb#y7<)V?u)BI5JxN}tI1QB+=~(xb-pG>W>O zpG@W_A(`%!SuxkEkU%SC`l3`A?_21`K&B4N7i-L+fdcFhKj!a08UHv4e_wTVA5Xw5 zXdk8232OgY81I}R%rSj^O6YK`oxC>q9;`-&pR|n%qMx5}Wa}OW`@b(H5f2=hL*%ic zXWbX75Ss@UI)39Y-o=b>OQb%6_DQx8eQ@-j`Q|hGnG6M#;k>iU7DIBoC*+^=p6;^T zhai6Rv;7Z`xXPinDi&Aiyotya+Npd!k|0_v?y9=FMnILGNv+(vJc-bGhU_wu1n>L$ zA>H&hA)d2;$WaaKhWHU5eOzpO2-ZIu{+ew5A`n0R_o}Dc=Hd5g>W1$1W`B^U#}aG~hF2k*`Y!ieBVhhov;08M-Ua!P zV%LbH+d-JGIII(UUbez~Q94MY_*Xc8hcB5VQiBI8;Jop*V*%lW6k)Vy^8}N#=P)A2 zE|zO7bPzSs82cW-s)#br7BF_dSVENNiwW&A z9qaO5LmbwtpYrYBAC80d&u0H!3s(jBJy?B0VU~DI7QMEMI|vyD{Xf4y>swqDLVuip zuW!m{({&+m?0XmqtgC9zEDIxFF|3BYMo@*x8ujGIK2Vcfo&;S4c literal 0 HcmV?d00001 diff --git a/tests/deplete_tests/test_utilities.py b/tests/deplete_tests/test_utilities.py new file mode 100644 index 0000000000..faa1a500b4 --- /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()