From 6c78d129c3db8b09b0e8b21cb3a4965460127666 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 12 Jul 2022 15:49:40 +0100 Subject: [PATCH 01/16] added bbox by @pshriwise --- openmc/universe.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/openmc/universe.py b/openmc/universe.py index f3f95aebb..bb1554de8 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -7,6 +7,7 @@ from pathlib import Path from tempfile import TemporaryDirectory from xml.etree import ElementTree as ET +import h5py import numpy as np import openmc @@ -630,6 +631,9 @@ class DAGMCUniverse(UniverseBase): auto_mat_ids : bool Set IDs automatically on initialization (True) or report overlaps in ID space between OpenMC and UWUW materials (False) + bounding_box : 2-tuple of numpy.array + Lower-left and upper-right coordinates of an axis-aligned bounding box + of the universe. """ def __init__(self, @@ -650,6 +654,15 @@ class DAGMCUniverse(UniverseBase): string += '{: <16}=\t{}\n'.format('\tFile', self.filename) return string + @property + def bounding_box(self): + dagmc_file = h5py.File(self.filename) + coords = dagmc_file['tstt']['nodes']['coordinates'][()] + lower_left_corner = coords.min(axis=0) + upper_right_corner = coords.max(axis=0) + bounding_box = (lower_left_corner, upper_right_corner) + return bounding_box + @property def filename(self): return self._filename From 48a3484c39ebd5737e1ced1a576160b0efb78c69 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 12 Jul 2022 16:36:43 +0100 Subject: [PATCH 02/16] context manager review suggestion from @pshriwise --- openmc/universe.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index bb1554de8..e7be46822 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -656,12 +656,12 @@ class DAGMCUniverse(UniverseBase): @property def bounding_box(self): - dagmc_file = h5py.File(self.filename) - coords = dagmc_file['tstt']['nodes']['coordinates'][()] - lower_left_corner = coords.min(axis=0) - upper_right_corner = coords.max(axis=0) - bounding_box = (lower_left_corner, upper_right_corner) - return bounding_box + with h5py.File(self.filename) as dagmc_file: + coords = dagmc_file['tstt']['nodes']['coordinates'][()] + lower_left_corner = coords.min(axis=0) + upper_right_corner = coords.max(axis=0) + bounding_box = (lower_left_corner, upper_right_corner) + return bounding_box @property def filename(self): From 8b06863b37ff3a258080a33058ceb9311f15dea3 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 12 Jul 2022 16:52:02 +0100 Subject: [PATCH 03/16] added dagmc bounding box test --- tests/regression_tests/dagmc/universes/test.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/regression_tests/dagmc/universes/test.py b/tests/regression_tests/dagmc/universes/test.py index 01963986c..6726b4b44 100644 --- a/tests/regression_tests/dagmc/universes/test.py +++ b/tests/regression_tests/dagmc/universes/test.py @@ -46,6 +46,11 @@ class DAGMCUniverseTest(PyAPITestHarness): # create the DAGMC universe pincell_univ = openmc.DAGMCUniverse(filename='dagmc.h5m', auto_geom_ids=True) + # checks that the bounding box is calculated correctly + bounding_box = pincell_univ.bounding_box + assert bounding_box[0].tolist() == [-25., -25., -25.] + assert bounding_box[1].tolist() == [25., 25., 25.] + # create a 2 x 2 lattice using the DAGMC pincell pitch = np.asarray((24.0, 24.0)) lattice = openmc.RectLattice() From 4784be1fc96d5b688681bedf22c752a55b719b0a Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 12 Jul 2022 21:03:24 +0100 Subject: [PATCH 04/16] review line removal from @paulromano Co-authored-by: Paul Romano --- openmc/universe.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index e7be46822..b65263ef9 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -660,8 +660,7 @@ class DAGMCUniverse(UniverseBase): coords = dagmc_file['tstt']['nodes']['coordinates'][()] lower_left_corner = coords.min(axis=0) upper_right_corner = coords.max(axis=0) - bounding_box = (lower_left_corner, upper_right_corner) - return bounding_box + return (lower_left_corner, upper_right_corner) @property def filename(self): From 39b008df26a449c737caa41238725eb4ed47fd04 Mon Sep 17 00:00:00 2001 From: yardasol Date: Tue, 12 Jul 2022 16:07:59 -0500 Subject: [PATCH 05/16] Remove __enter__, __exit__ methods from TransportOperator This commit removes the context manager for switching directories from the TransportOperator class, and moves this functionality to a module-level function, change_directory. --- openmc/deplete/abc.py | 45 +++++++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index a4f80e3d6..a7cde4afa 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -10,6 +10,7 @@ from collections.abc import Iterable, Callable from copy import deepcopy from inspect import signature from numbers import Real, Integral +from contextlib import contextmanager import os from pathlib import Path import sys @@ -56,6 +57,23 @@ except AttributeError: # Can't set __doc__ on properties on Python 3.4 pass +@contextmanager +def change_directory(output_dir): + """ + Helper function for managing the current directory. + + Parameters + ---------- + output_dir : pathlib.Path + Directory to switch to. + """ + _orig_dir = os.getcwd() + try: + output_dir.mkdir(exist_ok=True) + os.chdir(output_dir) + yield + finally: + os.chdir(_orig_dir) class TransportOperator(ABC): """Abstract class defining a transport operator @@ -87,9 +105,14 @@ class TransportOperator(ABC): Initial atom density [atoms/cm^3] to add for nuclides that are zero in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. + output_dir : pathlib.Path + Path to output directory to save results. prev_res : Results or None Results from a previous depletion calculation. ``None`` if no results are to be used. + chain : openmc.deplete.Chain + The depletion chain information necessary to form matrices and tallies. + """ def __init__(self, chain_file, fission_q=None, dilute_initial=1.0e3, prev_results=None): @@ -133,18 +156,6 @@ class TransportOperator(ABC): """ - def __enter__(self): - # Save current directory and move to specific output directory - self._orig_dir = os.getcwd() - self.output_dir.mkdir(exist_ok=True) - os.chdir(self.output_dir) - - return self.initial_condition() - - def __exit__(self, exc_type, exc_value, traceback): - self.finalize() - os.chdir(self._orig_dir) - @property def output_dir(self): return self._output_dir @@ -775,7 +786,8 @@ class Integrator(ABC): .. versionadded:: 0.13.1 """ - with self.operator as conc: + with change_directory(self.operator.output_dir): + conc = self.operator.initial_condition() t, self._i_res = self._get_start_data() for i, (dt, source_rate) in enumerate(self): @@ -814,6 +826,8 @@ class Integrator(ABC): source_rate, self._i_res + len(self), proc_time) self.operator.write_bos_data(len(self) + self._i_res) + self.operator.finalize() + @add_params class SIIntegrator(Integrator): @@ -931,7 +945,8 @@ class SIIntegrator(Integrator): .. versionadded:: 0.13.1 """ - with self.operator as conc: + with change_directory(self.operator.output_dir): + conc = self.operator.initial_condition() t, self._i_res = self._get_start_data() for i, (dt, p) in enumerate(self): @@ -967,6 +982,8 @@ class SIIntegrator(Integrator): p, self._i_res + len(self), proc_time) self.operator.write_bos_data(self._i_res + len(self)) + self.operator.finalize() + class DepSystemSolver(ABC): r"""Abstract class for solving depletion equations From f5db4e3afb41ff1b56639ccf157e0256b0b5a19f Mon Sep 17 00:00:00 2001 From: yardasol Date: Tue, 12 Jul 2022 18:30:29 -0500 Subject: [PATCH 06/16] move reuasble code to new parent class; 1st pass - new class OpenMCOperator to contain shared bits of code between the currently existing Operator class and the as-of-yet to be created FluxDepletionOperator class - Operator is now a subclass of OpenMCOperator - move _distribute method to OpenMCOperator - move _get_burnable_mats, _extract_number, _set_number_from_mat, _set_number_from_results to OpenMCOperator - split initial condition into openmc.lib dependent and non-dependent parts; dependent part goes to Operator, non-dependent part goes to OpenMCOperator - move _update_materials to OpenMCOperator - move _get_tally_nuclides to OpenMCOperator - rename _unpack_tallies_and_normalize to _calculate_reaction_rates; move to OpenMCOperator - move get_results_info to OpenMCOperator --- openmc/deplete/__init__.py | 1 + openmc/deplete/openmc_operator.py | 522 ++++++++++++++++++++++++++++++ openmc/deplete/operator.py | 331 +------------------ 3 files changed, 535 insertions(+), 319 deletions(-) create mode 100644 openmc/deplete/openmc_operator.py diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 5891555a2..95264308c 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -7,6 +7,7 @@ A depletion front-end tool. from .nuclide import * from .chain import * +from .openmc_operator import * from .operator import * from .reaction_rates import * from .atom_number import * diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py new file mode 100644 index 000000000..303a5208f --- /dev/null +++ b/openmc/deplete/openmc_operator.py @@ -0,0 +1,522 @@ +"""OpenMC transport operator + +This module implements a transport operator for OpenMC so that it can be used by +depletion integrators. The implementation makes use of the Python bindings to +OpenMC's C API so that reading tally results and updating material number +densities is all done in-memory instead of through the filesystem. + +""" + +import copy +from abc import abstractmethod +from collections import OrderedDict +import os +from warnings import warn + +import numpy as np + +import openmc +from openmc.exceptions import DataError +from openmc.mpi import comm +from .abc import TransportOperator, OperatorResult +from .atom_number import AtomNumber +from .reaction_rates import ReactionRates +from .results import Results +from .helpers import ( + DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, + FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, + SourceRateHelper, FluxCollapseHelper) + + +__all__ = ["OpenMCOperator", "OperatorResult"] + + +def _distribute(items): + """Distribute items across MPI communicator + + Parameters + ---------- + items : list + List of items of distribute + + Returns + ------- + list + Items assigned to process that called + + """ + min_size, extra = divmod(len(items), comm.size) + j = 0 + for i in range(comm.size): + chunk_size = min_size + int(i < extra) + if comm.rank == i: + return items[j:j + chunk_size] + j += chunk_size + +class OpenMCOperator(TransportOperator): + """OpenMC transport operator for depletion. + + Instances of this class can be used to perform depletion using OpenMC as the + transport operator. Normally, a user needn't call methods of this class + directly. Instead, an instance of this class is passed to an integrator + class, such as :class:`openmc.deplete.CECMIntegrator`. + + .. versionchanged:: 0.13.0 + The geometry and settings parameters have been replaced with a + model parameter that takes a :class:`~openmc.model.Model` object + + Parameters + ---------- + chain_file : str, optional + Path to the depletion chain XML file. Defaults to the file + listed under ``depletion_chain`` in + :envvar:`OPENMC_CROSS_SECTIONS` environment variable. + prev_results : Results, optional + Results from a previous depletion calculation. If this argument is + specified, the depletion calculation will start from the latest state + in the previous results. + fission_q : dict, optional + Dictionary of nuclides and their fission Q values [eV]. If not given, + values will be pulled from the ``chain_file``. Only applicable + if ``"normalization_mode" == "fission-q"`` + dilute_initial : float, optional + Initial atom density [atoms/cm^3] to add for nuclides that are zero + in initial condition to ensure they exist in the decay chain. + Only done for nuclides with reaction rates. + Defaults to 1.0e3. + + Attributes + ---------- + model : openmc.model.Model + OpenMC model object + geometry : openmc.Geometry + OpenMC geometry object + settings : openmc.Settings + OpenMC settings object + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the decay + chain. Only done for nuclides with reaction rates. + output_dir : pathlib.Path + Path to output directory to save results. + round_number : bool + Whether or not to round output to OpenMC to 8 digits. + Useful in testing, as OpenMC is incredibly sensitive to exact values. + number : openmc.deplete.AtomNumber + Total number of atoms in simulation. + nuclides_with_data : set of str + A set listing all unique nuclides available from cross_sections.xml. + chain : openmc.deplete.Chain + The depletion chain information necessary to form matrices and tallies. + reaction_rates : openmc.deplete.ReactionRates + Reaction rates from the last operator step. + burnable_mats : list of str + All burnable material IDs + heavy_metal : float + Initial heavy metal inventory [g] + local_mats : list of str + All burnable material IDs being managed by a single process + prev_res : Results or None + Results from a previous depletion calculation. ``None`` if no + results are to be used. + cleanup_when_done : bool + Whether to finalize and clear the shared library memory when the + depletion operation is complete. Defaults to clearing the library. + """ + + ## modify + def __init__(self, chain_file=None, fission_q=None, dilute_initial=0.0, prev_results=None): + + super().__init__(chain_file, fission_q, dilute_initial, prev_results) + + ## clear, but we must keep the method + def __call__(self, vec, source_rate): + """Runs a simulation. + + Simulation will abort under the following circumstances: + + 1) No energy is computed using OpenMC tallies. + + Parameters + ---------- + vec : list of numpy.ndarray + Total atoms to be used in function. + source_rate : float + Power in [W] or source rate in [neutron/sec] + + Returns + ------- + openmc.deplete.OperatorResult + Eigenvalue and reaction rates resulting from transport operator + + """ + + ## clear, but must keep the method + def write_bos_data(step): + """Write a state-point file with beginning of step data + + Parameters + ---------- + step : int + Current depletion step including restarts + """ + pass + + ## keep + def _get_burnable_mats(self): + """Determine depletable materials, volumes, and nuclides + + Returns + ------- + burnable_mats : list of str + List of burnable material IDs + volume : OrderedDict of str to float + Volume of each material in [cm^3] + nuclides : list of str + Nuclides in order of how they'll appear in the simulation. + + """ + + burnable_mats = set() + model_nuclides = set() + volume = OrderedDict() + + self.heavy_metal = 0.0 + + # Iterate once through the geometry to get dictionaries + for mat in self.materials: + for nuclide in mat.get_nuclides(): + model_nuclides.add(nuclide) + if mat.depletable: + burnable_mats.add(str(mat.id)) + if mat.volume is None: + raise RuntimeError("Volume not specified for depletable " + "material with ID={}.".format(mat.id)) + volume[str(mat.id)] = mat.volume + self.heavy_metal += mat.fissionable_mass + + # Make sure there are burnable materials + if not burnable_mats: + raise RuntimeError( + "No depletable materials were found in the model.") + + # Sort the sets + burnable_mats = sorted(burnable_mats, key=int) + model_nuclides = sorted(model_nuclides) + + # Construct a global nuclide dictionary, burned first + nuclides = list(self.chain.nuclide_dict) + for nuc in model_nuclides: + if nuc not in nuclides: + nuclides.append(nuc) + + return burnable_mats, volume, nuclides + + ## keep + def _extract_number(self, local_mats, volume, nuclides, prev_res=None): + """Construct AtomNumber using geometry + + Parameters + ---------- + local_mats : list of str + Material IDs to be managed by this process + volume : OrderedDict of str to float + Volumes for the above materials in [cm^3] + nuclides : list of str + Nuclides to be used in the simulation. + prev_res : Results, optional + Results from a previous depletion calculation + + """ + self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain)) + + if self.dilute_initial != 0.0: + for nuc in self._burnable_nucs: + self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) + + # Now extract and store the number densities + # From the geometry if no previous depletion results + if prev_res is None: + for mat in self.materials: + if str(mat.id) in local_mats: + self._set_number_from_mat(mat) + + # Else from previous depletion results + else: + for mat in self.materials: + if str(mat.id) in local_mats: + self._set_number_from_results(mat, prev_res) + + def _set_number_from_mat(self, mat): + """Extracts material and number densities from openmc.Material + + Parameters + ---------- + mat : openmc.Material + The material to read from + + """ + mat_id = str(mat.id) + + for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): + atom_per_cc = atom_per_bcm * 1.0e24 + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + + def _set_number_from_results(self, mat, prev_res): + """Extracts material nuclides and number densities. + + If the nuclide concentration's evolution is tracked, the densities come + from depletion results. Else, densities are extracted from the geometry + in the summary. + + Parameters + ---------- + mat : openmc.Material + The material to read from + prev_res : Results + Results from a previous depletion calculation + + """ + mat_id = str(mat.id) + + # Get nuclide lists from geometry and depletion results + depl_nuc = prev_res[-1].nuc_to_ind + geom_nuc_densities = mat.get_nuclide_atom_densities() + + # Merge lists of nuclides, with the same order for every calculation + geom_nuc_densities.update(depl_nuc) + + for nuclide, atom_per_bcm in geom_nuc_densities.items(): + if nuclide in depl_nuc: + concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] + volume = prev_res[-1].volume[mat_id] + atom_per_cc = concentration / volume + else: + atom_per_cc = atom_per_bcm * 1.0e24 + + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + + ## keep, will need to modify + def initial_condition(self, materials): + """Performs final setup and returns initial condition. + + Parameters + ---------- + materials : ??? + + Returns + ------- + list of numpy.ndarray + Total density for initial conditions. + """ + + self._rate_helper.generate_tallies(materials, self.chain.reactions) + self._normalization_helper.prepare( + self.chain.nuclides, self.reaction_rates.index_nuc) + # Tell fission yield helper what materials this process is + # responsible for + self._yield_helper.generate_tallies( + materials, tuple(sorted(self._mat_index_map.values()))) + + # Return number density vector + return list(self.number.get_mat_slice(np.s_[:])) + + ## keep? there is a non-removable part of this + ## that needs the openmc c executable + 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.materials: + nuclides = [] + densities = [] + for nuc in number_i.nuclides: + if nuc in self.nuclides_with_data: + 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.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 + + # Update densities on C API side + mat_internal = openmc.lib.materials[int(mat)] + mat_internal.set_densities(nuclides, densities) + + #TODO Update densities on the Python side, otherwise the + # summary.h5 file contains densities at the first time step + + ## keep, but rename and modify docstring + ## get_tally_nuclides + def _get_tally_nuclides(self): + """Determine nuclides that should be tallied for reaction rates. + + This method returns a list of all nuclides that have neutron data and + are listed in the depletion chain. Technically, we should tally nuclides + that may not appear in the depletion chain because we still need to get + the fission reaction rate for these nuclides in order to normalize + power, but that is left as a future exercise. + + Returns + ------- + list of str + Tally nuclides + + """ + nuc_set = set() + + # Create the set of all nuclides in the decay chain in materials marked + # for burning in which the number density is greater than zero. + for nuc in self.number.nuclides: + if nuc in self.nuclides_with_data: + if np.sum(self.number[:, nuc]) > 0.0: + nuc_set.add(nuc) + + # Communicate which nuclides have nonzeros to rank 0 + if comm.rank == 0: + for i in range(1, comm.size): + nuc_newset = comm.recv(source=i, tag=i) + nuc_set |= nuc_newset + else: + comm.send(nuc_set, dest=0, tag=comm.rank) + + if comm.rank == 0: + # Sort nuclides in the same order as self.number + nuc_list = [nuc for nuc in self.number.nuclides + if nuc in nuc_set] + else: + nuc_list = None + + # Store list of tally nuclides on each process + nuc_list = comm.bcast(nuc_list) + return [nuc for nuc in nuc_list if nuc in self.chain] + + ## modify + ## the tally-specific methods + ## for the normalization and fission helper classes + ## will already work with the current implementation. + ## For transport-less depletion, we'll need to write + ## a new ReactionRateHelper + def _calculate_reaction_rates(self, source_rate): + """Unpack tallies from OpenMC and return an operator result + + This method uses OpenMC's C API bindings to determine the k-effective + value and reaction rates from the simulation. The reaction rates are + normalized by a helper class depending on the method being used. + + Parameters + ---------- + source_rate : float + Power in [W] or source rate in [neutron/sec] + + Returns + ------- + rates : openmc.deplete.ReactionRates + Reaction rates for nuclides + + """ + rates = self.reaction_rates + rates.fill(0.0) + + # Extract tally bins + nuclides = self._rate_helper.nuclides + + # Form fast map + nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] + react_ind = [rates.index_rx[react] for react in self.chain.reactions] + + # Keep track of energy produced from all reactions in eV per source + # particle + self._normalization_helper.reset() + self._yield_helper.unpack() + + # Store fission yield dictionaries + fission_yields = [] + + # Create arrays to store fission Q values, reaction rates, and nuclide + # numbers, zeroed out in material iteration + number = np.empty(rates.n_nuc) + + fission_ind = rates.index_rx.get("fission") + + # Extract results + for i, mat in enumerate(self.local_mats): + # Get tally index + mat_index = self._mat_index_map[mat] + + # Zero out reaction rates and nuclide numbers + number.fill(0.0) + + # Get new number densities + for nuc, i_nuc_results in zip(nuclides, nuc_ind): + number[i_nuc_results] = self.number[mat, nuc] + + tally_rates = self._rate_helper.get_material_rates( + mat_index, nuc_ind, react_ind) + + # Compute fission yields for this material + fission_yields.append(self._yield_helper.weighted_yields(i)) + + # Accumulate energy from fission + if fission_ind is not None: + self._normalization_helper.update(tally_rates[:, fission_ind]) + + # Divide by total number and store + rates[i] = self._rate_helper.divide_by_adens(number) + + # Scale reaction rates to obtain units of reactions/sec + rates *= self._normalization_helper.factor(source_rate) + + # Store new fission yields on the chain + self.chain.fission_yields = fission_yields + + return rates + + @abstractmethod + def _get_nuclides_with_data(self): + """Find nuclides with cross section data""" + + ## Keep + def get_results_info(self): + """Returns volume list, material lists, and nuc lists. + + Returns + ------- + volume : dict of str float + Volumes corresponding to materials in full_burn_dict + nuc_list : list of str + A list of all nuclide names. Used for sorting the simulation. + burn_list : list of int + A list of all material IDs to be burned. Used for sorting the simulation. + full_burn_list : list + List of all burnable material IDs + + """ + nuc_list = self.number.burnable_nuclides + burn_list = self.local_mats + + volume = {} + for i, mat in enumerate(burn_list): + volume[mat] = self.number.volume[i] + + # Combine volume dictionaries across processes + volume_list = comm.allgather(volume) + volume = {k: v for d in volume_list for k, v in d.items()} + + return volume, nuc_list, burn_list, self.burnable_mats diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 7ef874e87..665de1999 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -24,6 +24,7 @@ from openmc.mpi import comm from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .chain import _find_chain_file +from .openmc_operator import OpenMCOperator, _distribute from .reaction_rates import ReactionRates from .results import Results from .helpers import ( @@ -34,30 +35,6 @@ from .helpers import ( __all__ = ["Operator", "OperatorResult"] - -def _distribute(items): - """Distribute items across MPI communicator - - Parameters - ---------- - items : list - List of items of distribute - - Returns - ------- - list - Items assigned to process that called - - """ - min_size, extra = divmod(len(items), comm.size) - j = 0 - for i in range(comm.size): - chunk_size = min_size + int(i < extra) - if comm.rank == i: - return items[j:j + chunk_size] - j += chunk_size - - def _find_cross_sections(model): """Determine cross sections to use for depletion""" if model.materials and model.materials.cross_sections is not None: @@ -74,7 +51,7 @@ def _find_cross_sections(model): return cross_sections -class Operator(TransportOperator): +class Operator(OpenMCOperator): """OpenMC transport operator for depletion. Instances of this class can be used to perform depletion using OpenMC as the @@ -267,7 +244,7 @@ class Operator(TransportOperator): # Clear out OpenMC, create task lists, distribute openmc.reset_auto_ids() - self.burnable_mats, volume, nuclides = self._get_burnable_mats() + self.burnable_mats, volume, nuclides = super()._get_burnable_mats() self.local_mats = _distribute(self.burnable_mats) # Generate map from local materials => material index @@ -298,7 +275,7 @@ class Operator(TransportOperator): if nuc.name in self.nuclides_with_data] # Extract number densities from the geometry / previous depletion run - self._extract_number(self.local_mats, volume, nuclides, self.prev_res) + super()._extract_number(self.local_mats, volume, nuclides, self.prev_res) # Create reaction rates array self.reaction_rates = ReactionRates( @@ -380,7 +357,7 @@ class Operator(TransportOperator): openmc.reset_auto_ids() # Update tally nuclides data in preparation for transport solve - nuclides = self._get_tally_nuclides() + nuclides = super()._get_tally_nuclides() self._rate_helper.nuclides = nuclides self._normalization_helper.nuclides = nuclides self._yield_helper.update_tally_nuclides(nuclides) @@ -390,7 +367,12 @@ class Operator(TransportOperator): openmc.lib.reset_timers() # Extract results - op_result = self._unpack_tallies_and_normalize(source_rate) + rates = self._calculate_reaction_rates(source_rate) + + # Get k and uncertainty + keff = ufloat(*openmc.lib.keff()) + + op_result = OperatorResult(keff, rates) return copy.deepcopy(op_result) @@ -434,138 +416,6 @@ class Operator(TransportOperator): cell.fill = [mat.clone() for i in range(cell.num_instances)] - def _get_burnable_mats(self): - """Determine depletable materials, volumes, and nuclides - - Returns - ------- - burnable_mats : list of str - List of burnable material IDs - volume : OrderedDict of str to float - Volume of each material in [cm^3] - nuclides : list of str - Nuclides in order of how they'll appear in the simulation. - - """ - - burnable_mats = set() - model_nuclides = set() - volume = OrderedDict() - - self.heavy_metal = 0.0 - - # Iterate once through the geometry to get dictionaries - for mat in self.materials: - for nuclide in mat.get_nuclides(): - model_nuclides.add(nuclide) - if mat.depletable: - burnable_mats.add(str(mat.id)) - if mat.volume is None: - raise RuntimeError("Volume not specified for depletable " - "material with ID={}.".format(mat.id)) - volume[str(mat.id)] = mat.volume - self.heavy_metal += mat.fissionable_mass - - # Make sure there are burnable materials - if not burnable_mats: - raise RuntimeError( - "No depletable materials were found in the model.") - - # Sort the sets - burnable_mats = sorted(burnable_mats, key=int) - model_nuclides = sorted(model_nuclides) - - # Construct a global nuclide dictionary, burned first - nuclides = list(self.chain.nuclide_dict) - for nuc in model_nuclides: - if nuc not in nuclides: - nuclides.append(nuc) - - return burnable_mats, volume, nuclides - - def _extract_number(self, local_mats, volume, nuclides, prev_res=None): - """Construct AtomNumber using geometry - - Parameters - ---------- - local_mats : list of str - Material IDs to be managed by this process - volume : OrderedDict of str to float - Volumes for the above materials in [cm^3] - nuclides : list of str - Nuclides to be used in the simulation. - prev_res : Results, optional - Results from a previous depletion calculation - - """ - self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain)) - - if self.dilute_initial != 0.0: - for nuc in self._burnable_nucs: - self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) - - # Now extract and store the number densities - # From the geometry if no previous depletion results - if prev_res is None: - for mat in self.materials: - if str(mat.id) in local_mats: - self._set_number_from_mat(mat) - - # Else from previous depletion results - else: - for mat in self.materials: - if str(mat.id) in local_mats: - self._set_number_from_results(mat, prev_res) - - def _set_number_from_mat(self, mat): - """Extracts material and number densities from openmc.Material - - Parameters - ---------- - mat : openmc.Material - The material to read from - - """ - mat_id = str(mat.id) - - for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): - atom_per_cc = atom_per_bcm * 1.0e24 - self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - - def _set_number_from_results(self, mat, prev_res): - """Extracts material nuclides and number densities. - - If the nuclide concentration's evolution is tracked, the densities come - from depletion results. Else, densities are extracted from the geometry - in the summary. - - Parameters - ---------- - mat : openmc.Material - The material to read from - prev_res : Results - Results from a previous depletion calculation - - """ - mat_id = str(mat.id) - - # Get nuclide lists from geometry and depletion results - depl_nuc = prev_res[-1].nuc_to_ind - geom_nuc_densities = mat.get_nuclide_atom_densities() - - # Merge lists of nuclides, with the same order for every calculation - geom_nuc_densities.update(depl_nuc) - - for nuclide, atom_per_bcm in geom_nuc_densities.items(): - if nuclide in depl_nuc: - concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] - volume = prev_res[-1].volume[mat_id] - atom_per_cc = concentration / volume - else: - atom_per_cc = atom_per_bcm * 1.0e24 - - self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - def initial_condition(self): """Performs final setup and returns initial condition. @@ -589,16 +439,8 @@ class Operator(TransportOperator): # Generate tallies in memory materials = [openmc.lib.materials[int(i)] for i in self.burnable_mats] - self._rate_helper.generate_tallies(materials, self.chain.reactions) - self._normalization_helper.prepare( - self.chain.nuclides, self.reaction_rates.index_nuc) - # Tell fission yield helper what materials this process is - # responsible for - self._yield_helper.generate_tallies( - materials, tuple(sorted(self._mat_index_map.values()))) - # Return number density vector - return list(self.number.get_mat_slice(np.s_[:])) + return super().initial_condition(materials) def finalize(self): """Finalize a depletion simulation and release resources.""" @@ -659,127 +501,6 @@ class Operator(TransportOperator): self.materials.export_to_xml() - def _get_tally_nuclides(self): - """Determine nuclides that should be tallied for reaction rates. - - This method returns a list of all nuclides that have neutron data and - are listed in the depletion chain. Technically, we should tally nuclides - that may not appear in the depletion chain because we still need to get - the fission reaction rate for these nuclides in order to normalize - power, but that is left as a future exercise. - - Returns - ------- - list of str - Tally nuclides - - """ - nuc_set = set() - - # Create the set of all nuclides in the decay chain in materials marked - # for burning in which the number density is greater than zero. - for nuc in self.number.nuclides: - if nuc in self.nuclides_with_data: - if np.sum(self.number[:, nuc]) > 0.0: - nuc_set.add(nuc) - - # Communicate which nuclides have nonzeros to rank 0 - if comm.rank == 0: - for i in range(1, comm.size): - nuc_newset = comm.recv(source=i, tag=i) - nuc_set |= nuc_newset - else: - comm.send(nuc_set, dest=0, tag=comm.rank) - - if comm.rank == 0: - # Sort nuclides in the same order as self.number - nuc_list = [nuc for nuc in self.number.nuclides - if nuc in nuc_set] - else: - nuc_list = None - - # Store list of tally nuclides on each process - nuc_list = comm.bcast(nuc_list) - return [nuc for nuc in nuc_list if nuc in self.chain] - - def _unpack_tallies_and_normalize(self, source_rate): - """Unpack tallies from OpenMC and return an operator result - - This method uses OpenMC's C API bindings to determine the k-effective - value and reaction rates from the simulation. The reaction rates are - normalized by a helper class depending on the method being used. - - Parameters - ---------- - source_rate : float - Power in [W] or source rate in [neutron/sec] - - Returns - ------- - openmc.deplete.OperatorResult - Eigenvalue and reaction rates resulting from transport operator - - """ - rates = self.reaction_rates - rates.fill(0.0) - - # Get k and uncertainty - keff = ufloat(*openmc.lib.keff()) - - # Extract tally bins - nuclides = self._rate_helper.nuclides - - # Form fast map - nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] - react_ind = [rates.index_rx[react] for react in self.chain.reactions] - - # Keep track of energy produced from all reactions in eV per source - # particle - self._normalization_helper.reset() - self._yield_helper.unpack() - - # Store fission yield dictionaries - fission_yields = [] - - # Create arrays to store fission Q values, reaction rates, and nuclide - # numbers, zeroed out in material iteration - number = np.empty(rates.n_nuc) - - fission_ind = rates.index_rx.get("fission") - - # Extract results - for i, mat in enumerate(self.local_mats): - # Get tally index - mat_index = self._mat_index_map[mat] - - # Zero out reaction rates and nuclide numbers - number.fill(0.0) - - # Get new number densities - for nuc, i_nuc_results in zip(nuclides, nuc_ind): - number[i_nuc_results] = self.number[mat, nuc] - - tally_rates = self._rate_helper.get_material_rates( - mat_index, nuc_ind, react_ind) - - # Compute fission yields for this material - fission_yields.append(self._yield_helper.weighted_yields(i)) - - # Accumulate energy from fission - if fission_ind is not None: - self._normalization_helper.update(tally_rates[:, fission_ind]) - - # Divide by total number and store - rates[i] = self._rate_helper.divide_by_adens(number) - - # Scale reaction rates to obtain units of reactions/sec - rates *= self._normalization_helper.factor(source_rate) - - # Store new fission yields on the chain - self.chain.fission_yields = fission_yields - - return OperatorResult(keff, rates) - def _get_nuclides_with_data(self, cross_sections): """Loads cross_sections.xml file to find nuclides with neutron data""" nuclides = set() @@ -792,31 +513,3 @@ class Operator(TransportOperator): nuclides.add(name) return nuclides - - def get_results_info(self): - """Returns volume list, material lists, and nuc lists. - - Returns - ------- - volume : dict of str float - Volumes corresponding to materials in full_burn_dict - nuc_list : list of str - A list of all nuclide names. Used for sorting the simulation. - burn_list : list of int - A list of all material IDs to be burned. Used for sorting the simulation. - full_burn_list : list - List of all burnable material IDs - - """ - nuc_list = self.number.burnable_nuclides - burn_list = self.local_mats - - volume = {} - for i, mat in enumerate(burn_list): - volume[mat] = self.number.volume[i] - - # Combine volume dictionaries across processes - volume_list = comm.allgather(volume) - volume = {k: v for d in volume_list for k, v in d.items()} - - return volume, nuc_list, burn_list, self.burnable_mats From 5d5a2d6fa5fa3c3f6f0fada60115f1e3e7add984 Mon Sep 17 00:00:00 2001 From: Olek <45364492+yardasol@users.noreply.github.com> Date: Wed, 13 Jul 2022 15:55:15 -0500 Subject: [PATCH 07/16] _orig_dir -> orig_dir Co-authored-by: Paul Romano --- openmc/deplete/abc.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index a7cde4afa..8bbe126df 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -67,13 +67,14 @@ def change_directory(output_dir): output_dir : pathlib.Path Directory to switch to. """ - _orig_dir = os.getcwd() + orig_dir = os.getcwd() try: output_dir.mkdir(exist_ok=True) os.chdir(output_dir) yield finally: - os.chdir(_orig_dir) + os.chdir(orig_dir) + class TransportOperator(ABC): """Abstract class defining a transport operator From 35ddf8b23c85e969b41e47caba65fa91099626a2 Mon Sep 17 00:00:00 2001 From: yardasol Date: Wed, 13 Jul 2022 17:17:26 -0500 Subject: [PATCH 08/16] move reusable code to new parent class; 2nd pass; some minor docstring edits --- openmc/deplete/openmc_operator.py | 167 +++++++++++++---------------- openmc/deplete/operator.py | 168 +++++++++++++----------------- 2 files changed, 143 insertions(+), 192 deletions(-) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 303a5208f..a2296d88a 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -1,9 +1,6 @@ """OpenMC transport operator -This module implements a transport operator for OpenMC so that it can be used by -depletion integrators. The implementation makes use of the Python bindings to -OpenMC's C API so that reading tally results and updating material number -densities is all done in-memory instead of through the filesystem. +This module implements functions used by both OpenMC transport operators as well as pure depletion operators. """ @@ -11,9 +8,9 @@ import copy from abc import abstractmethod from collections import OrderedDict import os -from warnings import warn import numpy as np +from uncertainties import ufloat import openmc from openmc.exceptions import DataError @@ -22,10 +19,6 @@ from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .reaction_rates import ReactionRates from .results import Results -from .helpers import ( - DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, - FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, - SourceRateHelper, FluxCollapseHelper) __all__ = ["OpenMCOperator", "OperatorResult"] @@ -124,45 +117,61 @@ class OpenMCOperator(TransportOperator): depletion operation is complete. Defaults to clearing the library. """ - ## modify - def __init__(self, chain_file=None, fission_q=None, dilute_initial=0.0, prev_results=None): + def __init__(self, materials=None, cross_sections=None, chain_file=None, prev_results=None, diff_burnable_mats=False, normalization_mode=None, fission_q=None, dilute_initial=0.0, fission_yield_mode=None, fission_yield_opts=None, + reaction_rate_mode=None, reaction_rate_opts=None, + reduce_chain=False, reduce_chain_level=None): super().__init__(chain_file, fission_q, dilute_initial, prev_results) + self.round_number=False + self.materials = materials + self.cross_sections = cross_sections - ## clear, but we must keep the method - def __call__(self, vec, source_rate): - """Runs a simulation. + # Reduce the chain to only those nuclides present + if reduce_chain: + init_nuclides = set() + for material in self.materials: + if not material.depletable: + continue + for name, _dens_percent, _dens_type in material.nuclides: + init_nuclides.add(name) - Simulation will abort under the following circumstances: + self.chain = self.chain.reduce(init_nuclides, reduce_chain_level) - 1) No energy is computed using OpenMC tallies. + if diff_burnable_mats: + self._differentiate_burnable_mats() - Parameters - ---------- - vec : list of numpy.ndarray - Total atoms to be used in function. - source_rate : float - Power in [W] or source rate in [neutron/sec] + # Determine which nuclides have cross section data + # This nuclides variables contains every nuclides + # for which there is an entry in the micro_xs parameter + openmc.reset_auto_ids() + self.burnable_mats, volumes, all_nuclides = self._get_burnable_mats() + self.local_mats = _distribute(self.burnable_mats) - Returns - ------- - openmc.deplete.OperatorResult - Eigenvalue and reaction rates resulting from transport operator + self._mat_index_map = { + lm: self.burnable_mats.index(lm) for lm in self.local_mats} - """ + if self.prev_res is not None: + self._load_previous_results() - ## clear, but must keep the method - def write_bos_data(step): - """Write a state-point file with beginning of step data - Parameters - ---------- - step : int - Current depletion step including restarts - """ - pass + self.nuclides_with_data = self._get_nuclides_with_data(self.cross_sections) + + # Select nuclides with data that are also in the chain + self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides + if nuc.name in self.nuclides_with_data] + + # Extract number densities from the geometry / previous depletion run + self._extract_number(self.local_mats, + volumes, + all_nuclides, + self.prev_res) + + # Create reaction rates array + self.reaction_rates = ReactionRates( + self.local_mats, self._burnable_nucs, self.chain.reactions) + + self._get_helper_classes(reaction_rate_mode, reaction_rate_opts, normalization_mode, fission_yield_mode, fission_yield_opts) - ## keep def _get_burnable_mats(self): """Determine depletable materials, volumes, and nuclides @@ -212,8 +221,7 @@ class OpenMCOperator(TransportOperator): return burnable_mats, volume, nuclides - ## keep - def _extract_number(self, local_mats, volume, nuclides, prev_res=None): + def _extract_number(self, local_mats, volume, all_nuclides, prev_res=None): """Construct AtomNumber using geometry Parameters @@ -222,13 +230,13 @@ class OpenMCOperator(TransportOperator): Material IDs to be managed by this process volume : OrderedDict of str to float Volumes for the above materials in [cm^3] - nuclides : list of str + all_nuclides : list of str Nuclides to be used in the simulation. prev_res : Results, optional Results from a previous depletion calculation """ - self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain)) + self.number = AtomNumber(local_mats, all_nuclides, volume, len(self.chain)) if self.dilute_initial != 0.0: for nuc in self._burnable_nucs: @@ -296,7 +304,6 @@ class OpenMCOperator(TransportOperator): self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - ## keep, will need to modify def initial_condition(self, materials): """Performs final setup and returns initial condition. @@ -321,50 +328,11 @@ class OpenMCOperator(TransportOperator): # Return number density vector return list(self.number.get_mat_slice(np.s_[:])) - ## keep? there is a non-removable part of this - ## that needs the openmc c executable + @abstractmethod 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.materials: - nuclides = [] - densities = [] - for nuc in number_i.nuclides: - if nuc in self.nuclides_with_data: - 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.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 - - # Update densities on C API side - mat_internal = openmc.lib.materials[int(mat)] - mat_internal.set_densities(nuclides, densities) - - #TODO Update densities on the Python side, otherwise the - # summary.h5 file contains densities at the first time step - - ## keep, but rename and modify docstring - ## get_tally_nuclides - def _get_tally_nuclides(self): + def _get_reaction_nuclides(self): """Determine nuclides that should be tallied for reaction rates. This method returns a list of all nuclides that have neutron data and @@ -407,12 +375,6 @@ class OpenMCOperator(TransportOperator): nuc_list = comm.bcast(nuc_list) return [nuc for nuc in nuc_list if nuc in self.chain] - ## modify - ## the tally-specific methods - ## for the normalization and fission helper classes - ## will already work with the current implementation. - ## For transport-less depletion, we'll need to write - ## a new ReactionRateHelper def _calculate_reaction_rates(self, source_rate): """Unpack tallies from OpenMC and return an operator result @@ -434,11 +396,11 @@ class OpenMCOperator(TransportOperator): rates = self.reaction_rates rates.fill(0.0) - # Extract tally bins - nuclides = self._rate_helper.nuclides + # Extract reaction nuclides + rxn_nuclides = self._rate_helper.nuclides # Form fast map - nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] + nuc_ind = [rates.index_nuc[nuc] for nuc in rxn_nuclides] react_ind = [rates.index_rx[react] for react in self.chain.reactions] # Keep track of energy produced from all reactions in eV per source @@ -464,7 +426,7 @@ class OpenMCOperator(TransportOperator): number.fill(0.0) # Get new number densities - for nuc, i_nuc_results in zip(nuclides, nuc_ind): + for nuc, i_nuc_results in zip(rxn_nuclides, nuc_ind): number[i_nuc_results] = self.number[mat, nuc] tally_rates = self._rate_helper.get_material_rates( @@ -488,11 +450,26 @@ class OpenMCOperator(TransportOperator): return rates + @abstractmethod - def _get_nuclides_with_data(self): + def _get_helper_classes(self, reaction_rate_mode, reaction_rate_opts, normalization_mode, fission_yield_mode, fission_yield_opts): + """Create the ``_rate_helper``, ``normalization_helper``, and + ``_yield_helper`` attributes""" + + @abstractmethod + def _load_previous_results(self): + """Load reuslts from a previous depletion simulation""" + + @abstractmethod + def _get_nuclides_with_data(self, cross_sections): """Find nuclides with cross section data""" - ## Keep + @abstractmethod + def _differentiate_burnable_mats(self): + """Assign distribmats for each burnable material + + """ + def get_results_info(self): """Returns volume list, material lists, and nuc lists. diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 665de1999..2770494ba 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -210,8 +210,6 @@ class Operator(OpenMCOperator): if fission_q is not None: warn("Fission Q dictionary will not be used") fission_q = None - super().__init__(chain_file, fission_q, dilute_initial, prev_results) - self.round_number = False self.model = model self.settings = model.settings self.geometry = model.geometry @@ -221,103 +219,15 @@ class Operator(OpenMCOperator): model.materials = openmc.Materials( model.geometry.get_all_materials().values() ) - self.materials = model.materials self.cleanup_when_done = True - # Reduce the chain before we create more materials - if reduce_chain: - all_isotopes = set() - for material in self.materials: - if not material.depletable: - continue - for name, _dens_percent, _dens_type in material.nuclides: - all_isotopes.add(name) - self.chain = self.chain.reduce(all_isotopes, reduce_chain_level) - - # Differentiate burnable materials with multiple instances - if diff_burnable_mats: - self._differentiate_burnable_mats() - self.materials = openmc.Materials( - model.geometry.get_all_materials().values() - ) - - # Clear out OpenMC, create task lists, distribute - openmc.reset_auto_ids() - self.burnable_mats, volume, nuclides = super()._get_burnable_mats() - self.local_mats = _distribute(self.burnable_mats) - - # Generate map from local materials => material index - self._mat_index_map = { - lm: self.burnable_mats.index(lm) for lm in self.local_mats} - - if self.prev_res is not None: - # Reload volumes into geometry - prev_results[-1].transfer_volumes(self.model) - - # Store previous results in operator - # Distribute reaction rates according to those tracked - # on this process - if comm.size == 1: - self.prev_res = prev_results - else: - self.prev_res = Results() - mat_indexes = _distribute(range(len(self.burnable_mats))) - for res_obj in prev_results: - new_res = res_obj.distribute(self.local_mats, mat_indexes) - self.prev_res.append(new_res) - - # Determine which nuclides have incident neutron data - self.nuclides_with_data = self._get_nuclides_with_data(cross_sections) - - # Select nuclides with data that are also in the chain - self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides - if nuc.name in self.nuclides_with_data] - - # Extract number densities from the geometry / previous depletion run - super()._extract_number(self.local_mats, volume, nuclides, self.prev_res) - - # Create reaction rates array - self.reaction_rates = ReactionRates( - self.local_mats, self._burnable_nucs, self.chain.reactions) - - # Get classes to assist working with tallies - if reaction_rate_mode == "direct": - self._rate_helper = DirectReactionRateHelper( - self.reaction_rates.n_nuc, self.reaction_rates.n_react) - elif reaction_rate_mode == "flux": - if reaction_rate_opts is None: - reaction_rate_opts = {} - - # Ensure energy group boundaries were specified - if 'energies' not in reaction_rate_opts: - raise ValueError( - "Energy group boundaries must be specified in the " - "reaction_rate_opts argument when reaction_rate_mode is" - "set to 'flux'.") - - self._rate_helper = FluxCollapseHelper( - self.reaction_rates.n_nuc, - self.reaction_rates.n_react, - **reaction_rate_opts - ) - else: - raise ValueError("Invalid reaction rate mode.") - - if normalization_mode == "fission-q": - self._normalization_helper = ChainFissionHelper() - elif normalization_mode == "energy-deposition": - score = "heating" if self.settings.photon_transport else "heating-local" - self._normalization_helper = EnergyScoreHelper(score) - else: - self._normalization_helper = SourceRateHelper() - - # Select and create fission yield helper - fission_helper = self._fission_helpers[fission_yield_mode] - fission_yield_opts = ( - {} if fission_yield_opts is None else fission_yield_opts) - self._yield_helper = fission_helper.from_operator( - self, **fission_yield_opts) + super().__init__(model.materials, cross_sections, chain_file, prev_results, + diff_burnable_mats, normalization_mode, + fission_q, dilute_initial, + fission_yield_mode, fission_yield_opts, + reaction_rate_mode, reaction_rate_opts, + reduce_chain, reduce_chain_level) def __call__(self, vec, source_rate): """Runs a simulation. @@ -357,11 +267,12 @@ class Operator(OpenMCOperator): openmc.reset_auto_ids() # Update tally nuclides data in preparation for transport solve - nuclides = super()._get_tally_nuclides() + nuclides = self._get_reaction_nuclides() self._rate_helper.nuclides = nuclides self._normalization_helper.nuclides = nuclides self._yield_helper.update_tally_nuclides(nuclides) + # Run OpenMC openmc.lib.run() openmc.lib.reset_timers() @@ -416,6 +327,10 @@ class Operator(OpenMCOperator): cell.fill = [mat.clone() for i in range(cell.num_instances)] + self.materials = openmc.Materials( + model.geometry.get_all_materials().values() + ) + def initial_condition(self): """Performs final setup and returns initial condition. @@ -513,3 +428,62 @@ class Operator(OpenMCOperator): nuclides.add(name) return nuclides + + def _load_previous_results(self): + """Load reuslts from a previous depletion simulation""" + # Reload volumes into geometry + prev_results[-1].transfer_volumes(self.model) + + # Store previous results in operator + # Distribute reaction rates according to those tracked + # on this process + if comm.size == 1: + self.prev_res = prev_results + else: + self.prev_res = Results() + mat_indexes = _distribute(range(len(self.burnable_mats))) + for res_obj in prev_results: + new_res = res_obj.distribute(self.local_mats, mat_indexes) + self.prev_res.append(new_res) + + def _get_helper_classes(self, reaction_rate_mode, reaction_rate_opts, normalization_mode, fission_yield_mode, fission_yield_opts): + """Create the ``_rate_helper``, ``normalization_helper``, and + ``_yield_helper`` attributes""" + # Get classes to assist working with tallies + if reaction_rate_mode == "direct": + self._rate_helper = DirectReactionRateHelper( + self.reaction_rates.n_nuc, self.reaction_rates.n_react) + elif reaction_rate_mode == "flux": + if reaction_rate_opts is None: + reaction_rate_opts = {} + + # Ensure energy group boundaries were specified + if 'energies' not in reaction_rate_opts: + raise ValueError( + "Energy group boundaries must be specified in the " + "reaction_rate_opts argument when reaction_rate_mode is" + "set to 'flux'.") + + self._rate_helper = FluxCollapseHelper( + self.reaction_rates.n_nuc, + self.reaction_rates.n_react, + **reaction_rate_opts + ) + else: + raise ValueError("Invalid reaction rate mode.") + + if normalization_mode == "fission-q": + self._normalization_helper = ChainFissionHelper() + elif normalization_mode == "energy-deposition": + score = "heating" if self.settings.photon_transport else "heating-local" + self._normalization_helper = EnergyScoreHelper(score) + else: + self._normalization_helper = SourceRateHelper() + + # Select and create fission yield helper + fission_helper = self._fission_helpers[fission_yield_mode] + fission_yield_opts = ( + {} if fission_yield_opts is None else fission_yield_opts) + self._yield_helper = fission_helper.from_operator( + self, **fission_yield_opts) + From 6ff164131c12db9c9b7c3d79a93095593e29c9d9 Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 14 Jul 2022 12:53:41 -0500 Subject: [PATCH 09/16] clean up docstrings, function ordering --- openmc/deplete/openmc_operator.py | 141 +++++++++---- openmc/deplete/operator.py | 332 +++++++++++++++--------------- 2 files changed, 264 insertions(+), 209 deletions(-) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index a2296d88a..0a64c6fee 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -7,7 +7,6 @@ This module implements functions used by both OpenMC transport operators as well import copy from abc import abstractmethod from collections import OrderedDict -import os import numpy as np from uncertainties import ufloat @@ -18,8 +17,6 @@ from openmc.mpi import comm from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .reaction_rates import ReactionRates -from .results import Results - __all__ = ["OpenMCOperator", "OperatorResult"] @@ -47,16 +44,11 @@ def _distribute(items): j += chunk_size class OpenMCOperator(TransportOperator): - """OpenMC transport operator for depletion. + """Abstrct class holding OpenMC-specific functions for running + depletion calculations. - Instances of this class can be used to perform depletion using OpenMC as the - transport operator. Normally, a user needn't call methods of this class - directly. Instead, an instance of this class is passed to an integrator - class, such as :class:`openmc.deplete.CECMIntegrator`. - - .. versionchanged:: 0.13.0 - The geometry and settings parameters have been replaced with a - model parameter that takes a :class:`~openmc.model.Model` object + Specific classes for running couples transport-depleton calculations or + depletion-only calculations are implemented as subclasses of OpenMCOperator Parameters ---------- @@ -68,6 +60,12 @@ class OpenMCOperator(TransportOperator): Results from a previous depletion calculation. If this argument is specified, the depletion calculation will start from the latest state in the previous results. + diff_burnable_mats : bool, optional + Whether to differentiate burnable materials with multiple instances. + Volumes are divided equally from the original material volume. + Default: False. + normalization_mode : str + Indicate how reaction rates should be normalized. fission_q : dict, optional Dictionary of nuclides and their fission Q values [eV]. If not given, values will be pulled from the ``chain_file``. Only applicable @@ -77,6 +75,30 @@ class OpenMCOperator(TransportOperator): in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. + fission_yield_mode : str + Key indicating what fission product yield scheme to use. + fission_yield_opts : dict of str to option, optional + Optional arguments to pass to the helper determined by + ``fission_yield_mode``. Will be passed directly on to the + helper. Passing a value of None will use the defaults for + the associated helper. + reaction_rate_mode : str, optional + Indicate how one-group reaction rates should be calculated. + reaction_rate_opts : dict, optional + Keyword arguments that are passed to the reaction rate helper class. + reduce_chain : bool, optional + If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the + depletion chain up to ``reduce_chain_level``. Default is False. + + .. versionadded:: 0.12 + reduce_chain_level : int, optional + Depth of the search when reducing the depletion chain. Only used + if ``reduce_chain`` evaluates to true. The default value of + ``None`` implies no limit on the depth. + + .. versionadded:: 0.12 + + Attributes ---------- @@ -112,9 +134,6 @@ class OpenMCOperator(TransportOperator): prev_res : Results or None Results from a previous depletion calculation. ``None`` if no results are to be used. - cleanup_when_done : bool - Whether to finalize and clear the shared library memory when the - depletion operation is complete. Defaults to clearing the library. """ def __init__(self, materials=None, cross_sections=None, chain_file=None, prev_results=None, diff_burnable_mats=False, normalization_mode=None, fission_q=None, dilute_initial=0.0, fission_yield_mode=None, fission_yield_opts=None, @@ -170,7 +189,11 @@ class OpenMCOperator(TransportOperator): self.reaction_rates = ReactionRates( self.local_mats, self._burnable_nucs, self.chain.reactions) - self._get_helper_classes(reaction_rate_mode, reaction_rate_opts, normalization_mode, fission_yield_mode, fission_yield_opts) + self._get_helper_classes(reaction_rate_mode, normalization_mode, fission_yield_mode, reaction_rate_opts, fission_yield_opts) + + @abstractmethod + def _differentiate_burnable_mats(self): + """Assign distribmats for each burnable material""" def _get_burnable_mats(self): """Determine depletable materials, volumes, and nuclides @@ -221,6 +244,15 @@ class OpenMCOperator(TransportOperator): return burnable_mats, volume, nuclides + + @abstractmethod + def _load_previous_results(self): + """Load reuslts from a previous depletion simulation""" + + @abstractmethod + def _get_nuclides_with_data(self, cross_sections): + """Find nuclides with cross section data""" + def _extract_number(self, local_mats, volume, all_nuclides, prev_res=None): """Construct AtomNumber using geometry @@ -304,12 +336,37 @@ class OpenMCOperator(TransportOperator): self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + @abstractmethod + def _get_helper_classes(self, reaction_rate_mode, normalization_mode, fission_yield_mode, reaction_rate_opts, fission_yield_opts): + """Create the ``_rate_helper``, ``_normalization_helper``, and + ``_yield_helper`` objects. + + Parameters + ---------- + reaction_rate_mode : str + Indicates the subclass of :class:`ReactionRateHelper` to + instantiate. + normalization_mode : str + Indicates the subclass of :class:`NormalizationHelper` to + instatiate. + fission_yield_mode : str + Indicates the subclass of :class:`FissionYieldHelper` to instatiate. + reaction_rate_opts : dict + Keyword arguments that are passed to the :class:`ReactionRateHelper` + subclass. + fission_yield_opts : dict + Keyword arguments that are passed to the :class:`FissionYieldHelper` + subclass. + + """ + def initial_condition(self, materials): """Performs final setup and returns initial condition. Parameters ---------- - materials : ??? + materials : list of str + list of material IDs Returns ------- @@ -328,6 +385,22 @@ class OpenMCOperator(TransportOperator): # Return number density vector return list(self.number.get_mat_slice(np.s_[:])) + def _update_materials_and_nuclides(self, vec): + """Update the number density, material compositions, and nuclide + lists in helper objects""" + # Update the number densities regardless of the source rate + self.number.set_density(vec) + self._update_materials() + + # Prevent OpenMC from complaining about re-creating tallies + openmc.reset_auto_ids() + + # Update tally nuclides data in preparation for transport solve + nuclides = self._get_reaction_nuclides() + self._rate_helper.nuclides = nuclides + self._normalization_helper.nuclides = nuclides + self._yield_helper.update_tally_nuclides(nuclides) + @abstractmethod def _update_materials(self): """Updates material compositions in OpenMC on all processes.""" @@ -335,16 +408,16 @@ class OpenMCOperator(TransportOperator): def _get_reaction_nuclides(self): """Determine nuclides that should be tallied for reaction rates. - This method returns a list of all nuclides that have neutron data and - are listed in the depletion chain. Technically, we should tally nuclides - that may not appear in the depletion chain because we still need to get - the fission reaction rate for these nuclides in order to normalize - power, but that is left as a future exercise. + This method returns a list of all nuclides that have cross section data + and are listed in the depletion chain. Technically, we should count + nuclides that may not appear in the depletion chain because we still + need to get the fission reaction rate for these nuclides in order to + normalize power, but that is left as a future exercise. Returns ------- list of str - Tally nuclides + Nuclides with reaction rates """ nuc_set = set() @@ -371,7 +444,7 @@ class OpenMCOperator(TransportOperator): else: nuc_list = None - # Store list of tally nuclides on each process + # Store list of nuclides on each process nuc_list = comm.bcast(nuc_list) return [nuc for nuc in nuc_list if nuc in self.chain] @@ -450,26 +523,6 @@ class OpenMCOperator(TransportOperator): return rates - - @abstractmethod - def _get_helper_classes(self, reaction_rate_mode, reaction_rate_opts, normalization_mode, fission_yield_mode, fission_yield_opts): - """Create the ``_rate_helper``, ``normalization_helper``, and - ``_yield_helper`` attributes""" - - @abstractmethod - def _load_previous_results(self): - """Load reuslts from a previous depletion simulation""" - - @abstractmethod - def _get_nuclides_with_data(self, cross_sections): - """Find nuclides with cross section data""" - - @abstractmethod - def _differentiate_burnable_mats(self): - """Assign distribmats for each burnable material - - """ - def get_results_info(self): """Returns volume list, material lists, and nuc lists. diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 2770494ba..30563623c 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -8,7 +8,6 @@ densities is all done in-memory instead of through the filesystem. """ import copy -from collections import OrderedDict import os from warnings import warn @@ -21,11 +20,9 @@ from openmc.data import DataLibrary from openmc.exceptions import DataError import openmc.lib from openmc.mpi import comm -from .abc import TransportOperator, OperatorResult -from .atom_number import AtomNumber +from .abc import OperatorResult from .chain import _find_chain_file from .openmc_operator import OpenMCOperator, _distribute -from .reaction_rates import ReactionRates from .results import Results from .helpers import ( DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, @@ -229,81 +226,8 @@ class Operator(OpenMCOperator): reaction_rate_mode, reaction_rate_opts, reduce_chain, reduce_chain_level) - def __call__(self, vec, source_rate): - """Runs a simulation. - - Simulation will abort under the following circumstances: - - 1) No energy is computed using OpenMC tallies. - - Parameters - ---------- - vec : list of numpy.ndarray - Total atoms to be used in function. - source_rate : float - Power in [W] or source rate in [neutron/sec] - - Returns - ------- - openmc.deplete.OperatorResult - Eigenvalue and reaction rates resulting from transport operator - - """ - # Reset results in OpenMC - openmc.lib.reset() - - # Update the number densities regardless of the source rate - self.number.set_density(vec) - self._update_materials() - - # If the source rate is zero, return zero reaction rates without running - # a transport solve - if source_rate == 0.0: - rates = self.reaction_rates.copy() - rates.fill(0.0) - return OperatorResult(ufloat(0.0, 0.0), rates) - - # Prevent OpenMC from complaining about re-creating tallies - openmc.reset_auto_ids() - - # Update tally nuclides data in preparation for transport solve - nuclides = self._get_reaction_nuclides() - self._rate_helper.nuclides = nuclides - self._normalization_helper.nuclides = nuclides - self._yield_helper.update_tally_nuclides(nuclides) - - - # Run OpenMC - openmc.lib.run() - openmc.lib.reset_timers() - - # Extract results - rates = self._calculate_reaction_rates(source_rate) - - # Get k and uncertainty - keff = ufloat(*openmc.lib.keff()) - - op_result = OperatorResult(keff, rates) - - return copy.deepcopy(op_result) - - @staticmethod - def write_bos_data(step): - """Write a state-point file with beginning of step data - - Parameters - ---------- - step : int - Current depletion step including restarts - """ - openmc.lib.statepoint_write( - "openmc_simulation_n{}.h5".format(step), - write_source=False) - def _differentiate_burnable_mats(self): - """Assign distribmats for each burnable material - - """ + """Assign distribmats for each burnable material""" # Count the number of instances for each cell and material self.geometry.determine_paths(instances_only=True) @@ -331,6 +255,97 @@ class Operator(OpenMCOperator): model.geometry.get_all_materials().values() ) + def _load_previous_results(self): + """Load results from a previous depletion simulation""" + # Reload volumes into geometry + prev_results[-1].transfer_volumes(self.model) + + # Store previous results in operator + # Distribute reaction rates according to those tracked + # on this process + if comm.size == 1: + self.prev_res = prev_results + else: + self.prev_res = Results() + mat_indexes = _distribute(range(len(self.burnable_mats))) + for res_obj in prev_results: + new_res = res_obj.distribute(self.local_mats, mat_indexes) + self.prev_res.append(new_res) + + def _get_nuclides_with_data(self, cross_sections): + """Loads cross_sections.xml file to find nuclides with neutron data""" + nuclides = set() + data_lib = DataLibrary.from_xml(cross_sections) + for library in data_lib.libraries: + if library['type'] != 'neutron': + continue + for name in library['materials']: + if name not in nuclides: + nuclides.add(name) + + return nuclides + + def _get_helper_classes(self, reaction_rate_mode, normalization_mode, fission_yield_mode, + reaction_rate_opts, fission_yield_opts): + """Create the ``_rate_helper``, ``_normalization_helper``, and + ``_yield_helper`` objects. + + Parameters + ---------- + reaction_rate_mode : str + Indicates the subclass of :class:`ReactionRateHelper` to + instantiate. + normalization_mode : str + Indicates the subclass of :class:`NormalizationHelper` to + instatiate. + fission_yield_mode : str + Indicates the subclass of :class:`FissionYieldHelper` to instatiate. + reaction_rate_opts : dict + Keyword arguments that are passed to the :class:`ReactionRateHelper` + subclass. + fission_yield_opts : dict + Keyword arguments that are passed to the :class:`FissionYieldHelper` + subclass. + + """ + # Get classes to assist working with tallies + if reaction_rate_mode == "direct": + self._rate_helper = DirectReactionRateHelper( + self.reaction_rates.n_nuc, self.reaction_rates.n_react) + elif reaction_rate_mode == "flux": + if reaction_rate_opts is None: + reaction_rate_opts = {} + + # Ensure energy group boundaries were specified + if 'energies' not in reaction_rate_opts: + raise ValueError( + "Energy group boundaries must be specified in the " + "reaction_rate_opts argument when reaction_rate_mode is" + "set to 'flux'.") + + self._rate_helper = FluxCollapseHelper( + self.reaction_rates.n_nuc, + self.reaction_rates.n_react, + **reaction_rate_opts + ) + else: + raise ValueError("Invalid reaction rate mode.") + + if normalization_mode == "fission-q": + self._normalization_helper = ChainFissionHelper() + elif normalization_mode == "energy-deposition": + score = "heating" if self.settings.photon_transport else "heating-local" + self._normalization_helper = EnergyScoreHelper(score) + else: + self._normalization_helper = SourceRateHelper() + + # Select and create fission yield helper + fission_helper = self._fission_helpers[fission_yield_mode] + fission_yield_opts = ( + {} if fission_yield_opts is None else fission_yield_opts) + self._yield_helper = fission_helper.from_operator( + self, **fission_yield_opts) + def initial_condition(self): """Performs final setup and returns initial condition. @@ -357,10 +372,66 @@ class Operator(OpenMCOperator): return super().initial_condition(materials) - def finalize(self): - """Finalize a depletion simulation and release resources.""" - if self.cleanup_when_done: - openmc.lib.finalize() + def _generate_materials_xml(self): + """Creates materials.xml from self.number. + + Due to uncertainty with how MPI interacts with OpenMC API, this + constructs the XML manually. The long term goal is to do this + through direct memory writing. + + """ + # Sort nuclides according to order in AtomNumber object + nuclides = list(self.number.nuclides) + for mat in self.materials: + mat._nuclides.sort(key=lambda x: nuclides.index(x[0])) + + self.materials.export_to_xml() + + def __call__(self, vec, source_rate): + """Runs a simulation. + + Simulation will abort under the following circumstances: + + 1) No energy is computed using OpenMC tallies. + + Parameters + ---------- + vec : list of numpy.ndarray + Total atoms to be used in function. + source_rate : float + Power in [W] or source rate in [neutron/sec] + + Returns + ------- + openmc.deplete.OperatorResult + Eigenvalue and reaction rates resulting from transport operator + + """ + # Reset results in OpenMC + openmc.lib.reset() + + self._update_materials_and_nuclides(vec) + + # If the source rate is zero, return zero reaction rates without running + # a transport solve + if source_rate == 0.0: + rates = self.reaction_rates.copy() + rates.fill(0.0) + return OperatorResult(ufloat(0.0, 0.0), rates) + + # Run OpenMC + openmc.lib.run() + openmc.lib.reset_timers() + + # Extract results + rates = self._calculate_reaction_rates(source_rate) + + # Get k and uncertainty + keff = ufloat(*openmc.lib.keff()) + + op_result = OperatorResult(keff, rates) + + return copy.deepcopy(op_result) def _update_materials(self): """Updates material compositions in OpenMC on all processes.""" @@ -401,89 +472,20 @@ class Operator(OpenMCOperator): #TODO Update densities on the Python side, otherwise the # summary.h5 file contains densities at the first time step - def _generate_materials_xml(self): - """Creates materials.xml from self.number. - - Due to uncertainty with how MPI interacts with OpenMC API, this - constructs the XML manually. The long term goal is to do this - through direct memory writing. + @staticmethod + def write_bos_data(step): + """Write a state-point file with beginning of step data + Parameters + ---------- + step : int + Current depletion step including restarts """ - # Sort nuclides according to order in AtomNumber object - nuclides = list(self.number.nuclides) - for mat in self.materials: - mat._nuclides.sort(key=lambda x: nuclides.index(x[0])) - - self.materials.export_to_xml() - - def _get_nuclides_with_data(self, cross_sections): - """Loads cross_sections.xml file to find nuclides with neutron data""" - nuclides = set() - data_lib = DataLibrary.from_xml(cross_sections) - for library in data_lib.libraries: - if library['type'] != 'neutron': - continue - for name in library['materials']: - if name not in nuclides: - nuclides.add(name) - - return nuclides - - def _load_previous_results(self): - """Load reuslts from a previous depletion simulation""" - # Reload volumes into geometry - prev_results[-1].transfer_volumes(self.model) - - # Store previous results in operator - # Distribute reaction rates according to those tracked - # on this process - if comm.size == 1: - self.prev_res = prev_results - else: - self.prev_res = Results() - mat_indexes = _distribute(range(len(self.burnable_mats))) - for res_obj in prev_results: - new_res = res_obj.distribute(self.local_mats, mat_indexes) - self.prev_res.append(new_res) - - def _get_helper_classes(self, reaction_rate_mode, reaction_rate_opts, normalization_mode, fission_yield_mode, fission_yield_opts): - """Create the ``_rate_helper``, ``normalization_helper``, and - ``_yield_helper`` attributes""" - # Get classes to assist working with tallies - if reaction_rate_mode == "direct": - self._rate_helper = DirectReactionRateHelper( - self.reaction_rates.n_nuc, self.reaction_rates.n_react) - elif reaction_rate_mode == "flux": - if reaction_rate_opts is None: - reaction_rate_opts = {} - - # Ensure energy group boundaries were specified - if 'energies' not in reaction_rate_opts: - raise ValueError( - "Energy group boundaries must be specified in the " - "reaction_rate_opts argument when reaction_rate_mode is" - "set to 'flux'.") - - self._rate_helper = FluxCollapseHelper( - self.reaction_rates.n_nuc, - self.reaction_rates.n_react, - **reaction_rate_opts - ) - else: - raise ValueError("Invalid reaction rate mode.") - - if normalization_mode == "fission-q": - self._normalization_helper = ChainFissionHelper() - elif normalization_mode == "energy-deposition": - score = "heating" if self.settings.photon_transport else "heating-local" - self._normalization_helper = EnergyScoreHelper(score) - else: - self._normalization_helper = SourceRateHelper() - - # Select and create fission yield helper - fission_helper = self._fission_helpers[fission_yield_mode] - fission_yield_opts = ( - {} if fission_yield_opts is None else fission_yield_opts) - self._yield_helper = fission_helper.from_operator( - self, **fission_yield_opts) + openmc.lib.statepoint_write( + "openmc_simulation_n{}.h5".format(step), + write_source=False) + def finalize(self): + """Finalize a depletion simulation and release resources.""" + if self.cleanup_when_done: + openmc.lib.finalize() From 0a05f9dcb4320ae0c7dbd88547ca2bf9b35b68f5 Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 14 Jul 2022 12:54:51 -0500 Subject: [PATCH 10/16] pep8 fixes --- openmc/deplete/openmc_operator.py | 49 ++++++++++++++++++++++------- openmc/deplete/operator.py | 52 ++++++++++++++++++++++--------- 2 files changed, 75 insertions(+), 26 deletions(-) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 0a64c6fee..2f1b93294 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -43,6 +43,7 @@ def _distribute(items): return items[j:j + chunk_size] j += chunk_size + class OpenMCOperator(TransportOperator): """Abstrct class holding OpenMC-specific functions for running depletion calculations. @@ -136,12 +137,25 @@ class OpenMCOperator(TransportOperator): results are to be used. """ - def __init__(self, materials=None, cross_sections=None, chain_file=None, prev_results=None, diff_burnable_mats=False, normalization_mode=None, fission_q=None, dilute_initial=0.0, fission_yield_mode=None, fission_yield_opts=None, - reaction_rate_mode=None, reaction_rate_opts=None, - reduce_chain=False, reduce_chain_level=None): + def __init__( + self, + materials=None, + cross_sections=None, + chain_file=None, + prev_results=None, + diff_burnable_mats=False, + normalization_mode=None, + fission_q=None, + dilute_initial=0.0, + fission_yield_mode=None, + fission_yield_opts=None, + reaction_rate_mode=None, + reaction_rate_opts=None, + reduce_chain=False, + reduce_chain_level=None): super().__init__(chain_file, fission_q, dilute_initial, prev_results) - self.round_number=False + self.round_number = False self.materials = materials self.cross_sections = cross_sections @@ -172,8 +186,8 @@ class OpenMCOperator(TransportOperator): if self.prev_res is not None: self._load_previous_results() - - self.nuclides_with_data = self._get_nuclides_with_data(self.cross_sections) + self.nuclides_with_data = self._get_nuclides_with_data( + self.cross_sections) # Select nuclides with data that are also in the chain self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides @@ -189,7 +203,12 @@ class OpenMCOperator(TransportOperator): self.reaction_rates = ReactionRates( self.local_mats, self._burnable_nucs, self.chain.reactions) - self._get_helper_classes(reaction_rate_mode, normalization_mode, fission_yield_mode, reaction_rate_opts, fission_yield_opts) + self._get_helper_classes( + reaction_rate_mode, + normalization_mode, + fission_yield_mode, + reaction_rate_opts, + fission_yield_opts) @abstractmethod def _differentiate_burnable_mats(self): @@ -244,7 +263,6 @@ class OpenMCOperator(TransportOperator): return burnable_mats, volume, nuclides - @abstractmethod def _load_previous_results(self): """Load reuslts from a previous depletion simulation""" @@ -268,11 +286,14 @@ class OpenMCOperator(TransportOperator): Results from a previous depletion calculation """ - self.number = AtomNumber(local_mats, all_nuclides, volume, len(self.chain)) + self.number = AtomNumber( + local_mats, all_nuclides, volume, len( + self.chain)) if self.dilute_initial != 0.0: for nuc in self._burnable_nucs: - self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) + self.number.set_atom_density( + np.s_[:], nuc, self.dilute_initial) # Now extract and store the number densities # From the geometry if no previous depletion results @@ -337,7 +358,13 @@ class OpenMCOperator(TransportOperator): self.number.set_atom_density(mat_id, nuclide, atom_per_cc) @abstractmethod - def _get_helper_classes(self, reaction_rate_mode, normalization_mode, fission_yield_mode, reaction_rate_opts, fission_yield_opts): + def _get_helper_classes( + self, + reaction_rate_mode, + normalization_mode, + fission_yield_mode, + reaction_rate_opts, + fission_yield_opts): """Create the ``_rate_helper``, ``_normalization_helper``, and ``_yield_helper`` objects. diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 30563623c..7149cddad 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -32,6 +32,7 @@ from .helpers import ( __all__ = ["Operator", "OperatorResult"] + def _find_cross_sections(model): """Determine cross sections to use for depletion""" if model.materials and model.materials.cross_sections is not None: @@ -219,12 +220,21 @@ class Operator(OpenMCOperator): self.cleanup_when_done = True - super().__init__(model.materials, cross_sections, chain_file, prev_results, - diff_burnable_mats, normalization_mode, - fission_q, dilute_initial, - fission_yield_mode, fission_yield_opts, - reaction_rate_mode, reaction_rate_opts, - reduce_chain, reduce_chain_level) + super().__init__( + model.materials, + cross_sections, + chain_file, + prev_results, + diff_burnable_mats, + normalization_mode, + fission_q, + dilute_initial, + fission_yield_mode, + fission_yield_opts, + reaction_rate_mode, + reaction_rate_opts, + reduce_chain, + reduce_chain_level) def _differentiate_burnable_mats(self): """Assign distribmats for each burnable material""" @@ -240,7 +250,7 @@ class Operator(OpenMCOperator): for mat in distribmats: if mat.volume is None: raise RuntimeError("Volume not specified for depletable " - "material with ID={}.".format(mat.id)) + "material with ID={}.".format(mat.id)) mat.volume /= mat.num_instances if distribmats: @@ -252,8 +262,8 @@ class Operator(OpenMCOperator): for i in range(cell.num_instances)] self.materials = openmc.Materials( - model.geometry.get_all_materials().values() - ) + model.geometry.get_all_materials().values() + ) def _load_previous_results(self): """Load results from a previous depletion simulation""" @@ -285,8 +295,13 @@ class Operator(OpenMCOperator): return nuclides - def _get_helper_classes(self, reaction_rate_mode, normalization_mode, fission_yield_mode, - reaction_rate_opts, fission_yield_opts): + def _get_helper_classes( + self, + reaction_rate_mode, + normalization_mode, + fission_yield_mode, + reaction_rate_opts, + fission_yield_opts): """Create the ``_rate_helper``, ``_normalization_helper``, and ``_yield_helper`` objects. @@ -459,17 +474,24 @@ class Operator(OpenMCOperator): densities.append(val) else: # Only output warnings if values are significantly - # negative. CRAM does not guarantee positive values. + # negative. CRAM does not guarantee positive + # values. if val < -1.0e-21: - print("WARNING: nuclide ", nuc, " in material ", mat, - " is negative (density = ", val, " at/barn-cm)") + print( + "WARNING: nuclide ", + nuc, + " in material ", + mat, + " is negative (density = ", + val, + " at/barn-cm)") number_i[mat, nuc] = 0.0 # Update densities on C API side mat_internal = openmc.lib.materials[int(mat)] mat_internal.set_densities(nuclides, densities) - #TODO Update densities on the Python side, otherwise the + # TODO Update densities on the Python side, otherwise the # summary.h5 file contains densities at the first time step @staticmethod From 39472e286c34f01e1379c457beacf9dc84cd3abc Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 14 Jul 2022 12:56:19 -0500 Subject: [PATCH 11/16] add OpenMCOperator to docs, fix referecne to TalliedFissionYieldHelper --- docs/source/pythonapi/deplete.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 9f7d8c447..11062cd87 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -204,6 +204,7 @@ prior to depleting materials :template: mycallable.rst abc.TransportOperator + openmc_operator.OpenMCOperator The following classes are abstract classes used to pass information from OpenMC simulations back on to the :class:`abc.TransportOperator` @@ -216,7 +217,7 @@ OpenMC simulations back on to the :class:`abc.TransportOperator` abc.NormalizationHelper abc.FissionYieldHelper abc.ReactionRateHelper - abc.TalliedFissionYieldHelper + helpers.TalliedFissionYieldHelper Custom integrators or depletion solvers can be developed by subclassing from the following abstract base classes: From 47344ca5cb18441e59910c734e2a7c5a88f6e7d1 Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 14 Jul 2022 14:22:33 -0500 Subject: [PATCH 12/16] spell checker typo fixes --- openmc/deplete/openmc_operator.py | 8 ++++---- openmc/deplete/operator.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 2f1b93294..6b7ff2920 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -45,10 +45,10 @@ def _distribute(items): class OpenMCOperator(TransportOperator): - """Abstrct class holding OpenMC-specific functions for running + """Abstract class holding OpenMC-specific functions for running depletion calculations. - Specific classes for running couples transport-depleton calculations or + Specific classes for running coupled transport-depletion calculations or depletion-only calculations are implemented as subclasses of OpenMCOperator Parameters @@ -265,7 +265,7 @@ class OpenMCOperator(TransportOperator): @abstractmethod def _load_previous_results(self): - """Load reuslts from a previous depletion simulation""" + """Load results from a previous depletion simulation""" @abstractmethod def _get_nuclides_with_data(self, cross_sections): @@ -375,7 +375,7 @@ class OpenMCOperator(TransportOperator): instantiate. normalization_mode : str Indicates the subclass of :class:`NormalizationHelper` to - instatiate. + instantiate. fission_yield_mode : str Indicates the subclass of :class:`FissionYieldHelper` to instatiate. reaction_rate_opts : dict diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 7149cddad..a86eb182f 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -314,7 +314,7 @@ class Operator(OpenMCOperator): Indicates the subclass of :class:`NormalizationHelper` to instatiate. fission_yield_mode : str - Indicates the subclass of :class:`FissionYieldHelper` to instatiate. + Indicates the subclass of :class:`FissionYieldHelper` to instantiate. reaction_rate_opts : dict Keyword arguments that are passed to the :class:`ReactionRateHelper` subclass. From 7b041114209c6073793ad8c33d87db7e07afda6b Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 14 Jul 2022 14:23:11 -0500 Subject: [PATCH 13/16] add small blurb about OpenMCOperator; move TalliedFissionYieldHelper to different section --- docs/source/pythonapi/deplete.rst | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 11062cd87..d09a665d8 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -188,6 +188,7 @@ total system energy. helpers.EnergyScoreHelper helpers.FissionYieldCutoffHelper helpers.FluxCollapseHelper + helpers.TalliedFissionYieldHelper Abstract Base Classes --------------------- @@ -204,7 +205,16 @@ prior to depleting materials :template: mycallable.rst abc.TransportOperator - openmc_operator.OpenMCOperator + +Methods common to OpenMC-specific implementations are stored in :class:`openmc_operator.OpenMCOperator` + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: mycallable.rst + + abc.TransportOperator + The following classes are abstract classes used to pass information from OpenMC simulations back on to the :class:`abc.TransportOperator` @@ -217,7 +227,6 @@ OpenMC simulations back on to the :class:`abc.TransportOperator` abc.NormalizationHelper abc.FissionYieldHelper abc.ReactionRateHelper - helpers.TalliedFissionYieldHelper Custom integrators or depletion solvers can be developed by subclassing from the following abstract base classes: From 80fbcdcdf380e583c346d140f0fa8724b526cfa9 Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 14 Jul 2022 15:14:35 -0500 Subject: [PATCH 14/16] move intermediate classes to their own section in the api docs --- docs/source/pythonapi/deplete.rst | 38 +++++++++++++++++++++++-------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index d09a665d8..e275f38f2 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -188,8 +188,36 @@ total system energy. helpers.EnergyScoreHelper helpers.FissionYieldCutoffHelper helpers.FluxCollapseHelper + + +Intermediate Classes +-------------------- + +Specific implementations of abstract base classes may utilize some of +the same methods and data structures. These methods and data are stored +in intermediate classes. + +Methods common to tally-based implementation of :class:`FissionYieldHelper` +are stored in :class:`helpers.TalliedFissionYieldHelper` + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + helpers.TalliedFissionYieldHelper +Methods common to OpenMC-specific implementations of :class:`TransportOperator` +are stored in :class:`openmc_operator.OpenMCOperator` + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: mycallable.rst + + openmc_operator.OpenMCOperator + + Abstract Base Classes --------------------- @@ -206,16 +234,6 @@ prior to depleting materials abc.TransportOperator -Methods common to OpenMC-specific implementations are stored in :class:`openmc_operator.OpenMCOperator` - -.. autosummary:: - :toctree: generated - :nosignatures: - :template: mycallable.rst - - abc.TransportOperator - - The following classes are abstract classes used to pass information from OpenMC simulations back on to the :class:`abc.TransportOperator` From 8bf4465b27ac32b2e1a7d152e93276620541d39d Mon Sep 17 00:00:00 2001 From: yardasol Date: Thu, 14 Jul 2022 16:50:05 -0500 Subject: [PATCH 15/16] docstring updates --- openmc/deplete/openmc_operator.py | 49 +++++++++++++++++++++---------- openmc/deplete/operator.py | 25 ++++++++++++++-- 2 files changed, 56 insertions(+), 18 deletions(-) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 6b7ff2920..047a7f967 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -53,6 +53,11 @@ class OpenMCOperator(TransportOperator): Parameters ---------- + materials : openmc.Materials + List of all materials in the model + cross_sections : str or pandas.DataFrame + Path to continuous energy cross section library, or object containing + one-group cross-sections. chain_file : str, optional Path to the depletion chain XML file. Defaults to the file listed under ``depletion_chain`` in @@ -68,9 +73,7 @@ class OpenMCOperator(TransportOperator): normalization_mode : str Indicate how reaction rates should be normalized. fission_q : dict, optional - Dictionary of nuclides and their fission Q values [eV]. If not given, - values will be pulled from the ``chain_file``. Only applicable - if ``"normalization_mode" == "fission-q"`` + Dictionary of nuclides and their fission Q values [eV]. dilute_initial : float, optional Initial atom density [atoms/cm^3] to add for nuclides that are zero in initial condition to ensure they exist in the decay chain. @@ -90,25 +93,16 @@ class OpenMCOperator(TransportOperator): reduce_chain : bool, optional If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the depletion chain up to ``reduce_chain_level``. Default is False. - - .. versionadded:: 0.12 reduce_chain_level : int, optional Depth of the search when reducing the depletion chain. Only used if ``reduce_chain`` evaluates to true. The default value of ``None`` implies no limit on the depth. - .. versionadded:: 0.12 - - Attributes ---------- - model : openmc.model.Model - OpenMC model object - geometry : openmc.Geometry - OpenMC geometry object - settings : openmc.Settings - OpenMC settings object + materials : openmc.Materials + All materials present in the model dilute_initial : float Initial atom density [atoms/cm^3] to add for nuclides that are zero in initial condition to ensure they exist in the decay @@ -135,6 +129,7 @@ class OpenMCOperator(TransportOperator): prev_res : Results or None Results from a previous depletion calculation. ``None`` if no results are to be used. + """ def __init__( @@ -269,7 +264,20 @@ class OpenMCOperator(TransportOperator): @abstractmethod def _get_nuclides_with_data(self, cross_sections): - """Find nuclides with cross section data""" + """Find nuclides with cross section data + + Parameters + ---------- + cross_sections : str or pandas.DataFrame + Path to continuous energy cross section library, or object + containing one-group cross-sections. + + Returns + ------- + nuclides : set of str + Set of nuclide names that have cross secton data + + """ def _extract_number(self, local_mats, volume, all_nuclides, prev_res=None): """Construct AtomNumber using geometry @@ -399,6 +407,7 @@ class OpenMCOperator(TransportOperator): ------- list of numpy.ndarray Total density for initial conditions. + """ self._rate_helper.generate_tallies(materials, self.chain.reactions) @@ -414,7 +423,15 @@ class OpenMCOperator(TransportOperator): def _update_materials_and_nuclides(self, vec): """Update the number density, material compositions, and nuclide - lists in helper objects""" + lists in helper objects + + Parameters + ---------- + vec : list of numpy.ndarray + Total atoms. + + """ + # Update the number densities regardless of the source rate self.number.set_density(vec) self._update_materials() diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index a86eb182f..d8b84f5a0 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -34,7 +34,14 @@ __all__ = ["Operator", "OperatorResult"] def _find_cross_sections(model): - """Determine cross sections to use for depletion""" + """Determine cross sections to use for depletion + + Parameters + ---------- + model : openmc.model.Model + Reactor model + + """ if model.materials and model.materials.cross_sections is not None: # Prefer info from Model class if available return model.materials.cross_sections @@ -283,7 +290,19 @@ class Operator(OpenMCOperator): self.prev_res.append(new_res) def _get_nuclides_with_data(self, cross_sections): - """Loads cross_sections.xml file to find nuclides with neutron data""" + """Loads cross_sections.xml file to find nuclides with neutron data + + Parameters + ---------- + cross_sections : str + Path to cross_sections.xml file + + Returns + ------- + nuclides : set of str + Set of nuclide names that have cross secton data + + """ nuclides = set() data_lib = DataLibrary.from_xml(cross_sections) for library in data_lib.libraries: @@ -368,6 +387,7 @@ class Operator(OpenMCOperator): ------- list of numpy.ndarray Total density for initial conditions. + """ # Create XML files @@ -502,6 +522,7 @@ class Operator(OpenMCOperator): ---------- step : int Current depletion step including restarts + """ openmc.lib.statepoint_write( "openmc_simulation_n{}.h5".format(step), From efa92ad865d7abee86bd21524164433d96957185 Mon Sep 17 00:00:00 2001 From: yardasol Date: Tue, 19 Jul 2022 10:05:47 -0500 Subject: [PATCH 16/16] Address pauromano's comments --- openmc/deplete/openmc_operator.py | 3 --- openmc/deplete/operator.py | 9 ++++----- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 047a7f967..e057b3467 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -4,15 +4,12 @@ This module implements functions used by both OpenMC transport operators as well """ -import copy from abc import abstractmethod from collections import OrderedDict import numpy as np -from uncertainties import ufloat import openmc -from openmc.exceptions import DataError from openmc.mpi import comm from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index d8b84f5a0..b3c0756b8 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -269,20 +269,19 @@ class Operator(OpenMCOperator): for i in range(cell.num_instances)] self.materials = openmc.Materials( - model.geometry.get_all_materials().values() + self.model.geometry.get_all_materials().values() ) def _load_previous_results(self): """Load results from a previous depletion simulation""" # Reload volumes into geometry - prev_results[-1].transfer_volumes(self.model) + self.prev_res[-1].transfer_volumes(self.model) # Store previous results in operator # Distribute reaction rates according to those tracked # on this process - if comm.size == 1: - self.prev_res = prev_results - else: + if comm.size != 1: + prev_results = self.prev_res self.prev_res = Results() mat_indexes = _distribute(range(len(self.burnable_mats))) for res_obj in prev_results: