From 608a1c3386db6af88d70d35cc2202d3955aea30b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 23 Jun 2026 02:00:47 -0500 Subject: [PATCH] Fix several issues related to independent operator depletion (#3977) --- openmc/deplete/independent_operator.py | 4 ++ openmc/deplete/microxs.py | 63 ++++++++++++++++++------ tests/unit_tests/test_deplete_microxs.py | 60 ++++++++++++++++++++++ 3 files changed, 112 insertions(+), 15 deletions(-) diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index c12863956b..bdfc32763b 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -339,8 +339,12 @@ class IndependentOperator(OpenMCOperator): for i_nuc in nuc_index: nuc = self.nuc_ind_map[i_nuc] + if nuc not in xs._index_nuc: + continue for i_rx in react_index: rx = self.rx_ind_map[i_rx] + if rx not in xs._index_rx: + continue # Determine reaction rate by multiplying xs in [b] by flux # in [n-cm/src] to give [(reactions/src)*b-cm/atom] diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 687cf646f2..42bb958caf 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -84,7 +84,12 @@ def get_microxs_and_flux( reactions listed in the depletion chain file are used. energies : iterable of float or str Energy group boundaries in [eV] or the name of the group structure. - If left as None energies will default to [0.0, 100e6] + If left as None, no energy filter is applied to the flux tally. When + `reaction_rate_mode` is "direct", these boundaries define the output + flux and microscopic cross section energy group structure. When + `reaction_rate_mode` is "flux", these boundaries define the multigroup + flux tally used to collapse continuous-energy cross sections; returned + fluxes and microscopic cross sections are one-group. reaction_rate_mode : {"direct", "flux"}, optional The "direct" method tallies reaction rates directly (per energy group). The "flux" method tallies a multigroup flux spectrum and then @@ -110,7 +115,9 @@ def get_microxs_and_flux( reaction_rate_opts : dict, optional When `reaction_rate_mode="flux"`, allows selecting a subset of nuclide/reaction pairs to be computed via direct reaction-rate tallies - (per energy group). Supported keys: "nuclides", "reactions". + over one energy bin spanning the full `energies` range. Supported keys: + "nuclides", "reactions". If "reactions" are specified without + "nuclides", all selected nuclides are used. Returns ------- @@ -139,10 +146,14 @@ def get_microxs_and_flux( nuclides = [nuc.name for nuc in chain.nuclides if nuc.name in nuclides_with_data] - # Set up the reaction rate and flux tallies + # Set up the reaction rate and flux tallies. When energies are omitted, no + # energy filter is needed for the transport calculation. A one-group energy + # range is still needed later if flux collapse is requested. + collapse_energies = energies if energies is None: - energies = [0.0, 100.0e6] - if isinstance(energies, str): + energy_filter = None + collapse_energies = [0.0, 100.0e6] + elif isinstance(energies, str): energy_filter = openmc.EnergyFilter.from_group_structure(energies) else: energy_filter = openmc.EnergyFilter(energies) @@ -172,8 +183,11 @@ def get_microxs_and_flux( rr_reactions = list(reactions) elif reaction_rate_mode == 'flux' and reaction_rate_opts: opts = reaction_rate_opts or {} - rr_nuclides = list(opts.get('nuclides', [])) rr_reactions = list(opts.get('reactions', [])) + if rr_reactions: + rr_nuclides = list(opts.get('nuclides', nuclides)) + else: + rr_nuclides = list(opts.get('nuclides', [])) # Keep only requested pairs within overall sets if rr_nuclides: rr_nuclides = [n for n in rr_nuclides if n in set(nuclides)] @@ -182,7 +196,7 @@ def get_microxs_and_flux( # Use 1-group energy filter for RR in flux mode has_rr = bool(rr_nuclides and rr_reactions) - if has_rr and reaction_rate_mode == 'flux': + if has_rr and reaction_rate_mode == 'flux' and energy_filter is not None: rr_energy_filter = openmc.EnergyFilter( [energy_filter.values[0], energy_filter.values[-1]]) else: @@ -194,14 +208,18 @@ def get_microxs_and_flux( model.tallies = [] for i, domain_filter in enumerate(domain_filters): flux_tally = openmc.Tally(name=f'MicroXS flux {i}') - flux_tally.filters = [domain_filter, energy_filter] + flux_tally.filters = [domain_filter] + if energy_filter is not None: + flux_tally.filters.append(energy_filter) flux_tally.scores = ['flux'] model.tallies.append(flux_tally) flux_tallies.append(flux_tally) if has_rr: rr_tally = openmc.Tally(name=f'MicroXS RR {i}') - rr_tally.filters = [domain_filter, rr_energy_filter] + rr_tally.filters = [domain_filter] + if rr_energy_filter is not None: + rr_tally.filters.append(rr_energy_filter) rr_tally.nuclides = rr_nuclides rr_tally.multiply_density = False rr_tally.scores = rr_reactions @@ -255,8 +273,12 @@ def get_microxs_and_flux( all_flux_arrays = [] for flux_tally in flux_tallies: # Get flux values and make energy groups last dimension - flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1) - flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups) + flux = flux_tally.get_reshaped_data() + if energy_filter is None: + flux = flux[..., np.newaxis] # (domains, 1, 1, groups) + else: + # (domains, groups, 1, 1) -> (domains, 1, 1, groups) + flux = np.moveaxis(flux, 1, -1) all_flux_arrays.append(flux) fluxes.extend(flux.squeeze((1, 2))) @@ -266,8 +288,15 @@ def get_microxs_and_flux( for flux_arr, rr_tally in zip(all_flux_arrays, rr_tallies): flux = flux_arr # Get reaction rates and make energy groups last dimension - reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions) - reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups) + reaction_rates = rr_tally.get_reshaped_data() + if rr_energy_filter is None: + # (domains, nuclides, reactions) -> + # (domains, nuclides, reactions, groups) + reaction_rates = reaction_rates[..., np.newaxis] + else: + # (domains, groups, nuclides, reactions) -> + # (domains, nuclides, reactions, groups) + reaction_rates = np.moveaxis(reaction_rates, 1, -1) # If RR is 1-group, sum flux over groups if reaction_rate_mode == "flux": @@ -279,16 +308,20 @@ def get_microxs_and_flux( direct_micros.extend( MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs) - # If using flux mode, compute flux-collapsed microscopic XS if reaction_rate_mode == 'flux': + # Compute flux-collapsed microscopic XS flux_micros = [MicroXS.from_multigroup_flux( - energies=energies, + energies=collapse_energies, multigroup_flux=flux_i, chain_file=chain_file, nuclides=nuclides, reactions=reactions ) for flux_i in fluxes] + # We need to return one-group fluxes to match the microscopic cross + # sections, which are always one-group by virtue of the collapse + fluxes = [flux.sum(keepdims=True) for flux in fluxes] + # Decide which micros to use and merge if needed if reaction_rate_mode == 'flux' and rr_tallies: micros = [m1.merge(m2) for m1, m2 in zip(flux_micros, direct_micros)] diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index 26529e6ce9..0b1937facd 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -179,6 +179,66 @@ def test_hybrid_tally_setup(): assert ef.values[0] == pytest.approx(energies[0]) assert ef.values[-1] == pytest.approx(energies[-1]) + +def _simple_model(): + model = openmc.Model() + mat = openmc.Material(components={'H1': 1.0, 'H2': 1.0}, + density=5.0, density_units='g/cm3') + sphere = openmc.Sphere(r=10.0, boundary_type='vacuum') + cell = openmc.Cell(region=-sphere, fill=mat) + model.geometry = openmc.Geometry([cell]) + model.settings.particles = 100 + model.settings.batches = 5 + model.settings.run_mode = 'fixed source' + return model, mat + + +def test_hybrid_tally_defaults_to_all_nuclides(run_in_tmpdir): + energies = [0., 0.625, 2.0e7] + kwargs = { + 'nuclides': ['H1', 'H2'], + 'reactions': ['(n,2n)', '(n,gamma)'], + 'energies': energies, + 'reaction_rate_mode': 'flux', + 'chain_file': CHAIN_FILE, + } + + model, mat = _simple_model() + default_fluxes, default_micros = get_microxs_and_flux( + model, [mat], reaction_rate_opts={'reactions': ['(n,2n)']}, **kwargs + ) + + model, mat = _simple_model() + explicit_fluxes, explicit_micros = get_microxs_and_flux( + model, [mat], + reaction_rate_opts={ + 'nuclides': ['H1', 'H2'], + 'reactions': ['(n,2n)'] + }, + **kwargs + ) + + np.testing.assert_allclose(default_fluxes[0], explicit_fluxes[0]) + np.testing.assert_allclose(default_micros[0].data, explicit_micros[0].data) + assert default_micros[0].nuclides == explicit_micros[0].nuclides + assert default_micros[0].reactions == explicit_micros[0].reactions + + +def test_flux_mode_returns_one_group_flux(run_in_tmpdir): + model, mat = _simple_model() + fluxes, micros = get_microxs_and_flux( + model, [mat], + nuclides=['H1'], + reactions=['(n,2n)'], + energies=[0., 0.625, 2.0e7], + reaction_rate_mode='flux', + chain_file=CHAIN_FILE, + ) + + assert fluxes[0].shape == (1,) + assert micros[0].data.shape == (1, 1, 1) + assert fluxes[0][0] > 0.0 + # --------------------------------------------------------------------------- # Tests for MicroXS.merge() # ---------------------------------------------------------------------------