diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 383268f8b..9f7d8c447 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -134,7 +134,7 @@ data, such as number densities and reaction rates for each material. OperatorResult ReactionRates Results - ResultsList + StepResult The following class and functions are used to solve the depletion equations, with :func:`cram.CRAM48` being the default. diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 78a12c4ca..57619611f 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -49,11 +49,11 @@ one of these functions along with the timesteps and power level:: The coupled transport-depletion problem is executed, and once it is done a ``depletion_results.h5`` file is written. The results can be analyzed using the -:class:`openmc.deplete.ResultsList` class. This class has methods that allow for +:class:`openmc.deplete.Results` class. This class has methods that allow for easy retrieval of k-effective, nuclide concentrations, and reaction rates over time:: - results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5") + results = openmc.deplete.Results("depletion_results.h5") time, keff = results.get_keff() Note that the coupling between the transport solver and the transmutation solver diff --git a/examples/pincell_depletion/restart_depletion.py b/examples/pincell_depletion/restart_depletion.py index 72209d2ab..95bbb9954 100644 --- a/examples/pincell_depletion/restart_depletion.py +++ b/examples/pincell_depletion/restart_depletion.py @@ -12,7 +12,7 @@ with openmc.StatePoint(statepoint) as sp: geometry = sp.summary.geometry # Load previous depletion results -previous_results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5") +previous_results = openmc.deplete.Results("depletion_results.h5") ############################################################################### # Transport calculation settings @@ -56,7 +56,7 @@ integrator.integrate() ############################################################################### # Open results file -results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5") +results = openmc.deplete.Results("depletion_results.h5") # Obtain K_eff as a function of time time, keff = results.get_keff(time_units='d') diff --git a/examples/pincell_depletion/run_depletion.py b/examples/pincell_depletion/run_depletion.py index 53b0614d1..ae069334e 100644 --- a/examples/pincell_depletion/run_depletion.py +++ b/examples/pincell_depletion/run_depletion.py @@ -101,7 +101,7 @@ integrator.integrate() ############################################################################### # Open results file -results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5") +results = openmc.deplete.Results("depletion_results.h5") # Obtain K_eff as a function of time time, keff = results.get_keff(time_units='d') diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index b8c1cdfff..5891555a2 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -10,8 +10,8 @@ from .chain import * from .operator import * from .reaction_rates import * from .atom_number import * +from .stepresult import * from .results import * -from .results_list import * from .integrators import * from . import abc from . import cram diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 0d1e8ad41..457c51a8a 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -22,9 +22,9 @@ from uncertainties import ufloat from openmc.lib import MaterialFilter, Tally from openmc.checkvalue import check_type, check_greater_than from openmc.mpi import comm -from .results import Results +from .stepresult import StepResult from .chain import Chain -from .results_list import ResultsList +from .results import Results from .pool import deplete @@ -79,7 +79,7 @@ class TransportOperator(ABC): in initial condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. - prev_results : ResultsList, optional + prev_results : Results, optional Results from a previous depletion calculation. Attributes @@ -88,7 +88,7 @@ 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. - prev_res : ResultsList or None + prev_res : Results or None Results from a previous depletion calculation. ``None`` if no results are to be used. """ @@ -102,7 +102,7 @@ class TransportOperator(ABC): if prev_results is None: self.prev_res = None else: - check_type("previous results", prev_results, ResultsList) + check_type("previous results", prev_results, Results) self.prev_res = prev_results @property @@ -719,6 +719,9 @@ class Integrator(ABC): elif unit.lower() == 'mwd/kg': watt_days_per_kg = 1e6*timestep kilograms = 1e-3*operator.heavy_metal + if rate == 0.0: + raise ValueError("Cannot specify a timestep in [MWd/kg] when" + " the power is zero.") days = watt_days_per_kg * kilograms / rate seconds.append(days*_SECONDS_PER_DAY) else: @@ -842,7 +845,7 @@ class Integrator(ABC): k = ufloat(res.k[0, 0], res.k[0, 1]) # Scale reaction rates by ratio of source rates - rates *= source_rate / res.source_rate[0] + rates *= source_rate / res.source_rate return bos_conc, OperatorResult(k, rates) def _get_start_data(self): @@ -889,7 +892,7 @@ class Integrator(ABC): # Remove actual EOS concentration for next step conc = conc_list.pop() - Results.save(self.operator, conc_list, res_list, [t, t + dt], + StepResult.save(self.operator, conc_list, res_list, [t, t + dt], source_rate, self._i_res + i, proc_time) t += dt @@ -901,7 +904,7 @@ class Integrator(ABC): if output and final_step: print(f"[openmc.deplete] t={t} (final operator evaluation)") res_list = [self.operator(conc, source_rate if final_step else 0.0)] - Results.save(self.operator, [conc], res_list, [t, t], + StepResult.save(self.operator, [conc], res_list, [t, t], source_rate, self._i_res + len(self), proc_time) self.operator.write_bos_data(len(self) + self._i_res) @@ -1048,13 +1051,13 @@ class SIIntegrator(Integrator): # Remove actual EOS concentration for next step conc = conc_list.pop() - Results.save(self.operator, conc_list, res_list, [t, t + dt], + StepResult.save(self.operator, conc_list, res_list, [t, t + dt], p, self._i_res + i, proc_time) t += dt # No final simulation for SIE, use last iteration results - Results.save(self.operator, [conc], [res_list[-1]], [t, t], + StepResult.save(self.operator, [conc], [res_list[-1]], [t, t], p, self._i_res + len(self), proc_time) self.operator.write_bos_data(self._i_res + len(self)) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 1f31a22b6..327c762e0 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -25,7 +25,7 @@ from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .chain import _find_chain_file from .reaction_rates import ReactionRates -from .results_list import ResultsList +from .results import Results from .helpers import ( DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, @@ -94,7 +94,7 @@ class Operator(TransportOperator): Path to the depletion chain XML file. Defaults to the file listed under ``depletion_chain`` in :envvar:`OPENMC_CROSS_SECTIONS` environment variable. - prev_results : ResultsList, optional + 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. @@ -192,7 +192,7 @@ class Operator(TransportOperator): Initial heavy metal inventory [g] local_mats : list of str All burnable material IDs being managed by a single process - prev_res : ResultsList or None + prev_res : Results or None Results from a previous depletion calculation. ``None`` if no results are to be used. cleanup_when_done : bool @@ -284,7 +284,7 @@ class Operator(TransportOperator): if comm.size == 1: self.prev_res = prev_results else: - self.prev_res = ResultsList() + 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) @@ -494,7 +494,7 @@ class Operator(TransportOperator): Volumes for the above materials in [cm^3] nuclides : list of str Nuclides to be used in the simulation. - prev_res : ResultsList, optional + prev_res : Results, optional Results from a previous depletion calculation """ @@ -543,7 +543,7 @@ class Operator(TransportOperator): ---------- mat : openmc.Material The material to read from - prev_res : ResultsList + prev_res : Results Results from a previous depletion calculation """ diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index 3916fb4d7..c238002fe 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -1,521 +1,416 @@ -"""The results module. - -Contains results generation and saving capabilities. -""" - -from collections import OrderedDict -import copy +import numbers +import bisect +import math +from warnings import warn import h5py import numpy as np -import openmc -from openmc.mpi import comm, MPI -from .reaction_rates import ReactionRates +from .stepresult import StepResult, VERSION_RESULTS +import openmc.checkvalue as cv +from openmc.data.library import DataLibrary +from openmc.material import Material, Materials +from openmc.exceptions import DataError -VERSION_RESULTS = (1, 1) +__all__ = ["Results", "ResultsList"] -__all__ = ["Results"] +def _get_time_as(seconds, units): + if units == "d": + return seconds / (60 * 60 * 24) + elif units == "h": + return seconds / (60 * 60) + elif units == "min": + return seconds / 60 + else: + return seconds -class Results: - """Output of a depletion run +class Results(list): + """Results from a depletion simulation - Attributes + The :class:`Results` class acts as a list that stores the results from + each depletion step and provides extra methods for interrogating these + results. + + Parameters ---------- - k : list of (float, float) - Eigenvalue and uncertainty for each substep. - time : list of float - Time at beginning, end of step, in seconds. - source_rate : float - Source rate during timestep in [W] or [neutron/sec] - n_mat : int - Number of mats. - n_nuc : int - Number of nuclides. - rates : list of ReactionRates - The reaction rates for each substep. - volume : OrderedDict of str to float - Dictionary mapping mat id to volume. - mat_to_ind : OrderedDict of str to int - A dictionary mapping mat ID as string to index. - nuc_to_ind : OrderedDict of str to int - A dictionary mapping nuclide name as string to index. - mat_to_hdf5_ind : OrderedDict of str to int - A dictionary mapping mat ID as string to global index. - n_hdf5_mats : int - Number of materials in entire geometry. - n_stages : int - Number of stages in simulation. - data : numpy.ndarray - Atom quantity, stored by stage, mat, then by nuclide. - proc_time : int - Average time spent depleting a material across all - materials and processes + filename : str + Path to depletion result file """ - def __init__(self): - self.k = None - self.time = None - self.source_rate = None - self.rates = None - self.volume = None - self.proc_time = None + def __init__(self, filename=None): + data = [] + if filename is not None: + with h5py.File(str(filename), "r") as fh: + cv.check_filetype_version(fh, 'depletion results', VERSION_RESULTS[0]) - self.mat_to_ind = None - self.nuc_to_ind = None - self.mat_to_hdf5_ind = None + # Get number of results stored + n = fh["number"][...].shape[0] - self.data = None + for i in range(n): + data.append(StepResult.from_hdf5(fh, i)) + super().__init__(data) - def __getitem__(self, pos): - """Retrieves an item from results. - Parameters - ---------- - pos : tuple - A three-length tuple containing a stage index, mat index and a nuc - index. All can be integers or slices. The second two can be - strings corresponding to their respective dictionary. - - Returns - ------- - float - The atoms for stage, mat, nuc - - """ - stage, mat, nuc = pos - if isinstance(mat, openmc.Material): - mat = str(mat.id) - if isinstance(mat, str): - mat = self.mat_to_ind[mat] - if isinstance(nuc, str): - nuc = self.nuc_to_ind[nuc] - - return self.data[stage, mat, nuc] - - def __setitem__(self, pos, val): - """Sets an item from results. - - Parameters - ---------- - pos : tuple - A three-length tuple containing a stage index, mat index and a nuc - index. All can be integers or slices. The second two can be - strings corresponding to their respective dictionary. - - val : float - The value to set data to. - - """ - stage, mat, nuc = pos - if isinstance(mat, str): - mat = self.mat_to_ind[mat] - if isinstance(nuc, str): - nuc = self.nuc_to_ind[nuc] - - self.data[stage, mat, nuc] = val - - @property - def n_mat(self): - return len(self.mat_to_ind) - - @property - def n_nuc(self): - return len(self.nuc_to_ind) - - @property - def n_hdf5_mats(self): - return len(self.mat_to_hdf5_ind) - - @property - def n_stages(self): - return self.data.shape[0] - - def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages): - """Allocates memory of Results. - - Parameters - ---------- - volume : dict of str float - Volumes corresponding to materials in full_burn_dict - nuc_list : list of str - A list of all nuclide names. Used for sorting the simulation. - burn_list : list of int - A list of all mat IDs to be burned. Used for sorting the simulation. - full_burn_list : list of str - List of all burnable material IDs - stages : int - Number of stages in simulation. - - """ - self.volume = copy.deepcopy(volume) - self.nuc_to_ind = {nuc: i for i, nuc in enumerate(nuc_list)} - self.mat_to_ind = {mat: i for i, mat in enumerate(burn_list)} - self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)} - - # Create storage array - self.data = np.zeros((stages, self.n_mat, self.n_nuc)) - - def distribute(self, local_materials, ranges): - """Create a new object containing data for distributed materials - - Parameters - ---------- - local_materials : iterable of str - Materials for this process - ranges : iterable of int - Slice-like object indicating indicies of ``local_materials`` - in the material dimension of :attr:`data` and each element - in :attr:`rates` - - Returns - ------- - Results - New results object - """ - new = Results() - new.volume = {lm: self.volume[lm] for lm in local_materials} - new.mat_to_ind = {mat: idx for (idx, mat) in enumerate(local_materials)} - - # Direct transfer - direct_attrs = ("time", "k", "source_rate", "nuc_to_ind", - "mat_to_hdf5_ind", "proc_time") - for attr in direct_attrs: - setattr(new, attr, getattr(self, attr)) - # Get applicable slice of data - new.data = self.data[:, ranges] - new.rates = [r[ranges] for r in self.rates] - return new - - def export_to_hdf5(self, filename, step): - """Export results to an HDF5 file + @classmethod + def from_hdf5(cls, filename): + """Load in depletion results from a previous file Parameters ---------- filename : str - The filename to write to - step : int - What step is this? + Path to depletion result file + + Returns + ------- + Results + New instance of depletion results """ - # Write new file if first time step, else add to existing file - kwargs = {'mode': "w" if step == 0 else "a"} + warn( + "The ResultsList.from_hdf5(...) method is no longer necessary and will " + "be removed in a future version of OpenMC. Use Results(...) instead.", + FutureWarning + ) + return cls(filename) - if h5py.get_config().mpi and comm.size > 1: - # Write results in parallel - kwargs['driver'] = 'mpio' - kwargs['comm'] = comm - with h5py.File(filename, **kwargs) as handle: - self._to_hdf5(handle, step, parallel=True) + def get_atoms(self, mat, nuc, nuc_units="atoms", time_units="s"): + """Get number of nuclides over time from a single material + + .. note:: + Initial values for some isotopes that do not appear in + initial concentrations may be non-zero, depending on the + value of the :attr:`openmc.deplete.Operator.dilute_initial` + attribute. The :class:`openmc.deplete.Operator` class adds isotopes + according to this setting, which can be set to zero. + + Parameters + ---------- + mat : openmc.Material, str + Material object or material id to evaluate + nuc : str + Nuclide name to evaluate + nuc_units : {"atoms", "atom/b-cm", "atom/cm3"}, optional + Units for the returned concentration. Default is ``"atoms"`` + + .. versionadded:: 0.12 + time_units : {"s", "min", "h", "d"}, optional + Units for the returned time array. Default is ``"s"`` to + return the value in seconds. + + .. versionadded:: 0.12 + + Returns + ------- + times : numpy.ndarray + Array of times in units of ``time_units`` + concentrations : numpy.ndarray + Concentration of specified nuclide in units of ``nuc_units`` + + """ + cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) + cv.check_value("nuc_units", nuc_units, + {"atoms", "atom/b-cm", "atom/cm3"}) + + if isinstance(mat, Material): + mat_id = str(mat.id) + elif isinstance(mat, str): + mat_id = mat else: - # Gather results at root process - all_results = comm.gather(self) + raise TypeError('mat should be of type openmc.Material or str') + times = np.empty_like(self, dtype=float) + concentrations = np.empty_like(self, dtype=float) - # Only root process writes results - if comm.rank == 0: - with h5py.File(filename, **kwargs) as handle: - for res in all_results: - res._to_hdf5(handle, step, parallel=False) + # Evaluate value in each region + for i, result in enumerate(self): + times[i] = result.time[0] + concentrations[i] = result[0, mat_id, nuc] - def _write_hdf5_metadata(self, handle): - """Writes result metadata in HDF5 file + # Unit conversions + times = _get_time_as(times, time_units) + if nuc_units != "atoms": + # Divide by volume to get density + concentrations /= self[0].volume[mat_id] + if nuc_units == "atom/b-cm": + # 1 barn = 1e-24 cm^2 + concentrations *= 1e-24 + + return times, concentrations + + def get_reaction_rate(self, mat, nuc, rx): + """Get reaction rate in a single material/nuclide over time + + .. note:: + + Initial values for some isotopes that do not appear in + initial concentrations may be non-zero, depending on the + value of :class:`openmc.deplete.Operator` ``dilute_initial`` + The :class:`openmc.deplete.Operator` adds isotopes according + to this setting, which can be set to zero. Parameters ---------- - handle : h5py.File or h5py.Group - An hdf5 file or group type to store this in. + mat : openmc.Material, str + Material object or material id to evaluate + nuc : str + Nuclide name to evaluate + rx : str + Reaction rate to evaluate + + Returns + ------- + times : numpy.ndarray + Array of times in [s] + rates : numpy.ndarray + Array of reaction rates """ - # Create and save the 5 dictionaries: - # quantities - # self.mat_to_ind -> self.volume (TODO: support for changing volumes) - # self.nuc_to_ind - # reactions - # self.rates[0].nuc_to_ind (can be different from above, above is superset) - # self.rates[0].react_to_ind - # these are shared by every step of the simulation, and should be deduplicated. + times = np.empty_like(self, dtype=float) + rates = np.empty_like(self, dtype=float) - # Store concentration mat and nuclide dictionaries (along with volumes) - - handle.attrs['version'] = np.array(VERSION_RESULTS) - handle.attrs['filetype'] = np.string_('depletion results') - - mat_list = sorted(self.mat_to_hdf5_ind, key=int) - nuc_list = sorted(self.nuc_to_ind) - rxn_list = sorted(self.rates[0].index_rx) - - n_mats = self.n_hdf5_mats - n_nuc_number = len(nuc_list) - n_nuc_rxn = len(self.rates[0].index_nuc) - n_rxn = len(rxn_list) - n_stages = self.n_stages - - mat_group = handle.create_group("materials") - - for mat in mat_list: - mat_single_group = mat_group.create_group(mat) - mat_single_group.attrs["index"] = self.mat_to_hdf5_ind[mat] - mat_single_group.attrs["volume"] = self.volume[mat] - - nuc_group = handle.create_group("nuclides") - - for nuc in nuc_list: - nuc_single_group = nuc_group.create_group(nuc) - nuc_single_group.attrs["atom number index"] = self.nuc_to_ind[nuc] - if nuc in self.rates[0].index_nuc: - nuc_single_group.attrs["reaction rate index"] = self.rates[0].index_nuc[nuc] - - rxn_group = handle.create_group("reactions") - - for rxn in rxn_list: - rxn_single_group = rxn_group.create_group(rxn) - rxn_single_group.attrs["index"] = self.rates[0].index_rx[rxn] - - # Construct array storage - - handle.create_dataset("number", (1, n_stages, n_mats, n_nuc_number), - maxshape=(None, n_stages, n_mats, n_nuc_number), - chunks=(1, 1, n_mats, n_nuc_number), - dtype='float64') - - handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn), - maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn), - chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn), - dtype='float64') - - handle.create_dataset("eigenvalues", (1, n_stages, 2), - maxshape=(None, n_stages, 2), dtype='float64') - - handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64') - - handle.create_dataset("source_rate", (1, n_stages), maxshape=(None, n_stages), - dtype='float64') - - handle.create_dataset( - "depletion time", (1,), maxshape=(None,), - dtype="float64") - - def _to_hdf5(self, handle, index, parallel=False): - """Converts results object into an hdf5 object. - - Parameters - ---------- - handle : h5py.File or h5py.Group - An HDF5 file or group type to store this in. - index : int - What step is this? - parallel : bool - Being called with parallel HDF5? - - """ - if "/number" not in handle: - if parallel: - comm.barrier() - self._write_hdf5_metadata(handle) - - if parallel: - comm.barrier() - - # Grab handles - number_dset = handle["/number"] - rxn_dset = handle["/reaction rates"] - eigenvalues_dset = handle["/eigenvalues"] - time_dset = handle["/time"] - source_rate_dset = handle["/source_rate"] - proc_time_dset = handle["/depletion time"] - - # Get number of results stored - number_shape = list(number_dset.shape) - number_results = number_shape[0] - - new_shape = index + 1 - - if number_results < new_shape: - # Extend first dimension by 1 - number_shape[0] = new_shape - number_dset.resize(number_shape) - - rxn_shape = list(rxn_dset.shape) - rxn_shape[0] = new_shape - rxn_dset.resize(rxn_shape) - - eigenvalues_shape = list(eigenvalues_dset.shape) - eigenvalues_shape[0] = new_shape - eigenvalues_dset.resize(eigenvalues_shape) - - time_shape = list(time_dset.shape) - time_shape[0] = new_shape - time_dset.resize(time_shape) - - source_rate_shape = list(source_rate_dset.shape) - source_rate_shape[0] = new_shape - source_rate_dset.resize(source_rate_shape) - - proc_shape = list(proc_time_dset.shape) - proc_shape[0] = new_shape - proc_time_dset.resize(proc_shape) - - # If nothing to write, just return - if len(self.mat_to_ind) == 0: - return - - # Add data - # Note, for the last step, self.n_stages = 1, even if n_stages != 1. - n_stages = self.n_stages - inds = [self.mat_to_hdf5_ind[mat] for mat in self.mat_to_ind] - low = min(inds) - high = max(inds) - for i in range(n_stages): - number_dset[index, i, low:high+1] = self.data[i] - rxn_dset[index, i, low:high+1] = self.rates[i] - if comm.rank == 0: - eigenvalues_dset[index, i] = self.k[i] - if comm.rank == 0: - time_dset[index] = self.time - source_rate_dset[index] = self.source_rate - if self.proc_time is not None: - proc_time_dset[index] = ( - self.proc_time / (comm.size * self.n_hdf5_mats) - ) - - @classmethod - def from_hdf5(cls, handle, step): - """Loads results object from HDF5. - - Parameters - ---------- - handle : h5py.File or h5py.Group - An HDF5 file or group type to load from. - step : int - Index for depletion step - """ - results = cls() - - # Grab handles - number_dset = handle["/number"] - eigenvalues_dset = handle["/eigenvalues"] - time_dset = handle["/time"] - if "source_rate" in handle: - source_rate_dset = handle["/source_rate"] + if isinstance(mat, Material): + mat_id = str(mat.id) + elif isinstance(mat, str): + mat_id = mat else: - # Older versions used "power" instead of "source_rate" - source_rate_dset = handle["/power"] + raise TypeError('mat should be of type openmc.Material or str') - results.data = number_dset[step, :, :, :] - results.k = eigenvalues_dset[step, :] - results.time = time_dset[step, :] - results.source_rate = source_rate_dset[step, :] + # Evaluate value in each region + for i, result in enumerate(self): + times[i] = result.time[0] + rates[i] = result.rates[0].get(mat_id, nuc, rx) * result[0, mat, nuc] - if "depletion time" in handle: - proc_time_dset = handle["/depletion time"] - if step < proc_time_dset.shape[0]: - results.proc_time = proc_time_dset[step] + return times, rates - if results.proc_time is None: - results.proc_time = np.array([np.nan]) + def get_keff(self, time_units='s'): + """Evaluates the eigenvalue from a results list. - # Reconstruct dictionaries - results.volume = OrderedDict() - results.mat_to_ind = OrderedDict() - results.nuc_to_ind = OrderedDict() - rxn_nuc_to_ind = OrderedDict() - rxn_to_ind = OrderedDict() - - for mat, mat_handle in handle["/materials"].items(): - vol = mat_handle.attrs["volume"] - ind = mat_handle.attrs["index"] - - results.volume[mat] = vol - results.mat_to_ind[mat] = ind - - for nuc, nuc_handle in handle["/nuclides"].items(): - ind_atom = nuc_handle.attrs["atom number index"] - results.nuc_to_ind[nuc] = ind_atom - - if "reaction rate index" in nuc_handle.attrs: - rxn_nuc_to_ind[nuc] = nuc_handle.attrs["reaction rate index"] - - for rxn, rxn_handle in handle["/reactions"].items(): - rxn_to_ind[rxn] = rxn_handle.attrs["index"] - - results.rates = [] - # Reconstruct reactions - for i in range(results.n_stages): - rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind, True) - - rate[:] = handle["/reaction rates"][step, i, :, :, :] - results.rates.append(rate) - - return results - - @staticmethod - def save(op, x, op_results, t, source_rate, step_ind, proc_time=None): - """Creates and writes depletion results to disk + .. versionadded:: 0.13.1 Parameters ---------- - op : openmc.deplete.TransportOperator - The operator used to generate these results. - x : list of list of numpy.array - The prior x vectors. Indexed [i][cell] using the above equation. - op_results : list of openmc.deplete.OperatorResult - Results of applying transport operator - t : list of float - Time indices. - source_rate : float - Source rate during time step in [W] or [neutron/sec] - step_ind : int - Step index. - proc_time : float or None - Total process time spent depleting materials. This may - be process-dependent and will be reduced across MPI - processes. + time_units : {"s", "d", "h", "min"}, optional + Desired units for the times array + + Returns + ------- + times : numpy.ndarray + Array of times in specified units + eigenvalues : numpy.ndarray + k-eigenvalue at each time. Column 0 + contains the eigenvalue, while column + 1 contains the associated uncertainty """ - # Get indexing terms - vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info() + cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) - stages = len(x) + times = np.empty_like(self, dtype=float) + eigenvalues = np.empty((len(self), 2), dtype=float) - # Create results - results = Results() - results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages) + # Get time/eigenvalue at each point + for i, result in enumerate(self): + times[i] = result.time[0] + eigenvalues[i] = result.k[0] - n_mat = len(burn_list) + # Convert time units if necessary + times = _get_time_as(times, time_units) + return times, eigenvalues - for i in range(stages): - for mat_i in range(n_mat): - results[i, mat_i, :] = x[i][mat_i] + def get_eigenvalue(self, time_units='s'): + warn("The get_eigenvalue(...) function has been renamed get_keff and " + "will be removed in a future version of OpenMC.", FutureWarning) + return self.get_keff(time_units) - results.k = [(r.k.nominal_value, r.k.std_dev) for r in op_results] - results.rates = [r.rates for r in op_results] - results.time = t - results.source_rate = source_rate - results.proc_time = proc_time - if results.proc_time is not None: - results.proc_time = comm.reduce(proc_time, op=MPI.SUM) + def get_depletion_time(self): + """Return an array of the average time to deplete a material - results.export_to_hdf5("depletion_results.h5", step_ind) + .. note:: + The return value will have one fewer values than several other + methods, such as :meth:`get_keff`, because no depletion is performed + at the final transport stage. - def transfer_volumes(self, model): - """Transfers volumes from depletion results to geometry - - Parameters - ---------- - model : OpenMC model to be used in a depletion restart - calculation + Returns + ------- + times : numpy.ndarray + Vector of average time to deplete a single material + across all processes and materials. """ - - if not model.materials: - materials = openmc.Materials( - model.geometry.get_all_materials().values() - ) + times = np.empty(len(self) - 1) + # Need special logic because the predictor + # writes EOS values for step i as BOS values + # for step i+1 + # The first proc_time may be zero + if self[0].proc_time > 0.0: + items = self[:-1] else: - materials = model.materials + items = self[1:] + for ix, res in enumerate(items): + times[ix] = res.proc_time + return times - for material in materials: - if material.depletable: - material.volume = self.volume[str(material.id)] + def get_times(self, time_units="d") -> np.ndarray: + """Return the points in time that define the depletion schedule + + .. versionadded:: 0.12.1 + + Parameters + ---------- + time_units : {"s", "d", "h", "min"}, optional + Return the vector in these units. Default is to + convert to days + + Returns + ------- + numpy.ndarray + 1-D vector of time points + + """ + cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) + + times = np.fromiter( + (r.time[0] for r in self), + dtype=self[0].time.dtype, + count=len(self), + ) + + return _get_time_as(times, time_units) + + def get_step_where( + self, time, time_units="d", atol=1e-6, rtol=1e-3 + ) -> int: + """Return the index closest to a given point in time + + In the event ``time`` lies exactly between two points, the + lower index will be returned. It is possible that the index + will be at most one past the point in time requested, but only + according to tolerances requested. + + Passing ``atol=math.inf`` and ``rtol=math.inf`` will return + the closest index to the requested point. + + .. versionadded:: 0.12.1 + + Parameters + ---------- + time : float + Desired point in time + time_units : {"s", "d", "min", "h"}, optional + Units on ``time``. Default: days + atol : float, optional + Absolute tolerance (in ``time_units``) if ``time`` is not + found. + rtol : float, optional + Relative tolerance if ``time`` is not found. + + Returns + ------- + int + + """ + cv.check_type("time", time, numbers.Real) + cv.check_type("atol", atol, numbers.Real) + cv.check_type("rtol", rtol, numbers.Real) + + times = self.get_times(time_units) + + if times[0] < time < times[-1]: + ix = bisect.bisect_left(times, time) + if ix == times.size: + ix -= 1 + # Bisection will place us either directly on the point + # or one-past the first value less than time + elif time - times[ix - 1] <= times[ix] - time: + ix -= 1 + elif times[0] >= time: + ix = 0 + elif time >= times[-1]: + ix = times.size - 1 + + if math.isclose(time, times[ix], rel_tol=rtol, abs_tol=atol): + return ix + + raise ValueError( + "A value of {} {} was not found given absolute and " + "relative tolerances {} and {}.".format( + time, time_units, atol, rtol) + ) + + def export_to_materials(self, burnup_index, nuc_with_data=None) -> Materials: + """Return openmc.Materials object based on results at a given step + + .. versionadded:: 0.12.1 + + Parameters + ---------- + burn_index : int + Index of burnup step to evaluate. See also: get_step_where for + obtaining burnup step indices from other data such as the time. + nuc_with_data : Iterable of str, optional + Nuclides to include in resulting materials. + This can be specified if not all nuclides appearing in + depletion results have associated neutron cross sections, and + as such cannot be used in subsequent transport calculations. + If not provided, nuclides from the cross_sections element of + materials.xml will be used. If that element is not present, + nuclides from OPENMC_CROSS_SECTIONS will be used. + + Returns + ------- + mat_file : Materials + A modified Materials instance containing depleted material data + and original isotopic compositions of non-depletable materials + """ + result = self[burnup_index] + + # Only materials found in the original materials.xml file will be + # updated. If for some reason you have modified OpenMC to produce + # new materials as depletion takes place, this method will not + # work as expected and leave out that material. + mat_file = Materials.from_xml("materials.xml") + + # Only nuclides with valid transport data will be written to + # the new materials XML file. The precedence of nuclides to select + # is first ones provided as a kwarg here, then ones specified + # in the materials.xml file if provided, then finally from + # the environment variable OPENMC_CROSS_SECTIONS. + if nuc_with_data: + cv.check_iterable_type('nuclide names', nuc_with_data, str) + available_cross_sections = nuc_with_data + else: + # select cross_sections.xml file to use + if mat_file.cross_sections: + this_library = DataLibrary.from_xml(path=mat_file.cross_sections) + else: + this_library = DataLibrary.from_xml() + + # Find neutron libraries we have access to + available_cross_sections = set() + for lib in this_library.libraries: + if lib['type'] == 'neutron': + available_cross_sections.update(lib['materials']) + if not available_cross_sections: + raise DataError('No neutron libraries found in cross_sections.xml') + + # Overwrite material definitions, if they can be found in the depletion + # results, and save them to the new depleted xml file. + for mat in mat_file: + mat_id = str(mat.id) + if mat_id in result.mat_to_ind: + mat.volume = result.volume[mat_id] + mat.set_density('sum') + for nuc in result.nuc_to_ind: + if nuc not in available_cross_sections: + continue + atoms = result[0, mat_id, nuc] + if atoms > 0.0: + atoms_per_barn_cm = 1e-24 * atoms / mat.volume + mat.remove_nuclide(nuc) # Replace if it's there + mat.add_nuclide(nuc, atoms_per_barn_cm) + + return mat_file + + +# Retain deprecated name for the time being +ResultsList = Results diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py deleted file mode 100644 index f2f2be0ab..000000000 --- a/openmc/deplete/results_list.py +++ /dev/null @@ -1,396 +0,0 @@ -import numbers -import bisect -import math -from warnings import warn - -import h5py -import numpy as np - -from .results import Results, VERSION_RESULTS -import openmc.checkvalue as cv -from openmc.data.library import DataLibrary -from openmc.material import Material, Materials -from openmc.exceptions import DataError, InvalidArgumentError - -__all__ = ["ResultsList"] - - -def _get_time_as(seconds, units): - if units == "d": - return seconds / (60 * 60 * 24) - elif units == "h": - return seconds / (60 * 60) - elif units == "min": - return seconds / 60 - else: - return seconds - - -class ResultsList(list): - """A list of openmc.deplete.Results objects - - It is recommended to use :meth:`from_hdf5` over - direct creation. - """ - - @classmethod - def from_hdf5(cls, filename): - """Load in depletion results from a previous file - - Parameters - ---------- - filename : str - Path to depletion result file - - Returns - ------- - new : ResultsList - New instance of depletion results - """ - with h5py.File(str(filename), "r") as fh: - cv.check_filetype_version(fh, 'depletion results', VERSION_RESULTS[0]) - new = cls() - - # Get number of results stored - n = fh["number"][...].shape[0] - - for i in range(n): - new.append(Results.from_hdf5(fh, i)) - return new - - def get_atoms(self, mat, nuc, nuc_units="atoms", time_units="s"): - """Get number of nuclides over time from a single material - - .. note:: - Initial values for some isotopes that do not appear in - initial concentrations may be non-zero, depending on the - value of :class:`openmc.deplete.Operator` ``dilute_initial``. - The :class:`openmc.deplete.Operator` adds isotopes according - to this setting, which can be set to zero. - - Parameters - ---------- - mat : openmc.Material, str - Material object or material id to evaluate - nuc : str - Nuclide name to evaluate - nuc_units : {"atoms", "atom/b-cm", "atom/cm3"}, optional - Units for the returned concentration. Default is ``"atoms"`` - - .. versionadded:: 0.12 - time_units : {"s", "min", "h", "d"}, optional - Units for the returned time array. Default is ``"s"`` to - return the value in seconds. - - .. versionadded:: 0.12 - - Returns - ------- - times : numpy.ndarray - Array of times in units of ``time_units`` - concentrations : numpy.ndarray - Concentration of specified nuclide in units of ``nuc_units`` - - """ - cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) - cv.check_value("nuc_units", nuc_units, - {"atoms", "atom/b-cm", "atom/cm3"}) - - if isinstance(mat, Material): - mat_id = str(mat.id) - elif isinstance(mat, str): - mat_id = mat - else: - raise TypeError('mat should be of type openmc.Material or str') - times = np.empty_like(self, dtype=float) - concentrations = np.empty_like(self, dtype=float) - - # Evaluate value in each region - for i, result in enumerate(self): - times[i] = result.time[0] - concentrations[i] = result[0, mat_id, nuc] - - # Unit conversions - times = _get_time_as(times, time_units) - if nuc_units != "atoms": - # Divide by volume to get density - concentrations /= self[0].volume[mat_id] - if nuc_units == "atom/b-cm": - # 1 barn = 1e-24 cm^2 - concentrations *= 1e-24 - - return times, concentrations - - def get_reaction_rate(self, mat, nuc, rx): - """Get reaction rate in a single material/nuclide over time - - .. note:: - - Initial values for some isotopes that do not appear in - initial concentrations may be non-zero, depending on the - value of :class:`openmc.deplete.Operator` ``dilute_initial`` - The :class:`openmc.deplete.Operator` adds isotopes according - to this setting, which can be set to zero. - - Parameters - ---------- - mat : openmc.Material, str - Material object or material id to evaluate - nuc : str - Nuclide name to evaluate - rx : str - Reaction rate to evaluate - - Returns - ------- - times : numpy.ndarray - Array of times in [s] - rates : numpy.ndarray - Array of reaction rates - - """ - times = np.empty_like(self, dtype=float) - rates = np.empty_like(self, dtype=float) - - if isinstance(mat, Material): - mat_id = str(mat.id) - elif isinstance(mat, str): - mat_id = mat - else: - raise TypeError('mat should be of type openmc.Material or str') - - # Evaluate value in each region - for i, result in enumerate(self): - times[i] = result.time[0] - rates[i] = result.rates[0].get(mat_id, nuc, rx) * result[0, mat, nuc] - - return times, rates - - def get_keff(self, time_units='s'): - """Evaluates the eigenvalue from a results list. - - Parameters - ---------- - time_units : {"s", "d", "h", "min"}, optional - Desired units for the times array - - Returns - ------- - times : numpy.ndarray - Array of times in specified units - eigenvalues : numpy.ndarray - k-eigenvalue at each time. Column 0 - contains the eigenvalue, while column - 1 contains the associated uncertainty - - """ - cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) - - times = np.empty_like(self, dtype=float) - eigenvalues = np.empty((len(self), 2), dtype=float) - - # Get time/eigenvalue at each point - for i, result in enumerate(self): - times[i] = result.time[0] - eigenvalues[i] = result.k[0] - - # Convert time units if necessary - times = _get_time_as(times, time_units) - return times, eigenvalues - - def get_eigenvalue(self, time_units='s'): - warn("The get_eigenvalue(...) function has been renamed get_keff and " - "will be removed in a future version of OpenMC.", FutureWarning) - return self.get_keff(time_units) - - - def get_depletion_time(self): - """Return an array of the average time to deplete a material - - .. note:: - The return value will have one fewer values than several other - methods, such as :meth:`get_keff`, because no depletion is performed - at the final transport stage. - - Returns - ------- - times : numpy.ndarray - Vector of average time to deplete a single material - across all processes and materials. - - """ - times = np.empty(len(self) - 1) - # Need special logic because the predictor - # writes EOS values for step i as BOS values - # for step i+1 - # The first proc_time may be zero - if self[0].proc_time > 0.0: - items = self[:-1] - else: - items = self[1:] - for ix, res in enumerate(items): - times[ix] = res.proc_time - return times - - def get_times(self, time_units="d") -> np.ndarray: - """Return the points in time that define the depletion schedule - - - .. versionadded:: 0.12.1 - - Parameters - ---------- - time_units : {"s", "d", "h", "min"}, optional - Return the vector in these units. Default is to - convert to days - - Returns - ------- - numpy.ndarray - 1-D vector of time points - - """ - cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) - - times = np.fromiter( - (r.time[0] for r in self), - dtype=self[0].time.dtype, - count=len(self), - ) - - return _get_time_as(times, time_units) - - def get_step_where( - self, time, time_units="d", atol=1e-6, rtol=1e-3 - ) -> int: - """Return the index closest to a given point in time - - In the event ``time`` lies exactly between two points, the - lower index will be returned. It is possible that the index - will be at most one past the point in time requested, but only - according to tolerances requested. - - Passing ``atol=math.inf`` and ``rtol=math.inf`` will return - the closest index to the requested point. - - - .. versionadded:: 0.12.1 - - Parameters - ---------- - time : float - Desired point in time - time_units : {"s", "d", "min", "h"}, optional - Units on ``time``. Default: days - atol : float, optional - Absolute tolerance (in ``time_units``) if ``time`` is not - found. - rtol : float, optional - Relative tolerance if ``time`` is not found. - - Returns - ------- - int - - """ - cv.check_type("time", time, numbers.Real) - cv.check_type("atol", atol, numbers.Real) - cv.check_type("rtol", rtol, numbers.Real) - - times = self.get_times(time_units) - - if times[0] < time < times[-1]: - ix = bisect.bisect_left(times, time) - if ix == times.size: - ix -= 1 - # Bisection will place us either directly on the point - # or one-past the first value less than time - elif time - times[ix - 1] <= times[ix] - time: - ix -= 1 - elif times[0] >= time: - ix = 0 - elif time >= times[-1]: - ix = times.size - 1 - - if math.isclose(time, times[ix], rel_tol=rtol, abs_tol=atol): - return ix - - raise ValueError( - "A value of {} {} was not found given absolute and " - "relative tolerances {} and {}.".format( - time, time_units, atol, rtol) - ) - - def export_to_materials(self, burnup_index, nuc_with_data=None) -> Materials: - """Return openmc.Materials object based on results at a given step - - .. versionadded:: 0.12.1 - - Parameters - ---------- - burn_index : int - Index of burnup step to evaluate. See also: get_step_where for - obtaining burnup step indices from other data such as the time. - nuc_with_data : Iterable of str, optional - Nuclides to include in resulting materials. - This can be specified if not all nuclides appearing in - depletion results have associated neutron cross sections, and - as such cannot be used in subsequent transport calculations. - If not provided, nuclides from the cross_sections element of - materials.xml will be used. If that element is not present, - nuclides from OPENMC_CROSS_SECTIONS will be used. - - Returns - ------- - mat_file : Materials - A modified Materials instance containing depleted material data - and original isotopic compositions of non-depletable materials - """ - result = self[burnup_index] - - # Only materials found in the original materials.xml file will be - # updated. If for some reason you have modified OpenMC to produce - # new materials as depletion takes place, this method will not - # work as expected and leave out that material. - mat_file = Materials.from_xml("materials.xml") - - # Only nuclides with valid transport data will be written to - # the new materials XML file. The precedence of nuclides to select - # is first ones provided as a kwarg here, then ones specified - # in the materials.xml file if provided, then finally from - # the environment variable OPENMC_CROSS_SECTIONS. - if nuc_with_data: - cv.check_iterable_type('nuclide names', nuc_with_data, str) - available_cross_sections = nuc_with_data - else: - # select cross_sections.xml file to use - if mat_file.cross_sections: - this_library = DataLibrary.from_xml(path=mat_file.cross_sections) - else: - this_library = DataLibrary.from_xml() - - # Find neutron libraries we have access to - available_cross_sections = set() - for lib in this_library.libraries: - if lib['type'] == 'neutron': - available_cross_sections.update(lib['materials']) - if not available_cross_sections: - raise DataError('No neutron libraries found in cross_sections.xml') - - # Overwrite material definitions, if they can be found in the depletion - # results, and save them to the new depleted xml file. - for mat in mat_file: - mat_id = str(mat.id) - if mat_id in result.mat_to_ind: - mat.volume = result.volume[mat_id] - mat.set_density('sum') - for nuc in result.nuc_to_ind: - if nuc not in available_cross_sections: - continue - atoms = result[0, mat_id, nuc] - if atoms > 0.0: - atoms_per_barn_cm = 1e-24 * atoms / mat.volume - mat.remove_nuclide(nuc) # Replace if it's there - mat.add_nuclide(nuc, atoms_per_barn_cm) - - return mat_file diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py new file mode 100644 index 000000000..6e2373c45 --- /dev/null +++ b/openmc/deplete/stepresult.py @@ -0,0 +1,527 @@ +"""The stepresult module. + +Contains capabilities for generating and saving results of a single depletion +timestep. +""" + +from collections import OrderedDict +import copy + +import h5py +import numpy as np + +import openmc +from openmc.mpi import comm, MPI +from .reaction_rates import ReactionRates + +VERSION_RESULTS = (1, 1) + + +__all__ = ["StepResult"] + + +class StepResult: + """Result of a single depletion timestep + + Attributes + ---------- + k : list of (float, float) + Eigenvalue and uncertainty for each substep. + time : list of float + Time at beginning, end of step, in seconds. + source_rate : float + Source rate during timestep in [W] or [neutron/sec] + n_mat : int + Number of mats. + n_nuc : int + Number of nuclides. + rates : list of ReactionRates + The reaction rates for each substep. + volume : OrderedDict of str to float + Dictionary mapping mat id to volume. + mat_to_ind : OrderedDict of str to int + A dictionary mapping mat ID as string to index. + nuc_to_ind : OrderedDict of str to int + A dictionary mapping nuclide name as string to index. + mat_to_hdf5_ind : OrderedDict of str to int + A dictionary mapping mat ID as string to global index. + n_hdf5_mats : int + Number of materials in entire geometry. + n_stages : int + Number of stages in simulation. + data : numpy.ndarray + Atom quantity, stored by stage, mat, then by nuclide. + proc_time : int + Average time spent depleting a material across all + materials and processes + + """ + def __init__(self): + self.k = None + self.time = None + self.source_rate = None + self.rates = None + self.volume = None + self.proc_time = None + + self.mat_to_ind = None + self.nuc_to_ind = None + self.mat_to_hdf5_ind = None + + self.data = None + + def __repr__(self): + t = self.time[0] + dt = self.time[1] - self.time[0] + return f"" + + def __getitem__(self, pos): + """Retrieves an item from results. + + Parameters + ---------- + pos : tuple + A three-length tuple containing a stage index, mat index and a nuc + index. All can be integers or slices. The second two can be + strings corresponding to their respective dictionary. + + Returns + ------- + float + The atoms for stage, mat, nuc + + """ + stage, mat, nuc = pos + if isinstance(mat, openmc.Material): + mat = str(mat.id) + if isinstance(mat, str): + mat = self.mat_to_ind[mat] + if isinstance(nuc, str): + nuc = self.nuc_to_ind[nuc] + + return self.data[stage, mat, nuc] + + def __setitem__(self, pos, val): + """Sets an item from results. + + Parameters + ---------- + pos : tuple + A three-length tuple containing a stage index, mat index and a nuc + index. All can be integers or slices. The second two can be + strings corresponding to their respective dictionary. + + val : float + The value to set data to. + + """ + stage, mat, nuc = pos + if isinstance(mat, str): + mat = self.mat_to_ind[mat] + if isinstance(nuc, str): + nuc = self.nuc_to_ind[nuc] + + self.data[stage, mat, nuc] = val + + @property + def n_mat(self): + return len(self.mat_to_ind) + + @property + def n_nuc(self): + return len(self.nuc_to_ind) + + @property + def n_hdf5_mats(self): + return len(self.mat_to_hdf5_ind) + + @property + def n_stages(self): + return self.data.shape[0] + + def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages): + """Allocate memory for depletion step data + + Parameters + ---------- + volume : dict of str float + Volumes corresponding to materials in full_burn_dict + nuc_list : list of str + A list of all nuclide names. Used for sorting the simulation. + burn_list : list of int + A list of all mat IDs to be burned. Used for sorting the simulation. + full_burn_list : list of str + List of all burnable material IDs + stages : int + Number of stages in simulation. + + """ + self.volume = copy.deepcopy(volume) + self.nuc_to_ind = {nuc: i for i, nuc in enumerate(nuc_list)} + self.mat_to_ind = {mat: i for i, mat in enumerate(burn_list)} + self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)} + + # Create storage array + self.data = np.zeros((stages, self.n_mat, self.n_nuc)) + + def distribute(self, local_materials, ranges): + """Create a new object containing data for distributed materials + + Parameters + ---------- + local_materials : iterable of str + Materials for this process + ranges : iterable of int + Slice-like object indicating indicies of ``local_materials`` + in the material dimension of :attr:`data` and each element + in :attr:`rates` + + Returns + ------- + StepResult + New results object + """ + new = StepResult() + new.volume = {lm: self.volume[lm] for lm in local_materials} + new.mat_to_ind = {mat: idx for (idx, mat) in enumerate(local_materials)} + + # Direct transfer + direct_attrs = ("time", "k", "source_rate", "nuc_to_ind", + "mat_to_hdf5_ind", "proc_time") + for attr in direct_attrs: + setattr(new, attr, getattr(self, attr)) + # Get applicable slice of data + new.data = self.data[:, ranges] + new.rates = [r[ranges] for r in self.rates] + return new + + def export_to_hdf5(self, filename, step): + """Export results to an HDF5 file + + Parameters + ---------- + filename : str + The filename to write to + step : int + What step is this? + + """ + # Write new file if first time step, else add to existing file + kwargs = {'mode': "w" if step == 0 else "a"} + + if h5py.get_config().mpi and comm.size > 1: + # Write results in parallel + kwargs['driver'] = 'mpio' + kwargs['comm'] = comm + with h5py.File(filename, **kwargs) as handle: + self._to_hdf5(handle, step, parallel=True) + else: + # Gather results at root process + all_results = comm.gather(self) + + # Only root process writes results + if comm.rank == 0: + with h5py.File(filename, **kwargs) as handle: + for res in all_results: + res._to_hdf5(handle, step, parallel=False) + + def _write_hdf5_metadata(self, handle): + """Writes result metadata in HDF5 file + + Parameters + ---------- + handle : h5py.File or h5py.Group + An hdf5 file or group type to store this in. + + """ + # Create and save the 5 dictionaries: + # quantities + # self.mat_to_ind -> self.volume (TODO: support for changing volumes) + # self.nuc_to_ind + # reactions + # self.rates[0].nuc_to_ind (can be different from above, above is superset) + # self.rates[0].react_to_ind + # these are shared by every step of the simulation, and should be deduplicated. + + # Store concentration mat and nuclide dictionaries (along with volumes) + + handle.attrs['version'] = np.array(VERSION_RESULTS) + handle.attrs['filetype'] = np.string_('depletion results') + + mat_list = sorted(self.mat_to_hdf5_ind, key=int) + nuc_list = sorted(self.nuc_to_ind) + rxn_list = sorted(self.rates[0].index_rx) + + n_mats = self.n_hdf5_mats + n_nuc_number = len(nuc_list) + n_nuc_rxn = len(self.rates[0].index_nuc) + n_rxn = len(rxn_list) + n_stages = self.n_stages + + mat_group = handle.create_group("materials") + + for mat in mat_list: + mat_single_group = mat_group.create_group(mat) + mat_single_group.attrs["index"] = self.mat_to_hdf5_ind[mat] + mat_single_group.attrs["volume"] = self.volume[mat] + + nuc_group = handle.create_group("nuclides") + + for nuc in nuc_list: + nuc_single_group = nuc_group.create_group(nuc) + nuc_single_group.attrs["atom number index"] = self.nuc_to_ind[nuc] + if nuc in self.rates[0].index_nuc: + nuc_single_group.attrs["reaction rate index"] = self.rates[0].index_nuc[nuc] + + rxn_group = handle.create_group("reactions") + + for rxn in rxn_list: + rxn_single_group = rxn_group.create_group(rxn) + rxn_single_group.attrs["index"] = self.rates[0].index_rx[rxn] + + # Construct array storage + + handle.create_dataset("number", (1, n_stages, n_mats, n_nuc_number), + maxshape=(None, n_stages, n_mats, n_nuc_number), + chunks=(1, 1, n_mats, n_nuc_number), + dtype='float64') + + handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn), + maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn), + chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn), + dtype='float64') + + handle.create_dataset("eigenvalues", (1, n_stages, 2), + maxshape=(None, n_stages, 2), dtype='float64') + + handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64') + + handle.create_dataset("source_rate", (1, n_stages), maxshape=(None, n_stages), + dtype='float64') + + handle.create_dataset( + "depletion time", (1,), maxshape=(None,), + dtype="float64") + + def _to_hdf5(self, handle, index, parallel=False): + """Converts results object into an hdf5 object. + + Parameters + ---------- + handle : h5py.File or h5py.Group + An HDF5 file or group type to store this in. + index : int + What step is this? + parallel : bool + Being called with parallel HDF5? + + """ + if "/number" not in handle: + if parallel: + comm.barrier() + self._write_hdf5_metadata(handle) + + if parallel: + comm.barrier() + + # Grab handles + number_dset = handle["/number"] + rxn_dset = handle["/reaction rates"] + eigenvalues_dset = handle["/eigenvalues"] + time_dset = handle["/time"] + source_rate_dset = handle["/source_rate"] + proc_time_dset = handle["/depletion time"] + + # Get number of results stored + number_shape = list(number_dset.shape) + number_results = number_shape[0] + + new_shape = index + 1 + + if number_results < new_shape: + # Extend first dimension by 1 + number_shape[0] = new_shape + number_dset.resize(number_shape) + + rxn_shape = list(rxn_dset.shape) + rxn_shape[0] = new_shape + rxn_dset.resize(rxn_shape) + + eigenvalues_shape = list(eigenvalues_dset.shape) + eigenvalues_shape[0] = new_shape + eigenvalues_dset.resize(eigenvalues_shape) + + time_shape = list(time_dset.shape) + time_shape[0] = new_shape + time_dset.resize(time_shape) + + source_rate_shape = list(source_rate_dset.shape) + source_rate_shape[0] = new_shape + source_rate_dset.resize(source_rate_shape) + + proc_shape = list(proc_time_dset.shape) + proc_shape[0] = new_shape + proc_time_dset.resize(proc_shape) + + # If nothing to write, just return + if len(self.mat_to_ind) == 0: + return + + # Add data + # Note, for the last step, self.n_stages = 1, even if n_stages != 1. + n_stages = self.n_stages + inds = [self.mat_to_hdf5_ind[mat] for mat in self.mat_to_ind] + low = min(inds) + high = max(inds) + for i in range(n_stages): + number_dset[index, i, low:high+1] = self.data[i] + rxn_dset[index, i, low:high+1] = self.rates[i] + if comm.rank == 0: + eigenvalues_dset[index, i] = self.k[i] + if comm.rank == 0: + time_dset[index] = self.time + source_rate_dset[index] = self.source_rate + if self.proc_time is not None: + proc_time_dset[index] = ( + self.proc_time / (comm.size * self.n_hdf5_mats) + ) + + @classmethod + def from_hdf5(cls, handle, step): + """Loads results object from HDF5. + + Parameters + ---------- + handle : h5py.File or h5py.Group + An HDF5 file or group type to load from. + step : int + Index for depletion step + """ + results = cls() + + # Grab handles + number_dset = handle["/number"] + eigenvalues_dset = handle["/eigenvalues"] + time_dset = handle["/time"] + if "source_rate" in handle: + source_rate_dset = handle["/source_rate"] + else: + # Older versions used "power" instead of "source_rate" + source_rate_dset = handle["/power"] + + results.data = number_dset[step, :, :, :] + results.k = eigenvalues_dset[step, :] + results.time = time_dset[step, :] + results.source_rate = source_rate_dset[step, 0] + + if "depletion time" in handle: + proc_time_dset = handle["/depletion time"] + if step < proc_time_dset.shape[0]: + results.proc_time = proc_time_dset[step] + + if results.proc_time is None: + results.proc_time = np.array([np.nan]) + + # Reconstruct dictionaries + results.volume = OrderedDict() + results.mat_to_ind = OrderedDict() + results.nuc_to_ind = OrderedDict() + rxn_nuc_to_ind = OrderedDict() + rxn_to_ind = OrderedDict() + + for mat, mat_handle in handle["/materials"].items(): + vol = mat_handle.attrs["volume"] + ind = mat_handle.attrs["index"] + + results.volume[mat] = vol + results.mat_to_ind[mat] = ind + + for nuc, nuc_handle in handle["/nuclides"].items(): + ind_atom = nuc_handle.attrs["atom number index"] + results.nuc_to_ind[nuc] = ind_atom + + if "reaction rate index" in nuc_handle.attrs: + rxn_nuc_to_ind[nuc] = nuc_handle.attrs["reaction rate index"] + + for rxn, rxn_handle in handle["/reactions"].items(): + rxn_to_ind[rxn] = rxn_handle.attrs["index"] + + results.rates = [] + # Reconstruct reactions + for i in range(results.n_stages): + rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind, True) + + rate[:] = handle["/reaction rates"][step, i, :, :, :] + results.rates.append(rate) + + return results + + @staticmethod + def save(op, x, op_results, t, source_rate, step_ind, proc_time=None): + """Creates and writes depletion results to disk + + Parameters + ---------- + op : openmc.deplete.TransportOperator + The operator used to generate these results. + x : list of list of numpy.array + The prior x vectors. Indexed [i][cell] using the above equation. + op_results : list of openmc.deplete.OperatorResult + Results of applying transport operator + t : list of float + Time indices. + source_rate : float + Source rate during time step in [W] or [neutron/sec] + step_ind : int + Step index. + proc_time : float or None + Total process time spent depleting materials. This may + be process-dependent and will be reduced across MPI + processes. + + """ + # Get indexing terms + vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info() + + stages = len(x) + + # Create results + results = StepResult() + results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages) + + n_mat = len(burn_list) + + for i in range(stages): + for mat_i in range(n_mat): + results[i, mat_i, :] = x[i][mat_i] + + results.k = [(r.k.nominal_value, r.k.std_dev) for r in op_results] + results.rates = [r.rates for r in op_results] + results.time = t + results.source_rate = source_rate + results.proc_time = proc_time + if results.proc_time is not None: + results.proc_time = comm.reduce(proc_time, op=MPI.SUM) + + results.export_to_hdf5("depletion_results.h5", step_ind) + + def transfer_volumes(self, model): + """Transfers volumes from depletion results to geometry + + Parameters + ---------- + model : OpenMC model to be used in a depletion restart + calculation + + """ + + if not model.materials: + materials = openmc.Materials( + model.geometry.get_all_materials().values() + ) + else: + materials = model.materials + + for material in materials: + if material.depletable: + material.volume = self.volume[str(material.id)] diff --git a/tests/regression_tests/deplete/test.py b/tests/regression_tests/deplete/test.py index 4f0e738a3..6c431fb33 100644 --- a/tests/regression_tests/deplete/test.py +++ b/tests/regression_tests/deplete/test.py @@ -77,8 +77,8 @@ def test_full(run_in_tmpdir, problem, multiproc): return # Load the reference/test results - res_test = openmc.deplete.ResultsList.from_hdf5(path_test) - res_ref = openmc.deplete.ResultsList.from_hdf5(path_reference) + res_test = openmc.deplete.Results(path_test) + res_ref = openmc.deplete.Results(path_reference) # Assert same mats for mat in res_ref[0].mat_to_ind: @@ -139,7 +139,7 @@ def test_depletion_results_to_material(run_in_tmpdir, problem): """Checks openmc.Materials objects can be created from depletion results""" # Load the reference/test results path_reference = Path(__file__).with_name('test_reference.h5') - res_ref = openmc.deplete.ResultsList.from_hdf5(path_reference) + res_ref = openmc.deplete.Results(path_reference) # Firstly need to export materials.xml file for the initial simulation state geometry, lower_left, upper_right = problem diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index 2f02ec4b8..0b82a5fbc 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -107,7 +107,7 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts integrator.integrate() # Get resulting number of atoms - results = openmc.deplete.ResultsList.from_hdf5('depletion_results.h5') + results = openmc.deplete.Results('depletion_results.h5') _, atoms = results.get_atoms(w, "W186") assert atoms[0] == pytest.approx(n0) @@ -155,7 +155,7 @@ def test_decay(run_in_tmpdir): integrator.integrate() # Get resulting number of atoms - results = openmc.deplete.ResultsList.from_hdf5('depletion_results.h5') + results = openmc.deplete.Results('depletion_results.h5') _, atoms = results.get_atoms(mat, "Sr89") # Ensure density goes down by a factor of 2 after each half-life diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index a08d8738c..cb6374152 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -16,7 +16,7 @@ import pytest from openmc.mpi import comm from openmc.deplete import ( - ReactionRates, Results, ResultsList, OperatorResult, PredictorIntegrator, + ReactionRates, StepResult, Results, OperatorResult, PredictorIntegrator, CECMIntegrator, CF4Integrator, CELIIntegrator, EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator, cram) @@ -99,11 +99,11 @@ def test_results_save(run_in_tmpdir): for k, rates in zip(eigvl1, rate1)] op_result2 = [OperatorResult(ufloat(*k), rates) for k, rates in zip(eigvl2, rate2)] - Results.save(op, x1, op_result1, t1, 0, 0) - Results.save(op, x2, op_result2, t2, 0, 1) + StepResult.save(op, x1, op_result1, t1, 0, 0) + StepResult.save(op, x2, op_result2, t2, 0, 1) # Load the files - res = ResultsList.from_hdf5("depletion_results.h5") + res = Results("depletion_results.h5") for i in range(stages): for mat_i, mat in enumerate(burn_list): @@ -176,8 +176,7 @@ def test_integrator(run_in_tmpdir, scheme): # get expected results - res = ResultsList.from_hdf5( - operator.output_dir / "depletion_results.h5") + res = Results(operator.output_dir / "depletion_results.h5") t1, y1 = res.get_atoms("1", "1") t2, y2 = res.get_atoms("1", "2") diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py index 82ef05ac7..e8bfc062a 100644 --- a/tests/unit_tests/test_deplete_restart.py +++ b/tests/unit_tests/test_deplete_restart.py @@ -24,8 +24,7 @@ def test_restart_predictor_cecm(run_in_tmpdir): openmc.deplete.PredictorIntegrator(op, dt, power).integrate() # Load the files - prev_res = openmc.deplete.ResultsList.from_hdf5( - op.output_dir / "depletion_results.h5") + prev_res = openmc.deplete.Results(op.output_dir / "depletion_results.h5") # Re-create depletion operator and load previous results op = dummy_operator.DummyOperator(prev_res) @@ -51,8 +50,7 @@ def test_restart_cecm_predictor(run_in_tmpdir): cecm.integrate() # Load the files - prev_res = openmc.deplete.ResultsList.from_hdf5( - op.output_dir / "depletion_results.h5") + prev_res = openmc.deplete.Results(op.output_dir / "depletion_results.h5") # Re-create depletion operator and load previous results op = dummy_operator.DummyOperator(prev_res) @@ -75,7 +73,7 @@ def test_restart(run_in_tmpdir, scheme): bundle.solver(operator, [0.75], 1.0).integrate() # restart - prev_res = openmc.deplete.ResultsList.from_hdf5( + prev_res = openmc.deplete.Results( operator.output_dir / "depletion_results.h5") operator = dummy_operator.DummyOperator(prev_res) @@ -84,7 +82,7 @@ def test_restart(run_in_tmpdir, scheme): # compare results - results = openmc.deplete.ResultsList.from_hdf5( + results = openmc.deplete.Results( operator.output_dir / "depletion_results.h5") _t, y1 = results.get_atoms("1", "1") diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index 03e7c3eb0..6ce15601f 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -1,4 +1,4 @@ -"""Tests the ResultsList class""" +"""Tests the Results class""" from pathlib import Path from math import inf @@ -13,7 +13,7 @@ def res(): """Load the reference results""" filename = (Path(__file__).parents[1] / 'regression_tests' / 'deplete' / 'test_reference.h5') - return openmc.deplete.ResultsList.from_hdf5(filename) + return openmc.deplete.Results(filename) def test_get_atoms(res): @@ -72,9 +72,9 @@ def test_get_keff(res): @pytest.mark.parametrize("unit", ("s", "d", "min", "h")) def test_get_steps(unit): - # Make a ResultsList full of near-empty Result instances + # Make a Results full of near-empty Result instances # Just fill out a time schedule - results = openmc.deplete.ResultsList() + results = openmc.deplete.Results() # Time in units of unit times = np.linspace(0, 100, num=5) if unit == "d": @@ -87,7 +87,7 @@ def test_get_steps(unit): conversion_to_seconds = 1 for ix in range(times.size): - res = openmc.deplete.Results() + res = openmc.deplete.StepResult() res.time = times[ix:ix + 1] * conversion_to_seconds results.append(res) diff --git a/tests/unit_tests/test_transfer_volumes.py b/tests/unit_tests/test_transfer_volumes.py index e1f3c1053..b768bc88a 100644 --- a/tests/unit_tests/test_transfer_volumes.py +++ b/tests/unit_tests/test_transfer_volumes.py @@ -2,7 +2,7 @@ from pytest import approx import openmc -from openmc.deplete import PredictorIntegrator, ResultsList +from openmc.deplete import PredictorIntegrator, Results from tests import dummy_operator @@ -19,7 +19,7 @@ def test_transfer_volumes(run_in_tmpdir): PredictorIntegrator(op, dt, power).integrate() # Load the files - res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5") + res = openmc.deplete.Results(op.output_dir / "depletion_results.h5") # Create a dictionary of volumes to transfer res[0].volume['1'] = 1.5