diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 32b449a13..1f4321642 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -190,6 +190,35 @@ total system energy. helpers.FissionYieldCutoffHelper helpers.FluxCollapseHelper + +Intermediate Classes +-------------------- + +Specific implementations of abstract base classes may utilize some of +the same methods and data structures. These methods and data are stored +in intermediate classes. + +Methods common to tally-based implementation of :class:`FissionYieldHelper` +are stored in :class:`helpers.TalliedFissionYieldHelper` + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + helpers.TalliedFissionYieldHelper + +Methods common to OpenMC-specific implementations of :class:`TransportOperator` +are stored in :class:`openmc_operator.OpenMCOperator` + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: mycallable.rst + + openmc_operator.OpenMCOperator + + Abstract Base Classes --------------------- @@ -217,7 +246,6 @@ OpenMC simulations back on to the :class:`abc.TransportOperator` abc.NormalizationHelper abc.FissionYieldHelper abc.ReactionRateHelper - abc.TalliedFissionYieldHelper Custom integrators or depletion solvers can be developed by subclassing from the following abstract base classes: diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 753167741..deddf1100 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 .flux_operator import * from .reaction_rates import * diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index a4f80e3d6..8bbe126df 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,24 @@ except AttributeError: # Can't set __doc__ on properties on Python 3.4 pass +@contextmanager +def change_directory(output_dir): + """ + Helper function for managing the current directory. + + Parameters + ---------- + output_dir : pathlib.Path + Directory to switch to. + """ + orig_dir = os.getcwd() + try: + output_dir.mkdir(exist_ok=True) + os.chdir(output_dir) + yield + finally: + os.chdir(orig_dir) + class TransportOperator(ABC): """Abstract class defining a transport operator @@ -87,9 +106,14 @@ class TransportOperator(ABC): Initial atom density [atoms/cm^3] to add for nuclides that are zero in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. + output_dir : pathlib.Path + Path to output directory to save results. prev_res : Results or None Results from a previous depletion calculation. ``None`` if no results are to be used. + chain : openmc.deplete.Chain + The depletion chain information necessary to form matrices and tallies. + """ def __init__(self, chain_file, fission_q=None, dilute_initial=1.0e3, prev_results=None): @@ -133,18 +157,6 @@ class TransportOperator(ABC): """ - def __enter__(self): - # Save current directory and move to specific output directory - self._orig_dir = os.getcwd() - self.output_dir.mkdir(exist_ok=True) - os.chdir(self.output_dir) - - return self.initial_condition() - - def __exit__(self, exc_type, exc_value, traceback): - self.finalize() - os.chdir(self._orig_dir) - @property def output_dir(self): return self._output_dir @@ -775,7 +787,8 @@ class Integrator(ABC): .. versionadded:: 0.13.1 """ - with self.operator as conc: + with change_directory(self.operator.output_dir): + conc = self.operator.initial_condition() t, self._i_res = self._get_start_data() for i, (dt, source_rate) in enumerate(self): @@ -814,6 +827,8 @@ class Integrator(ABC): source_rate, self._i_res + len(self), proc_time) self.operator.write_bos_data(len(self) + self._i_res) + self.operator.finalize() + @add_params class SIIntegrator(Integrator): @@ -931,7 +946,8 @@ class SIIntegrator(Integrator): .. versionadded:: 0.13.1 """ - with self.operator as conc: + with change_directory(self.operator.output_dir): + conc = self.operator.initial_condition() t, self._i_res = self._get_start_data() for i, (dt, p) in enumerate(self): @@ -967,6 +983,8 @@ class SIIntegrator(Integrator): p, self._i_res + len(self), proc_time) self.operator.write_bos_data(self._i_res + len(self)) + self.operator.finalize() + class DepSystemSolver(ABC): r"""Abstract class for solving depletion equations diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py new file mode 100644 index 000000000..e057b3467 --- /dev/null +++ b/openmc/deplete/openmc_operator.py @@ -0,0 +1,593 @@ +"""OpenMC transport operator + +This module implements functions used by both OpenMC transport operators as well as pure depletion operators. + +""" + +from abc import abstractmethod +from collections import OrderedDict + +import numpy as np + +import openmc +from openmc.mpi import comm +from .abc import TransportOperator, OperatorResult +from .atom_number import AtomNumber +from .reaction_rates import ReactionRates + +__all__ = ["OpenMCOperator", "OperatorResult"] + + +def _distribute(items): + """Distribute items across MPI communicator + + Parameters + ---------- + items : list + List of items of distribute + + Returns + ------- + list + Items assigned to process that called + + """ + min_size, extra = divmod(len(items), comm.size) + j = 0 + for i in range(comm.size): + chunk_size = min_size + int(i < extra) + if comm.rank == i: + return items[j:j + chunk_size] + j += chunk_size + + +class OpenMCOperator(TransportOperator): + """Abstract class holding OpenMC-specific functions for running + depletion calculations. + + Specific classes for running coupled transport-depletion calculations or + depletion-only calculations are implemented as subclasses of OpenMCOperator + + Parameters + ---------- + materials : openmc.Materials + List of all materials in the model + cross_sections : str or pandas.DataFrame + Path to continuous energy cross section library, or object containing + one-group cross-sections. + chain_file : str, optional + Path to the depletion chain XML file. Defaults to the file + listed under ``depletion_chain`` in + :envvar:`OPENMC_CROSS_SECTIONS` environment variable. + prev_results : Results, optional + Results from a previous depletion calculation. If this argument is + specified, the depletion calculation will start from the latest state + in the previous results. + diff_burnable_mats : bool, optional + Whether to differentiate burnable materials with multiple instances. + Volumes are divided equally from the original material volume. + Default: False. + normalization_mode : str + Indicate how reaction rates should be normalized. + fission_q : dict, optional + Dictionary of nuclides and their fission Q values [eV]. + dilute_initial : float, optional + Initial atom density [atoms/cm^3] to add for nuclides that are zero + in initial condition to ensure they exist in the decay chain. + Only done for nuclides with reaction rates. + Defaults to 1.0e3. + fission_yield_mode : str + Key indicating what fission product yield scheme to use. + fission_yield_opts : dict of str to option, optional + Optional arguments to pass to the helper determined by + ``fission_yield_mode``. Will be passed directly on to the + helper. Passing a value of None will use the defaults for + the associated helper. + reaction_rate_mode : str, optional + Indicate how one-group reaction rates should be calculated. + reaction_rate_opts : dict, optional + Keyword arguments that are passed to the reaction rate helper class. + reduce_chain : bool, optional + If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the + depletion chain up to ``reduce_chain_level``. Default is False. + reduce_chain_level : int, optional + Depth of the search when reducing the depletion chain. Only used + if ``reduce_chain`` evaluates to true. The default value of + ``None`` implies no limit on the depth. + + + Attributes + ---------- + materials : openmc.Materials + All materials present in the model + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the decay + chain. Only done for nuclides with reaction rates. + output_dir : pathlib.Path + Path to output directory to save results. + round_number : bool + Whether or not to round output to OpenMC to 8 digits. + Useful in testing, as OpenMC is incredibly sensitive to exact values. + number : openmc.deplete.AtomNumber + Total number of atoms in simulation. + nuclides_with_data : set of str + A set listing all unique nuclides available from cross_sections.xml. + chain : openmc.deplete.Chain + The depletion chain information necessary to form matrices and tallies. + reaction_rates : openmc.deplete.ReactionRates + Reaction rates from the last operator step. + burnable_mats : list of str + All burnable material IDs + heavy_metal : float + Initial heavy metal inventory [g] + local_mats : list of str + All burnable material IDs being managed by a single process + prev_res : Results or None + Results from a previous depletion calculation. ``None`` if no + results are to be used. + + """ + + def __init__( + self, + materials=None, + cross_sections=None, + chain_file=None, + prev_results=None, + diff_burnable_mats=False, + normalization_mode=None, + fission_q=None, + dilute_initial=0.0, + fission_yield_mode=None, + fission_yield_opts=None, + reaction_rate_mode=None, + reaction_rate_opts=None, + reduce_chain=False, + reduce_chain_level=None): + + super().__init__(chain_file, fission_q, dilute_initial, prev_results) + self.round_number = False + self.materials = materials + self.cross_sections = cross_sections + + # Reduce the chain to only those nuclides present + if reduce_chain: + init_nuclides = set() + for material in self.materials: + if not material.depletable: + continue + for name, _dens_percent, _dens_type in material.nuclides: + init_nuclides.add(name) + + self.chain = self.chain.reduce(init_nuclides, reduce_chain_level) + + if diff_burnable_mats: + self._differentiate_burnable_mats() + + # Determine which nuclides have cross section data + # This nuclides variables contains every nuclides + # for which there is an entry in the micro_xs parameter + openmc.reset_auto_ids() + self.burnable_mats, volumes, all_nuclides = self._get_burnable_mats() + self.local_mats = _distribute(self.burnable_mats) + + self._mat_index_map = { + lm: self.burnable_mats.index(lm) for lm in self.local_mats} + + if self.prev_res is not None: + self._load_previous_results() + + self.nuclides_with_data = self._get_nuclides_with_data( + self.cross_sections) + + # Select nuclides with data that are also in the chain + self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides + if nuc.name in self.nuclides_with_data] + + # Extract number densities from the geometry / previous depletion run + self._extract_number(self.local_mats, + volumes, + all_nuclides, + self.prev_res) + + # Create reaction rates array + self.reaction_rates = ReactionRates( + self.local_mats, self._burnable_nucs, self.chain.reactions) + + self._get_helper_classes( + reaction_rate_mode, + normalization_mode, + fission_yield_mode, + reaction_rate_opts, + fission_yield_opts) + + @abstractmethod + def _differentiate_burnable_mats(self): + """Assign distribmats for each burnable material""" + + def _get_burnable_mats(self): + """Determine depletable materials, volumes, and nuclides + + Returns + ------- + burnable_mats : list of str + List of burnable material IDs + volume : OrderedDict of str to float + Volume of each material in [cm^3] + nuclides : list of str + Nuclides in order of how they'll appear in the simulation. + + """ + + burnable_mats = set() + model_nuclides = set() + volume = OrderedDict() + + self.heavy_metal = 0.0 + + # Iterate once through the geometry to get dictionaries + for mat in self.materials: + for nuclide in mat.get_nuclides(): + model_nuclides.add(nuclide) + if mat.depletable: + burnable_mats.add(str(mat.id)) + if mat.volume is None: + raise RuntimeError("Volume not specified for depletable " + "material with ID={}.".format(mat.id)) + volume[str(mat.id)] = mat.volume + self.heavy_metal += mat.fissionable_mass + + # Make sure there are burnable materials + if not burnable_mats: + raise RuntimeError( + "No depletable materials were found in the model.") + + # Sort the sets + burnable_mats = sorted(burnable_mats, key=int) + model_nuclides = sorted(model_nuclides) + + # Construct a global nuclide dictionary, burned first + nuclides = list(self.chain.nuclide_dict) + for nuc in model_nuclides: + if nuc not in nuclides: + nuclides.append(nuc) + + return burnable_mats, volume, nuclides + + @abstractmethod + def _load_previous_results(self): + """Load results from a previous depletion simulation""" + + @abstractmethod + def _get_nuclides_with_data(self, cross_sections): + """Find nuclides with cross section data + + Parameters + ---------- + cross_sections : str or pandas.DataFrame + Path to continuous energy cross section library, or object + containing one-group cross-sections. + + Returns + ------- + nuclides : set of str + Set of nuclide names that have cross secton data + + """ + + def _extract_number(self, local_mats, volume, all_nuclides, prev_res=None): + """Construct AtomNumber using geometry + + Parameters + ---------- + local_mats : list of str + Material IDs to be managed by this process + volume : OrderedDict of str to float + Volumes for the above materials in [cm^3] + all_nuclides : list of str + Nuclides to be used in the simulation. + prev_res : Results, optional + Results from a previous depletion calculation + + """ + self.number = AtomNumber( + local_mats, all_nuclides, volume, len( + self.chain)) + + if self.dilute_initial != 0.0: + for nuc in self._burnable_nucs: + self.number.set_atom_density( + np.s_[:], nuc, self.dilute_initial) + + # Now extract and store the number densities + # From the geometry if no previous depletion results + if prev_res is None: + for mat in self.materials: + if str(mat.id) in local_mats: + self._set_number_from_mat(mat) + + # Else from previous depletion results + else: + for mat in self.materials: + if str(mat.id) in local_mats: + self._set_number_from_results(mat, prev_res) + + def _set_number_from_mat(self, mat): + """Extracts material and number densities from openmc.Material + + Parameters + ---------- + mat : openmc.Material + The material to read from + + """ + mat_id = str(mat.id) + + for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): + atom_per_cc = atom_per_bcm * 1.0e24 + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + + def _set_number_from_results(self, mat, prev_res): + """Extracts material nuclides and number densities. + + If the nuclide concentration's evolution is tracked, the densities come + from depletion results. Else, densities are extracted from the geometry + in the summary. + + Parameters + ---------- + mat : openmc.Material + The material to read from + prev_res : Results + Results from a previous depletion calculation + + """ + mat_id = str(mat.id) + + # Get nuclide lists from geometry and depletion results + depl_nuc = prev_res[-1].nuc_to_ind + geom_nuc_densities = mat.get_nuclide_atom_densities() + + # Merge lists of nuclides, with the same order for every calculation + geom_nuc_densities.update(depl_nuc) + + for nuclide, atom_per_bcm in geom_nuc_densities.items(): + if nuclide in depl_nuc: + concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] + volume = prev_res[-1].volume[mat_id] + atom_per_cc = concentration / volume + else: + atom_per_cc = atom_per_bcm * 1.0e24 + + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + + @abstractmethod + def _get_helper_classes( + self, + reaction_rate_mode, + normalization_mode, + fission_yield_mode, + reaction_rate_opts, + fission_yield_opts): + """Create the ``_rate_helper``, ``_normalization_helper``, and + ``_yield_helper`` objects. + + Parameters + ---------- + reaction_rate_mode : str + Indicates the subclass of :class:`ReactionRateHelper` to + instantiate. + normalization_mode : str + Indicates the subclass of :class:`NormalizationHelper` to + instantiate. + fission_yield_mode : str + Indicates the subclass of :class:`FissionYieldHelper` to instatiate. + reaction_rate_opts : dict + Keyword arguments that are passed to the :class:`ReactionRateHelper` + subclass. + fission_yield_opts : dict + Keyword arguments that are passed to the :class:`FissionYieldHelper` + subclass. + + """ + + def initial_condition(self, materials): + """Performs final setup and returns initial condition. + + Parameters + ---------- + materials : list of str + list of material IDs + + Returns + ------- + list of numpy.ndarray + Total density for initial conditions. + + """ + + self._rate_helper.generate_tallies(materials, self.chain.reactions) + self._normalization_helper.prepare( + self.chain.nuclides, self.reaction_rates.index_nuc) + # Tell fission yield helper what materials this process is + # responsible for + self._yield_helper.generate_tallies( + materials, tuple(sorted(self._mat_index_map.values()))) + + # Return number density vector + return list(self.number.get_mat_slice(np.s_[:])) + + def _update_materials_and_nuclides(self, vec): + """Update the number density, material compositions, and nuclide + lists in helper objects + + Parameters + ---------- + vec : list of numpy.ndarray + Total atoms. + + """ + + # Update the number densities regardless of the source rate + self.number.set_density(vec) + self._update_materials() + + # Prevent OpenMC from complaining about re-creating tallies + openmc.reset_auto_ids() + + # Update tally nuclides data in preparation for transport solve + nuclides = self._get_reaction_nuclides() + self._rate_helper.nuclides = nuclides + self._normalization_helper.nuclides = nuclides + self._yield_helper.update_tally_nuclides(nuclides) + + @abstractmethod + def _update_materials(self): + """Updates material compositions in OpenMC on all processes.""" + + def _get_reaction_nuclides(self): + """Determine nuclides that should be tallied for reaction rates. + + This method returns a list of all nuclides that have cross section data + and are listed in the depletion chain. Technically, we should count + nuclides that may not appear in the depletion chain because we still + need to get the fission reaction rate for these nuclides in order to + normalize power, but that is left as a future exercise. + + Returns + ------- + list of str + Nuclides with reaction rates + + """ + nuc_set = set() + + # Create the set of all nuclides in the decay chain in materials marked + # for burning in which the number density is greater than zero. + for nuc in self.number.nuclides: + if nuc in self.nuclides_with_data: + if np.sum(self.number[:, nuc]) > 0.0: + nuc_set.add(nuc) + + # Communicate which nuclides have nonzeros to rank 0 + if comm.rank == 0: + for i in range(1, comm.size): + nuc_newset = comm.recv(source=i, tag=i) + nuc_set |= nuc_newset + else: + comm.send(nuc_set, dest=0, tag=comm.rank) + + if comm.rank == 0: + # Sort nuclides in the same order as self.number + nuc_list = [nuc for nuc in self.number.nuclides + if nuc in nuc_set] + else: + nuc_list = None + + # Store list of nuclides on each process + nuc_list = comm.bcast(nuc_list) + return [nuc for nuc in nuc_list if nuc in self.chain] + + def _calculate_reaction_rates(self, source_rate): + """Unpack tallies from OpenMC and return an operator result + + This method uses OpenMC's C API bindings to determine the k-effective + value and reaction rates from the simulation. The reaction rates are + normalized by a helper class depending on the method being used. + + Parameters + ---------- + source_rate : float + Power in [W] or source rate in [neutron/sec] + + Returns + ------- + rates : openmc.deplete.ReactionRates + Reaction rates for nuclides + + """ + rates = self.reaction_rates + rates.fill(0.0) + + # Extract reaction nuclides + rxn_nuclides = self._rate_helper.nuclides + + # Form fast map + nuc_ind = [rates.index_nuc[nuc] for nuc in rxn_nuclides] + react_ind = [rates.index_rx[react] for react in self.chain.reactions] + + # Keep track of energy produced from all reactions in eV per source + # particle + self._normalization_helper.reset() + self._yield_helper.unpack() + + # Store fission yield dictionaries + fission_yields = [] + + # Create arrays to store fission Q values, reaction rates, and nuclide + # numbers, zeroed out in material iteration + number = np.empty(rates.n_nuc) + + fission_ind = rates.index_rx.get("fission") + + # Extract results + for i, mat in enumerate(self.local_mats): + # Get tally index + mat_index = self._mat_index_map[mat] + + # Zero out reaction rates and nuclide numbers + number.fill(0.0) + + # Get new number densities + for nuc, i_nuc_results in zip(rxn_nuclides, nuc_ind): + number[i_nuc_results] = self.number[mat, nuc] + + tally_rates = self._rate_helper.get_material_rates( + mat_index, nuc_ind, react_ind) + + # Compute fission yields for this material + fission_yields.append(self._yield_helper.weighted_yields(i)) + + # Accumulate energy from fission + if fission_ind is not None: + self._normalization_helper.update(tally_rates[:, fission_ind]) + + # Divide by total number and store + rates[i] = self._rate_helper.divide_by_adens(number) + + # Scale reaction rates to obtain units of reactions/sec + rates *= self._normalization_helper.factor(source_rate) + + # Store new fission yields on the chain + self.chain.fission_yields = fission_yields + + return rates + + def get_results_info(self): + """Returns volume list, material lists, and nuc lists. + + Returns + ------- + volume : dict of str float + Volumes corresponding to materials in full_burn_dict + nuc_list : list of str + A list of all nuclide names. Used for sorting the simulation. + burn_list : list of int + A list of all material IDs to be burned. Used for sorting the simulation. + full_burn_list : list + List of all burnable material IDs + + """ + nuc_list = self.number.burnable_nuclides + burn_list = self.local_mats + + volume = {} + for i, mat in enumerate(burn_list): + volume[mat] = self.number.volume[i] + + # Combine volume dictionaries across processes + volume_list = comm.allgather(volume) + volume = {k: v for d in volume_list for k, v in d.items()} + + return volume, nuc_list, burn_list, self.burnable_mats diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 7ef874e87..b3c0756b8 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,10 +20,9 @@ from openmc.data import DataLibrary from openmc.exceptions import DataError import openmc.lib from openmc.mpi import comm -from .abc import TransportOperator, OperatorResult -from .atom_number import AtomNumber +from .abc import OperatorResult from .chain import _find_chain_file -from .reaction_rates import ReactionRates +from .openmc_operator import OpenMCOperator, _distribute from .results import Results from .helpers import ( DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, @@ -35,31 +33,15 @@ from .helpers import ( __all__ = ["Operator", "OperatorResult"] -def _distribute(items): - """Distribute items across MPI communicator +def _find_cross_sections(model): + """Determine cross sections to use for depletion Parameters ---------- - items : list - List of items of distribute - - Returns - ------- - list - Items assigned to process that called + model : openmc.model.Model + Reactor model """ - min_size, extra = divmod(len(items), comm.size) - j = 0 - for i in range(comm.size): - chunk_size = min_size + int(i < extra) - if comm.rank == i: - return items[j:j + chunk_size] - j += chunk_size - - -def _find_cross_sections(model): - """Determine cross sections to use for depletion""" if model.materials and model.materials.cross_sections is not None: # Prefer info from Model class if available return model.materials.cross_sections @@ -74,7 +56,7 @@ def _find_cross_sections(model): return cross_sections -class Operator(TransportOperator): +class Operator(OpenMCOperator): """OpenMC transport operator for depletion. Instances of this class can be used to perform depletion using OpenMC as the @@ -233,8 +215,6 @@ class Operator(TransportOperator): if fission_q is not None: warn("Fission Q dictionary will not be used") fission_q = None - super().__init__(chain_file, fission_q, dilute_initial, prev_results) - self.round_number = False self.model = model self.settings = model.settings self.geometry = model.geometry @@ -244,66 +224,123 @@ class Operator(TransportOperator): model.materials = openmc.Materials( model.geometry.get_all_materials().values() ) - self.materials = model.materials self.cleanup_when_done = True - # Reduce the chain before we create more materials - if reduce_chain: - all_isotopes = set() - for material in self.materials: - if not material.depletable: - continue - for name, _dens_percent, _dens_type in material.nuclides: - all_isotopes.add(name) - self.chain = self.chain.reduce(all_isotopes, reduce_chain_level) + super().__init__( + model.materials, + cross_sections, + chain_file, + prev_results, + diff_burnable_mats, + normalization_mode, + fission_q, + dilute_initial, + fission_yield_mode, + fission_yield_opts, + reaction_rate_mode, + reaction_rate_opts, + reduce_chain, + reduce_chain_level) - # Differentiate burnable materials with multiple instances - if diff_burnable_mats: - self._differentiate_burnable_mats() - self.materials = openmc.Materials( - model.geometry.get_all_materials().values() - ) + def _differentiate_burnable_mats(self): + """Assign distribmats for each burnable material""" - # Clear out OpenMC, create task lists, distribute - openmc.reset_auto_ids() - self.burnable_mats, volume, nuclides = self._get_burnable_mats() - self.local_mats = _distribute(self.burnable_mats) + # Count the number of instances for each cell and material + self.geometry.determine_paths(instances_only=True) - # Generate map from local materials => material index - self._mat_index_map = { - lm: self.burnable_mats.index(lm) for lm in self.local_mats} + # Extract all burnable materials which have multiple instances + distribmats = set( + [mat for mat in self.materials + if mat.depletable and mat.num_instances > 1]) - if self.prev_res is not None: - # Reload volumes into geometry - prev_results[-1].transfer_volumes(self.model) + for mat in distribmats: + if mat.volume is None: + raise RuntimeError("Volume not specified for depletable " + "material with ID={}.".format(mat.id)) + mat.volume /= mat.num_instances - # Store previous results in operator - # Distribute reaction rates according to those tracked - # on this process - if comm.size == 1: - self.prev_res = prev_results - else: - self.prev_res = Results() - mat_indexes = _distribute(range(len(self.burnable_mats))) - for res_obj in prev_results: - new_res = res_obj.distribute(self.local_mats, mat_indexes) - self.prev_res.append(new_res) + if distribmats: + # Assign distribmats to cells + for cell in self.geometry.get_all_material_cells().values(): + if cell.fill in distribmats: + mat = cell.fill + cell.fill = [mat.clone() + for i in range(cell.num_instances)] - # Determine which nuclides have incident neutron data - self.nuclides_with_data = self._get_nuclides_with_data(cross_sections) + self.materials = openmc.Materials( + self.model.geometry.get_all_materials().values() + ) - # Select nuclides with data that are also in the chain - self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides - if nuc.name in self.nuclides_with_data] + def _load_previous_results(self): + """Load results from a previous depletion simulation""" + # Reload volumes into geometry + self.prev_res[-1].transfer_volumes(self.model) - # Extract number densities from the geometry / previous depletion run - self._extract_number(self.local_mats, volume, nuclides, self.prev_res) + # Store previous results in operator + # Distribute reaction rates according to those tracked + # on this process + if comm.size != 1: + prev_results = self.prev_res + self.prev_res = Results() + mat_indexes = _distribute(range(len(self.burnable_mats))) + for res_obj in prev_results: + new_res = res_obj.distribute(self.local_mats, mat_indexes) + self.prev_res.append(new_res) - # Create reaction rates array - self.reaction_rates = ReactionRates( - self.local_mats, self._burnable_nucs, self.chain.reactions) + def _get_nuclides_with_data(self, cross_sections): + """Loads cross_sections.xml file to find nuclides with neutron data + Parameters + ---------- + cross_sections : str + Path to cross_sections.xml file + + Returns + ------- + nuclides : set of str + Set of nuclide names that have cross secton data + + """ + nuclides = set() + data_lib = DataLibrary.from_xml(cross_sections) + for library in data_lib.libraries: + if library['type'] != 'neutron': + continue + for name in library['materials']: + if name not in nuclides: + nuclides.add(name) + + return nuclides + + def _get_helper_classes( + self, + reaction_rate_mode, + normalization_mode, + fission_yield_mode, + reaction_rate_opts, + fission_yield_opts): + """Create the ``_rate_helper``, ``_normalization_helper``, and + ``_yield_helper`` objects. + + Parameters + ---------- + reaction_rate_mode : str + Indicates the subclass of :class:`ReactionRateHelper` to + instantiate. + normalization_mode : str + Indicates the subclass of :class:`NormalizationHelper` to + instatiate. + fission_yield_mode : str + Indicates the subclass of :class:`FissionYieldHelper` to instantiate. + reaction_rate_opts : dict + Keyword arguments that are passed to the :class:`ReactionRateHelper` + subclass. + fission_yield_opts : dict + Keyword arguments that are passed to the :class:`FissionYieldHelper` + subclass. + + """ # Get classes to assist working with tallies if reaction_rate_mode == "direct": self._rate_helper = DirectReactionRateHelper( @@ -342,6 +379,48 @@ class Operator(TransportOperator): self._yield_helper = fission_helper.from_operator( self, **fission_yield_opts) + def initial_condition(self): + """Performs final setup and returns initial condition. + + Returns + ------- + list of numpy.ndarray + Total density for initial conditions. + + """ + + # Create XML files + if comm.rank == 0: + self.geometry.export_to_xml() + self.settings.export_to_xml() + self._generate_materials_xml() + + # Initialize OpenMC library + comm.barrier() + if not openmc.lib.is_initialized: + openmc.lib.init(intracomm=comm) + + # Generate tallies in memory + materials = [openmc.lib.materials[int(i)] + for i in self.burnable_mats] + + return super().initial_condition(materials) + + def _generate_materials_xml(self): + """Creates materials.xml from self.number. + + Due to uncertainty with how MPI interacts with OpenMC API, this + constructs the XML manually. The long term goal is to do this + through direct memory writing. + + """ + # Sort nuclides according to order in AtomNumber object + nuclides = list(self.number.nuclides) + for mat in self.materials: + mat._nuclides.sort(key=lambda x: nuclides.index(x[0])) + + self.materials.export_to_xml() + def __call__(self, vec, source_rate): """Runs a simulation. @@ -365,9 +444,7 @@ class Operator(TransportOperator): # Reset results in OpenMC openmc.lib.reset() - # Update the number densities regardless of the source rate - self.number.set_density(vec) - self._update_materials() + self._update_materials_and_nuclides(vec) # If the source rate is zero, return zero reaction rates without running # a transport solve @@ -376,235 +453,20 @@ class Operator(TransportOperator): rates.fill(0.0) return OperatorResult(ufloat(0.0, 0.0), rates) - # Prevent OpenMC from complaining about re-creating tallies - openmc.reset_auto_ids() - - # Update tally nuclides data in preparation for transport solve - nuclides = self._get_tally_nuclides() - self._rate_helper.nuclides = nuclides - self._normalization_helper.nuclides = nuclides - self._yield_helper.update_tally_nuclides(nuclides) - # Run OpenMC openmc.lib.run() openmc.lib.reset_timers() # Extract results - op_result = self._unpack_tallies_and_normalize(source_rate) + rates = self._calculate_reaction_rates(source_rate) + + # Get k and uncertainty + keff = ufloat(*openmc.lib.keff()) + + op_result = OperatorResult(keff, rates) return copy.deepcopy(op_result) - @staticmethod - def write_bos_data(step): - """Write a state-point file with beginning of step data - - Parameters - ---------- - step : int - Current depletion step including restarts - """ - openmc.lib.statepoint_write( - "openmc_simulation_n{}.h5".format(step), - write_source=False) - - def _differentiate_burnable_mats(self): - """Assign distribmats for each burnable material - - """ - - # Count the number of instances for each cell and material - self.geometry.determine_paths(instances_only=True) - - # Extract all burnable materials which have multiple instances - distribmats = set( - [mat for mat in self.materials - if mat.depletable and mat.num_instances > 1]) - - for mat in distribmats: - if mat.volume is None: - raise RuntimeError("Volume not specified for depletable " - "material with ID={}.".format(mat.id)) - mat.volume /= mat.num_instances - - if distribmats: - # Assign distribmats to cells - for cell in self.geometry.get_all_material_cells().values(): - if cell.fill in distribmats: - mat = cell.fill - cell.fill = [mat.clone() - for i in range(cell.num_instances)] - - def _get_burnable_mats(self): - """Determine depletable materials, volumes, and nuclides - - Returns - ------- - burnable_mats : list of str - List of burnable material IDs - volume : OrderedDict of str to float - Volume of each material in [cm^3] - nuclides : list of str - Nuclides in order of how they'll appear in the simulation. - - """ - - burnable_mats = set() - model_nuclides = set() - volume = OrderedDict() - - self.heavy_metal = 0.0 - - # Iterate once through the geometry to get dictionaries - for mat in self.materials: - for nuclide in mat.get_nuclides(): - model_nuclides.add(nuclide) - if mat.depletable: - burnable_mats.add(str(mat.id)) - if mat.volume is None: - raise RuntimeError("Volume not specified for depletable " - "material with ID={}.".format(mat.id)) - volume[str(mat.id)] = mat.volume - self.heavy_metal += mat.fissionable_mass - - # Make sure there are burnable materials - if not burnable_mats: - raise RuntimeError( - "No depletable materials were found in the model.") - - # Sort the sets - burnable_mats = sorted(burnable_mats, key=int) - model_nuclides = sorted(model_nuclides) - - # Construct a global nuclide dictionary, burned first - nuclides = list(self.chain.nuclide_dict) - for nuc in model_nuclides: - if nuc not in nuclides: - nuclides.append(nuc) - - return burnable_mats, volume, nuclides - - def _extract_number(self, local_mats, volume, nuclides, prev_res=None): - """Construct AtomNumber using geometry - - Parameters - ---------- - local_mats : list of str - Material IDs to be managed by this process - volume : OrderedDict of str to float - Volumes for the above materials in [cm^3] - nuclides : list of str - Nuclides to be used in the simulation. - prev_res : Results, optional - Results from a previous depletion calculation - - """ - self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain)) - - if self.dilute_initial != 0.0: - for nuc in self._burnable_nucs: - self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) - - # Now extract and store the number densities - # From the geometry if no previous depletion results - if prev_res is None: - for mat in self.materials: - if str(mat.id) in local_mats: - self._set_number_from_mat(mat) - - # Else from previous depletion results - else: - for mat in self.materials: - if str(mat.id) in local_mats: - self._set_number_from_results(mat, prev_res) - - def _set_number_from_mat(self, mat): - """Extracts material and number densities from openmc.Material - - Parameters - ---------- - mat : openmc.Material - The material to read from - - """ - mat_id = str(mat.id) - - for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): - atom_per_cc = atom_per_bcm * 1.0e24 - self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - - def _set_number_from_results(self, mat, prev_res): - """Extracts material nuclides and number densities. - - If the nuclide concentration's evolution is tracked, the densities come - from depletion results. Else, densities are extracted from the geometry - in the summary. - - Parameters - ---------- - mat : openmc.Material - The material to read from - prev_res : Results - Results from a previous depletion calculation - - """ - mat_id = str(mat.id) - - # Get nuclide lists from geometry and depletion results - depl_nuc = prev_res[-1].nuc_to_ind - geom_nuc_densities = mat.get_nuclide_atom_densities() - - # Merge lists of nuclides, with the same order for every calculation - geom_nuc_densities.update(depl_nuc) - - for nuclide, atom_per_bcm in geom_nuc_densities.items(): - if nuclide in depl_nuc: - concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] - volume = prev_res[-1].volume[mat_id] - atom_per_cc = concentration / volume - else: - atom_per_cc = atom_per_bcm * 1.0e24 - - self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - - def initial_condition(self): - """Performs final setup and returns initial condition. - - Returns - ------- - list of numpy.ndarray - Total density for initial conditions. - """ - - # Create XML files - if comm.rank == 0: - self.geometry.export_to_xml() - self.settings.export_to_xml() - self._generate_materials_xml() - - # Initialize OpenMC library - comm.barrier() - if not openmc.lib.is_initialized: - openmc.lib.init(intracomm=comm) - - # Generate tallies in memory - materials = [openmc.lib.materials[int(i)] - for i in self.burnable_mats] - self._rate_helper.generate_tallies(materials, self.chain.reactions) - self._normalization_helper.prepare( - self.chain.nuclides, self.reaction_rates.index_nuc) - # Tell fission yield helper what materials this process is - # responsible for - self._yield_helper.generate_tallies( - materials, tuple(sorted(self._mat_index_map.values()))) - - # Return number density vector - return list(self.number.get_mat_slice(np.s_[:])) - - def finalize(self): - """Finalize a depletion simulation and release resources.""" - if self.cleanup_when_done: - openmc.lib.finalize() - def _update_materials(self): """Updates material compositions in OpenMC on all processes.""" @@ -631,192 +493,41 @@ class Operator(TransportOperator): densities.append(val) else: # Only output warnings if values are significantly - # negative. CRAM does not guarantee positive values. + # negative. CRAM does not guarantee positive + # values. if val < -1.0e-21: - print("WARNING: nuclide ", nuc, " in material ", mat, - " is negative (density = ", val, " at/barn-cm)") + print( + "WARNING: nuclide ", + nuc, + " in material ", + mat, + " is negative (density = ", + val, + " at/barn-cm)") number_i[mat, nuc] = 0.0 # Update densities on C API side mat_internal = openmc.lib.materials[int(mat)] mat_internal.set_densities(nuclides, densities) - #TODO Update densities on the Python side, otherwise the + # TODO Update densities on the Python side, otherwise the # summary.h5 file contains densities at the first time step - def _generate_materials_xml(self): - """Creates materials.xml from self.number. - - Due to uncertainty with how MPI interacts with OpenMC API, this - constructs the XML manually. The long term goal is to do this - through direct memory writing. - - """ - # Sort nuclides according to order in AtomNumber object - nuclides = list(self.number.nuclides) - for mat in self.materials: - mat._nuclides.sort(key=lambda x: nuclides.index(x[0])) - - self.materials.export_to_xml() - - def _get_tally_nuclides(self): - """Determine nuclides that should be tallied for reaction rates. - - This method returns a list of all nuclides that have neutron data and - are listed in the depletion chain. Technically, we should tally nuclides - that may not appear in the depletion chain because we still need to get - the fission reaction rate for these nuclides in order to normalize - power, but that is left as a future exercise. - - Returns - ------- - list of str - Tally nuclides - - """ - nuc_set = set() - - # Create the set of all nuclides in the decay chain in materials marked - # for burning in which the number density is greater than zero. - for nuc in self.number.nuclides: - if nuc in self.nuclides_with_data: - if np.sum(self.number[:, nuc]) > 0.0: - nuc_set.add(nuc) - - # Communicate which nuclides have nonzeros to rank 0 - if comm.rank == 0: - for i in range(1, comm.size): - nuc_newset = comm.recv(source=i, tag=i) - nuc_set |= nuc_newset - else: - comm.send(nuc_set, dest=0, tag=comm.rank) - - if comm.rank == 0: - # Sort nuclides in the same order as self.number - nuc_list = [nuc for nuc in self.number.nuclides - if nuc in nuc_set] - else: - nuc_list = None - - # Store list of tally nuclides on each process - nuc_list = comm.bcast(nuc_list) - return [nuc for nuc in nuc_list if nuc in self.chain] - - def _unpack_tallies_and_normalize(self, source_rate): - """Unpack tallies from OpenMC and return an operator result - - This method uses OpenMC's C API bindings to determine the k-effective - value and reaction rates from the simulation. The reaction rates are - normalized by a helper class depending on the method being used. + @staticmethod + def write_bos_data(step): + """Write a state-point file with beginning of step data Parameters ---------- - source_rate : float - Power in [W] or source rate in [neutron/sec] - - Returns - ------- - openmc.deplete.OperatorResult - Eigenvalue and reaction rates resulting from transport operator + step : int + Current depletion step including restarts """ - rates = self.reaction_rates - rates.fill(0.0) + openmc.lib.statepoint_write( + "openmc_simulation_n{}.h5".format(step), + write_source=False) - # Get k and uncertainty - keff = ufloat(*openmc.lib.keff()) - - # Extract tally bins - nuclides = self._rate_helper.nuclides - - # Form fast map - nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] - react_ind = [rates.index_rx[react] for react in self.chain.reactions] - - # Keep track of energy produced from all reactions in eV per source - # particle - self._normalization_helper.reset() - self._yield_helper.unpack() - - # Store fission yield dictionaries - fission_yields = [] - - # Create arrays to store fission Q values, reaction rates, and nuclide - # numbers, zeroed out in material iteration - number = np.empty(rates.n_nuc) - - fission_ind = rates.index_rx.get("fission") - - # Extract results - for i, mat in enumerate(self.local_mats): - # Get tally index - mat_index = self._mat_index_map[mat] - - # Zero out reaction rates and nuclide numbers - number.fill(0.0) - - # Get new number densities - for nuc, i_nuc_results in zip(nuclides, nuc_ind): - number[i_nuc_results] = self.number[mat, nuc] - - tally_rates = self._rate_helper.get_material_rates( - mat_index, nuc_ind, react_ind) - - # Compute fission yields for this material - fission_yields.append(self._yield_helper.weighted_yields(i)) - - # Accumulate energy from fission - if fission_ind is not None: - self._normalization_helper.update(tally_rates[:, fission_ind]) - - # Divide by total number and store - rates[i] = self._rate_helper.divide_by_adens(number) - - # Scale reaction rates to obtain units of reactions/sec - rates *= self._normalization_helper.factor(source_rate) - - # Store new fission yields on the chain - self.chain.fission_yields = fission_yields - - return OperatorResult(keff, rates) - - def _get_nuclides_with_data(self, cross_sections): - """Loads cross_sections.xml file to find nuclides with neutron data""" - nuclides = set() - data_lib = DataLibrary.from_xml(cross_sections) - for library in data_lib.libraries: - if library['type'] != 'neutron': - continue - for name in library['materials']: - if name not in nuclides: - nuclides.add(name) - - return nuclides - - def get_results_info(self): - """Returns volume list, material lists, and nuc lists. - - Returns - ------- - volume : dict of str float - Volumes corresponding to materials in full_burn_dict - nuc_list : list of str - A list of all nuclide names. Used for sorting the simulation. - burn_list : list of int - A list of all material IDs to be burned. Used for sorting the simulation. - full_burn_list : list - List of all burnable material IDs - - """ - nuc_list = self.number.burnable_nuclides - burn_list = self.local_mats - - volume = {} - for i, mat in enumerate(burn_list): - volume[mat] = self.number.volume[i] - - # Combine volume dictionaries across processes - volume_list = comm.allgather(volume) - volume = {k: v for d in volume_list for k, v in d.items()} - - return volume, nuc_list, burn_list, self.burnable_mats + def finalize(self): + """Finalize a depletion simulation and release resources.""" + if self.cleanup_when_done: + openmc.lib.finalize() diff --git a/openmc/universe.py b/openmc/universe.py index f3f95aebb..b65263ef9 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,14 @@ class DAGMCUniverse(UniverseBase): string += '{: <16}=\t{}\n'.format('\tFile', self.filename) return string + @property + def bounding_box(self): + with h5py.File(self.filename) as dagmc_file: + coords = dagmc_file['tstt']['nodes']['coordinates'][()] + lower_left_corner = coords.min(axis=0) + upper_right_corner = coords.max(axis=0) + return (lower_left_corner, upper_right_corner) + @property def filename(self): return self._filename 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()