diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index c2d35565a8a..ddd60ae22d7 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -35,3 +35,4 @@ from .function import * from .effective_dose.dose import dose_coefficients +from .mass_energy_absorption import mu_en_coefficients diff --git a/openmc/data/data.py b/openmc/data/data.py index 5ecadd37be0..757540f983c 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -274,6 +274,7 @@ # Unit conversions EV_PER_MEV = 1.0e6 JOULE_PER_EV = 1.602176634e-19 +BARN_PER_CM_SQ = 1.0e24 # Avogadro's constant AVOGADRO = 6.02214076e23 diff --git a/openmc/data/endf.py b/openmc/data/endf.py index eca37446933..d93a022365c 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -59,7 +59,9 @@ 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/mass_energy_absorption.py b/openmc/data/mass_energy_absorption.py new file mode 100644 index 00000000000..e0a2e5692a6 --- /dev/null +++ b/openmc/data/mass_energy_absorption.py @@ -0,0 +1,93 @@ +import numpy as np + +import openmc.checkvalue as cv +from openmc.data import EV_PER_MEV + +# Embedded NIST-126 data +# Air (Dry Near Sea Level) — NIST Standard Reference Database 126 Table 4 (doi: 10.18434/T4D01F) +# Columns: Energy (MeV), μen/ρ (cm^2/g) +_NIST126_AIR = np.array( + [ + [1.00000e-03, 3.599e03], + [1.50000e-03, 1.188e03], + [2.00000e-03, 5.262e02], + [3.00000e-03, 1.614e02], + [3.20290e-03, 1.330e02], + [3.20290e-03, 1.460e02], + [4.00000e-03, 7.636e01], + [5.00000e-03, 3.931e01], + [6.00000e-03, 2.270e01], + [8.00000e-03, 9.446e00], + [1.00000e-02, 4.742e00], + [1.50000e-02, 1.334e00], + [2.00000e-02, 5.389e-01], + [3.00000e-02, 1.537e-01], + [4.00000e-02, 6.833e-02], + [5.00000e-02, 4.098e-02], + [6.00000e-02, 3.041e-02], + [8.00000e-02, 2.407e-02], + [1.00000e-01, 2.325e-02], + [1.50000e-01, 2.496e-02], + [2.00000e-01, 2.672e-02], + [3.00000e-01, 2.872e-02], + [4.00000e-01, 2.949e-02], + [5.00000e-01, 2.966e-02], + [6.00000e-01, 2.953e-02], + [8.00000e-01, 2.882e-02], + [1.00000e00, 2.789e-02], + [1.25000e00, 2.666e-02], + [1.50000e00, 2.547e-02], + [2.00000e00, 2.345e-02], + [3.00000e00, 2.057e-02], + [4.00000e00, 1.870e-02], + [5.00000e00, 1.740e-02], + [6.00000e00, 1.647e-02], + [8.00000e00, 1.525e-02], + [1.00000e01, 1.450e-02], + [1.50000e01, 1.353e-02], + [2.00000e01, 1.311e-02], + ], + dtype=float, +) + +# Registry of embedded tables: (data_source, material) -> ndarray +# Table shape: (N, 2) with columns [Energy (MeV), μen/ρ (cm^2/g)] +_MUEN_TABLES = { + ("nist126", "air"): _NIST126_AIR, +} + + +def mu_en_coefficients( + material: str, data_source: str = "nist126" +) -> tuple[np.ndarray, np.ndarray]: + """Return tabulated mass energy-absorption coefficients. + + Parameters + ---------- + material : {'air'} + Material compound for which to load coefficients. + data_source : {'nist126'} + Source library. + + Returns + ------- + energy : numpy.ndarray + Energies [eV] + mu_en_coeffs : numpy.ndarray + Mass energy-absorption coefficients [cm^2/g] + """ + cv.check_value("material", material, {"air"}) + cv.check_value("data_source", data_source, {"nist126"}) + + key = (data_source, material) + if key not in _MUEN_TABLES: + available = sorted({m for (ds, m) in _MUEN_TABLES.keys() if ds == data_source}) + raise ValueError( + f"'{material}' has no embedded μen/ρ table for data source {data_source}. " + f"Available materials for {data_source}: {available}" + ) + + data = _MUEN_TABLES[key] + energy = data[:, 0].copy() * EV_PER_MEV # MeV -> eV + mu_en_coeffs = data[:, 1].copy() + return energy, mu_en_coeffs diff --git a/openmc/data/photon.py b/openmc/data/photon.py index bc21b2e56a5..43aacf962b8 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 Evaluation, SUM_RULES, 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/material.py b/openmc/material.py index 735a0574326..87686db141d 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -2,6 +2,7 @@ from collections import defaultdict, namedtuple, Counter from collections.abc import Iterable from copy import deepcopy +from functools import reduce from numbers import Real from pathlib import Path import re @@ -22,8 +23,10 @@ from .utility_funcs import input_path from . import waste from openmc.checkvalue import PathLike -from openmc.stats import Univariate, Discrete, Mixture -from openmc.data.data import _get_element_symbol +from openmc.stats import Univariate, Discrete, Mixture, Tabular +from openmc.data.data import _get_element_symbol, BARN_PER_CM_SQ, JOULE_PER_EV +from openmc.data.function import Combination, Tabulated1D +from openmc.data import mu_en_coefficients, dose_coefficients # Units for density supported by OpenMC @@ -409,6 +412,234 @@ def get_decay_photon_energy( return combined + + + def get_photon_contact_dose_rate(self, dose_quantity:str = "absorbed-air", build_up:float = 2.0, by_nuclide: bool = False) -> float | dict[str, float]: + """Compute the photon contact dose rate (CDR) produced by radioactive decay + of the material. + + The contact dose rate is calculated from decay photon energy spectra for + each nuclide in the material, combined with photon mass attenuation data + for the material and the appropriate response function for the dose quantity. + A slab-geometry approximation and a photon build-up factor are used. + + Absorbed-air dose: + The approach follows the FISPACT-II manual (UKAEA-CCFE-RE(21)02 - May 2021). + Appendix C.7.1. + This method integrates over the photon energy: + + (B/2) * (mu_en_air(E) / mu_material(E)) * E * S(E) + + Effective dose: + The approach uses ICRP-116 effective dose coefficients to convert the photon + fluence due to decay photons to effective dose. + This method integrates over the photon energy: + + (B/2) * (h_e(E) / mu_material(E)) * S(E) + + where: + - mu_en_air(E) is the air mass energy-absorption coefficient, + - mu_material(E) is the photon mass attenuation coefficient of the material, + - S(E) is the photon emission spectrum per atom, + - h_e(E) is the ICRP-116 effective dose coefficient, + - B is the build-up factor, + - E is the photon energy. + + Parameters + ---------- + dose_quantity : {'absorbed-air', 'effective'}, optional + Specifies the dose quantity to be calculated. + The only supported options are 'absorbed-air' which implements the methodology + from FISPACT-II, and 'effective' which uses ICRP-116 effective dose coefficients. + build_up : float, optional. The default value is 2.0 as suggested in the FISPACT-II + manual. + by_nuclide : bool, optional + Specifies if the cdr should be returned for the material as a + whole or per nuclide. Default is False. + + Limitations + ---------- + This method does not implement correction from Bremsstrahlung particles which can be + relevant at close distances. + In addition, it computes the gamma contact dose rate only for the unstable nuclides + for which the radiation source specification is present in the chain file. + + Returns + ------- + cdr : float or dict[str, float] + Contact Dose Rate due to decay photons. + 'absorbed-air': returns the absorbed dose in air [Gy/hr]. + 'effective': returns the effective dose [Sv/hr]. + """ + + cv.check_type("by_nuclide", by_nuclide, bool) + cv.check_type("dose_quantity", dose_quantity, str) + cv.check_value("dose_quantity", dose_quantity, ['absorbed-air', 'effective']) + cv.check_type("build_up", build_up, Real) + cv.check_greater_than("build_up", build_up, 0.0) + + # photon linear attenuation distribution as a function of energy + # distribution values in [cm-1] + from openmc.plotter import _calculate_cexs_elem_mat + + mu_e_vals, cexs = _calculate_cexs_elem_mat( + this=self, + types=[501], + incident_particle="photon" + ) + mu_y_vals = np.array(cexs[0]) # total mass attenuation coeffs + if mu_y_vals is None: + raise ValueError("Cannot compute photon mass attenuation for material") + linear_attenuation_dist = Tabulated1D(mu_e_vals, mu_y_vals, breakpoints=[len(mu_e_vals)], interpolation=[5]) + + # CDR computation + cdr = {} + + geometry_factor_slab = 0.5 + + # ancillary conversion factors for clarity + seconds_per_hour = 3600.0 + grams_per_kg = 1000.0 + sv_per_psv = 1e-12 + + if dose_quantity == 'absorbed-air': + + # mu_en/ rho for air distribution, [eV, cm2/g] + response_f_x, response_f_y = mu_en_coefficients("air", data_source="nist126") + + # converts [eV cm2 barns-1 g-1 s-1] to [Gy hr-1] + multiplier = ( + build_up + * geometry_factor_slab + * seconds_per_hour + * grams_per_kg + * BARN_PER_CM_SQ + * JOULE_PER_EV + ) + + elif dose_quantity == 'effective': + + # effective dose as a function of photon fluence [pSv cm2] + response_f_x, response_f_y = dose_coefficients("photon", geometry='AP', data_source='icrp116') + + # converts [pSv cm2 barns-1 s-1] to [Sv hr-1] + multiplier = ( + build_up + * geometry_factor_slab + * seconds_per_hour + * sv_per_psv + * BARN_PER_CM_SQ + ) + + response_f = Tabulated1D(response_f_x, response_f_y, breakpoints=[len(response_f_x)], interpolation=[5]) + + for nuc, nuc_atoms_per_bcm in self.get_nuclide_atom_densities().items(): + + cdr_nuc = 0.0 + + photon_source_per_atom = openmc.data.decay_photon_energy(nuc) + + # nuclides with no contribution + if photon_source_per_atom is None or nuc_atoms_per_bcm <= 0.0: + cdr[nuc] = 0.0 + continue + + if isinstance(photon_source_per_atom, (Discrete, Tabular)): + e_vals = np.array(photon_source_per_atom.x) + p_vals = np.array(photon_source_per_atom.p) + + e_lists = [e_vals, response_f.x] + e_lists.append(mu_e_vals) + + # clip distributions for values outside the tabulated values + left_bound = max(a.min() for a in e_lists) + right_bound = min(a.max() for a in e_lists) + + mask = (e_vals >= left_bound) & (e_vals <= right_bound) + e_vals = e_vals[mask] + p_vals = p_vals[mask] + + else: + raise ValueError( + f"Unknown decay photon energy data type for nuclide {nuc}" + f"value returned: {type(photon_source_per_atom)}" + ) + + if isinstance(photon_source_per_atom, Discrete): + mu_vals = np.array(linear_attenuation_dist(e_vals)) + if np.any(mu_vals <= 0.0): + zero_vals = e_vals[mu_vals <= 0.0] + raise ValueError( + f"Mass attenuation coefficient <= 0 at energies: {zero_vals}" + ) + if dose_quantity == 'absorbed-air': + # units [eV cm3 g-1 atoms-1 s-1] + cdr_nuc += np.sum((response_f(e_vals) / mu_vals) * p_vals * e_vals) + elif dose_quantity == 'effective': + # units [pSv cm3 atoms-1 s-1] + cdr_nuc += np.sum((response_f(e_vals) / mu_vals) * p_vals) + + + elif isinstance(photon_source_per_atom, Tabular): + + # generate the tabulated1D functions + e_p_dist = Tabulated1D( + e_vals, p_vals, breakpoints=[len(e_vals)], interpolation=[1] + ) + + + # limit the computation to the tabulated mu_en_air range + e_union = reduce(np.union1d, e_lists) + e_union = e_union[(e_union >= left_bound) & (e_union <= right_bound)] + if len(e_union) < 2: + raise ValueError("Not enough overlapping energy points to compute CDR") + + # check for negative denominator valuenters + mu_vals_check = np.array(linear_attenuation_dist(e_union)) + if np.any(mu_vals_check <= 0.0): + zero_vals = e_union[mu_vals_check <= 0.0] + raise ValueError( + f"Mass attenuation coefficient <= 0 at energies: {zero_vals}" + ) + + if dose_quantity == 'absorbed-air': + # units [eV cm3 g-1 atoms-1 s-1] + e_e_dist = Tabulated1D( + e_vals, e_vals, breakpoints=[len(e_vals)], interpolation=[2] + ) + integrand_operator = Combination( + functions=[response_f, e_p_dist, e_e_dist, linear_attenuation_dist], + operations=[np.multiply, np.multiply, np.divide], + ) + elif dose_quantity == 'effective': + # units [pSv cm3 atoms-1 s-1] + integrand_operator = Combination( + functions=[response_f, e_p_dist, linear_attenuation_dist], + operations=[np.multiply, np.divide], + ) + + + y_evaluated = integrand_operator(e_union) + + integrand_function = Tabulated1D( + e_union, y_evaluated, breakpoints=[len(e_union)], interpolation=[2] + ) + + cdr_nuc += integrand_function.integral()[-1] + + + # units air-absorbed dose [eV cm2 barns-1 g-1 s-1] + # units effective dose [pSv cm2 barns-1 s-1] + cdr_nuc *= nuc_atoms_per_bcm + + # units air-absorbed dose [Gy hr-1] + # units effective dose [Sv hr-1] + cdr_nuc *= multiplier + + cdr[nuc] = cdr_nuc + + return cdr if by_nuclide else sum(cdr.values()) + @classmethod def from_hdf5(cls, group: h5py.Group) -> Material: """Create material from HDF5 group 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_data_mu_en_coefficients.py b/tests/unit_tests/test_data_mu_en_coefficients.py new file mode 100644 index 00000000000..cf4173f3265 --- /dev/null +++ b/tests/unit_tests/test_data_mu_en_coefficients.py @@ -0,0 +1,18 @@ +from pytest import approx, raises + +from openmc.data import mu_en_coefficients + + +def test_mu_en_coefficients(): + # Spot checks on values from NIST tables + energy, mu_en = mu_en_coefficients("air") + assert energy[0] == approx(1e3) + assert mu_en[0] == approx(3.599e3) + assert energy[-1] == approx(2e7) + assert mu_en[-1] == approx(1.311e-2) + + # Invalid particle/geometry should raise an exception + with raises(ValueError): + mu_en_coefficients("pasta") + with raises(ValueError): + mu_en_coefficients("air", data_source="nist000") diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 9aca8b59656..89f64f1c933 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -610,6 +610,7 @@ def test_mesh_get_homogenized_materials(): @pytest.fixture def sphere_model(): + openmc.reset_auto_ids() # Model with three materials separated by planes x=0 and z=0 mats = [] for i in range(3): diff --git a/tests/unit_tests/test_plotter.py b/tests/unit_tests/test_plotter.py index 0220cec3e71..aaf131bd9c3 100644 --- a/tests/unit_tests/test_plotter.py +++ b/tests/unit_tests/test_plotter.py @@ -1,9 +1,10 @@ import numpy as np -import openmc import pytest +import openmc + -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def test_mat(): mat_1 = openmc.Material() mat_1.add_element("H", 4.0, "ao") @@ -35,9 +36,7 @@ def test_calculate_cexs_elem_mat_sab(test_mat): @pytest.mark.parametrize("this", ["Li", "Li6"]) def test_calculate_cexs_with_nuclide_and_element(this): # single type (reaction) - energy_grid, data = openmc.plotter.calculate_cexs( - this=this, types=[205] - ) + energy_grid, data = openmc.plotter.calculate_cexs(this=this, types=[205]) assert isinstance(energy_grid, np.ndarray) assert isinstance(data, np.ndarray) @@ -46,9 +45,7 @@ def test_calculate_cexs_with_nuclide_and_element(this): assert len(data[0]) == len(energy_grid) # two types (reactions) - energy_grid, data = openmc.plotter.calculate_cexs( - this=this, types=[2, "elastic"] - ) + energy_grid, data = openmc.plotter.calculate_cexs(this=this, types=[2, "elastic"]) assert isinstance(energy_grid, np.ndarray) assert isinstance(data, np.ndarray) @@ -61,9 +58,7 @@ def test_calculate_cexs_with_nuclide_and_element(this): def test_calculate_cexs_with_materials(test_mat): - energy_grid, data = openmc.plotter.calculate_cexs( - this=test_mat, types=[205] - ) + energy_grid, data = openmc.plotter.calculate_cexs(this=test_mat, types=[205]) assert isinstance(energy_grid, np.ndarray) assert isinstance(data, np.ndarray) @@ -75,48 +70,55 @@ def test_calculate_cexs_with_materials(test_mat): @pytest.mark.parametrize("this", ["Be", "Be9"]) def test_plot_xs(this): from matplotlib.figure import Figure - assert isinstance(openmc.plot_xs({this: ['total', 'elastic', 16, '(n,2n)']}), Figure) + + assert isinstance( + openmc.plot_xs({this: ["total", "elastic", 16, "(n,2n)"]}), Figure + ) def test_plot_xs_mat(test_mat): from matplotlib.figure import Figure - assert isinstance(openmc.plot_xs({test_mat: ['total']}), Figure) + + assert isinstance(openmc.plot_xs({test_mat: ["total"]}), Figure) @pytest.mark.parametrize("units", ["eV", "keV", "MeV"]) def test_plot_xs_energy_axis(units): - plot = openmc.plot_xs({'Be9': ['(n,2n)']}, energy_axis_units=units) + plot = openmc.plot_xs({"Be9": ["(n,2n)"]}, energy_axis_units=units) axis_text = plot.get_axes()[0].get_xaxis().get_label().get_text() - assert axis_text == f'Energy [{units}]' + assert axis_text == f"Energy [{units}]" def test_plot_axes_labels(): # just nuclides axis_label = openmc.plotter._get_yaxis_label( reactions={ - 'Li6': [205], - 'Li7': [205], - }, divisor_types=False + "Li6": [205], + "Li7": [205], + }, + divisor_types=False, ) - assert axis_label == 'Microscopic Cross Section [b]' + assert axis_label == "Microscopic Cross Section [b]" # just elements axis_label = openmc.plotter._get_yaxis_label( reactions={ - 'Li': [205], - 'Be': [16], - }, divisor_types=False + "Li": [205], + "Be": [16], + }, + divisor_types=False, ) - assert axis_label == 'Microscopic Cross Section [b]' + assert axis_label == "Microscopic Cross Section [b]" # mixed nuclide and element axis_label = openmc.plotter._get_yaxis_label( reactions={ - 'Li': [205], - 'Li7': [205], - }, divisor_types=False + "Li": [205], + "Li7": [205], + }, + divisor_types=False, ) - assert axis_label == 'Microscopic Cross Section [b]' + assert axis_label == "Microscopic Cross Section [b]" axis_label = openmc.plotter._get_yaxis_label( reactions={ @@ -135,49 +137,233 @@ def test_plot_axes_labels(): # just materials mat1 = openmc.Material() - mat1.add_nuclide('Fe56', 1) - mat1.set_density('g/cm3', 1) + mat1.add_nuclide("Fe56", 1) + mat1.set_density("g/cm3", 1) mat2 = openmc.Material() - mat2.add_element('Fe', 1) - mat2.add_nuclide('Fe55', 1) - mat2.set_density('g/cm3', 1) + mat2.add_element("Fe", 1) + mat2.add_nuclide("Fe55", 1) + mat2.set_density("g/cm3", 1) axis_label = openmc.plotter._get_yaxis_label( reactions={ mat1: [205], mat2: [16], - }, divisor_types=False + }, + divisor_types=False, ) - assert axis_label == 'Macroscopic Cross Section [1/cm]' + assert axis_label == "Macroscopic Cross Section [1/cm]" # mixed materials and nuclides with pytest.raises(TypeError): openmc.plotter._get_yaxis_label( - reactions={'Li6': [205], mat2: [16]}, - divisor_types=False + reactions={"Li6": [205], mat2: [16]}, divisor_types=False ) # mixed materials and elements with pytest.raises(TypeError): openmc.plotter._get_yaxis_label( - reactions={'Li': [205], mat2: [16]}, - divisor_types=False + reactions={"Li": [205], mat2: [16]}, divisor_types=False ) def test_get_title(): - title = openmc.plotter._get_title(reactions={'Li': [205]}) - assert title == 'Cross Section Plot For Li' - title = openmc.plotter._get_title(reactions={'Li6': [205]}) - assert title == 'Cross Section Plot For Li6' - title = openmc.plotter._get_title(reactions={ - 'Li6': [205], - 'Li7': [205] - }) - assert title == 'Cross Section Plot' + title = openmc.plotter._get_title(reactions={"Li": [205]}) + assert title == "Cross Section Plot For Li" + title = openmc.plotter._get_title(reactions={"Li6": [205]}) + assert title == "Cross Section Plot For Li6" + title = openmc.plotter._get_title(reactions={"Li6": [205], "Li7": [205]}) + assert title == "Cross Section Plot" mat1 = openmc.Material() - mat1.add_nuclide('Fe56', 1) - mat1.set_density('g/cm3', 1) - mat1.name = 'my_mat' + mat1.add_nuclide("Fe56", 1) + mat1.set_density("g/cm3", 1) + mat1.name = "my_mat" title = openmc.plotter._get_title(reactions={mat1: [205]}) - assert title == 'Cross Section Plot For my_mat' + 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 _any_photon_mt(element_symbol="C", cross_sections=None): + """Pick an MT that actually exists in the configured photon 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"]) + return next(iter(inc.reactions.keys())) + + +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)