diff --git a/docs/source/methods/depletion.rst b/docs/source/methods/depletion.rst index edcf2c3f534..aaf0ff228b3 100644 --- a/docs/source/methods/depletion.rst +++ b/docs/source/methods/depletion.rst @@ -176,14 +176,16 @@ only solving the transport equation to estimate transmutation reaction rates sections (in the case of transport-independent depletion), but also a series of choices about what data to include. In OpenMC, the burnup matrix is constructed based on data inside of a *depletion chain* file, which includes fundamental -data gathered from ENDF incident neutron, decay, and fission product yield -sublibraries. For each nuclide, this file includes: +data gathered from ENDF incident neutron, decay, fission product yield and +spontaneous fission yield sublibraries. For each nuclide, this file includes: - What transmutation reactions are possible, their Q values, and their products; - If a nuclide is not stable, what decay modes are possible, their branching - ratios, and their products; and + ratios, and their products; - If a nuclide is fissionable, the fission products yields at any number of - incident neutron energies. + incident neutron energies; and +- If a nuclide can undergo spontanous fission, the products yields for + spontaneous fission. Transmutation Reactions ----------------------- diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index f112cf8ccff..ac00679338f 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -135,6 +135,19 @@ The :class:`Chain` class uses information from the following module variable: :type: dict +It contains the following switch: + +.. data:: chain.INCLUDE_SPONT_FISSION + + Boolean switch to enable or disable the use of spontaneous fission + product yields when solving the Bateman equations. The default is to + include spontaneous fission product yields if they are present in + the Chain object. + + :type: bool + + + The following classes are used during a depletion simulation and store auxiliary data, such as number densities and reaction rates for each material. diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 261900ce61a..8fc269cb221 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -107,14 +107,15 @@ Caveats .. _energy-deposition: -Energy Deposition +Energy Deposition and Spontaneous Fission ~~~~~~~~~~~~~~~~~ The default energy deposition mode, ``"fission-q"``, instructs the :class:`~openmc.deplete.CoupledOperator` to normalize reaction rates using the product of fission reaction rates and fission Q values taken from the depletion chain. This approach does not consider indirect contributions to energy -deposition, such as neutron heating and energy from secondary photons. In doing +deposition, such as neutron heating and energy from secondary photons, nor does it +consider energy released in spontaneous fission. In doing this, the energy deposited during a transport calculation will be lower than expected. This causes the reaction rates to be over-adjusted to hit the user-specific power, or power density, leading to an over-depletion of burnable @@ -151,6 +152,13 @@ These modified heating libraries can be generated by running the latest version of :meth:`openmc.data.IncidentNeutron.from_njoy()`, and will eventually be bundled into the distributed libraries. +By default, depletion calculations include fission yields from spontaneous fission +decay. Such decay should result in the emission of neutrons but a spontaneous fission +source term is not included in the Boltzmann equation in current version of OpenMC. +In practice, this should be a minor issue as neutron-induced fission is much more +important than spontaneous fission for the sourcing of neutrons. + + Local Spectra and Repeated Materials ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 8e24f716bc2..be56bcd073d 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -25,6 +25,11 @@ import openmc.data +# Configurable switch that enables / disables the inclusion +# of spontantous fission yields when calculating the matrix in +# the Batemans equation. +INCLUDE_SPONT_FISSION = True + # tuple of (possible MT values, secondaries) ReactionInfo = namedtuple('ReactionInfo', ('mts', 'secondaries')) @@ -228,12 +233,61 @@ def replace_missing_fpy(actinide, fpy_data, decay_data): return 'U235' +def replace_missing_sfy(actinide, sfy_data, decay_data): + """Replace missing fission product yields + + Parameters + ---------- + actinide : str + Name of actinide missing SFY data + sfy_data : dict + Dictionary of SFY data + decay_data : dict + Dictionary of decay data + + Returns + ------- + str + Actinide that can be used as replacement for SFY purposes + + """ + + # Check if metastable state has data (e.g., Am242m) + Z, A, m = zam(actinide) + if m == 0: + metastable = gnds_name(Z, A, 1) + if metastable in sfy_data: + return metastable + + # Try increasing Z, holding N constant + isotone = actinide + while isotone in decay_data: + Z += 1 + A += 1 + isotone = gnds_name(Z, A, 0) + if isotone in sfy_data: + return isotone + + # Try decreasing Z, holding N constant + isotone = actinide + while isotone in decay_data: + Z -= 1 + A -= 1 + isotone = gnds_name(Z, A, 0) + if isotone in sfy_data: + return isotone + + # If all else fails, use U238 yields + return 'U238' + + class Chain: """Full representation of a depletion chain. A depletion chain can be created by using the :meth:`from_endf` method which requires a list of ENDF incident neutron, decay, and neutron fission product - yield sublibrary files. The depletion chain used during a depletion + yield sublibrary files and optionally spontanteous fission product yield sublibrary files. + The depletion chain used during a depletion simulation is indicated by either an argument to :class:`openmc.deplete.CoupledOperator` or :class:`openmc.deplete.IndependentOperator`, or through @@ -260,6 +314,9 @@ class Chain: Otherwise, an entry can be added for each material to be burned. Ordering should be identical to how the operator orders reaction rates for burnable materials. + spont_fission_yields: None or iterable of dict + List of effective spontaneous fission yields. + The dictionary takes the same form as fission_yields. """ def __init__(self): @@ -267,6 +324,7 @@ def __init__(self): self.reactions = [] self.nuclide_dict = {} self._fission_yields = None + self._spont_fission_yields = None def __contains__(self, nuclide): return nuclide in self.nuclide_dict @@ -308,7 +366,7 @@ def add_nuclide(self, nuclide: Nuclide): self.reactions.append(rx.type) @classmethod - def from_endf(cls, decay_files, fpy_files, neutron_files, + def from_endf(cls, decay_files, fpy_files, neutron_files, sfy_files=[], reactions=('(n,2n)', '(n,3n)', '(n,4n)', '(n,gamma)', '(n,p)', '(n,a)'), progress=True ): @@ -327,6 +385,8 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, List of ENDF neutron-induced fission product yield sub-library files neutron_files : list of str or openmc.data.endf.Evaluation List of ENDF neutron reaction sub-library files + sfy_files : list of str or openmc.data.endf.Evaluation + List of ENDF spontaneous fission product yield sub-library files reactions : iterable of str, optional Transmutation reactions to include in the depletion chain, e.g., `["(n,2n)", "(n,gamma)"]`. Note that fission is always included if @@ -388,12 +448,22 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, data = openmc.data.FissionProductYields(f) fpy_data[data.nuclide['name']] = data + if progress: + print('Processing spontaneous fission product yield sub-library files...') + sfy_data = {} + for f in sfy_files: + data = openmc.data.FissionProductYields(f) + sfy_data[data.nuclide['name']] = data + if progress: print('Creating depletion_chain...') missing_daughter = [] missing_rx_product = [] missing_fpy = [] + missing_sfy = [] missing_fp = [] + missing_sfp = [] + sfp_without_fp = [] chain = cls() for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)): @@ -494,14 +564,51 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, else: nuclide._fpy = replace_missing_fpy(parent, fpy_data, decay_data) missing_fpy.append((parent, nuclide._fpy)) + + # Spontaneous fission yield data + if 'sf' in [decay_mode.type for decay_mode in nuclide.decay_modes]: + if parent in sfy_data: + sfy = sfy_data[parent] + + yield_data = {} + for yield_table in sfy.independent: + yield_replace = 0.0 + yields = defaultdict(float) + for product, y in yield_table.items(): + # Handle fission products that have no decay data + if product not in decay_data: + daughter = replace_missing(product, decay_data) + product = daughter + yield_replace += y.nominal_value + + yields[product] += y.nominal_value + + if yield_replace > 0.0: + missing_sfp.append((parent, yield_replace)) + yield_data[float('inf')] = yields + + nuclide.spont_yield_data = FissionYieldDistribution(yield_data) + else: + nuclide._sfy = replace_missing_sfy(parent, sfy_data, decay_data) + missing_sfy.append((parent, nuclide._sfy)) + + if not fissionable: + sfp_without_fp.append(parent) + # Add nuclide to chain chain.add_nuclide(nuclide) + # Replace missing FPY data for nuclide in chain.nuclides: if hasattr(nuclide, '_fpy'): nuclide.yield_data = chain[nuclide._fpy].yield_data + + # Replace missing SFY data + for nuclide in chain.nuclides: + if hasattr(nuclide, '_sfy'): + nuclide.spont_yield_data = chain[nuclide._sfy].spont_yield_data # Display warnings if missing_daughter: @@ -522,11 +629,30 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, print(f' {parent}, replaced with {replacement}') print('') + if missing_sfy: + print("The following nuclides with spontaneous fission " + "have no spontaneous fission product yields:") + for parent, replacement in missing_sfy: + print(f' {parent}, replaced with {replacement}') + print('') + if missing_fp: print('The following nuclides have fission products with no decay data:') for vals in missing_fp: print(' {}, E={} eV (total yield={})'.format(*vals)) + if missing_sfp: + print("The following nuclides have spontaneous fission" + "products with no decay data:") + for vals in missing_sfp: + print(' {}, (total yield={})'.format(*vals)) + + if sfp_without_fp: + print("The following nuclides have spontaneous fission products" + " but no induced fission products:") + for parent in sfp_without_fp: + print(f' {parent}') + return chain @classmethod @@ -602,6 +728,25 @@ def get_default_fission_yields(self): out[nuc.name] = dict(yield_obj) return out + def get_spont_fission_yields(self): + """Return spontaneous fission yields. + + Returns + ------- + spont_fission_yields : dict + Dictionary of ``{parent: {product: f_yield}}`` + where ``parent`` and ``product`` are both string + names of nuclides with yield data and ``f_yield`` + is a float for the fission yield. + """ + out = defaultdict(dict) + for nuc in self.nuclides: + if nuc.spont_yield_data is None: + continue + yield_obj = nuc.spont_yield_data[float('inf')] + out[nuc.name] = dict(yield_obj) + return out + def form_matrix(self, rates, fission_yields=None): """Forms depletion matrix. @@ -633,6 +778,8 @@ def form_matrix(self, rates, fission_yields=None): if fission_yields is None: fission_yields = self.get_default_fission_yields() + spont_fission_yields = self.get_spont_fission_yields() + for i, nuc in enumerate(self.nuclides): # Loss from radioactive decay if nuc.half_life is not None: @@ -647,6 +794,15 @@ def form_matrix(self, rates, fission_yields=None): # Allow for total annihilation for debug purposes if branch_val != 0.0: + if decay_type == "sf": + if INCLUDE_SPONT_FISSION and nuc.spont_yield_data is not None: + print("HERE FOR SPONT FISS") + for product, y in spont_fission_yields[nuc.name].items(): + yield_val = y * branch_val + if yield_val != 0.0: + k = self.nuclide_dict[product] + matrix[k, i] += yield_val + if target is not None: k = self.nuclide_dict[target] matrix[k, i] += branch_val @@ -1028,6 +1184,12 @@ def fission_yields(self): self._fission_yields = [self.get_default_fission_yields()] return self._fission_yields + @property + def spont_fission_yields(self): + if self._spont_fission_yields is None: + self._spont_fission_yields = [self.get_spont_fission_yields()] + return self._spont_fission_yields + @fission_yields.setter def fission_yields(self, yields): _invalidate_chain_cache(self) @@ -1037,6 +1199,15 @@ def fission_yields(self, yields): check_type("fission_yields", yields, Iterable, Mapping) self._fission_yields = yields + @spont_fission_yields.setter + def spont_fission_yields(self, spont_yields): + _invalidate_chain_cache(self) + if spont_yields is not None: + if isinstance(spont_yields, Mapping): + spont_yields = [spont_yields] + check_type("spont_fission_yields", spont_yields, Iterable, Mapping) + self._spont_fission_yields = spont_yields + def validate(self, strict=True, quiet=False, tolerance=1e-4): """Search for possible inconsistencies @@ -1046,6 +1217,8 @@ def validate(self, strict=True, quiet=False, tolerance=1e-4): ratios equal about one? 2) For fission reactions, does the sum of fission yield fractions equal about two? + 3) For spontanous fission, does the sum of fission yield + fractions equal about two? Parameters ---------- @@ -1155,13 +1328,15 @@ def reduce(self, initial_isotopes, level=None): new_nuclide.sources = previous.sources.copy() if hasattr(previous, '_fpy'): new_nuclide._fpy = previous._fpy + if hasattr(previous, '_sfy'): + new_nuclide._sfy = previous._sfy for mode in previous.decay_modes: if mode.target in all_isotopes: new_nuclide.add_decay_mode(*mode) else: new_nuclide.add_decay_mode(mode.type, None, mode.branching_ratio) - + for rx in previous.reactions: if rx.target in all_isotopes: new_nuclide.add_reaction(*rx) @@ -1174,6 +1349,12 @@ def reduce(self, initial_isotopes, level=None): else: new_nuclide.add_reaction(rx.type, None, rx.Q, rx.branching_ratio) + if previous.spont_yield_data is not None: + new_sf_yields = new_nuclide.spont_yield_data = ( + previous.spont_yield_data.restrict_products(name_sort)) + else: + new_nuclide.spont_yield_data = None + new_chain.add_nuclide(new_nuclide) # Doesn't appear that the ordering matters for the reactions, @@ -1229,8 +1410,8 @@ def _follow(self, isotopes, level): continue next_iso.add(product) - if nuclide.yield_data is not None: - for product in nuclide.yield_data.products: + if nuclide.yield_data is not None or nuclide.spont_yield_data is not None: + for product in nuclide.yield_data.products + nuclide.spont_yield_data.products: if (product in next_iso or product in found or product in isotopes): continue diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 60e3e5317f1..1b6a64a6bba 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -105,6 +105,10 @@ class Nuclide: yield_data : FissionYieldDistribution or None Fission product yields at tabulated energies for this nuclide. Can be treated as a nested dictionary ``{energy: {product: yield}}`` + spont_yield_data : FissionYieldDistribution or None + Spontaneous fission product yields. To keep the same structure as for + yield_data this is treated as a nested dictionary + ``{energy: {product: yield}}`` with one energy value = Inf. yield_energies : tuple of float or None Energies at which fission product yields exist """ @@ -127,6 +131,9 @@ def __init__(self, name=None): # Neutron fission yields, if present self._yield_data = None + # Spontanteous fission yields, if present + self._spont_yield_data = None + def __repr__(self): n_modes, n_rx = self.n_decay_modes, self.n_reaction_paths return f"" @@ -145,6 +152,12 @@ def yield_data(self): return None return self._yield_data + @property + def spont_yield_data(self): + if self._spont_yield_data is None: + return None + return self._spont_yield_data + @yield_data.setter def yield_data(self, fission_yields): if fission_yields is None: @@ -156,6 +169,17 @@ def yield_data(self, fission_yields): else: self._yield_data = FissionYieldDistribution(fission_yields) + @spont_yield_data.setter + def spont_yield_data(self, fission_yields): + if fission_yields is None: + self._spont_yield_data = None + else: + check_type("fission_yields", fission_yields, Mapping) + if isinstance(fission_yields, FissionYieldDistribution): + self._spont_yield_data = fission_yields + else: + self._spont_yield_data = FissionYieldDistribution(fission_yields) + @property def yield_energies(self): if self._yield_data is None: @@ -287,6 +311,27 @@ def from_xml(cls, element, root=None, fission_q=None): nuc.yield_data = FissionYieldDistribution.from_xml_element(fpy_elem) + + sfy_elem = element.find('spont_fission_yields') + if sfy_elem is not None: + # Check for use of SFY from other nuclide + parent = sfy_elem.get('parent') + if parent is not None: + assert root is not None + sfy_elem = root.find( + f'.//nuclide[@name="{parent}"]/spont_fission_yields' + ) + if sfy_elem is None: + raise ValueError( + "Spontanteous fission product yields for {0} borrow" + " from {1}, but {1} is not present in the chain file" + " or has no yields.".format( + nuc.name, parent + )) + nuc._sfy = parent + + nuc.spont_yield_data = FissionYieldDistribution.from_xml_element(sfy_elem) + return nuc def to_xml_element(self): @@ -340,6 +385,15 @@ def to_xml_element(self): energy_elem.text = ' '.join(str(E) for E in self.yield_energies) self.yield_data.to_xml_element(fpy_elem) + if self.spont_yield_data: + sfy_elem = ET.SubElement(elem, 'spont_fission_yields') + + if hasattr(self, '_sfy'): + # Check for link to other nuclide data + sfy_elem.set('parent', self._sfy) + else: + self.spont_yield_data.to_xml_element(sfy_elem) + return elem def validate(self, strict=True, quiet=False, tolerance=1e-4): @@ -351,6 +405,8 @@ def validate(self, strict=True, quiet=False, tolerance=1e-4): does the sum of branching ratios equal about one? 2) for fission reactions, does the sum of fission yield fractions equal about two? + 3) For spontanous fission, does the sum of fission yield + fractions equal about two? Parameters ---------- @@ -438,6 +494,24 @@ def validate(self, strict=True, quiet=False, tolerance=1e-4): warn(msg) valid = False + if self.spont_yield_data: + for _, fission_yield in self.spont_yield_data.items(): + sum_yield = fission_yield.yields.sum() + stat = 2.0 - tolerance <= sum_yield <= 2.0 + tolerance + if stat: + continue + msg = msg_func( + name=self.name, actual=sum_yield, + expected=2.0, tol=tolerance, + prop=f"spontaneous fission yields") + if strict: + raise ValueError(msg) + elif quiet: + return False + warn(msg) + valid = False + + return valid