diff --git a/openmc/deplete/coupled_operator.py b/openmc/deplete/coupled_operator.py index 34bb28b491b..486dde38f4e 100644 --- a/openmc/deplete/coupled_operator.py +++ b/openmc/deplete/coupled_operator.py @@ -26,7 +26,9 @@ from .results import Results from .helpers import ( DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, - FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, + FissionYieldCutoffHelper, AveragedFissionYieldHelper, + LogLinInterpolateFissionYieldHelper, LinLinInterpolateFissionYieldHelper, + EnergyScoreHelper, SourceRateHelper, FluxCollapseHelper) @@ -123,13 +125,15 @@ class CoupledOperator(OpenMCOperator): 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"`` - fission_yield_mode : {"constant", "cutoff", "average"} + fission_yield_mode : {"constant", "cutoff", "average", "loglin", "linlin"} Key indicating what fission product yield scheme to use. The key determines what fission energy helper is used: * "constant": :class:`~openmc.deplete.helpers.ConstantFissionYieldHelper` * "cutoff": :class:`~openmc.deplete.helpers.FissionYieldCutoffHelper` * "average": :class:`~openmc.deplete.helpers.AveragedFissionYieldHelper` + * "loglin": :class:`~openmc.deplete.helpers.LogLinInterpolateFissionYieldHelper` + * "linlin": :class:`~openmc.deplete.helpers.LinLinInterpolateFissionYieldHelper` The documentation on these classes describe their methodology and differences. Default: ``"constant"`` @@ -201,6 +205,8 @@ class CoupledOperator(OpenMCOperator): "average": AveragedFissionYieldHelper, "constant": ConstantFissionYieldHelper, "cutoff": FissionYieldCutoffHelper, + "loglin": LogLinInterpolateFissionYieldHelper, + "linlin": LinLinInterpolateFissionYieldHelper, } def __init__(self, model, chain_file=None, prev_results=None, @@ -469,6 +475,12 @@ def _update_materials(self): if nuc in self.nuclides_with_data: val = 1.0e-24 * number_i.get_atom_density(mat, nuc) + # Guard against non-finite values that can arise when + # depletion solver linear systems become ill-conditioned. + if not np.isfinite(val): + number_i[mat, nuc] = 0.0 + continue + # If nuclide is zero, do not add to the problem. if val > 0.0: if self.round_number: diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 14780d61a20..a5e32ba045f 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -9,7 +9,7 @@ from numbers import Real import sys -from numpy import dot, zeros, newaxis, asarray +from numpy import dot, zeros, newaxis, asarray, isfinite, clip from openmc.mpi import comm from openmc.checkvalue import check_type, check_greater_than @@ -21,10 +21,12 @@ ReactionRateHelper, NormalizationHelper, FissionYieldHelper) __all__ = ( - "DirectReactionRateHelper", "ChainFissionHelper", "EnergyScoreHelper" + "DirectReactionRateHelper", "ChainFissionHelper", "EnergyScoreHelper", "SourceRateHelper", "TalliedFissionYieldHelper", "ConstantFissionYieldHelper", "FissionYieldCutoffHelper", - "AveragedFissionYieldHelper", "FluxCollapseHelper") + "AveragedFissionYieldHelper", "LogLinInterpolateFissionYieldHelper", + "LinLinInterpolateFissionYieldHelper", + "FluxCollapseHelper") class TalliedFissionYieldHelper(FissionYieldHelper): @@ -1019,3 +1021,306 @@ def from_operator(cls, operator, **kwargs): AveragedFissionYieldHelper """ return cls(operator.chain.nuclides) + +class LogLinInterpolateFissionYieldHelper(TalliedFissionYieldHelper): + r"""Helper that computes fission yields from log-lin weighted groups. + + Parameters + ---------- + chain_nuclides : iterable of openmc.deplete.Nuclide + Nuclides tracked in the depletion chain. All nuclides are + not required to have fission yield data. + energy_bins : iterable of float, optional + Energy points [eV] that define ``len(energy_bins)`` log-lin + groups. Must be strictly increasing and positive. + Defaults to ``(0.025, 500e3, 14e6)``. + + Attributes + ---------- + constant_yields : collections.defaultdict + Fission yields for all nuclides that only have one set of + fission yield data. Dictionary of form ``{str: {str: float}}`` + representing yields for ``{parent: {product: yield}}``. Default + return object is an empty dictionary + results : None or numpy.ndarray + If tallies have been generated and unpacked, then the array will + have shape ``(n_mats, n_groups, n_tnucs)``, where ``n_mats`` is the + number of materials where fission reactions were tallied, + ``n_groups = len(energy_bins)``, and ``n_tnucs`` is the number + of nuclides with multiple sets of fission yields. + Data are normalized group fractions corresponding to ``y1..yN``. + + Notes + ----- + Fission events are partitioned into ``len(energy_bins)`` group + fractions using the provided energy points. Adjacent groups share each + interval and are blended using log-linear interpolation in incident energy. + For energies below the first point, all contribution is assigned to the + first group. For energies above the last point, all contribution is + assigned to the last group. + + Group fractions are then used to blend the corresponding yield libraries + for each nuclide. + + Implementation overview + + The helper creates one denominator tally and ``N`` weighted numerator + tallies, where ``N = len(energy_bins)``: + + 1. Denominator tally + ``self._fission_rate_tally`` is given an + :class:`openmc.lib.EnergyFilter` over ``[0, self._upper_energy]``. This + stores total fission scores used for normalization. + 2. Group filters + A bank of ``N`` :class:`openmc.lib.EnergyFunctionFilter` objects is + built. Each one is configured as: + + - ``filt = EnergyFunctionFilter()`` + - ``filt.set_data(x_nodes, y_nodes)`` + - ``filt.interpolation = 'log-linear'`` (when supported by wrapper) + + Here ``x_nodes = [e0, *energy_bins, e_top]`` and each + group's ``y_nodes`` determines its weight-vs-energy shape. + 3. Weighted tallies + For each group filter, a tally with ``score='fission'`` is created in + ``self._weighted_tallies``. + 4. Normalized fractions + In :meth:`unpack`, each weighted tally is divided by the denominator + tally to produce per-group fractions ``w_g``. + + Depletion-chain combination + + For each nuclide, the final effective fission-yield vector is computed from + group-wise yield vectors in the depletion chain: + + .. math:: + + \vec{f} = \sum_g \vec{f}_g\, w_g + + where ``\vec{f}_g`` is the chain fission-yield vector selected for group + ``g``, and ``w_g`` is the normalized group fraction from tallies. + + Example with bins ``(0.025, 500e3, 14e6)`` + + This gives 3 groups ``(y1, y2, y3)`` with + ``x_nodes = [e0, 0.025, 500e3, 14e6, e_top]`` and: + + - Group 1 ``y_nodes``: ``[1, 1, 0, 0, 0]`` + - Group 2 ``y_nodes``: ``[0, 0, 1, 0, 0]`` + - Group 3 ``y_nodes``: ``[0, 0, 0, 1, 1]`` + + For a fission event at ``E = 1000 eV`` (between ``0.025`` and ``500e3``), + only groups 1 and 2 are active. Their log-lin interval weights are: + + .. math:: + + w_2(E) = \frac{\ln E - \ln E_1}{\ln E_2 - \ln E_1}, \qquad + w_1(E) = 1 - w_2(E) + + with ``E_1 = 0.025``, ``E_2 = 500e3`` and ``E = 1000`` giving + ``w_1 \approx 0.3697`` and ``w_2 \approx 0.6303``. + + For unit event score before unpacking: + + - ``T1 -> T1 + 0.3697`` + - ``T2 -> T2 + 0.6303`` + - ``T3 -> T3 + 0.0`` + - ``D -> D + 1.0`` + + So the event contributes to the chain-weighted yield as + ``\vec{f} \approx 0.3697\,\vec{f}_1 + 0.6303\,\vec{f}_2`` for this + interval. + + Examples + -------- + operator=openmc.deplete.CoupledOperator(model,chain, fission_yield_mode='loglin', fission_yield_opts={'energy_bins':(0.025, 500.0e3, 14.0e6)}) + + + """ + + def __init__(self, chain_nuclides, energy_bins=(0.025, 500.0e3, 14.0e6)): + super().__init__(chain_nuclides) + self._weighted_tallies = [] + + energy_bins = asarray(energy_bins).ravel() + if energy_bins.size < 1: + raise ValueError("energy_bins must have at least one entry") + for i, ene in enumerate(energy_bins): + check_greater_than(f"energy_bins[{i}]", ene, 0.0, equality=False) + if i > 0: + check_greater_than( + f"energy_bins[{i}]", ene, energy_bins[i - 1], equality=False) + check_greater_than( + "upper tally energy", self._upper_energy, energy_bins[-1], equality=False) + self._energy_bins = energy_bins + self._n_groups = len(self._energy_bins) + + # Store one yield set per group for each nuclide. Out-of-range incident + # energies are clamped by the group filters to the first/last group. + self._group_yields = {} + convert_to_constant = set() + for name, nuc in self._chain_nuclides.items(): + energies = nuc.yield_energies + yields = nuc.yield_data + + group_energies = [ + min(energies, key=lambda e: abs(e - boundary)) + for boundary in self._energy_bins + ] + + if len(set(group_energies)) == 1: + self._constant_yields[name] = yields[group_energies[0]] + convert_to_constant.add(name) + continue + + self._group_yields[name] = tuple(yields[e] for e in group_energies) + + for name in convert_to_constant: + self._chain_nuclides.pop(name) + self._chain_set = set(self._chain_nuclides) | set(self._constant_yields) + + def _set_loglin_interpolation(self, filt): + try: + filt.interpolation = 'log-linear' + except AttributeError: + # Older Python wrapper in some versions does not expose this + # property; C++ defaults to linear-linear in this case. + pass + + def generate_tallies(self, materials, mat_indexes): + """Construct weighted fission tallies for arbitrary log-lin groups. + + Detailed algorithm notes and worked examples are documented in the + class docstring. + """ + super().generate_tallies(materials, mat_indexes) + fission_tally = self._fission_rate_tally + filters = fission_tally.filters + + # Restrict denominator to the same energy range used by weighting + # functions so normalized group fractions remain consistent. + ene_filter = EnergyFilter([0.0, self._upper_energy]) + fission_tally.filters = filters + [ene_filter] + + e0 = min(1.0e-12, self._energy_bins[0] * 1.0e-6) + e_top = self._upper_energy + + # Clamp outside [min_bin, max_bin] by keeping first/last weights at 1. + x_nodes = [e0, *self._energy_bins, e_top] + + # Build one weighted tally filter per user-specified energy bin. + filters_by_group = [] + for i_group in range(self._n_groups): + y_nodes = [0.0] * len(x_nodes) + if i_group == 0: + y_nodes[0] = 1.0 + + # Anchor each boundary to one group index so adjacent intervals + # naturally blend between neighboring groups. + for i_boundary, _ in enumerate(self._energy_bins): + if i_group == i_boundary: + y_nodes[i_boundary + 1] = 1.0 + + # Final group receives full weight above the highest boundary. + if i_group == self._n_groups - 1: + y_nodes[-1] = 1.0 + + filt = EnergyFunctionFilter() + filt.set_data(tuple(x_nodes), tuple(y_nodes)) + self._set_loglin_interpolation(filt) + filters_by_group.append(filt) + + self._weighted_tallies = [] + for filt in filters_by_group: + weighted_tally = Tally() + weighted_tally.writable = False + weighted_tally.scores = ['fission'] + weighted_tally.filters = filters + [filt] + self._weighted_tallies.append(weighted_tally) + + def update_tally_nuclides(self, nuclides): + """Tally nuclides with non-zero density and multiple yields.""" + tally_nucs = super().update_tally_nuclides(nuclides) + for weighted_tally in self._weighted_tallies: + weighted_tally.nuclides = tally_nucs + return tally_nucs + + def unpack(self): + """Unpack tallies and populate normalized group fractions.""" + if not self._tally_nucs or self._local_indexes.size == 0: + self.results = None + return + + fission_results = self._fission_rate_tally.mean[self._local_indexes] + n_mats, n_nucs = fission_results.shape + self.results = zeros((n_mats, self._n_groups, n_nucs)) + + for i, weighted_tally in enumerate(self._weighted_tallies): + self.results[:, i, :] = weighted_tally.mean[self._local_indexes] + + nz_mat, nz_nuc = fission_results.nonzero() + self.results[nz_mat, :, nz_nuc] /= fission_results[nz_mat, newaxis, nz_nuc] + + def weighted_yields(self, local_mat_index): + """Return weighted fission yields for one material. + + Group fractions from :attr:`results` are applied to per-group yields + selected for each nuclide. + """ + if not self._tally_nucs: + return self.constant_yields + + mat_yields = defaultdict(dict) + group_fracs = self.results[local_mat_index] + for i_nuc, nuc in enumerate(self._tally_nucs): + group_yields = self._group_yields[nuc.name] + group_weights = asarray(group_fracs[:, i_nuc], dtype=float) + + # Low-statistics runs can produce non-finite or slightly negative + # normalized fractions. Keep physically meaningful non-negative + # weights and renormalize before combining yield libraries. + if not isfinite(group_weights).all(): + group_weights[:] = 0.0 + group_weights[0] = 1.0 + else: + group_weights = clip(group_weights, 0.0, None) + weight_sum = group_weights.sum() + if weight_sum > 0.0: + group_weights /= weight_sum + else: + group_weights[:] = 0.0 + group_weights[0] = 1.0 + + # FissionYield objects cannot be summed with the default integer + # start value used by sum(), so accumulate explicitly. + weighted = None + for y, w in zip(group_yields, group_weights): + term = y * w + weighted = term if weighted is None else weighted + term + mat_yields[nuc.name] = weighted + + mat_yields.update(self.constant_yields) + return mat_yields + + @classmethod + def from_operator(cls, operator, **kwargs): + """Return a new helper with data from an operator.""" + return cls(operator.chain.nuclides, **kwargs) + + +class LinLinInterpolateFissionYieldHelper(LogLinInterpolateFissionYieldHelper): + r"""Helper that computes fission yields from lin-lin weighted groups. + + This helper uses the same tally construction and yield-combination logic as + :class:`LogLinInterpolateFissionYieldHelper`, but the + :class:`openmc.lib.EnergyFunctionFilter` interpolation is set to + ``'linear-linear'`` when supported by the Python wrapper. + """ + + def _set_loglin_interpolation(self, filt): + try: + filt.interpolation = 'linear-linear' + except AttributeError: + # Older Python wrapper in some versions does not expose this + # property; C++ defaults to linear-linear in this case. + pass \ No newline at end of file diff --git a/tests/unit_tests/test_deplete_fission_yields.py b/tests/unit_tests/test_deplete_fission_yields.py index 1937e61e333..e17e0082784 100644 --- a/tests/unit_tests/test_deplete_fission_yields.py +++ b/tests/unit_tests/test_deplete_fission_yields.py @@ -12,7 +12,8 @@ from openmc.deplete.nuclide import Nuclide, FissionYieldDistribution from openmc.deplete.helpers import ( FissionYieldCutoffHelper, ConstantFissionYieldHelper, - AveragedFissionYieldHelper) + AveragedFissionYieldHelper, LogLinInterpolateFissionYieldHelper, + LinLinInterpolateFissionYieldHelper) @pytest.fixture(scope="module") @@ -289,3 +290,76 @@ def interp_average_yields(nuc, avg_energy): assert thermal_E < avg_energy < fast_E split = (avg_energy - thermal_E)/(fast_E - thermal_E) return yields[thermal_E]*(1 - split) + yields[fast_E]*split + + +@pytest.mark.parametrize( + "helper_cls, expected_split", + ( + ( + LogLinInterpolateFissionYieldHelper, + lambda E, E1, E2: (np.log(E) - np.log(E1)) / (np.log(E2) - np.log(E1)), + ), + ( + LinLinInterpolateFissionYieldHelper, + lambda E, E1, E2: (E - E1) / (E2 - E1), + ), + ), +) +def test_interpolate_helpers_single_event(materials, nuclide_bundle, helper_cls, + expected_split): + """Emulate one fission event at 1000 eV and verify group weighting.""" + energy_bins = (0.025, 500.0e3, 14.0e6) + event_energy = 1.0e3 + helper = helper_cls(nuclide_bundle, energy_bins=energy_bins) + helper.generate_tallies(materials, [0]) + + tallied_nucs = helper.update_tally_nuclides([n.name for n in nuclide_bundle]) + assert tallied_nucs == ["Pu239", "U235"] + + # Emulate denominator tally with exactly one fission event in material 0. + fission_data = proxy_tally_data(helper._fission_rate_tally, fill=0.0) + fission_data[0] = 1.0 + helper._fission_rate_tally = Mock() + helper._fission_rate_tally.mean = fission_data + + # Event at 1000 eV lies between first two bins, so only groups 1 and 2 + # contribute. Group 3 is zero by construction. + split = expected_split(event_energy, energy_bins[0], energy_bins[1]) + w1 = 1.0 - split + w2 = split + w3 = 0.0 + + weighted_vals = (w1, w2, w3) + for i_group, weighted_tally in enumerate(helper._weighted_tallies): + data = proxy_tally_data(weighted_tally, fill=0.0) + data[0] = weighted_vals[i_group] + helper._weighted_tallies[i_group] = Mock() + helper._weighted_tallies[i_group].mean = data + + helper.unpack() + + expected_results = np.empty((1, 3, len(tallied_nucs))) + expected_results[:, 0] = w1 + expected_results[:, 1] = w2 + expected_results[:, 2] = w3 + assert helper.results == pytest.approx(expected_results) + + actual_yields = helper.weighted_yields(0) + + # U238 has one yield set and is therefore constant. + assert actual_yields["U238"] == nuclide_bundle.u238.yield_data[5e5] + + nuc_by_name = {n.name: n for n in nuclide_bundle} + for nuc in tallied_nucs: + nuclide = nuc_by_name[nuc] + group_energies = [ + min(nuclide.yield_energies, + key=lambda e: abs(e - boundary)) + for boundary in energy_bins + ] + expected = ( + nuclide.yield_data[group_energies[0]] * w1 + + nuclide.yield_data[group_energies[1]] * w2 + + nuclide.yield_data[group_energies[2]] * w3 + ) + assert actual_yields[nuc] == expected