From f5db4e3afb41ff1b56639ccf157e0256b0b5a19f Mon Sep 17 00:00:00 2001 From: yardasol Date: Tue, 12 Jul 2022 18:30:29 -0500 Subject: [PATCH] move reuasble code to new parent class; 1st pass - new class OpenMCOperator to contain shared bits of code between the currently existing Operator class and the as-of-yet to be created FluxDepletionOperator class - Operator is now a subclass of OpenMCOperator - move _distribute method to OpenMCOperator - move _get_burnable_mats, _extract_number, _set_number_from_mat, _set_number_from_results to OpenMCOperator - split initial condition into openmc.lib dependent and non-dependent parts; dependent part goes to Operator, non-dependent part goes to OpenMCOperator - move _update_materials to OpenMCOperator - move _get_tally_nuclides to OpenMCOperator - rename _unpack_tallies_and_normalize to _calculate_reaction_rates; move to OpenMCOperator - move get_results_info to OpenMCOperator --- openmc/deplete/__init__.py | 1 + openmc/deplete/openmc_operator.py | 522 ++++++++++++++++++++++++++++++ openmc/deplete/operator.py | 331 +------------------ 3 files changed, 535 insertions(+), 319 deletions(-) create mode 100644 openmc/deplete/openmc_operator.py diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 5891555a2..95264308c 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -7,6 +7,7 @@ A depletion front-end tool. from .nuclide import * from .chain import * +from .openmc_operator import * from .operator import * from .reaction_rates import * from .atom_number import * diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py new file mode 100644 index 000000000..303a5208f --- /dev/null +++ b/openmc/deplete/openmc_operator.py @@ -0,0 +1,522 @@ +"""OpenMC transport operator + +This module implements a transport operator for OpenMC so that it can be used by +depletion integrators. The implementation makes use of the Python bindings to +OpenMC's C API so that reading tally results and updating material number +densities is all done in-memory instead of through the filesystem. + +""" + +import copy +from abc import abstractmethod +from collections import OrderedDict +import os +from warnings import warn + +import numpy as np + +import openmc +from openmc.exceptions import DataError +from openmc.mpi import comm +from .abc import TransportOperator, OperatorResult +from .atom_number import AtomNumber +from .reaction_rates import ReactionRates +from .results import Results +from .helpers import ( + DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, + FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, + SourceRateHelper, FluxCollapseHelper) + + +__all__ = ["OpenMCOperator", "OperatorResult"] + + +def _distribute(items): + """Distribute items across MPI communicator + + Parameters + ---------- + items : list + List of items of distribute + + Returns + ------- + list + Items assigned to process that called + + """ + min_size, extra = divmod(len(items), comm.size) + j = 0 + for i in range(comm.size): + chunk_size = min_size + int(i < extra) + if comm.rank == i: + return items[j:j + chunk_size] + j += chunk_size + +class OpenMCOperator(TransportOperator): + """OpenMC transport operator for depletion. + + Instances of this class can be used to perform depletion using OpenMC as the + transport operator. Normally, a user needn't call methods of this class + directly. Instead, an instance of this class is passed to an integrator + class, such as :class:`openmc.deplete.CECMIntegrator`. + + .. versionchanged:: 0.13.0 + The geometry and settings parameters have been replaced with a + model parameter that takes a :class:`~openmc.model.Model` object + + Parameters + ---------- + chain_file : str, optional + Path to the depletion chain XML file. Defaults to the file + listed under ``depletion_chain`` in + :envvar:`OPENMC_CROSS_SECTIONS` environment variable. + prev_results : Results, optional + Results from a previous depletion calculation. If this argument is + specified, the depletion calculation will start from the latest state + in the previous results. + fission_q : dict, optional + Dictionary of nuclides and their fission Q values [eV]. If not given, + values will be pulled from the ``chain_file``. Only applicable + if ``"normalization_mode" == "fission-q"`` + dilute_initial : float, optional + Initial atom density [atoms/cm^3] to add for nuclides that are zero + in initial condition to ensure they exist in the decay chain. + Only done for nuclides with reaction rates. + Defaults to 1.0e3. + + Attributes + ---------- + model : openmc.model.Model + OpenMC model object + geometry : openmc.Geometry + OpenMC geometry object + settings : openmc.Settings + OpenMC settings object + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the decay + chain. Only done for nuclides with reaction rates. + output_dir : pathlib.Path + Path to output directory to save results. + round_number : bool + Whether or not to round output to OpenMC to 8 digits. + Useful in testing, as OpenMC is incredibly sensitive to exact values. + number : openmc.deplete.AtomNumber + Total number of atoms in simulation. + nuclides_with_data : set of str + A set listing all unique nuclides available from cross_sections.xml. + chain : openmc.deplete.Chain + The depletion chain information necessary to form matrices and tallies. + reaction_rates : openmc.deplete.ReactionRates + Reaction rates from the last operator step. + burnable_mats : list of str + All burnable material IDs + heavy_metal : float + Initial heavy metal inventory [g] + local_mats : list of str + All burnable material IDs being managed by a single process + prev_res : Results or None + Results from a previous depletion calculation. ``None`` if no + results are to be used. + cleanup_when_done : bool + Whether to finalize and clear the shared library memory when the + depletion operation is complete. Defaults to clearing the library. + """ + + ## modify + def __init__(self, chain_file=None, fission_q=None, dilute_initial=0.0, prev_results=None): + + super().__init__(chain_file, fission_q, dilute_initial, prev_results) + + ## clear, but we must keep the method + def __call__(self, vec, source_rate): + """Runs a simulation. + + Simulation will abort under the following circumstances: + + 1) No energy is computed using OpenMC tallies. + + Parameters + ---------- + vec : list of numpy.ndarray + Total atoms to be used in function. + source_rate : float + Power in [W] or source rate in [neutron/sec] + + Returns + ------- + openmc.deplete.OperatorResult + Eigenvalue and reaction rates resulting from transport operator + + """ + + ## clear, but must keep the method + def write_bos_data(step): + """Write a state-point file with beginning of step data + + Parameters + ---------- + step : int + Current depletion step including restarts + """ + pass + + ## keep + def _get_burnable_mats(self): + """Determine depletable materials, volumes, and nuclides + + Returns + ------- + burnable_mats : list of str + List of burnable material IDs + volume : OrderedDict of str to float + Volume of each material in [cm^3] + nuclides : list of str + Nuclides in order of how they'll appear in the simulation. + + """ + + burnable_mats = set() + model_nuclides = set() + volume = OrderedDict() + + self.heavy_metal = 0.0 + + # Iterate once through the geometry to get dictionaries + for mat in self.materials: + for nuclide in mat.get_nuclides(): + model_nuclides.add(nuclide) + if mat.depletable: + burnable_mats.add(str(mat.id)) + if mat.volume is None: + raise RuntimeError("Volume not specified for depletable " + "material with ID={}.".format(mat.id)) + volume[str(mat.id)] = mat.volume + self.heavy_metal += mat.fissionable_mass + + # Make sure there are burnable materials + if not burnable_mats: + raise RuntimeError( + "No depletable materials were found in the model.") + + # Sort the sets + burnable_mats = sorted(burnable_mats, key=int) + model_nuclides = sorted(model_nuclides) + + # Construct a global nuclide dictionary, burned first + nuclides = list(self.chain.nuclide_dict) + for nuc in model_nuclides: + if nuc not in nuclides: + nuclides.append(nuc) + + return burnable_mats, volume, nuclides + + ## keep + def _extract_number(self, local_mats, volume, nuclides, prev_res=None): + """Construct AtomNumber using geometry + + Parameters + ---------- + local_mats : list of str + Material IDs to be managed by this process + volume : OrderedDict of str to float + Volumes for the above materials in [cm^3] + nuclides : list of str + Nuclides to be used in the simulation. + prev_res : Results, optional + Results from a previous depletion calculation + + """ + self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain)) + + if self.dilute_initial != 0.0: + for nuc in self._burnable_nucs: + self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) + + # Now extract and store the number densities + # From the geometry if no previous depletion results + if prev_res is None: + for mat in self.materials: + if str(mat.id) in local_mats: + self._set_number_from_mat(mat) + + # Else from previous depletion results + else: + for mat in self.materials: + if str(mat.id) in local_mats: + self._set_number_from_results(mat, prev_res) + + def _set_number_from_mat(self, mat): + """Extracts material and number densities from openmc.Material + + Parameters + ---------- + mat : openmc.Material + The material to read from + + """ + mat_id = str(mat.id) + + for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): + atom_per_cc = atom_per_bcm * 1.0e24 + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + + def _set_number_from_results(self, mat, prev_res): + """Extracts material nuclides and number densities. + + If the nuclide concentration's evolution is tracked, the densities come + from depletion results. Else, densities are extracted from the geometry + in the summary. + + Parameters + ---------- + mat : openmc.Material + The material to read from + prev_res : Results + Results from a previous depletion calculation + + """ + mat_id = str(mat.id) + + # Get nuclide lists from geometry and depletion results + depl_nuc = prev_res[-1].nuc_to_ind + geom_nuc_densities = mat.get_nuclide_atom_densities() + + # Merge lists of nuclides, with the same order for every calculation + geom_nuc_densities.update(depl_nuc) + + for nuclide, atom_per_bcm in geom_nuc_densities.items(): + if nuclide in depl_nuc: + concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] + volume = prev_res[-1].volume[mat_id] + atom_per_cc = concentration / volume + else: + atom_per_cc = atom_per_bcm * 1.0e24 + + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + + ## keep, will need to modify + def initial_condition(self, materials): + """Performs final setup and returns initial condition. + + Parameters + ---------- + materials : ??? + + Returns + ------- + list of numpy.ndarray + Total density for initial conditions. + """ + + self._rate_helper.generate_tallies(materials, self.chain.reactions) + self._normalization_helper.prepare( + self.chain.nuclides, self.reaction_rates.index_nuc) + # Tell fission yield helper what materials this process is + # responsible for + self._yield_helper.generate_tallies( + materials, tuple(sorted(self._mat_index_map.values()))) + + # Return number density vector + return list(self.number.get_mat_slice(np.s_[:])) + + ## keep? there is a non-removable part of this + ## that needs the openmc c executable + def _update_materials(self): + """Updates material compositions in OpenMC on all processes.""" + + for rank in range(comm.size): + number_i = comm.bcast(self.number, root=rank) + + for mat in number_i.materials: + nuclides = [] + densities = [] + for nuc in number_i.nuclides: + if nuc in self.nuclides_with_data: + val = 1.0e-24 * number_i.get_atom_density(mat, nuc) + + # If nuclide is zero, do not add to the problem. + if val > 0.0: + if self.round_number: + val_magnitude = np.floor(np.log10(val)) + val_scaled = val / 10**val_magnitude + val_round = round(val_scaled, 8) + + val = val_round * 10**val_magnitude + + nuclides.append(nuc) + densities.append(val) + else: + # Only output warnings if values are significantly + # negative. CRAM does not guarantee positive values. + if val < -1.0e-21: + print("WARNING: nuclide ", nuc, " in material ", mat, + " is negative (density = ", val, " at/barn-cm)") + number_i[mat, nuc] = 0.0 + + # Update densities on C API side + mat_internal = openmc.lib.materials[int(mat)] + mat_internal.set_densities(nuclides, densities) + + #TODO Update densities on the Python side, otherwise the + # summary.h5 file contains densities at the first time step + + ## keep, but rename and modify docstring + ## get_tally_nuclides + def _get_tally_nuclides(self): + """Determine nuclides that should be tallied for reaction rates. + + This method returns a list of all nuclides that have neutron data and + are listed in the depletion chain. Technically, we should tally nuclides + that may not appear in the depletion chain because we still need to get + the fission reaction rate for these nuclides in order to normalize + power, but that is left as a future exercise. + + Returns + ------- + list of str + Tally nuclides + + """ + nuc_set = set() + + # Create the set of all nuclides in the decay chain in materials marked + # for burning in which the number density is greater than zero. + for nuc in self.number.nuclides: + if nuc in self.nuclides_with_data: + if np.sum(self.number[:, nuc]) > 0.0: + nuc_set.add(nuc) + + # Communicate which nuclides have nonzeros to rank 0 + if comm.rank == 0: + for i in range(1, comm.size): + nuc_newset = comm.recv(source=i, tag=i) + nuc_set |= nuc_newset + else: + comm.send(nuc_set, dest=0, tag=comm.rank) + + if comm.rank == 0: + # Sort nuclides in the same order as self.number + nuc_list = [nuc for nuc in self.number.nuclides + if nuc in nuc_set] + else: + nuc_list = None + + # Store list of tally nuclides on each process + nuc_list = comm.bcast(nuc_list) + return [nuc for nuc in nuc_list if nuc in self.chain] + + ## modify + ## the tally-specific methods + ## for the normalization and fission helper classes + ## will already work with the current implementation. + ## For transport-less depletion, we'll need to write + ## a new ReactionRateHelper + def _calculate_reaction_rates(self, source_rate): + """Unpack tallies from OpenMC and return an operator result + + This method uses OpenMC's C API bindings to determine the k-effective + value and reaction rates from the simulation. The reaction rates are + normalized by a helper class depending on the method being used. + + Parameters + ---------- + source_rate : float + Power in [W] or source rate in [neutron/sec] + + Returns + ------- + rates : openmc.deplete.ReactionRates + Reaction rates for nuclides + + """ + rates = self.reaction_rates + rates.fill(0.0) + + # Extract tally bins + nuclides = self._rate_helper.nuclides + + # Form fast map + nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] + react_ind = [rates.index_rx[react] for react in self.chain.reactions] + + # Keep track of energy produced from all reactions in eV per source + # particle + self._normalization_helper.reset() + self._yield_helper.unpack() + + # Store fission yield dictionaries + fission_yields = [] + + # Create arrays to store fission Q values, reaction rates, and nuclide + # numbers, zeroed out in material iteration + number = np.empty(rates.n_nuc) + + fission_ind = rates.index_rx.get("fission") + + # Extract results + for i, mat in enumerate(self.local_mats): + # Get tally index + mat_index = self._mat_index_map[mat] + + # Zero out reaction rates and nuclide numbers + number.fill(0.0) + + # Get new number densities + for nuc, i_nuc_results in zip(nuclides, nuc_ind): + number[i_nuc_results] = self.number[mat, nuc] + + tally_rates = self._rate_helper.get_material_rates( + mat_index, nuc_ind, react_ind) + + # Compute fission yields for this material + fission_yields.append(self._yield_helper.weighted_yields(i)) + + # Accumulate energy from fission + if fission_ind is not None: + self._normalization_helper.update(tally_rates[:, fission_ind]) + + # Divide by total number and store + rates[i] = self._rate_helper.divide_by_adens(number) + + # Scale reaction rates to obtain units of reactions/sec + rates *= self._normalization_helper.factor(source_rate) + + # Store new fission yields on the chain + self.chain.fission_yields = fission_yields + + return rates + + @abstractmethod + def _get_nuclides_with_data(self): + """Find nuclides with cross section data""" + + ## Keep + def get_results_info(self): + """Returns volume list, material lists, and nuc lists. + + Returns + ------- + volume : dict of str float + Volumes corresponding to materials in full_burn_dict + nuc_list : list of str + A list of all nuclide names. Used for sorting the simulation. + burn_list : list of int + A list of all material IDs to be burned. Used for sorting the simulation. + full_burn_list : list + List of all burnable material IDs + + """ + nuc_list = self.number.burnable_nuclides + burn_list = self.local_mats + + volume = {} + for i, mat in enumerate(burn_list): + volume[mat] = self.number.volume[i] + + # Combine volume dictionaries across processes + volume_list = comm.allgather(volume) + volume = {k: v for d in volume_list for k, v in d.items()} + + return volume, nuc_list, burn_list, self.burnable_mats diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 7ef874e87..665de1999 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -24,6 +24,7 @@ from openmc.mpi import comm from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .chain import _find_chain_file +from .openmc_operator import OpenMCOperator, _distribute from .reaction_rates import ReactionRates from .results import Results from .helpers import ( @@ -34,30 +35,6 @@ from .helpers import ( __all__ = ["Operator", "OperatorResult"] - -def _distribute(items): - """Distribute items across MPI communicator - - Parameters - ---------- - items : list - List of items of distribute - - Returns - ------- - list - Items assigned to process that called - - """ - min_size, extra = divmod(len(items), comm.size) - j = 0 - for i in range(comm.size): - chunk_size = min_size + int(i < extra) - if comm.rank == i: - return items[j:j + chunk_size] - j += chunk_size - - def _find_cross_sections(model): """Determine cross sections to use for depletion""" if model.materials and model.materials.cross_sections is not None: @@ -74,7 +51,7 @@ def _find_cross_sections(model): return cross_sections -class Operator(TransportOperator): +class Operator(OpenMCOperator): """OpenMC transport operator for depletion. Instances of this class can be used to perform depletion using OpenMC as the @@ -267,7 +244,7 @@ class Operator(TransportOperator): # Clear out OpenMC, create task lists, distribute openmc.reset_auto_ids() - self.burnable_mats, volume, nuclides = self._get_burnable_mats() + self.burnable_mats, volume, nuclides = super()._get_burnable_mats() self.local_mats = _distribute(self.burnable_mats) # Generate map from local materials => material index @@ -298,7 +275,7 @@ class Operator(TransportOperator): if nuc.name in self.nuclides_with_data] # Extract number densities from the geometry / previous depletion run - self._extract_number(self.local_mats, volume, nuclides, self.prev_res) + super()._extract_number(self.local_mats, volume, nuclides, self.prev_res) # Create reaction rates array self.reaction_rates = ReactionRates( @@ -380,7 +357,7 @@ class Operator(TransportOperator): openmc.reset_auto_ids() # Update tally nuclides data in preparation for transport solve - nuclides = self._get_tally_nuclides() + nuclides = super()._get_tally_nuclides() self._rate_helper.nuclides = nuclides self._normalization_helper.nuclides = nuclides self._yield_helper.update_tally_nuclides(nuclides) @@ -390,7 +367,12 @@ class Operator(TransportOperator): openmc.lib.reset_timers() # Extract results - op_result = self._unpack_tallies_and_normalize(source_rate) + rates = self._calculate_reaction_rates(source_rate) + + # Get k and uncertainty + keff = ufloat(*openmc.lib.keff()) + + op_result = OperatorResult(keff, rates) return copy.deepcopy(op_result) @@ -434,138 +416,6 @@ class Operator(TransportOperator): cell.fill = [mat.clone() for i in range(cell.num_instances)] - def _get_burnable_mats(self): - """Determine depletable materials, volumes, and nuclides - - Returns - ------- - burnable_mats : list of str - List of burnable material IDs - volume : OrderedDict of str to float - Volume of each material in [cm^3] - nuclides : list of str - Nuclides in order of how they'll appear in the simulation. - - """ - - burnable_mats = set() - model_nuclides = set() - volume = OrderedDict() - - self.heavy_metal = 0.0 - - # Iterate once through the geometry to get dictionaries - for mat in self.materials: - for nuclide in mat.get_nuclides(): - model_nuclides.add(nuclide) - if mat.depletable: - burnable_mats.add(str(mat.id)) - if mat.volume is None: - raise RuntimeError("Volume not specified for depletable " - "material with ID={}.".format(mat.id)) - volume[str(mat.id)] = mat.volume - self.heavy_metal += mat.fissionable_mass - - # Make sure there are burnable materials - if not burnable_mats: - raise RuntimeError( - "No depletable materials were found in the model.") - - # Sort the sets - burnable_mats = sorted(burnable_mats, key=int) - model_nuclides = sorted(model_nuclides) - - # Construct a global nuclide dictionary, burned first - nuclides = list(self.chain.nuclide_dict) - for nuc in model_nuclides: - if nuc not in nuclides: - nuclides.append(nuc) - - return burnable_mats, volume, nuclides - - def _extract_number(self, local_mats, volume, nuclides, prev_res=None): - """Construct AtomNumber using geometry - - Parameters - ---------- - local_mats : list of str - Material IDs to be managed by this process - volume : OrderedDict of str to float - Volumes for the above materials in [cm^3] - nuclides : list of str - Nuclides to be used in the simulation. - prev_res : Results, optional - Results from a previous depletion calculation - - """ - self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain)) - - if self.dilute_initial != 0.0: - for nuc in self._burnable_nucs: - self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) - - # Now extract and store the number densities - # From the geometry if no previous depletion results - if prev_res is None: - for mat in self.materials: - if str(mat.id) in local_mats: - self._set_number_from_mat(mat) - - # Else from previous depletion results - else: - for mat in self.materials: - if str(mat.id) in local_mats: - self._set_number_from_results(mat, prev_res) - - def _set_number_from_mat(self, mat): - """Extracts material and number densities from openmc.Material - - Parameters - ---------- - mat : openmc.Material - The material to read from - - """ - mat_id = str(mat.id) - - for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): - atom_per_cc = atom_per_bcm * 1.0e24 - self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - - def _set_number_from_results(self, mat, prev_res): - """Extracts material nuclides and number densities. - - If the nuclide concentration's evolution is tracked, the densities come - from depletion results. Else, densities are extracted from the geometry - in the summary. - - Parameters - ---------- - mat : openmc.Material - The material to read from - prev_res : Results - Results from a previous depletion calculation - - """ - mat_id = str(mat.id) - - # Get nuclide lists from geometry and depletion results - depl_nuc = prev_res[-1].nuc_to_ind - geom_nuc_densities = mat.get_nuclide_atom_densities() - - # Merge lists of nuclides, with the same order for every calculation - geom_nuc_densities.update(depl_nuc) - - for nuclide, atom_per_bcm in geom_nuc_densities.items(): - if nuclide in depl_nuc: - concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] - volume = prev_res[-1].volume[mat_id] - atom_per_cc = concentration / volume - else: - atom_per_cc = atom_per_bcm * 1.0e24 - - self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - def initial_condition(self): """Performs final setup and returns initial condition. @@ -589,16 +439,8 @@ class Operator(TransportOperator): # Generate tallies in memory materials = [openmc.lib.materials[int(i)] for i in self.burnable_mats] - self._rate_helper.generate_tallies(materials, self.chain.reactions) - self._normalization_helper.prepare( - self.chain.nuclides, self.reaction_rates.index_nuc) - # Tell fission yield helper what materials this process is - # responsible for - self._yield_helper.generate_tallies( - materials, tuple(sorted(self._mat_index_map.values()))) - # Return number density vector - return list(self.number.get_mat_slice(np.s_[:])) + return super().initial_condition(materials) def finalize(self): """Finalize a depletion simulation and release resources.""" @@ -659,127 +501,6 @@ class Operator(TransportOperator): self.materials.export_to_xml() - def _get_tally_nuclides(self): - """Determine nuclides that should be tallied for reaction rates. - - This method returns a list of all nuclides that have neutron data and - are listed in the depletion chain. Technically, we should tally nuclides - that may not appear in the depletion chain because we still need to get - the fission reaction rate for these nuclides in order to normalize - power, but that is left as a future exercise. - - Returns - ------- - list of str - Tally nuclides - - """ - nuc_set = set() - - # Create the set of all nuclides in the decay chain in materials marked - # for burning in which the number density is greater than zero. - for nuc in self.number.nuclides: - if nuc in self.nuclides_with_data: - if np.sum(self.number[:, nuc]) > 0.0: - nuc_set.add(nuc) - - # Communicate which nuclides have nonzeros to rank 0 - if comm.rank == 0: - for i in range(1, comm.size): - nuc_newset = comm.recv(source=i, tag=i) - nuc_set |= nuc_newset - else: - comm.send(nuc_set, dest=0, tag=comm.rank) - - if comm.rank == 0: - # Sort nuclides in the same order as self.number - nuc_list = [nuc for nuc in self.number.nuclides - if nuc in nuc_set] - else: - nuc_list = None - - # Store list of tally nuclides on each process - nuc_list = comm.bcast(nuc_list) - return [nuc for nuc in nuc_list if nuc in self.chain] - - def _unpack_tallies_and_normalize(self, source_rate): - """Unpack tallies from OpenMC and return an operator result - - This method uses OpenMC's C API bindings to determine the k-effective - value and reaction rates from the simulation. The reaction rates are - normalized by a helper class depending on the method being used. - - Parameters - ---------- - source_rate : float - Power in [W] or source rate in [neutron/sec] - - Returns - ------- - openmc.deplete.OperatorResult - Eigenvalue and reaction rates resulting from transport operator - - """ - rates = self.reaction_rates - rates.fill(0.0) - - # Get k and uncertainty - keff = ufloat(*openmc.lib.keff()) - - # Extract tally bins - nuclides = self._rate_helper.nuclides - - # Form fast map - nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] - react_ind = [rates.index_rx[react] for react in self.chain.reactions] - - # Keep track of energy produced from all reactions in eV per source - # particle - self._normalization_helper.reset() - self._yield_helper.unpack() - - # Store fission yield dictionaries - fission_yields = [] - - # Create arrays to store fission Q values, reaction rates, and nuclide - # numbers, zeroed out in material iteration - number = np.empty(rates.n_nuc) - - fission_ind = rates.index_rx.get("fission") - - # Extract results - for i, mat in enumerate(self.local_mats): - # Get tally index - mat_index = self._mat_index_map[mat] - - # Zero out reaction rates and nuclide numbers - number.fill(0.0) - - # Get new number densities - for nuc, i_nuc_results in zip(nuclides, nuc_ind): - number[i_nuc_results] = self.number[mat, nuc] - - tally_rates = self._rate_helper.get_material_rates( - mat_index, nuc_ind, react_ind) - - # Compute fission yields for this material - fission_yields.append(self._yield_helper.weighted_yields(i)) - - # Accumulate energy from fission - if fission_ind is not None: - self._normalization_helper.update(tally_rates[:, fission_ind]) - - # Divide by total number and store - rates[i] = self._rate_helper.divide_by_adens(number) - - # Scale reaction rates to obtain units of reactions/sec - rates *= self._normalization_helper.factor(source_rate) - - # Store new fission yields on the chain - self.chain.fission_yields = fission_yields - - return OperatorResult(keff, rates) - def _get_nuclides_with_data(self, cross_sections): """Loads cross_sections.xml file to find nuclides with neutron data""" nuclides = set() @@ -792,31 +513,3 @@ class Operator(TransportOperator): nuclides.add(name) return nuclides - - def get_results_info(self): - """Returns volume list, material lists, and nuc lists. - - Returns - ------- - volume : dict of str float - Volumes corresponding to materials in full_burn_dict - nuc_list : list of str - A list of all nuclide names. Used for sorting the simulation. - burn_list : list of int - A list of all material IDs to be burned. Used for sorting the simulation. - full_burn_list : list - List of all burnable material IDs - - """ - nuc_list = self.number.burnable_nuclides - burn_list = self.local_mats - - volume = {} - for i, mat in enumerate(burn_list): - volume[mat] = self.number.volume[i] - - # Combine volume dictionaries across processes - volume_list = comm.allgather(volume) - volume = {k: v for d in volume_list for k, v in d.items()} - - return volume, nuc_list, burn_list, self.burnable_mats