Source code for openmc.deplete.results

"""The results module.

Contains results generation and saving capabilities.

from collections import OrderedDict
import copy
from warnings import warn

import numpy as np
import h5py

from . import comm, have_mpi, MPI
from .reaction_rates import ReactionRates


__all__ = ["Results"]

[docs]class Results(object): """Output of a depletion run Attributes ---------- k : list of (float, float) Eigenvalue and uncertainty for each substep. time : list of float Time at beginning, end of step, in seconds. power : float Power during time step, in Watts n_mat : int Number of mats. n_nuc : int Number of nuclides. rates : list of ReactionRates The reaction rates for each substep. volume : OrderedDict of int to float Dictionary mapping mat id to volume. mat_to_ind : OrderedDict of str to int A dictionary mapping mat ID as string to index. nuc_to_ind : OrderedDict of str to int A dictionary mapping nuclide name as string to index. mat_to_hdf5_ind : OrderedDict of str to int A dictionary mapping mat ID as string to global index. n_hdf5_mats : int Number of materials in entire geometry. n_stages : int Number of stages in simulation. data : numpy.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.power = 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 = None def __getitem__(self, pos): """Retrieves an item from results. Parameters ---------- pos : tuple A three-length tuple containing a stage index, mat index and a nuc index. All can be integers or slices. The second two can be strings corresponding to their respective dictionary. Returns ------- float The atoms for stage, mat, nuc """ stage, mat, nuc = pos if isinstance(mat, str): mat = self.mat_to_ind[mat] if isinstance(nuc, str): nuc = self.nuc_to_ind[nuc] return[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][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[0]
[docs] 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 = np.zeros((stages, self.n_mat, self.n_nuc))
[docs] 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 = dict(zip( local_materials, range(len(local_materials)))) # Direct transfer direct_attrs = ("time", "k", "power", "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 =[:, ranges] new.rates = [r[ranges] for r in self.rates] return new
[docs] 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? """ if have_mpi and h5py.get_config().mpi: kwargs = {'driver': 'mpio', 'comm': comm} else: kwargs = {} # Write new file if first time step, else add to existing file kwargs['mode'] = "w" if step == 0 else "a" with h5py.File(filename, **kwargs) as handle: self._to_hdf5(handle, step)
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("power", (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): """Converts results object into an hdf5 object. Parameters ---------- handle : h5py.File or h5py.Group An HDF5 file or group type to store this in. index : int What step is this? """ if "/number" not in handle: comm.barrier() self._write_hdf5_metadata(handle) comm.barrier() # Grab handles number_dset = handle["/number"] rxn_dset = handle["/reaction rates"] eigenvalues_dset = handle["/eigenvalues"] time_dset = handle["/time"] power_dset = handle["/power"] 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) power_shape = list(power_dset.shape) power_shape[0] = new_shape power_dset.resize(power_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, :] =[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 power_dset[index, :] = self.power if self.proc_time is not None: proc_time_dset[index] = ( self.proc_time / (comm.size * self.n_hdf5_mats) )
[docs] @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"] power_dset = handle["/power"] = number_dset[step, :, :, :] results.k = eigenvalues_dset[step, :] results.time = time_dset[step, :] results.power = power_dset[step, :] 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
[docs] @staticmethod def save(op, x, op_results, t, power, 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. power : float Power during time step 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 = Results() 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.power = power 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)
[docs] def transfer_volumes(self, geometry): """Transfers volumes from depletion results to geometry Parameters ---------- geometry : OpenMC geometry to be used in a depletion restart calculation """ for cell in geometry.get_all_material_cells().values(): for material in cell.get_all_materials().values(): if material.depletable: material.volume = self.volume[str(]