From 44756bc2f8c2e1f08ad068cfe56cbb739d8dcbbd Mon Sep 17 00:00:00 2001 From: marco-de-pietri Date: Thu, 29 Jan 2026 16:01:39 -0500 Subject: [PATCH 1/2] initial pr split --- openmc/data/endf.py | 5 +- openmc/data/photon.py | 26 ++- openmc/plotter.py | 317 +++++++++++++++++-------------- tests/unit_tests/test_mesh.py | 1 + tests/unit_tests/test_plotter.py | 172 +++++++++++++++++ 5 files changed, 381 insertions(+), 140 deletions(-) diff --git a/openmc/data/endf.py b/openmc/data/endf.py index eca37446933..d9d337f55f6 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -59,7 +59,10 @@ 104: list(range(650, 700)), 105: list(range(700, 750)), 106: list(range(750, 800)), - 107: list(range(800, 850))} + 107: list(range(800, 850)), + 501: [502, 504, 516, 522], + 516: [515, 517] + } ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]) ?(\d+)') diff --git a/openmc/data/photon.py b/openmc/data/photon.py index bc21b2e56a5..49b3bfba28e 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -15,7 +15,7 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR from .ace import Table, get_metadata, get_table from .data import ATOMIC_SYMBOL, EV_PER_MEV -from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record +from .endf import SUM_RULES, Evaluation, get_head_record, get_tab1_record, get_list_record from .function import Tabulated1D @@ -487,6 +487,30 @@ def atomic_relaxation(self, atomic_relaxation): def name(self): return ATOMIC_SYMBOL[self.atomic_number] + def get_reaction_components(self, mt): + """Determine what reactions make up redundant reaction. + + Parameters + ---------- + mt : int + ENDF MT number of the reaction to find components of. + + Returns + ------- + mts : list of int + ENDF MT numbers of reactions that make up the redundant reaction and + have cross sections provided. + + """ + mts = [] + if mt in SUM_RULES: + for mt_i in SUM_RULES[mt]: + mts += self.get_reaction_components(mt_i) + if mts: + return mts + else: + return [mt] if mt in self else [] + @classmethod def from_ace(cls, ace_or_filename): """Generate incident photon data from an ACE table diff --git a/openmc/plotter.py b/openmc/plotter.py index abd8ab6dd4b..06cb10f17ba 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -128,6 +128,7 @@ def plot_xs( temperature: float = 294.0, axis: "plt.Axes" | None = None, sab_name: str | None = None, + incident_particle: str = 'neutron', ce_cross_sections: str | None = None, mg_cross_sections: str | None = None, enrichment: float | None = None, @@ -215,11 +216,11 @@ def plot_xs( if plot_CE: cv.check_type("this", this, (str, openmc.Material)) # Calculate for the CE cross sections - E, data = calculate_cexs(this, types, temperature, sab_name, + E, data = calculate_cexs(this, types, incident_particle,temperature, sab_name, ce_cross_sections, enrichment) if divisor_types: cv.check_length('divisor types', divisor_types, len(types)) - Ediv, data_div = calculate_cexs(this, divisor_types, temperature, + Ediv, data_div = calculate_cexs(this, divisor_types, incident_particle, temperature, sab_name, ce_cross_sections, enrichment) @@ -288,7 +289,7 @@ def plot_xs( return fig -def calculate_cexs(this, types, temperature=294., sab_name=None, +def calculate_cexs(this, types, incident_particle='neutron', temperature=294., sab_name=None, cross_sections=None, enrichment=None, ncrystal_cfg=None): """Calculates continuous-energy cross sections of a requested type. @@ -299,6 +300,9 @@ def calculate_cexs(this, types, temperature=294., sab_name=None, str types : Iterable of values of PLOT_TYPES The type of cross sections to calculate + incident_particle : str + The incident particle used to fetch the appropriate library. + Can be only 'neutron' or 'photon'. temperature : float, optional Temperature in Kelvin to plot. If not specified, a default temperature of 294K will be plotted. Note that the nearest @@ -327,6 +331,7 @@ def calculate_cexs(this, types, temperature=294., sab_name=None, # Check types cv.check_type('this', this, (str, openmc.Material)) cv.check_type('temperature', temperature, Real) + cv.check_value("incident particle", incident_particle, ['neutron', 'photon']) if sab_name: cv.check_type('sab_name', sab_name, str) if enrichment: @@ -335,11 +340,11 @@ def calculate_cexs(this, types, temperature=294., sab_name=None, if isinstance(this, str): if this in ELEMENT_NAMES: energy_grid, data = _calculate_cexs_elem_mat( - this, types, temperature, cross_sections, sab_name, enrichment + this, types, incident_particle, temperature, cross_sections, sab_name, enrichment ) else: energy_grid, xs = _calculate_cexs_nuclide( - this, types, temperature, sab_name, cross_sections, + this, types, incident_particle, temperature, sab_name, cross_sections, ncrystal_cfg ) @@ -350,13 +355,13 @@ def calculate_cexs(this, types, temperature=294., sab_name=None, for line in range(len(types)): data[line, :] = xs[line](energy_grid) else: - energy_grid, data = _calculate_cexs_elem_mat(this, types, temperature, + energy_grid, data = _calculate_cexs_elem_mat(this, types, incident_particle, temperature, cross_sections) return energy_grid, data -def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, +def _calculate_cexs_nuclide(this, types, incident_particle='neutron', temperature=294., sab_name=None, cross_sections=None, ncrystal_cfg=None): """Calculates continuous-energy cross sections of a requested type. @@ -369,6 +374,9 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, in openmc.PLOT_TYPES or keys from openmc.data.REACTION_MT which correspond to a reaction description e.g '(n,2n)' or integers which correspond to reaction channel (MT) numbers. + incident_particle : str + The incident particle used to fetch the appropriate library. + Can be only 'neutron' or 'photon'. temperature : float, optional Temperature in Kelvin to plot. If not specified, a default temperature of 294K will be plotted. Note that the nearest @@ -392,6 +400,19 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, # Load the library library = openmc.data.DataLibrary.from_xml(cross_sections) + if incident_particle == 'photon': + try: + z = openmc.data.zam(this)[0] + nuclide = openmc.data.ATOMIC_SYMBOL[z] + except (ValueError, KeyError, TypeError): + if this not in openmc.data.ELEMENT_SYMBOL.values(): + raise ValueError(f"Element '{this}' not found in ELEMENT_SYMBOL.") + nuclide = this + else: + nuclide = this + lib = library.get_by_material(nuclide, data_type=incident_particle) + if lib is None: + raise ValueError(this + " not in library") # Convert temperature to format needed for access in the library strT = f"{int(round(temperature))}K" @@ -400,8 +421,8 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, # Now we can create the data sets to be plotted energy_grid = [] xs = [] - lib = library.get_by_material(this) - if lib is not None: + + if incident_particle == 'neutron': nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) # Obtain the nearest temperature if strT in nuc.temperatures: @@ -442,133 +463,147 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, inelastic = sab.inelastic.xs[sabT] grid = np.union1d(grid, inelastic.x) if inelastic.x[-1] > sab_Emax: - sab_Emax = inelastic.x[-1] + sab_Emax = inelastic.x[-1] sab_funcs.append(inelastic) energy_grid = grid else: energy_grid = nuc.energy[nucT] - - # Parse the types - mts = [] - ops = [] - yields = [] - for line in types: - if line in PLOT_TYPES: - tmp_mts = [mtj for mti in PLOT_TYPES_MT[line] for mtj in - nuc.get_reaction_components(mti)] - mts.append(tmp_mts) - if line.startswith('nu'): - yields.append(True) - else: - yields.append(False) - if XI_MT in tmp_mts: - ops.append((np.add,) * (len(tmp_mts) - 2) + (np.multiply,)) - else: - ops.append((np.add,) * (len(tmp_mts) - 1)) - elif line in openmc.data.REACTION_MT: - mt_number = openmc.data.REACTION_MT[line] - cv.check_type('MT in types', mt_number, Integral) - cv.check_greater_than('MT in types', mt_number, 0) - tmp_mts = nuc.get_reaction_components(mt_number) - mts.append(tmp_mts) - ops.append((np.add,) * (len(tmp_mts) - 1)) - yields.append(False) - elif isinstance(line, int): - # Not a built-in type, we have to parse it ourselves - cv.check_type('MT in types', line, Integral) - cv.check_greater_than('MT in types', line, 0) - tmp_mts = nuc.get_reaction_components(line) - mts.append(tmp_mts) - ops.append((np.add,) * (len(tmp_mts) - 1)) + if incident_particle == 'photon': + nuc = openmc.data.IncidentPhoton.from_hdf5(lib['path']) + if any(type(line) is not int for line in types): + raise TypeError("Photon cross sections can only be requested " + "with integer MT numbers.") + + + # Parse the types + mts = [] + ops = [] + yields = [] + for line in types: + if line in PLOT_TYPES: + tmp_mts = [mtj for mti in PLOT_TYPES_MT[line] for mtj in + nuc.get_reaction_components(mti)] + mts.append(tmp_mts) + if line.startswith('nu'): + yields.append(True) + else: yields.append(False) + if XI_MT in tmp_mts: + ops.append((np.add,) * (len(tmp_mts) - 2) + (np.multiply,)) else: - raise TypeError("Invalid type", line) - - for i, mt_set in enumerate(mts): - # Get the reaction xs data from the nuclide - funcs = [] - op = ops[i] - for mt in mt_set: - if mt == 2: - if sab_name: - # Then we need to do a piece-wise function of - # The S(a,b) and non-thermal data - sab_sum = openmc.data.Sum(sab_funcs) - pw_funcs = openmc.data.Regions1D( - [sab_sum, nuc[mt].xs[nucT]], - [sab_Emax]) - funcs.append(pw_funcs) - elif ncrystal_cfg: - import NCrystal - nc_scatter = NCrystal.createScatter(ncrystal_cfg) - nc_func = nc_scatter.xsect - nc_emax = 5 # eV # this should be obtained from NCRYSTAL_MAX_ENERGY - energy_grid = np.union1d(np.geomspace(min(energy_grid), - 1.1*nc_emax, - 1000),energy_grid) # NCrystal does not have - # an intrinsic energy grid - pw_funcs = openmc.data.Regions1D( - [nc_func, nuc[mt].xs[nucT]], - [nc_emax]) - funcs.append(pw_funcs) + ops.append((np.add,) * (len(tmp_mts) - 1)) + elif line in openmc.data.REACTION_MT: + mt_number = openmc.data.REACTION_MT[line] + cv.check_type('MT in types', mt_number, Integral) + cv.check_greater_than('MT in types', mt_number, 0) + tmp_mts = nuc.get_reaction_components(mt_number) + mts.append(tmp_mts) + ops.append((np.add,) * (len(tmp_mts) - 1)) + yields.append(False) + elif isinstance(line, int): + # Not a built-in type, we have to parse it ourselves + cv.check_type('MT in types', line, Integral) + cv.check_greater_than('MT in types', line, 0) + tmp_mts = nuc.get_reaction_components(line) + mts.append(tmp_mts) + ops.append((np.add,) * (len(tmp_mts) - 1)) + yields.append(False) + else: + raise TypeError("Invalid type", line) + + for i, mt_set in enumerate(mts): + # Get the reaction xs data from the nuclide + funcs = [] + op = ops[i] + for mt in mt_set: + if mt == 2: + if sab_name: + # Then we need to do a piece-wise function of + # The S(a,b) and non-thermal data + sab_sum = openmc.data.Sum(sab_funcs) + pw_funcs = openmc.data.Regions1D( + [sab_sum, nuc[mt].xs[nucT]], + [sab_Emax]) + funcs.append(pw_funcs) + elif ncrystal_cfg: + import NCrystal + nc_scatter = NCrystal.createScatter(ncrystal_cfg) + nc_func = nc_scatter.xsect + nc_emax = 5 # eV # this should be obtained from NCRYSTAL_MAX_ENERGY + energy_grid = np.union1d(np.geomspace(min(energy_grid), + 1.1*nc_emax, + 1000),energy_grid) # NCrystal does not have + # an intrinsic energy grid + pw_funcs = openmc.data.Regions1D( + [nc_func, nuc[mt].xs[nucT]], + [nc_emax]) + funcs.append(pw_funcs) + else: + funcs.append(nuc[mt].xs[nucT]) + elif mt in nuc: + if yields[i]: + # Get the total yield first if available. This will be + # used primarily for fission. + for prod in chain(nuc[mt].products, + nuc[mt].derived_products): + if prod.particle == 'neutron' and \ + prod.emission_mode == 'total': + func = openmc.data.Combination( + [nuc[mt].xs[nucT], prod.yield_], + [np.multiply]) + funcs.append(func) + break else: - funcs.append(nuc[mt].xs[nucT]) - elif mt in nuc: - if yields[i]: - # Get the total yield first if available. This will be - # used primarily for fission. + # Total doesn't exist so we have to create from + # prompt and delayed. This is used for scatter + # multiplication. + func = None for prod in chain(nuc[mt].products, - nuc[mt].derived_products): + nuc[mt].derived_products): if prod.particle == 'neutron' and \ - prod.emission_mode == 'total': - func = openmc.data.Combination( - [nuc[mt].xs[nucT], prod.yield_], - [np.multiply]) - funcs.append(func) - break + prod.emission_mode != 'total': + if func: + func = openmc.data.Combination( + [prod.yield_, func], [np.add]) + else: + func = prod.yield_ + if func: + funcs.append(openmc.data.Combination( + [func, nuc[mt].xs[nucT]], [np.multiply])) else: - # Total doesn't exist so we have to create from - # prompt and delayed. This is used for scatter - # multiplication. - func = None - for prod in chain(nuc[mt].products, - nuc[mt].derived_products): - if prod.particle == 'neutron' and \ - prod.emission_mode != 'total': - if func: - func = openmc.data.Combination( - [prod.yield_, func], [np.add]) - else: - func = prod.yield_ - if func: - funcs.append(openmc.data.Combination( - [func, nuc[mt].xs[nucT]], [np.multiply])) - else: - # If func is still None, then there were no - # products. In that case, assume the yield is - # one as its not provided for some summed - # reactions like MT=4 - funcs.append(nuc[mt].xs[nucT]) - else: - funcs.append(nuc[mt].xs[nucT]) - elif mt == UNITY_MT: - funcs.append(lambda x: 1.) - elif mt == XI_MT: - awr = nuc.atomic_weight_ratio - alpha = ((awr - 1.) / (awr + 1.))**2 - xi = 1. + alpha * np.log(alpha) / (1. - alpha) - funcs.append(lambda x: xi) + # If func is still None, then there were no + # products. In that case, assume the yield is + # one as its not provided for some summed + # reactions like MT=4 + funcs.append(nuc[mt].xs[nucT]) else: - funcs.append(lambda x: 0.) - funcs = funcs if funcs else [lambda x: 0.] - xs.append(openmc.data.Combination(funcs, op)) - else: - raise ValueError(this + " not in library") + # general MT that can called with + # photons or neutrons + if incident_particle == 'photon': + temp_xs = nuc[mt].xs + energy_grid = np.union1d(energy_grid, temp_xs.x) + if incident_particle == 'neutron': + temp_xs = nuc[mt].xs[nucT] + funcs.append(temp_xs) + elif mt == UNITY_MT: + funcs.append(lambda x: 1.) + elif mt == XI_MT: + awr = nuc.atomic_weight_ratio + alpha = ((awr - 1.) / (awr + 1.))**2 + xi = 1. + alpha * np.log(alpha) / (1. - alpha) + funcs.append(lambda x: xi) + else: + funcs.append(lambda x: 0.) + funcs = funcs if funcs else [lambda x: 0.] + xs.append(openmc.data.Combination(funcs, op)) + + if len(energy_grid) == 0: + energy_grid = np.array([_MIN_E, _MAX_E], dtype=float) return energy_grid, xs -def _calculate_cexs_elem_mat(this, types, temperature=294., +def _calculate_cexs_elem_mat(this, types, incident_particle='neutron', temperature=294., cross_sections=None, sab_name=None, enrichment=None): """Calculates continuous-energy cross sections of a requested type. @@ -579,6 +614,9 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., Object to source data from. Element can be input as str types : Iterable of values of PLOT_TYPES The type of cross sections to calculate + incident_particle : str + The incident particle used to fetch the appropriate library. + Can be only 'neutron' or 'photon'. temperature : float, optional Temperature in Kelvin to plot. If not specified, a default temperature of 294K will be plotted. Note that the nearest @@ -602,6 +640,8 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., """ + cv.check_value("incident particle", incident_particle, ['neutron', 'photon']) + if isinstance(this, openmc.Material): if this.temperature is not None: T = this.temperature @@ -632,22 +672,23 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., # with a common nuclides format between openmc.Material and Elements nuclides = {nuclide[0]: nuclide[0] for nuclide in nuclides} - # Identify the nuclides which have S(a,b) data sabs = {} - for nuclide in nuclides.items(): - sabs[nuclide[0]] = None - if isinstance(this, openmc.Material): - for sab_name, _ in this._sab: - sab = openmc.data.ThermalScattering.from_hdf5( - library.get_by_material(sab_name, data_type='thermal')['path']) - for nuc in sab.nuclides: - sabs[nuc] = sab_name - else: - if sab_name: - sab = openmc.data.ThermalScattering.from_hdf5( - library.get_by_material(sab_name, data_type='thermal')['path']) - for nuc in sab.nuclides: - sabs[nuc] = sab_name + if incident_particle == 'neutron': + # Identify the nuclides which have S(a,b) data + for nuclide in nuclides.items(): + sabs[nuclide[0]] = None + if isinstance(this, openmc.Material): + for mat_sab_name, _ in this._sab: + sab = openmc.data.ThermalScattering.from_hdf5( + library.get_by_material(mat_sab_name, data_type='thermal')['path']) + for nuc in sab.nuclides: + sabs[nuc] = mat_sab_name + else: + if sab_name: + sab = openmc.data.ThermalScattering.from_hdf5( + library.get_by_material(sab_name, data_type='thermal')['path']) + for nuc in sab.nuclides: + sabs[nuc] = sab_name # Now we can create the data sets to be plotted xs = {} @@ -655,8 +696,8 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., for nuclide in nuclides.items(): name = nuclide[0] nuc = nuclide[1] - sab_name = sabs[name] - temp_E, temp_xs = calculate_cexs(nuc, types, T, sab_name, cross_sections, + nuc_sab_name = sabs.get(name) + temp_E, temp_xs = calculate_cexs(nuc, types, incident_particle, T, nuc_sab_name, cross_sections, ncrystal_cfg=ncrystal_cfg ) E.append(temp_E) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index c5855a7b051..70136cd06f8 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -613,6 +613,7 @@ def test_mesh_get_homogenized_materials(): @pytest.fixture def sphere_model(): # Model with three materials separated by planes x=0 and z=0 + openmc.reset_auto_ids() mats = [] for i in range(3): mat = openmc.Material() diff --git a/tests/unit_tests/test_plotter.py b/tests/unit_tests/test_plotter.py index 0220cec3e71..4a821d7ea6e 100644 --- a/tests/unit_tests/test_plotter.py +++ b/tests/unit_tests/test_plotter.py @@ -181,3 +181,175 @@ def test_get_title(): mat1.name = 'my_mat' title = openmc.plotter._get_title(reactions={mat1: [205]}) assert title == 'Cross Section Plot For my_mat' + + +def _any_photon_mt(element_symbol, cross_sections=None): + """Return a photon MT that is guaranteed to exist for the given element + in the configured cross sections library. + """ + if cross_sections is None: + cross_sections = openmc.config.get("cross_sections") + + library = openmc.data.DataLibrary.from_xml(cross_sections) + lib = library.get_by_material(element_symbol, data_type="photon") + if lib is None: + raise RuntimeError(f"No photon library entry found for {element_symbol}") + + inc = openmc.data.IncidentPhoton.from_hdf5(lib["path"]) + # `reactions` is a dict keyed by MT + return next(iter(inc.reactions.keys())) + +@pytest.mark.parametrize("this", ["Be", "Be9"]) +def test_calculate_cexs_photon_with_element_and_nuclide(this): + mt = _any_photon_mt("Be") + + # Use a common photoatomic MT (total) and verify basic shape/types + energy_grid, data = openmc.plotter.calculate_cexs( + this=this, types=[mt], incident_particle="photon" + ) + + assert isinstance(energy_grid, np.ndarray) + assert isinstance(data, np.ndarray) + assert len(energy_grid) > 1 + assert len(data) == 1 + assert len(data[0]) == len(energy_grid) + + +def test_calculate_cexs_photon_requires_integer_mts(): + # Photon cross sections can only be requested with integer MT numbers + with pytest.raises(TypeError): + openmc.plotter.calculate_cexs( + this="Be", types=["total"], incident_particle="photon" + ) + + with pytest.raises(TypeError): + openmc.plotter.calculate_cexs( + this="Be", types=[502, "elastic"], incident_particle="photon" + ) + + +def test_calculate_cexs_photon_with_material(): + mat = openmc.Material() + mat.add_element("Be", 1.0, "ao") + mat.set_density("g/cm3", 1.85) + + mt = _any_photon_mt("Be") + + energy_grid, data = openmc.plotter.calculate_cexs( + this=mat, types=[mt], incident_particle="photon" + ) + + assert isinstance(energy_grid, np.ndarray) + assert isinstance(data, np.ndarray) + assert len(energy_grid) > 1 + assert len(data) == 1 + assert len(data[0]) == len(energy_grid) + +def test_calculate_cexs_photon_material_element_vs_explicit_natural_abundance(): + mt = _any_photon_mt("C") + + # Material 1: defined as a single element (uses natural abundance implicitly) + mat_elem = openmc.Material() + mat_elem.add_element("C", 1.0, "ao") + mat_elem.set_density("g/cm3", 1.0) + + # Material 2: defined by explicitly specifying natural isotopic abundance + # (values are standard natural abundances for carbon) + mat_iso = openmc.Material() + mat_iso.add_nuclide("C12", 0.988922, "ao") + mat_iso.add_nuclide("C13", 0.011078, "ao") + mat_iso.set_density("g/cm3", 1.0) + + E1, xs1 = openmc.plotter.calculate_cexs( + this=mat_elem, types=[mt], incident_particle="photon" + ) + E2, xs2 = openmc.plotter.calculate_cexs( + this=mat_iso, types=[mt], incident_particle="photon" + ) + + assert isinstance(E1, np.ndarray) + assert isinstance(E2, np.ndarray) + assert isinstance(xs1, np.ndarray) + assert isinstance(xs2, np.ndarray) + + assert len(E1) > 1 + assert len(E2) > 1 + assert len(xs1) == 1 + assert len(xs2) == 1 + assert len(xs1[0]) == len(E1) + assert len(xs2[0]) == len(E2) + + # For photon data, isotopes map to the same element library, so these should match. + assert np.array_equal(E1, E2) + assert np.allclose(xs1[0], xs2[0], rtol=1e-12, atol=0.0) + + +def test_calculate_cexs_photon_missing_mt_fallback(): + # Use an MT that should never exist in photon data + energy_grid, data = openmc.plotter.calculate_cexs( + this="Be", types=[9999], incident_particle="photon" + ) + + assert isinstance(energy_grid, np.ndarray) + assert isinstance(data, np.ndarray) + assert np.allclose(energy_grid, [openmc.plotter._MIN_E, openmc.plotter._MAX_E]) + assert data.shape == (1, 2) + assert np.allclose(data, 0.0) + + +def test_calculate_cexs_photon_total_attenuation_reference_values(): + """Check total photon interaction XS for Pb and V at two reference energies. + + Total interaction is approximated by summing MTs: 502, 504, 515, 517, 522. + Reference mass attenuation data from NIST. + """ + openmc.reset_auto_ids() + + # Total interaction channels (library must contain these to run the test) + types = [501] + energies = np.array([1.0e5, 1.0e6]) # eV + # data from https://physics.nist.gov/PhysRefData/XrayMassCoef/ElemTab/z23.html + v_density = 6.11 # g/cm3 + v_ref = np.array( + [ + 2.877e-01, + 5.794e-02, + ] + ) # cm2/g + v_expected = v_ref * v_density + + # data from https://physics.nist.gov/PhysRefData/XrayMassCoef/ElemTab/z82.html + pb_density = 11.35 # g/cm3 + pb_ref = np.array( + [ + 5.549e00, + 7.102e-02, + ] + ) # cm2/g + pb_expected = pb_ref * pb_density + + def _run_element(symbol: str, density: float): + mat = openmc.Material() + mat.add_element(symbol, 1.0) + mat.set_density("g/cm3", density) + + # Compute macroscopic total XS for the material + e_grid, data = openmc.plotter.calculate_cexs( + this=mat, types=types, incident_particle="photon" + ) + xs_grid = data.sum(axis=0) + tab = openmc.data.Tabulated1D(e_grid, xs_grid, [len(e_grid)], [5]) + xs_mat_eval = tab(energies) + + return xs_mat_eval + + try: + pb_vals = _run_element("Pb", pb_density) + v_vals = _run_element("V", v_density) + except Exception: + pytest.skip( + "Pb or V photon data / required MTs not available in cross section library." + ) + + assert np.allclose(pb_vals, pb_expected, rtol=1e-2, atol=1e-8) + assert np.allclose(v_vals, v_expected, rtol=1e-2, atol=1e-8) From 69ea7015d4066f49cd36119dae4af06eb638f567 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 30 Jan 2026 10:24:33 +0200 Subject: [PATCH 2/2] add some fixes --- openmc/data/endf.py | 3 ++- openmc/plotter.py | 2 +- tests/unit_tests/test_plotter.py | 12 ++++-------- 3 files changed, 7 insertions(+), 10 deletions(-) diff --git a/openmc/data/endf.py b/openmc/data/endf.py index d9d337f55f6..7f92cd3e1d7 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -61,7 +61,8 @@ 106: list(range(750, 800)), 107: list(range(800, 850)), 501: [502, 504, 516, 522], - 516: [515, 517] + 516: [515, 517], + 522: list(range(534, 573)), } ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]) ?(\d+)') diff --git a/openmc/plotter.py b/openmc/plotter.py index 06cb10f17ba..c7f894b5efe 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -468,7 +468,7 @@ def _calculate_cexs_nuclide(this, types, incident_particle='neutron', temperatur energy_grid = grid else: energy_grid = nuc.energy[nucT] - if incident_particle == 'photon': + elif incident_particle == 'photon': nuc = openmc.data.IncidentPhoton.from_hdf5(lib['path']) if any(type(line) is not int for line in types): raise TypeError("Photon cross sections can only be requested " diff --git a/tests/unit_tests/test_plotter.py b/tests/unit_tests/test_plotter.py index 4a821d7ea6e..db85752ffa3 100644 --- a/tests/unit_tests/test_plotter.py +++ b/tests/unit_tests/test_plotter.py @@ -343,13 +343,9 @@ def _run_element(symbol: str, density: float): return xs_mat_eval - try: - pb_vals = _run_element("Pb", pb_density) - v_vals = _run_element("V", v_density) - except Exception: - pytest.skip( - "Pb or V photon data / required MTs not available in cross section library." - ) - + + pb_vals = _run_element("Pb", pb_density) + v_vals = _run_element("V", v_density) + assert np.allclose(pb_vals, pb_expected, rtol=1e-2, atol=1e-8) assert np.allclose(v_vals, v_expected, rtol=1e-2, atol=1e-8)