From eed5dae4f32e325ba22bedf6069a6d4b23299a83 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Aug 2026 08:55:14 -0500 Subject: [PATCH 01/10] Initial support for photon MGXS --- docs/source/pythonapi/mgxs.rst | 1 + openmc/mgxs/library.py | 104 ++++++++++++- openmc/mgxs/mgxs.py | 137 +++++++++++++++++- openmc/mgxs_library.py | 73 +++++++++- .../mgxs_library_condense/test.py | 5 +- .../mgxs_library_distribcell/test.py | 5 +- .../mgxs_library_hdf5/test.py | 5 +- .../mgxs_library_mesh/test.py | 5 +- .../mgxs_library_no_nuclides/test.py | 5 +- .../mgxs_library_nuclides/test.py | 5 +- .../mgxs_library_specific_nuclides/test.py | 5 +- tests/unit_tests/test_photon_mgxs.py | 115 +++++++++++++++ 12 files changed, 443 insertions(+), 22 deletions(-) create mode 100644 tests/unit_tests/test_photon_mgxs.py diff --git a/docs/source/pythonapi/mgxs.rst b/docs/source/pythonapi/mgxs.rst index 4141aa0a09a..bda07adf35d 100644 --- a/docs/source/pythonapi/mgxs.rst +++ b/docs/source/pythonapi/mgxs.rst @@ -51,6 +51,7 @@ Multi-group Cross Sections openmc.mgxs.KappaFissionXS openmc.mgxs.MultiplicityMatrixXS openmc.mgxs.NuFissionMatrixXS + openmc.mgxs.PhotonProductionMatrixXS openmc.mgxs.ReducedAbsorptionXS openmc.mgxs.ScatterXS openmc.mgxs.ScatterMatrixXS diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index faa83c0481f..19da8d7cad5 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -35,6 +35,9 @@ class Library: A geometry which has been initialized with a root universe by_nuclide : bool If true, computes cross sections for each nuclide in each domain + particle_type : {'neutron', 'photon'}, optional + Particle type for which cross sections are computed. If not specified, + tallies are not filtered by particle type. mgxs_types : Iterable of str The types of cross sections in the library (e.g., ['total', 'scatter']) name : str, optional @@ -47,6 +50,8 @@ class Library: An geometry which has been initialized with a root universe by_nuclide : bool If true, computes cross sections for each nuclide in each domain + particle_type : openmc.ParticleType or None + Particle type for which cross sections are computed mgxs_types : Iterable of str The types of cross sections in the library (e.g., ['total', 'scatter']) domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} @@ -102,11 +107,12 @@ class Library: """ def __init__(self, geometry, by_nuclide=False, - mgxs_types=None, name=''): + mgxs_types=None, name='', particle_type=None): self._name = '' self._geometry = None self._by_nuclide = None + self._particle_type = None self._mgxs_types = [] self._domain_type = None self._domains = 'all' @@ -129,6 +135,8 @@ def __init__(self, geometry, by_nuclide=False, self.name = name self.geometry = geometry self.by_nuclide = by_nuclide + if particle_type is not None: + self.particle_type = particle_type if mgxs_types is not None: self.mgxs_types = mgxs_types @@ -142,6 +150,7 @@ def __deepcopy__(self, memo): clone._name = self.name clone._geometry = self.geometry clone._by_nuclide = self.by_nuclide + clone._particle_type = self.particle_type clone._mgxs_types = self.mgxs_types clone._domain_type = self.domain_type clone._domains = copy.deepcopy(self.domains) @@ -203,7 +212,13 @@ def mgxs_types(self, mgxs_types): openmc.mgxs.ARBITRARY_VECTOR_TYPES + \ openmc.mgxs.ARBITRARY_MATRIX_TYPES if mgxs_types == 'all': - self._mgxs_types = all_mgxs_types + if self.particle_type == openmc.ParticleType.PHOTON: + self._mgxs_types = ( + 'total', 'absorption', 'photon-production matrix') + else: + self._mgxs_types = tuple( + xs for xs in all_mgxs_types + if xs != 'photon-production matrix') else: cv.check_iterable_type('mgxs_types', mgxs_types, str) for mgxs_type in mgxs_types: @@ -224,6 +239,18 @@ def by_nuclide(self, by_nuclide): self._by_nuclide = by_nuclide + @property + def particle_type(self): + return self._particle_type + + @particle_type.setter + def particle_type(self, particle_type): + particle_type = openmc.ParticleType(particle_type) + cv.check_value('particle type', particle_type, + (openmc.ParticleType.NEUTRON, + openmc.ParticleType.PHOTON)) + self._particle_type = particle_type + @property def domain_type(self): return self._domain_type @@ -513,13 +540,15 @@ def build_library(self): else: mgxs = openmc.mgxs.MGXS.get_mgxs( mgxs_type, name=self.name, num_polar=self.num_polar, - num_azimuthal=self.num_azimuthal) + num_azimuthal=self.num_azimuthal, + particle_type=self.particle_type) mgxs.domain = domain mgxs.domain_type = self.domain_type mgxs.energy_groups = self.energy_groups mgxs.by_nuclide = self.by_nuclide - if self.estimator is not None: + if self.estimator is not None and not isinstance( + mgxs, openmc.mgxs.PhotonProductionMatrixXS): mgxs.estimator = self.estimator if mgxs_type in openmc.mgxs.MDGXS_TYPES: @@ -1194,14 +1223,27 @@ def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', xsdata.set_decay_rate_mgxs(mymgxs, temperature=temperature, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) + # Photon production includes both the surviving primary photon and all + # banked secondary photons. Its multiplicity is therefore already + # folded into the matrix values. + if 'photon-production matrix' in self.mgxs_types: + production = self.get_mgxs( + domain, 'photon-production matrix') + xsdata.set_photon_production_mgxs( + production, temperature=temperature, xs_type=xs_type, + subdomain=subdomain) + using_photon_production = True + using_multiplicity = False + # If multiplicity matrix is available, prefer that - if 'multiplicity matrix' in self.mgxs_types: + elif 'multiplicity matrix' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'multiplicity matrix') xsdata.set_multiplicity_matrix_mgxs(mymgxs, temperature=temperature, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True + using_photon_production = False # multiplicity will fall back to using scatter and nu-scatter elif 'scatter matrix' in self.mgxs_types and \ @@ -1214,6 +1256,7 @@ def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True + using_photon_production = False # multiplicity will fall back to using scatter and nu-scatter elif 'consistent scatter matrix' in self.mgxs_types and \ @@ -1227,11 +1270,15 @@ def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True + using_photon_production = False else: using_multiplicity = False + using_photon_production = False - if using_multiplicity: + if using_photon_production: + pass + elif using_multiplicity: if 'nu-scatter matrix' in self.mgxs_types: nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix') else: @@ -1369,7 +1416,8 @@ def create_mg_library(self, xs_type='macro', xsdata_names=None, # Initialize file mgxs_file = openmc.MGXSLibrary( - self.energy_groups, num_delayed_groups=self.num_delayed_groups) + self.energy_groups, num_delayed_groups=self.num_delayed_groups, + particle_type=self.particle_type) if self.domain_type == 'mesh': # Create the xsdata objects and add to the mgxs_file @@ -1585,6 +1633,48 @@ def check_library_for_openmc_mgxs(self): error_flag = False + if 'photon-production matrix' in self.mgxs_types: + photon_mgxs_types = { + 'total', 'absorption', 'photon-production matrix'} + unsupported = set(self.mgxs_types) - photon_mgxs_types + if unsupported: + warn('Photon MGXS libraries do not support the following ' + f'MGXS types: {sorted(unsupported)}.') + error_flag = True + if self.particle_type != openmc.ParticleType.PHOTON: + warn('A photon production matrix requires particle_type to ' + 'be "photon".') + error_flag = True + if self.by_nuclide: + warn('Photon production cannot be tallied by nuclide.') + error_flag = True + if self.correction is not None: + warn('Photon MGXS libraries do not support a transport ' + 'correction; correction must be None.') + error_flag = True + if self.num_polar != 1 or self.num_azimuthal != 1: + warn('Photon production only supports an isotropic ' + 'representation.') + error_flag = True + if self.estimator not in (None, 'tracklength'): + warn('Photon MGXS libraries require the default tracklength ' + 'estimator for flux-weighted cross sections.') + error_flag = True + for mgxs_type in ('total', 'absorption'): + if mgxs_type not in self.mgxs_types: + warn(f'A "{mgxs_type}" MGXS type is required for a ' + 'photon MGXS library.') + error_flag = True + if error_flag: + raise ValueError('Invalid photon MGXS configuration ' + 'encountered.') + return + + if self.particle_type == openmc.ParticleType.PHOTON: + warn('A "photon-production matrix" MGXS type is required for a ' + 'photon MGXS library.') + raise ValueError('Invalid photon MGXS configuration encountered.') + # if correction is 'P0', then transport must be provided # otherwise total must be provided if self.correction == 'P0': diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 7f0200c0976..d415e2939a0 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -30,6 +30,7 @@ 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', + 'photon-production matrix', 'multiplicity matrix', 'nu-fission matrix', 'scatter probability matrix', @@ -177,6 +178,9 @@ class MGXS: The energy group structure for energy condensation by_nuclide : bool If true, computes cross sections for each nuclide in domain + particle_type : {'neutron', 'photon'}, optional + Particle type for which cross sections are computed. If not specified, + tallies are not filtered by particle type. name : str, optional Name of the multi-group cross section. Used as a label to identify tallies in OpenMC 'tallies.xml' file. @@ -195,6 +199,8 @@ class MGXS: Reaction type (e.g., 'total', 'nu-fission', etc.) by_nuclide : bool If true, computes cross sections for each nuclide in domain + particle_type : openmc.ParticleType or None + Particle type for which cross sections are computed domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh Domain for spatial homogenization domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'} @@ -261,10 +267,11 @@ class MGXS: def __init__(self, domain=None, domain_type=None, energy_groups=None, by_nuclide=False, name='', num_polar=1, - num_azimuthal=1): + num_azimuthal=1, particle_type=None): self._name = '' self._rxn_type = None self._by_nuclide = None + self._particle_type = None self._nuclides = None self._estimator = 'tracklength' self._domain = None @@ -284,6 +291,8 @@ def __init__(self, domain=None, domain_type=None, self.name = name self.by_nuclide = by_nuclide + if particle_type is not None: + self.particle_type = particle_type if domain_type is not None: self.domain_type = domain_type @@ -306,6 +315,7 @@ def __deepcopy__(self, memo): clone._name = self.name clone._rxn_type = self.rxn_type clone._by_nuclide = self.by_nuclide + clone._particle_type = self.particle_type clone._nuclides = copy.deepcopy(self._nuclides, memo) clone._domain = self.domain clone._domain_type = self.domain_type @@ -472,6 +482,18 @@ def by_nuclide(self, by_nuclide): cv.check_type('by_nuclide', by_nuclide, bool) self._by_nuclide = by_nuclide + @property + def particle_type(self): + return self._particle_type + + @particle_type.setter + def particle_type(self, particle_type): + particle_type = openmc.ParticleType(particle_type) + cv.check_value('particle type', particle_type, + (openmc.ParticleType.NEUTRON, + openmc.ParticleType.PHOTON)) + self._particle_type = particle_type + @property def domain(self): return self._domain @@ -611,6 +633,10 @@ def tallies(self): for add_filter in filters: self._tallies[key].filters.append(add_filter) + if self.particle_type is not None: + self._tallies[key].filters.append( + openmc.ParticleFilter(self.particle_type)) + # If this is a by-nuclide cross-section, add nuclides to Tally if self.by_nuclide and score != 'flux': self._tallies[key].nuclides += self.get_nuclides() @@ -719,7 +745,7 @@ def mgxs_type(self): @staticmethod def get_mgxs(mgxs_type, domain=None, domain_type=None, energy_groups=None, by_nuclide=False, name='', num_polar=1, - num_azimuthal=1): + num_azimuthal=1, particle_type=None): """Return a MGXS subclass object for some energy group structure within some spatial domain for some reaction type. @@ -753,6 +779,9 @@ def get_mgxs(mgxs_type, domain=None, domain_type=None, num_azimuthal : Integral, optional Number of equi-width azimuthal angles for angle discretization; defaults to no discretization + particle_type : {'neutron', 'photon'}, optional + Particle type for which cross sections are computed. If not + specified, tallies are not filtered by particle type. Returns ------- @@ -792,6 +821,10 @@ def get_mgxs(mgxs_type, domain=None, domain_type=None, mgxs = ScatterMatrixXS(domain, domain_type, energy_groups) elif mgxs_type == 'nu-scatter matrix': mgxs = ScatterMatrixXS(domain, domain_type, energy_groups, nu=True) + elif mgxs_type == 'photon-production matrix': + mgxs = PhotonProductionMatrixXS( + domain, domain_type, energy_groups, by_nuclide, name, + num_polar, num_azimuthal) elif mgxs_type == 'multiplicity matrix': mgxs = MultiplicityMatrixXS(domain, domain_type, energy_groups) elif mgxs_type == 'scatter probability matrix': @@ -835,6 +868,8 @@ def get_mgxs(mgxs_type, domain=None, domain_type=None, mgxs.name = name mgxs.num_polar = num_polar mgxs.num_azimuthal = num_azimuthal + if particle_type is not None: + mgxs.particle_type = particle_type return mgxs def get_nuclides(self): @@ -4923,6 +4958,104 @@ def print_groups_and_histogram(avg_xs, err_xs, num_groups, print(string) +@add_params +class PhotonProductionMatrixXS(MatrixMGXS): + r"""A photon production matrix multigroup cross section. + + This matrix includes both the photon that survives a coherent or incoherent + scattering event and secondary photons banked during photon interactions. + The latter include photons from atomic relaxation, thick-target + bremsstrahlung, and positron annihilation. Since each banked photon is + scored with its statistical weight, photon multiplicity is included directly + in the production matrix. + + Photon production is only available as a macroscopic, isotropic cross + section. Per-nuclide production cannot be determined because secondary + photons are tallied from the collision bank rather than a reaction score. + """ + + def __init__(self, domain=None, domain_type=None, energy_groups=None, + by_nuclide=False, name='', num_polar=1, num_azimuthal=1): + if by_nuclide: + raise ValueError('Photon production cannot be tallied by nuclide') + if num_polar != 1 or num_azimuthal != 1: + raise ValueError('Photon production only supports an isotropic ' + 'representation') + super().__init__(domain, domain_type, energy_groups, False, name) + self._rxn_type = 'photon-production' + self._mgxs_type = 'photon-production matrix' + self._particle_type = openmc.ParticleType.PHOTON + self._valid_estimators = ['analog'] + + @property + def by_nuclide(self): + return self._by_nuclide + + @by_nuclide.setter + def by_nuclide(self, by_nuclide): + cv.check_type('by_nuclide', by_nuclide, bool) + if by_nuclide: + raise ValueError('Photon production cannot be tallied by nuclide') + self._by_nuclide = False + + @property + def particle_type(self): + return self._particle_type + + @particle_type.setter + def particle_type(self, particle_type): + particle_type = openmc.ParticleType(particle_type) + if particle_type != openmc.ParticleType.PHOTON: + raise ValueError('Photon production requires photon tallies') + self._particle_type = particle_type + + @property + def scores(self): + return ['flux', 'scatter', 'events'] + + @property + def tally_keys(self): + return ['flux', 'primary photon production', 'secondary photon production'] + + @property + def filters(self): + group_edges = self.energy_groups.group_edges + energy = openmc.EnergyFilter(group_edges) + energyout = openmc.EnergyoutFilter(group_edges) + production = openmc.ParticleProductionFilter( + 'photon', group_edges) + return [[energy], [energy, energyout], [energy, production]] + + @property + def estimator(self): + return ['tracklength', 'analog', 'analog'] + + @property + def rxn_rate_tally(self): + if self._rxn_rate_tally is None: + primary = copy.deepcopy( + self.tallies['primary photon production']) + secondary = copy.deepcopy( + self.tallies['secondary photon production']) + + # Give the secondary tally the same outgoing-energy filter as the + # primary tally so that the two production contributions can be + # combined with standard tally arithmetic. + energyout = copy.deepcopy( + primary.find_filter(openmc.EnergyoutFilter)) + secondary.filters = [ + energyout + if isinstance(f, openmc.ParticleProductionFilter) else f + for f in secondary.filters + ] + primary._scores = ['photon-production'] + secondary._scores = ['photon-production'] + self._rxn_rate_tally = primary + secondary + self._rxn_rate_tally.sparse = self.sparse + + return self._rxn_rate_tally + + @add_params class MultiplicityMatrixXS(MatrixMGXS): r"""The scattering multiplicity matrix. diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index ab9b58b7a3a..776aa8b5bec 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1634,6 +1634,49 @@ def set_scatter_matrix_mgxs(self, scatter, temperature=ROOM_TEMPERATURE_KELVIN, scatter.get_xs(nuclides=nuclide, xs_type=xs_type, subdomains=subdomain) + def set_photon_production_mgxs( + self, production, temperature=ROOM_TEMPERATURE_KELVIN, + xs_type='macro', subdomain=None): + """Set the scattering data from a photon production matrix MGXS. + + The production matrix is written as the zeroth-order scattering + matrix. Its multiplicity is already included in the matrix values, so + no separate multiplicity matrix is written. + + Parameters + ---------- + production : openmc.mgxs.PhotonProductionMatrixXS + Photon production matrix to use + temperature : float, optional + Temperature of the cross section data in kelvin + xs_type : {'macro'}, optional + Cross section type. Photon production is only available as a + macroscopic cross section. + subdomain : iterable of int, optional + Mesh subdomain for which data are requested + """ + + check_type('production', production, + openmc.mgxs.PhotonProductionMatrixXS) + check_value('energy_groups', production.energy_groups, + [self.energy_groups]) + check_value('domain_type', production.domain_type, + openmc.mgxs.DOMAIN_TYPES) + check_value('xs_type', xs_type, ['macro']) + self._check_temperature(temperature) + + if self.representation != REPRESENTATION_ISOTROPIC: + raise ValueError('Photon production only supports an isotropic ' + 'representation') + + self.scatter_format = SCATTER_LEGENDRE + self.order = 0 + i = self._temperature_index(temperature) + self._scatter_matrix[i] = np.zeros( + self.xs_shapes["[G][G'][Order]"]) + self._scatter_matrix[i][:, :, 0] = production.get_xs( + xs_type=xs_type, subdomains=subdomain) + def set_multiplicity_matrix_mgxs(self, nuscatter, scatter=None, temperature=ROOM_TEMPERATURE_KELVIN, nuclide='total', xs_type='macro', subdomain=None): @@ -2364,6 +2407,8 @@ class MGXSLibrary: Energy group structure num_delayed_groups : int Num delayed groups + particle_type : {'neutron', 'photon'}, optional + Particle type represented by the library Attributes ---------- @@ -2371,13 +2416,19 @@ class MGXSLibrary: Energy group structure. num_delayed_groups : int Num delayed groups + particle_type : openmc.ParticleType or None + Particle type represented by the library xsdatas : Iterable of openmc.XSdata Iterable of multi-Group cross section data objects """ - def __init__(self, energy_groups, num_delayed_groups=0): + def __init__(self, energy_groups, num_delayed_groups=0, + particle_type=None): self.energy_groups = energy_groups self.num_delayed_groups = num_delayed_groups + self._particle_type = None + if particle_type is not None: + self.particle_type = particle_type self._xsdatas = [] def __deepcopy__(self, memo): @@ -2388,6 +2439,7 @@ def __deepcopy__(self, memo): clone = type(self).__new__(type(self)) clone._energy_groups = copy.deepcopy(self.energy_groups, memo) clone._num_delayed_groups = self.num_delayed_groups + clone._particle_type = self.particle_type clone._xsdatas = copy.deepcopy(self.xsdatas, memo) memo[id(self)] = clone @@ -2420,6 +2472,18 @@ def num_delayed_groups(self, num_delayed_groups): openmc.mgxs.MAX_DELAYED_GROUPS, equality=True) self._num_delayed_groups = num_delayed_groups + @property + def particle_type(self): + return self._particle_type + + @particle_type.setter + def particle_type(self, particle_type): + particle_type = openmc.ParticleType(particle_type) + check_value('particle type', particle_type, + (openmc.ParticleType.NEUTRON, + openmc.ParticleType.PHOTON)) + self._particle_type = particle_type + @property def xsdatas(self): return self._xsdatas @@ -2592,6 +2656,8 @@ def export_to_hdf5(self, filename='mgxs.h5', libver='earliest'): file.attrs['energy_groups'] = self.energy_groups.num_groups file.attrs['delayed_groups'] = self.num_delayed_groups file.attrs['group structure'] = self.energy_groups.group_edges + if self.particle_type is not None: + file.attrs['particle_type'] = np.bytes_(str(self.particle_type)) for xsdata in self._xsdatas: xsdata.to_hdf5(file) @@ -2633,7 +2699,10 @@ def from_hdf5(cls, filename=None): group_structure = file.attrs['group structure'] num_delayed_groups = file.attrs['delayed_groups'] energy_groups = openmc.mgxs.EnergyGroups(group_structure) - data = cls(energy_groups, num_delayed_groups) + particle_type = file.attrs.get('particle_type') + if isinstance(particle_type, bytes): + particle_type = particle_type.decode() + data = cls(energy_groups, num_delayed_groups, particle_type) for group_name, group in file.items(): data.add_xsdata(openmc.XSdata.from_hdf5(group, group_name, diff --git a/tests/regression_tests/mgxs_library_condense/test.py b/tests/regression_tests/mgxs_library_condense/test.py index bbc4c11bfa9..c9eae94e6c5 100644 --- a/tests/regression_tests/mgxs_library_condense/test.py +++ b/tests/regression_tests/mgxs_library_condense/test.py @@ -19,7 +19,10 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = False # Test all MGXS types - self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ + neutron_mgxs_types = tuple( + xs for xs in openmc.mgxs.MGXS_TYPES + if xs != 'photon-production matrix') + self.mgxs_lib.mgxs_types = neutron_mgxs_types + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups self.mgxs_lib.num_delayed_groups = 6 diff --git a/tests/regression_tests/mgxs_library_distribcell/test.py b/tests/regression_tests/mgxs_library_distribcell/test.py index 464b309c002..c0df9fce8d9 100644 --- a/tests/regression_tests/mgxs_library_distribcell/test.py +++ b/tests/regression_tests/mgxs_library_distribcell/test.py @@ -23,8 +23,9 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = False # Test all relevant MGXS types - relevant_MGXS_TYPES = [item for item in openmc.mgxs.MGXS_TYPES - if item != 'current'] + relevant_MGXS_TYPES = [ + item for item in openmc.mgxs.MGXS_TYPES + if item not in ('current', 'photon-production matrix')] self.mgxs_lib.mgxs_types = tuple(relevant_MGXS_TYPES) + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups diff --git a/tests/regression_tests/mgxs_library_hdf5/test.py b/tests/regression_tests/mgxs_library_hdf5/test.py index 4fb4bf09369..709512cce2c 100644 --- a/tests/regression_tests/mgxs_library_hdf5/test.py +++ b/tests/regression_tests/mgxs_library_hdf5/test.py @@ -23,7 +23,10 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = False # Test all MGXS types - self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ + neutron_mgxs_types = tuple( + xs for xs in openmc.mgxs.MGXS_TYPES + if xs != 'photon-production matrix') + self.mgxs_lib.mgxs_types = neutron_mgxs_types + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups self.mgxs_lib.num_delayed_groups = 6 diff --git a/tests/regression_tests/mgxs_library_mesh/test.py b/tests/regression_tests/mgxs_library_mesh/test.py index c1a5980b5d1..3d220ee8c54 100644 --- a/tests/regression_tests/mgxs_library_mesh/test.py +++ b/tests/regression_tests/mgxs_library_mesh/test.py @@ -40,7 +40,10 @@ def model(): model.mgxs_lib.by_nuclide = False # Test all MGXS types - model.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + openmc.mgxs.MDGXS_TYPES + neutron_mgxs_types = tuple( + xs for xs in openmc.mgxs.MGXS_TYPES + if xs != 'photon-production matrix') + model.mgxs_lib.mgxs_types = neutron_mgxs_types + openmc.mgxs.MDGXS_TYPES model.mgxs_lib.energy_groups = energy_groups model.mgxs_lib.num_delayed_groups = 6 model.mgxs_lib.correction = None # Avoid warning about P0 correction diff --git a/tests/regression_tests/mgxs_library_no_nuclides/test.py b/tests/regression_tests/mgxs_library_no_nuclides/test.py index a02086af3ec..f7588bafc78 100644 --- a/tests/regression_tests/mgxs_library_no_nuclides/test.py +++ b/tests/regression_tests/mgxs_library_no_nuclides/test.py @@ -20,8 +20,9 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = False # Test relevant MGXS types - relevant_MGXS_TYPES = [item for item in openmc.mgxs.MGXS_TYPES - if item != 'current'] + relevant_MGXS_TYPES = [ + item for item in openmc.mgxs.MGXS_TYPES + if item not in ('current', 'photon-production matrix')] # Add in a subset of openmc.mgxs.ARBITRARY_VECTOR_TYPES and # openmc.mgxs.ARBITRARY_MATRIX_TYPES so we can see the code works, # but not use too much resources diff --git a/tests/regression_tests/mgxs_library_nuclides/test.py b/tests/regression_tests/mgxs_library_nuclides/test.py index a10070358ad..e11bb995e27 100644 --- a/tests/regression_tests/mgxs_library_nuclides/test.py +++ b/tests/regression_tests/mgxs_library_nuclides/test.py @@ -19,8 +19,9 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = True # Test relevant MGXS types - relevant_MGXS_TYPES = [item for item in openmc.mgxs.MGXS_TYPES - if item != 'current'] + relevant_MGXS_TYPES = [ + item for item in openmc.mgxs.MGXS_TYPES + if item not in ('current', 'photon-production matrix')] # Add in a subset of openmc.mgxs.ARBITRARY_VECTOR_TYPES and # openmc.mgxs.ARBITRARY_MATRIX_TYPES so we can see the code works, # but not use too much resources diff --git a/tests/regression_tests/mgxs_library_specific_nuclides/test.py b/tests/regression_tests/mgxs_library_specific_nuclides/test.py index 0ccbb83bdbe..688c0660c84 100644 --- a/tests/regression_tests/mgxs_library_specific_nuclides/test.py +++ b/tests/regression_tests/mgxs_library_specific_nuclides/test.py @@ -19,8 +19,9 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = True # Test relevant MGXS types - relevant_MGXS_TYPES = [item for item in openmc.mgxs.MGXS_TYPES - if item != 'current'] + relevant_MGXS_TYPES = [ + item for item in openmc.mgxs.MGXS_TYPES + if item not in ('current', 'photon-production matrix')] # Add in a subset of openmc.mgxs.ARBITRARY_VECTOR_TYPES and # openmc.mgxs.ARBITRARY_MATRIX_TYPES so we can see the code works, # but not use too much resources diff --git a/tests/unit_tests/test_photon_mgxs.py b/tests/unit_tests/test_photon_mgxs.py new file mode 100644 index 00000000000..08a8eccc07d --- /dev/null +++ b/tests/unit_tests/test_photon_mgxs.py @@ -0,0 +1,115 @@ +import h5py +import numpy as np +import pytest + +import openmc + + +def test_photon_production_matrix_tallies(): + material = openmc.Material() + groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) + production = openmc.mgxs.PhotonProductionMatrixXS( + material, 'material', groups) + + assert production.particle_type == openmc.ParticleType.PHOTON + assert production.estimator == ['tracklength', 'analog', 'analog'] + assert production.tally_keys == [ + 'flux', 'primary photon production', 'secondary photon production'] + + tallies = production.tallies + assert tallies['primary photon production'].scores == ['scatter'] + assert tallies['secondary photon production'].scores == ['events'] + assert tallies['secondary photon production'].contains_filter( + openmc.ParticleProductionFilter) + for tally in tallies.values(): + particle_filter = tally.find_filter(openmc.ParticleFilter) + assert particle_filter.bins == ['photon'] + + +def test_photon_production_matrix_constraints(): + with pytest.raises(ValueError, match='by nuclide'): + openmc.mgxs.MGXS.get_mgxs( + 'photon-production matrix', by_nuclide=True) + with pytest.raises(ValueError, match='isotropic'): + openmc.mgxs.MGXS.get_mgxs( + 'photon-production matrix', num_polar=2) + with pytest.raises(ValueError, match='photon tallies'): + openmc.mgxs.MGXS.get_mgxs( + 'photon-production matrix', particle_type='neutron') + + +def test_photon_production_matrix_combines_primary_and_secondary(): + material = openmc.Material() + groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) + production = openmc.mgxs.PhotonProductionMatrixXS( + material, 'material', groups) + + for tally in production.tallies.values(): + values = np.arange(1, tally.num_filter_bins + 1, dtype=float) + values.shape = (tally.num_filter_bins, 1, 1) + tally._sum = 2.0 * values + tally._sum_sq = 2.0 * values**2 + tally._num_realizations = 2 + tally._sp_filename = 'statepoint.h5' + tally._results_read = True + + np.testing.assert_allclose( + production.rxn_rate_tally.mean.ravel(), [2.0, 4.0, 6.0, 8.0]) + + +def test_set_photon_production_mgxs(monkeypatch): + material = openmc.Material() + groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) + production = openmc.mgxs.PhotonProductionMatrixXS( + material, 'material', groups) + values = np.array([[0.1, 0.2], [0.3, 0.4]]) + monkeypatch.setattr(production, 'get_xs', lambda **kwargs: values) + + xsdata = openmc.XSdata('photon', groups) + xsdata.order = 0 + absorption = np.array([0.5, 0.6]) + xsdata.set_absorption(absorption) + xsdata.set_photon_production_mgxs(production) + + assert xsdata.scatter_format == 'legendre' + assert xsdata.order == 0 + np.testing.assert_allclose(xsdata._scatter_matrix[0][:, :, 0], values) + np.testing.assert_allclose(xsdata._absorption[0], absorption) + assert xsdata._multiplicity_matrix[0] is None + + +def test_photon_mgxs_library(tmp_path): + material = openmc.Material() + geometry = openmc.Geometry([openmc.Cell(fill=material)]) + groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) + library = openmc.mgxs.Library( + geometry, mgxs_types=[ + 'total', 'absorption', 'photon-production matrix'], + particle_type='photon') + library.domain_type = 'material' + library.energy_groups = groups + library.correction = None + library.build_library() + library.check_library_for_openmc_mgxs() + + for mgxs in library.all_mgxs[material.id].values(): + assert mgxs.particle_type == openmc.ParticleType.PHOTON + + mg_library = openmc.MGXSLibrary(groups, particle_type='photon') + path = tmp_path / 'mgxs.h5' + mg_library.export_to_hdf5(path) + with h5py.File(path) as h5file: + assert h5file.attrs['particle_type'] == b'photon' + assert openmc.MGXSLibrary.from_hdf5(path).particle_type == \ + openmc.ParticleType.PHOTON + + +def test_all_mgxs_types_respects_particle_type(): + geometry = openmc.Geometry([openmc.Cell()]) + neutron_library = openmc.mgxs.Library(geometry, mgxs_types='all') + photon_library = openmc.mgxs.Library( + geometry, mgxs_types='all', particle_type='photon') + + assert 'photon-production matrix' not in neutron_library.mgxs_types + assert photon_library.mgxs_types == ( + 'total', 'absorption', 'photon-production matrix') From 39e4fb26061f592e1830a7a93f3b5286e15d1f82 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 8 Aug 2026 21:35:18 -0500 Subject: [PATCH 02/10] Add photon MGXS regression test --- .../regression_tests/mgxs_photon/__init__.py | 0 .../mgxs_photon/inputs_true.dat | 95 ++++++++++++++++ .../mgxs_photon/results_true.dat | 8 ++ tests/regression_tests/mgxs_photon/test.py | 102 ++++++++++++++++++ 4 files changed, 205 insertions(+) create mode 100644 tests/regression_tests/mgxs_photon/__init__.py create mode 100644 tests/regression_tests/mgxs_photon/inputs_true.dat create mode 100644 tests/regression_tests/mgxs_photon/results_true.dat create mode 100644 tests/regression_tests/mgxs_photon/test.py diff --git a/tests/regression_tests/mgxs_photon/__init__.py b/tests/regression_tests/mgxs_photon/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/mgxs_photon/inputs_true.dat b/tests/regression_tests/mgxs_photon/inputs_true.dat new file mode 100644 index 00000000000..aac04060864 --- /dev/null +++ b/tests/regression_tests/mgxs_photon/inputs_true.dat @@ -0,0 +1,95 @@ + + + + + + + + + + + + + + + + + fixed source + 2000 + 2 + + + 0.0 0.0 0.0 + + + 1000000.0 1.0 + + + ttb + true + true + + 1000.0 + + + + + 1 + + + 1000.0 100000.0 500000.0 1100000.0 + + + photon + + + 1000.0 100000.0 500000.0 1100000.0 + + + photon + 1000.0 100000.0 500000.0 1100000.0 + + + 1 2 3 + total + flux + tracklength + + + 1 2 3 + total + total + tracklength + + + 1 2 3 + total + flux + tracklength + + + 1 2 3 + total + absorption + tracklength + + + 1 2 3 + total + flux + tracklength + + + 1 2 11 3 + total + scatter + analog + + + 1 2 12 3 + total + events + analog + + + diff --git a/tests/regression_tests/mgxs_photon/results_true.dat b/tests/regression_tests/mgxs_photon/results_true.dat new file mode 100644 index 00000000000..0a84849600c --- /dev/null +++ b/tests/regression_tests/mgxs_photon/results_true.dat @@ -0,0 +1,8 @@ +primary photons: 6.55500000e-01 +secondary photons: 1.25150000e+00 +production matrix: +[[ 0.33790683 0.34765112 0.79557413] + [ 0. 1.17578035 5.49497349] + [ 0. 0. 19.5842782 ]] +absorption: +[ 0.23475329 3.62317561 76.43950604] diff --git a/tests/regression_tests/mgxs_photon/test.py b/tests/regression_tests/mgxs_photon/test.py new file mode 100644 index 00000000000..b7ea928c7c7 --- /dev/null +++ b/tests/regression_tests/mgxs_photon/test.py @@ -0,0 +1,102 @@ +import hashlib + +import numpy as np +import openmc + +from tests.testing_harness import PyAPITestHarness + + +class PhotonMGXSTestHarness(PyAPITestHarness): + """Run photon transport and verify photon MGXS post-processing.""" + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + self.material = openmc.Material() + self.material.set_density('g/cm3', 11.35) + self.material.add_element('Pb', 1.0) + self._model.materials = openmc.Materials([self.material]) + + sphere = openmc.Sphere(r=1.0, boundary_type='vacuum') + cell = openmc.Cell(fill=self.material, region=-sphere) + self._model.geometry = openmc.Geometry([cell]) + + self._model.settings.run_mode = 'fixed source' + self._model.settings.particles = 2000 + self._model.settings.batches = 2 + self._model.settings.photon_transport = True + self._model.settings.atomic_relaxation = True + self._model.settings.electron_treatment = 'ttb' + self._model.settings.cutoff = {'energy_photon': 1000.0} + self._model.settings.source = openmc.IndependentSource( + particle='photon', + space=openmc.stats.Point((0.0, 0.0, 0.0)), + energy=openmc.stats.Discrete([1.0e6], [1.0])) + + groups = openmc.mgxs.EnergyGroups( + group_edges=[1.0e3, 1.0e5, 5.0e5, 1.1e6]) + self.mgxs_lib = openmc.mgxs.Library( + self._model.geometry, + mgxs_types=['total', 'absorption', 'photon-production matrix'], + particle_type='photon') + self.mgxs_lib.energy_groups = groups + self.mgxs_lib.correction = None + self.mgxs_lib.domain_type = 'material' + self.mgxs_lib.build_library() + self.mgxs_lib.add_to_tallies(self._model.tallies, merge=False) + + def _get_results(self, hash_output=False): + with openmc.StatePoint(self._sp_name) as statepoint: + self.mgxs_lib.load_from_statepoint(statepoint) + + production = self.mgxs_lib.get_mgxs( + self.material, 'photon-production matrix') + absorption = self.mgxs_lib.get_mgxs(self.material, 'absorption') + + primary = production.tallies[ + 'primary photon production'].mean.sum() + secondary = production.tallies[ + 'secondary photon production'].mean.sum() + if secondary <= 0.0: + raise AssertionError( + 'Photon transport did not score any secondary photons') + + production_xs = production.get_xs() + absorption_xs = absorption.get_xs() + + # Exercise the XSdata conversion used by subsequent MG/RR + # workflows. + mg_library = self.mgxs_lib.create_mg_library() + xsdata = mg_library.xsdatas[0] + if xsdata.multiplicity_matrix[0] is not None: + raise AssertionError( + 'Photon production should not create a multiplicity ' + 'matrix') + if not np.allclose( + xsdata.scatter_matrix[0][:, :, 0], production_xs): + raise AssertionError('Photon production matrix was not ' + 'exported') + if not np.allclose(xsdata.absorption[0], absorption_xs): + raise AssertionError('Photon absorption was changed during ' + 'export') + + output = [ + f'primary photons: {primary:.8e}', + f'secondary photons: {secondary:.8e}', + 'production matrix:', + np.array2string(production_xs, precision=8), + 'absorption:', + np.array2string(absorption_xs, precision=8), + ] + output = '\n'.join(output) + '\n' + + if hash_output: + digest = hashlib.sha512(output.encode('utf-8')) + output = digest.hexdigest() + return output + + +def test_photon_mgxs(): + harness = PhotonMGXSTestHarness( + 'statepoint.2.h5', model=openmc.Model()) + harness.main() From 8d695ce382d990df631112a3a1f6104beb72f8ff Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 9 Aug 2026 16:47:34 -0500 Subject: [PATCH 03/10] Use nu-scatter formalism for photon production --- openmc/mgxs/library.py | 43 ++++---------- openmc/mgxs/mgxs.py | 22 +++++-- openmc/mgxs_library.py | 46 +-------------- .../mgxs_library_condense/test.py | 5 +- .../mgxs_library_distribcell/test.py | 5 +- .../mgxs_library_hdf5/test.py | 5 +- .../mgxs_library_mesh/test.py | 5 +- .../mgxs_library_no_nuclides/test.py | 5 +- .../mgxs_library_nuclides/test.py | 5 +- .../mgxs_library_specific_nuclides/test.py | 5 +- tests/regression_tests/mgxs_photon/test.py | 59 +++++++++++-------- tests/unit_tests/test_photon_mgxs.py | 20 +++---- 12 files changed, 86 insertions(+), 139 deletions(-) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 19da8d7cad5..6338c7d732d 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -214,11 +214,9 @@ def mgxs_types(self, mgxs_types): if mgxs_types == 'all': if self.particle_type == openmc.ParticleType.PHOTON: self._mgxs_types = ( - 'total', 'absorption', 'photon-production matrix') + 'total', 'absorption', 'nu-scatter matrix') else: - self._mgxs_types = tuple( - xs for xs in all_mgxs_types - if xs != 'photon-production matrix') + self._mgxs_types = all_mgxs_types else: cv.check_iterable_type('mgxs_types', mgxs_types, str) for mgxs_type in mgxs_types: @@ -1223,16 +1221,14 @@ def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', xsdata.set_decay_rate_mgxs(mymgxs, temperature=temperature, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) - # Photon production includes both the surviving primary photon and all - # banked secondary photons. Its multiplicity is therefore already - # folded into the matrix values. - if 'photon-production matrix' in self.mgxs_types: - production = self.get_mgxs( - domain, 'photon-production matrix') - xsdata.set_photon_production_mgxs( + # Photon nu-scatter includes both the surviving primary photon and all + # banked secondary photons, so its multiplicity is already folded into + # the matrix values. + if self.particle_type == openmc.ParticleType.PHOTON: + production = self.get_mgxs(domain, 'nu-scatter matrix') + xsdata.set_scatter_matrix_mgxs( production, temperature=temperature, xs_type=xs_type, - subdomain=subdomain) - using_photon_production = True + nuclide=[nuclide], subdomain=subdomain) using_multiplicity = False # If multiplicity matrix is available, prefer that @@ -1243,7 +1239,6 @@ def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True - using_photon_production = False # multiplicity will fall back to using scatter and nu-scatter elif 'scatter matrix' in self.mgxs_types and \ @@ -1256,7 +1251,6 @@ def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True - using_photon_production = False # multiplicity will fall back to using scatter and nu-scatter elif 'consistent scatter matrix' in self.mgxs_types and \ @@ -1270,13 +1264,11 @@ def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', nuclide=[nuclide], subdomain=subdomain) using_multiplicity = True - using_photon_production = False else: using_multiplicity = False - using_photon_production = False - if using_photon_production: + if self.particle_type == openmc.ParticleType.PHOTON: pass elif using_multiplicity: if 'nu-scatter matrix' in self.mgxs_types: @@ -1633,18 +1625,14 @@ def check_library_for_openmc_mgxs(self): error_flag = False - if 'photon-production matrix' in self.mgxs_types: + if self.particle_type == openmc.ParticleType.PHOTON: photon_mgxs_types = { - 'total', 'absorption', 'photon-production matrix'} + 'total', 'absorption', 'nu-scatter matrix'} unsupported = set(self.mgxs_types) - photon_mgxs_types if unsupported: warn('Photon MGXS libraries do not support the following ' f'MGXS types: {sorted(unsupported)}.') error_flag = True - if self.particle_type != openmc.ParticleType.PHOTON: - warn('A photon production matrix requires particle_type to ' - 'be "photon".') - error_flag = True if self.by_nuclide: warn('Photon production cannot be tallied by nuclide.') error_flag = True @@ -1660,7 +1648,7 @@ def check_library_for_openmc_mgxs(self): warn('Photon MGXS libraries require the default tracklength ' 'estimator for flux-weighted cross sections.') error_flag = True - for mgxs_type in ('total', 'absorption'): + for mgxs_type in ('total', 'absorption', 'nu-scatter matrix'): if mgxs_type not in self.mgxs_types: warn(f'A "{mgxs_type}" MGXS type is required for a ' 'photon MGXS library.') @@ -1670,11 +1658,6 @@ def check_library_for_openmc_mgxs(self): 'encountered.') return - if self.particle_type == openmc.ParticleType.PHOTON: - warn('A "photon-production matrix" MGXS type is required for a ' - 'photon MGXS library.') - raise ValueError('Invalid photon MGXS configuration encountered.') - # if correction is 'P0', then transport must be provided # otherwise total must be provided if self.correction == 'P0': diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index d415e2939a0..b4eeae63045 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -30,7 +30,6 @@ 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', - 'photon-production matrix', 'multiplicity matrix', 'nu-fission matrix', 'scatter probability matrix', @@ -819,12 +818,15 @@ def get_mgxs(mgxs_type, domain=None, domain_type=None, mgxs = ScatterXS(domain, domain_type, energy_groups, nu=True) elif mgxs_type == 'scatter matrix': mgxs = ScatterMatrixXS(domain, domain_type, energy_groups) - elif mgxs_type == 'nu-scatter matrix': - mgxs = ScatterMatrixXS(domain, domain_type, energy_groups, nu=True) - elif mgxs_type == 'photon-production matrix': + elif mgxs_type == 'nu-scatter matrix' and \ + particle_type is not None and \ + openmc.ParticleType(particle_type) == \ + openmc.ParticleType.PHOTON: mgxs = PhotonProductionMatrixXS( domain, domain_type, energy_groups, by_nuclide, name, num_polar, num_azimuthal) + elif mgxs_type == 'nu-scatter matrix': + mgxs = ScatterMatrixXS(domain, domain_type, energy_groups, nu=True) elif mgxs_type == 'multiplicity matrix': mgxs = MultiplicityMatrixXS(domain, domain_type, energy_groups) elif mgxs_type == 'scatter probability matrix': @@ -4982,8 +4984,8 @@ def __init__(self, domain=None, domain_type=None, energy_groups=None, raise ValueError('Photon production only supports an isotropic ' 'representation') super().__init__(domain, domain_type, energy_groups, False, name) - self._rxn_type = 'photon-production' - self._mgxs_type = 'photon-production matrix' + self._rxn_type = 'nu-scatter' + self._mgxs_type = 'nu-scatter matrix' self._particle_type = openmc.ParticleType.PHOTON self._valid_estimators = ['analog'] @@ -5030,6 +5032,14 @@ def filters(self): def estimator(self): return ['tracklength', 'analog', 'analog'] + @property + def scatter_format(self): + return SCATTER_LEGENDRE + + @property + def legendre_order(self): + return 0 + @property def rxn_rate_tally(self): if self._rxn_rate_tally is None: diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 776aa8b5bec..7d3c3d4151b 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1586,7 +1586,8 @@ def set_scatter_matrix_mgxs(self, scatter, temperature=ROOM_TEMPERATURE_KELVIN, """ - check_type('scatter', scatter, openmc.mgxs.ScatterMatrixXS) + check_type('scatter', scatter, (openmc.mgxs.ScatterMatrixXS, + openmc.mgxs.PhotonProductionMatrixXS)) check_value('energy_groups', scatter.energy_groups, [self.energy_groups]) check_value('domain_type', scatter.domain_type, @@ -1634,49 +1635,6 @@ def set_scatter_matrix_mgxs(self, scatter, temperature=ROOM_TEMPERATURE_KELVIN, scatter.get_xs(nuclides=nuclide, xs_type=xs_type, subdomains=subdomain) - def set_photon_production_mgxs( - self, production, temperature=ROOM_TEMPERATURE_KELVIN, - xs_type='macro', subdomain=None): - """Set the scattering data from a photon production matrix MGXS. - - The production matrix is written as the zeroth-order scattering - matrix. Its multiplicity is already included in the matrix values, so - no separate multiplicity matrix is written. - - Parameters - ---------- - production : openmc.mgxs.PhotonProductionMatrixXS - Photon production matrix to use - temperature : float, optional - Temperature of the cross section data in kelvin - xs_type : {'macro'}, optional - Cross section type. Photon production is only available as a - macroscopic cross section. - subdomain : iterable of int, optional - Mesh subdomain for which data are requested - """ - - check_type('production', production, - openmc.mgxs.PhotonProductionMatrixXS) - check_value('energy_groups', production.energy_groups, - [self.energy_groups]) - check_value('domain_type', production.domain_type, - openmc.mgxs.DOMAIN_TYPES) - check_value('xs_type', xs_type, ['macro']) - self._check_temperature(temperature) - - if self.representation != REPRESENTATION_ISOTROPIC: - raise ValueError('Photon production only supports an isotropic ' - 'representation') - - self.scatter_format = SCATTER_LEGENDRE - self.order = 0 - i = self._temperature_index(temperature) - self._scatter_matrix[i] = np.zeros( - self.xs_shapes["[G][G'][Order]"]) - self._scatter_matrix[i][:, :, 0] = production.get_xs( - xs_type=xs_type, subdomains=subdomain) - def set_multiplicity_matrix_mgxs(self, nuscatter, scatter=None, temperature=ROOM_TEMPERATURE_KELVIN, nuclide='total', xs_type='macro', subdomain=None): diff --git a/tests/regression_tests/mgxs_library_condense/test.py b/tests/regression_tests/mgxs_library_condense/test.py index c9eae94e6c5..bbc4c11bfa9 100644 --- a/tests/regression_tests/mgxs_library_condense/test.py +++ b/tests/regression_tests/mgxs_library_condense/test.py @@ -19,10 +19,7 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = False # Test all MGXS types - neutron_mgxs_types = tuple( - xs for xs in openmc.mgxs.MGXS_TYPES - if xs != 'photon-production matrix') - self.mgxs_lib.mgxs_types = neutron_mgxs_types + \ + self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups self.mgxs_lib.num_delayed_groups = 6 diff --git a/tests/regression_tests/mgxs_library_distribcell/test.py b/tests/regression_tests/mgxs_library_distribcell/test.py index c0df9fce8d9..464b309c002 100644 --- a/tests/regression_tests/mgxs_library_distribcell/test.py +++ b/tests/regression_tests/mgxs_library_distribcell/test.py @@ -23,9 +23,8 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = False # Test all relevant MGXS types - relevant_MGXS_TYPES = [ - item for item in openmc.mgxs.MGXS_TYPES - if item not in ('current', 'photon-production matrix')] + relevant_MGXS_TYPES = [item for item in openmc.mgxs.MGXS_TYPES + if item != 'current'] self.mgxs_lib.mgxs_types = tuple(relevant_MGXS_TYPES) + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups diff --git a/tests/regression_tests/mgxs_library_hdf5/test.py b/tests/regression_tests/mgxs_library_hdf5/test.py index 709512cce2c..4fb4bf09369 100644 --- a/tests/regression_tests/mgxs_library_hdf5/test.py +++ b/tests/regression_tests/mgxs_library_hdf5/test.py @@ -23,10 +23,7 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = False # Test all MGXS types - neutron_mgxs_types = tuple( - xs for xs in openmc.mgxs.MGXS_TYPES - if xs != 'photon-production matrix') - self.mgxs_lib.mgxs_types = neutron_mgxs_types + \ + self.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + \ openmc.mgxs.MDGXS_TYPES self.mgxs_lib.energy_groups = energy_groups self.mgxs_lib.num_delayed_groups = 6 diff --git a/tests/regression_tests/mgxs_library_mesh/test.py b/tests/regression_tests/mgxs_library_mesh/test.py index 3d220ee8c54..c1a5980b5d1 100644 --- a/tests/regression_tests/mgxs_library_mesh/test.py +++ b/tests/regression_tests/mgxs_library_mesh/test.py @@ -40,10 +40,7 @@ def model(): model.mgxs_lib.by_nuclide = False # Test all MGXS types - neutron_mgxs_types = tuple( - xs for xs in openmc.mgxs.MGXS_TYPES - if xs != 'photon-production matrix') - model.mgxs_lib.mgxs_types = neutron_mgxs_types + openmc.mgxs.MDGXS_TYPES + model.mgxs_lib.mgxs_types = openmc.mgxs.MGXS_TYPES + openmc.mgxs.MDGXS_TYPES model.mgxs_lib.energy_groups = energy_groups model.mgxs_lib.num_delayed_groups = 6 model.mgxs_lib.correction = None # Avoid warning about P0 correction diff --git a/tests/regression_tests/mgxs_library_no_nuclides/test.py b/tests/regression_tests/mgxs_library_no_nuclides/test.py index f7588bafc78..a02086af3ec 100644 --- a/tests/regression_tests/mgxs_library_no_nuclides/test.py +++ b/tests/regression_tests/mgxs_library_no_nuclides/test.py @@ -20,9 +20,8 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = False # Test relevant MGXS types - relevant_MGXS_TYPES = [ - item for item in openmc.mgxs.MGXS_TYPES - if item not in ('current', 'photon-production matrix')] + relevant_MGXS_TYPES = [item for item in openmc.mgxs.MGXS_TYPES + if item != 'current'] # Add in a subset of openmc.mgxs.ARBITRARY_VECTOR_TYPES and # openmc.mgxs.ARBITRARY_MATRIX_TYPES so we can see the code works, # but not use too much resources diff --git a/tests/regression_tests/mgxs_library_nuclides/test.py b/tests/regression_tests/mgxs_library_nuclides/test.py index e11bb995e27..a10070358ad 100644 --- a/tests/regression_tests/mgxs_library_nuclides/test.py +++ b/tests/regression_tests/mgxs_library_nuclides/test.py @@ -19,9 +19,8 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = True # Test relevant MGXS types - relevant_MGXS_TYPES = [ - item for item in openmc.mgxs.MGXS_TYPES - if item not in ('current', 'photon-production matrix')] + relevant_MGXS_TYPES = [item for item in openmc.mgxs.MGXS_TYPES + if item != 'current'] # Add in a subset of openmc.mgxs.ARBITRARY_VECTOR_TYPES and # openmc.mgxs.ARBITRARY_MATRIX_TYPES so we can see the code works, # but not use too much resources diff --git a/tests/regression_tests/mgxs_library_specific_nuclides/test.py b/tests/regression_tests/mgxs_library_specific_nuclides/test.py index 688c0660c84..0ccbb83bdbe 100644 --- a/tests/regression_tests/mgxs_library_specific_nuclides/test.py +++ b/tests/regression_tests/mgxs_library_specific_nuclides/test.py @@ -19,9 +19,8 @@ def __init__(self, *args, **kwargs): self.mgxs_lib.by_nuclide = True # Test relevant MGXS types - relevant_MGXS_TYPES = [ - item for item in openmc.mgxs.MGXS_TYPES - if item not in ('current', 'photon-production matrix')] + relevant_MGXS_TYPES = [item for item in openmc.mgxs.MGXS_TYPES + if item != 'current'] # Add in a subset of openmc.mgxs.ARBITRARY_VECTOR_TYPES and # openmc.mgxs.ARBITRARY_MATRIX_TYPES so we can see the code works, # but not use too much resources diff --git a/tests/regression_tests/mgxs_photon/test.py b/tests/regression_tests/mgxs_photon/test.py index b7ea928c7c7..d52dc1e9723 100644 --- a/tests/regression_tests/mgxs_photon/test.py +++ b/tests/regression_tests/mgxs_photon/test.py @@ -2,42 +2,52 @@ import numpy as np import openmc +import pytest from tests.testing_harness import PyAPITestHarness +@pytest.fixture +def model(): + model = openmc.Model() + + material = openmc.Material() + material.set_density('g/cm3', 11.35) + material.add_element('Pb', 1.0) + model.materials.append(material) + + sphere = openmc.Sphere(r=1.0, boundary_type='vacuum') + cell = openmc.Cell(fill=material, region=-sphere) + model.geometry = openmc.Geometry([cell]) + + model.settings.run_mode = 'fixed source' + model.settings.particles = 2000 + model.settings.batches = 2 + model.settings.photon_transport = True + model.settings.atomic_relaxation = True + model.settings.electron_treatment = 'ttb' + model.settings.cutoff = {'energy_photon': 1000.0} + model.settings.source = openmc.IndependentSource( + particle='photon', + space=openmc.stats.Point((0.0, 0.0, 0.0)), + energy=openmc.stats.Discrete([1.0e6], [1.0])) + + return model + + class PhotonMGXSTestHarness(PyAPITestHarness): """Run photon transport and verify photon MGXS post-processing.""" def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) - self.material = openmc.Material() - self.material.set_density('g/cm3', 11.35) - self.material.add_element('Pb', 1.0) - self._model.materials = openmc.Materials([self.material]) - - sphere = openmc.Sphere(r=1.0, boundary_type='vacuum') - cell = openmc.Cell(fill=self.material, region=-sphere) - self._model.geometry = openmc.Geometry([cell]) - - self._model.settings.run_mode = 'fixed source' - self._model.settings.particles = 2000 - self._model.settings.batches = 2 - self._model.settings.photon_transport = True - self._model.settings.atomic_relaxation = True - self._model.settings.electron_treatment = 'ttb' - self._model.settings.cutoff = {'energy_photon': 1000.0} - self._model.settings.source = openmc.IndependentSource( - particle='photon', - space=openmc.stats.Point((0.0, 0.0, 0.0)), - energy=openmc.stats.Discrete([1.0e6], [1.0])) + self.material = self._model.materials[0] groups = openmc.mgxs.EnergyGroups( group_edges=[1.0e3, 1.0e5, 5.0e5, 1.1e6]) self.mgxs_lib = openmc.mgxs.Library( self._model.geometry, - mgxs_types=['total', 'absorption', 'photon-production matrix'], + mgxs_types=['total', 'absorption', 'nu-scatter matrix'], particle_type='photon') self.mgxs_lib.energy_groups = groups self.mgxs_lib.correction = None @@ -50,7 +60,7 @@ def _get_results(self, hash_output=False): self.mgxs_lib.load_from_statepoint(statepoint) production = self.mgxs_lib.get_mgxs( - self.material, 'photon-production matrix') + self.material, 'nu-scatter matrix') absorption = self.mgxs_lib.get_mgxs(self.material, 'absorption') primary = production.tallies[ @@ -96,7 +106,6 @@ def _get_results(self, hash_output=False): return output -def test_photon_mgxs(): - harness = PhotonMGXSTestHarness( - 'statepoint.2.h5', model=openmc.Model()) +def test_photon_mgxs(model): + harness = PhotonMGXSTestHarness('statepoint.2.h5', model) harness.main() diff --git a/tests/unit_tests/test_photon_mgxs.py b/tests/unit_tests/test_photon_mgxs.py index 08a8eccc07d..175ac245422 100644 --- a/tests/unit_tests/test_photon_mgxs.py +++ b/tests/unit_tests/test_photon_mgxs.py @@ -29,13 +29,14 @@ def test_photon_production_matrix_tallies(): def test_photon_production_matrix_constraints(): with pytest.raises(ValueError, match='by nuclide'): openmc.mgxs.MGXS.get_mgxs( - 'photon-production matrix', by_nuclide=True) + 'nu-scatter matrix', by_nuclide=True, particle_type='photon') with pytest.raises(ValueError, match='isotropic'): openmc.mgxs.MGXS.get_mgxs( - 'photon-production matrix', num_polar=2) - with pytest.raises(ValueError, match='photon tallies'): - openmc.mgxs.MGXS.get_mgxs( - 'photon-production matrix', particle_type='neutron') + 'nu-scatter matrix', num_polar=2, particle_type='photon') + + neutron_production = openmc.mgxs.MGXS.get_mgxs( + 'nu-scatter matrix', particle_type='neutron') + assert isinstance(neutron_production, openmc.mgxs.ScatterMatrixXS) def test_photon_production_matrix_combines_primary_and_secondary(): @@ -57,7 +58,7 @@ def test_photon_production_matrix_combines_primary_and_secondary(): production.rxn_rate_tally.mean.ravel(), [2.0, 4.0, 6.0, 8.0]) -def test_set_photon_production_mgxs(monkeypatch): +def test_set_photon_nu_scatter_mgxs(monkeypatch): material = openmc.Material() groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) production = openmc.mgxs.PhotonProductionMatrixXS( @@ -69,7 +70,7 @@ def test_set_photon_production_mgxs(monkeypatch): xsdata.order = 0 absorption = np.array([0.5, 0.6]) xsdata.set_absorption(absorption) - xsdata.set_photon_production_mgxs(production) + xsdata.set_scatter_matrix_mgxs(production) assert xsdata.scatter_format == 'legendre' assert xsdata.order == 0 @@ -84,7 +85,7 @@ def test_photon_mgxs_library(tmp_path): groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) library = openmc.mgxs.Library( geometry, mgxs_types=[ - 'total', 'absorption', 'photon-production matrix'], + 'total', 'absorption', 'nu-scatter matrix'], particle_type='photon') library.domain_type = 'material' library.energy_groups = groups @@ -110,6 +111,5 @@ def test_all_mgxs_types_respects_particle_type(): photon_library = openmc.mgxs.Library( geometry, mgxs_types='all', particle_type='photon') - assert 'photon-production matrix' not in neutron_library.mgxs_types assert photon_library.mgxs_types == ( - 'total', 'absorption', 'photon-production matrix') + 'total', 'absorption', 'nu-scatter matrix') From e9dd574b3ce70efc625bc9ec944840480e418e55 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 16 Aug 2026 16:12:12 -0500 Subject: [PATCH 04/10] Simplify mgxs_photon test --- tests/regression_tests/mgxs_photon/inputs_true.dat | 8 -------- tests/regression_tests/mgxs_photon/test.py | 6 +----- 2 files changed, 1 insertion(+), 13 deletions(-) diff --git a/tests/regression_tests/mgxs_photon/inputs_true.dat b/tests/regression_tests/mgxs_photon/inputs_true.dat index aac04060864..130deaa0698 100644 --- a/tests/regression_tests/mgxs_photon/inputs_true.dat +++ b/tests/regression_tests/mgxs_photon/inputs_true.dat @@ -18,19 +18,11 @@ 2000 2 - - 0.0 0.0 0.0 - 1000000.0 1.0 - ttb - true true - - 1000.0 - diff --git a/tests/regression_tests/mgxs_photon/test.py b/tests/regression_tests/mgxs_photon/test.py index d52dc1e9723..e2d93a7e86c 100644 --- a/tests/regression_tests/mgxs_photon/test.py +++ b/tests/regression_tests/mgxs_photon/test.py @@ -24,13 +24,9 @@ def model(): model.settings.particles = 2000 model.settings.batches = 2 model.settings.photon_transport = True - model.settings.atomic_relaxation = True - model.settings.electron_treatment = 'ttb' - model.settings.cutoff = {'energy_photon': 1000.0} model.settings.source = openmc.IndependentSource( particle='photon', - space=openmc.stats.Point((0.0, 0.0, 0.0)), - energy=openmc.stats.Discrete([1.0e6], [1.0])) + energy=openmc.stats.delta_function(1.0e6)) return model From 2f954c1201c0782ce56ee94c432a09e41426e978 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 16 Aug 2026 23:49:46 -0500 Subject: [PATCH 05/10] Simplify unit tests --- tests/unit_tests/test_photon_mgxs.py | 60 ---------------------------- 1 file changed, 60 deletions(-) diff --git a/tests/unit_tests/test_photon_mgxs.py b/tests/unit_tests/test_photon_mgxs.py index 175ac245422..3bdd79c764f 100644 --- a/tests/unit_tests/test_photon_mgxs.py +++ b/tests/unit_tests/test_photon_mgxs.py @@ -1,7 +1,3 @@ -import h5py -import numpy as np -import pytest - import openmc @@ -26,59 +22,6 @@ def test_photon_production_matrix_tallies(): assert particle_filter.bins == ['photon'] -def test_photon_production_matrix_constraints(): - with pytest.raises(ValueError, match='by nuclide'): - openmc.mgxs.MGXS.get_mgxs( - 'nu-scatter matrix', by_nuclide=True, particle_type='photon') - with pytest.raises(ValueError, match='isotropic'): - openmc.mgxs.MGXS.get_mgxs( - 'nu-scatter matrix', num_polar=2, particle_type='photon') - - neutron_production = openmc.mgxs.MGXS.get_mgxs( - 'nu-scatter matrix', particle_type='neutron') - assert isinstance(neutron_production, openmc.mgxs.ScatterMatrixXS) - - -def test_photon_production_matrix_combines_primary_and_secondary(): - material = openmc.Material() - groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) - production = openmc.mgxs.PhotonProductionMatrixXS( - material, 'material', groups) - - for tally in production.tallies.values(): - values = np.arange(1, tally.num_filter_bins + 1, dtype=float) - values.shape = (tally.num_filter_bins, 1, 1) - tally._sum = 2.0 * values - tally._sum_sq = 2.0 * values**2 - tally._num_realizations = 2 - tally._sp_filename = 'statepoint.h5' - tally._results_read = True - - np.testing.assert_allclose( - production.rxn_rate_tally.mean.ravel(), [2.0, 4.0, 6.0, 8.0]) - - -def test_set_photon_nu_scatter_mgxs(monkeypatch): - material = openmc.Material() - groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) - production = openmc.mgxs.PhotonProductionMatrixXS( - material, 'material', groups) - values = np.array([[0.1, 0.2], [0.3, 0.4]]) - monkeypatch.setattr(production, 'get_xs', lambda **kwargs: values) - - xsdata = openmc.XSdata('photon', groups) - xsdata.order = 0 - absorption = np.array([0.5, 0.6]) - xsdata.set_absorption(absorption) - xsdata.set_scatter_matrix_mgxs(production) - - assert xsdata.scatter_format == 'legendre' - assert xsdata.order == 0 - np.testing.assert_allclose(xsdata._scatter_matrix[0][:, :, 0], values) - np.testing.assert_allclose(xsdata._absorption[0], absorption) - assert xsdata._multiplicity_matrix[0] is None - - def test_photon_mgxs_library(tmp_path): material = openmc.Material() geometry = openmc.Geometry([openmc.Cell(fill=material)]) @@ -99,15 +42,12 @@ def test_photon_mgxs_library(tmp_path): mg_library = openmc.MGXSLibrary(groups, particle_type='photon') path = tmp_path / 'mgxs.h5' mg_library.export_to_hdf5(path) - with h5py.File(path) as h5file: - assert h5file.attrs['particle_type'] == b'photon' assert openmc.MGXSLibrary.from_hdf5(path).particle_type == \ openmc.ParticleType.PHOTON def test_all_mgxs_types_respects_particle_type(): geometry = openmc.Geometry([openmc.Cell()]) - neutron_library = openmc.mgxs.Library(geometry, mgxs_types='all') photon_library = openmc.mgxs.Library( geometry, mgxs_types='all', particle_type='photon') From 33c379303e1a0cb2b2e91e42b0608afa1846e153 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 23 Aug 2026 09:24:59 -0500 Subject: [PATCH 06/10] Fix slicing and condensation on photon production MGXS --- openmc/mgxs/library.py | 3 ++ openmc/mgxs/mgxs.py | 56 ++++++++++++++++++++++++++++ openmc/tallies.py | 2 +- tests/unit_tests/test_photon_mgxs.py | 49 +++++++++++++++++++++++- 4 files changed, 108 insertions(+), 2 deletions(-) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 6338c7d732d..08cf5350d90 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -168,6 +168,7 @@ def __deepcopy__(self, memo): clone._sp_filename = self._sp_filename clone._keff = self._keff clone._sparse = self.sparse + clone._estimator = self.estimator clone._all_mgxs = {} for domain in self.domains: @@ -248,6 +249,8 @@ def particle_type(self, particle_type): (openmc.ParticleType.NEUTRON, openmc.ParticleType.PHOTON)) self._particle_type = particle_type + if particle_type == openmc.ParticleType.PHOTON: + self._correction = None @property def domain_type(self): diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index b4eeae63045..e7ba529ecf7 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -5040,6 +5040,62 @@ def scatter_format(self): def legendre_order(self): return 0 + def _with_energyout_filter(self): + """Return a copy using an outgoing-energy production filter.""" + production = copy.deepcopy(self) + secondary = production.tallies['secondary photon production'] + if not secondary.contains_filter(openmc.ParticleProductionFilter): + return production + particle_production = secondary.find_filter( + openmc.ParticleProductionFilter) + energyout = openmc.EnergyoutFilter(particle_production.energies) + secondary.filters = [ + energyout + if isinstance(f, openmc.ParticleProductionFilter) else f + for f in secondary.filters + ] + return production + + def get_condensed_xs(self, coarse_groups): + """Construct an energy-condensed photon production matrix. + + Parameters + ---------- + coarse_groups : openmc.mgxs.EnergyGroups + Coarse energy group structure + + Returns + ------- + openmc.mgxs.PhotonProductionMatrixXS + Photon production matrix condensed to the coarse group structure + + """ + production = self._with_energyout_filter() + return MGXS.get_condensed_xs(production, coarse_groups) + + def get_slice(self, nuclides=[], in_groups=[], out_groups=[]): + """Build a sliced photon production matrix. + + Parameters + ---------- + nuclides : list of str + Nuclides to include; photon production only supports macroscopic + data + in_groups : list of int + Incoming energy groups to include + out_groups : list of int + Outgoing energy groups to include + + Returns + ------- + openmc.mgxs.PhotonProductionMatrixXS + Sliced photon production matrix + + """ + production = self._with_energyout_filter() + return MatrixMGXS.get_slice( + production, nuclides, in_groups, out_groups) + @property def rxn_rate_tally(self): if self._rxn_rate_tally is None: diff --git a/openmc/tallies.py b/openmc/tallies.py index a503d04bc91..4393dc90b71 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3383,7 +3383,7 @@ def get_slice(self, scores=[], filters=[], filter_bins=[], nuclides=[], # Replace existing filter with new one for j, test_filter in enumerate(new_tally.filters): - if isinstance(test_filter, filter_type): + if type(test_filter) is filter_type: new_tally.filters[j] = new_filter # If original tally was sparse, sparsify the sliced tally diff --git a/tests/unit_tests/test_photon_mgxs.py b/tests/unit_tests/test_photon_mgxs.py index 3bdd79c764f..a1f1e3eb0cc 100644 --- a/tests/unit_tests/test_photon_mgxs.py +++ b/tests/unit_tests/test_photon_mgxs.py @@ -1,6 +1,16 @@ +import copy + +import numpy as np import openmc +def _add_tally_data(mgxs): + for tally in mgxs.tallies.values(): + tally._derived = True + tally._mean = np.ones(tally.shape) + tally._std_dev = np.zeros(tally.shape) + + def test_photon_production_matrix_tallies(): material = openmc.Material() groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) @@ -32,12 +42,23 @@ def test_photon_mgxs_library(tmp_path): particle_type='photon') library.domain_type = 'material' library.energy_groups = groups - library.correction = None library.build_library() library.check_library_for_openmc_mgxs() + assert library.correction is None + + clone = copy.deepcopy(library) + clone.check_library_for_openmc_mgxs() + assert clone.estimator is None for mgxs in library.all_mgxs[material.id].values(): assert mgxs.particle_type == openmc.ParticleType.PHOTON + _add_tally_data(mgxs) + + library._sp_filename = 'statepoint.h5' + coarse_groups = openmc.mgxs.EnergyGroups([1.0, 100.0]) + condensed = library.get_condensed_library(coarse_groups) + xsdata = condensed.create_mg_library().xsdatas[0] + assert xsdata.scatter_matrix[0].shape == (1, 1, 1) mg_library = openmc.MGXSLibrary(groups, particle_type='photon') path = tmp_path / 'mgxs.h5' @@ -53,3 +74,29 @@ def test_all_mgxs_types_respects_particle_type(): assert photon_library.mgxs_types == ( 'total', 'absorption', 'nu-scatter matrix') + assert photon_library.correction is None + + +def test_photon_production_matrix_group_transformations(): + material = openmc.Material() + groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0, 1000.0]) + production = openmc.mgxs.PhotonProductionMatrixXS( + material, 'material', groups) + _add_tally_data(production) + + coarse_groups = openmc.mgxs.EnergyGroups([1.0, 100.0, 1000.0]) + condensed = production.get_condensed_xs(coarse_groups) + assert condensed.get_xs().shape == (2, 2) + assert condensed.tallies['secondary photon production'].contains_filter( + openmc.EnergyoutFilter) + + sliced = production.get_slice(in_groups=[1, 2], out_groups=[1, 2]) + assert sliced.get_xs().shape == (2, 2) + assert sliced.tallies['secondary photon production'].contains_filter( + openmc.EnergyoutFilter) + assert production.tallies[ + 'secondary photon production'].contains_filter( + openmc.ParticleProductionFilter) + + resliced = sliced.get_slice(in_groups=[1], out_groups=[1]) + assert resliced.get_xs().shape == (1, 1) From d02c18115995c0992b30c66f391e98e4bc85e704 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 27 Aug 2026 01:03:10 -0500 Subject: [PATCH 07/10] Rename photon production matrix MGXS --- docs/source/pythonapi/mgxs.rst | 2 +- openmc/mgxs/library.py | 2 +- openmc/mgxs/mgxs.py | 14 +++++++------- openmc/mgxs_library.py | 2 +- tests/unit_tests/test_photon_mgxs.py | 4 ++-- 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/docs/source/pythonapi/mgxs.rst b/docs/source/pythonapi/mgxs.rst index bda07adf35d..10bd54ae280 100644 --- a/docs/source/pythonapi/mgxs.rst +++ b/docs/source/pythonapi/mgxs.rst @@ -51,7 +51,7 @@ Multi-group Cross Sections openmc.mgxs.KappaFissionXS openmc.mgxs.MultiplicityMatrixXS openmc.mgxs.NuFissionMatrixXS - openmc.mgxs.PhotonProductionMatrixXS + openmc.mgxs.PhotonTransferMatrixXS openmc.mgxs.ReducedAbsorptionXS openmc.mgxs.ScatterXS openmc.mgxs.ScatterMatrixXS diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 08cf5350d90..71929d74309 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -549,7 +549,7 @@ def build_library(self): mgxs.energy_groups = self.energy_groups mgxs.by_nuclide = self.by_nuclide if self.estimator is not None and not isinstance( - mgxs, openmc.mgxs.PhotonProductionMatrixXS): + mgxs, openmc.mgxs.PhotonTransferMatrixXS): mgxs.estimator = self.estimator if mgxs_type in openmc.mgxs.MDGXS_TYPES: diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index e7ba529ecf7..3079f8e22ab 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -822,7 +822,7 @@ def get_mgxs(mgxs_type, domain=None, domain_type=None, particle_type is not None and \ openmc.ParticleType(particle_type) == \ openmc.ParticleType.PHOTON: - mgxs = PhotonProductionMatrixXS( + mgxs = PhotonTransferMatrixXS( domain, domain_type, energy_groups, by_nuclide, name, num_polar, num_azimuthal) elif mgxs_type == 'nu-scatter matrix': @@ -4961,8 +4961,8 @@ def print_groups_and_histogram(avg_xs, err_xs, num_groups, @add_params -class PhotonProductionMatrixXS(MatrixMGXS): - r"""A photon production matrix multigroup cross section. +class PhotonTransferMatrixXS(MatrixMGXS): + r"""A photon transfer matrix multigroup cross section. This matrix includes both the photon that survives a coherent or incoherent scattering event and secondary photons banked during photon interactions. @@ -5066,8 +5066,8 @@ def get_condensed_xs(self, coarse_groups): Returns ------- - openmc.mgxs.PhotonProductionMatrixXS - Photon production matrix condensed to the coarse group structure + openmc.mgxs.PhotonTransferMatrixXS + Photon transfer matrix condensed to the coarse group structure """ production = self._with_energyout_filter() @@ -5088,8 +5088,8 @@ def get_slice(self, nuclides=[], in_groups=[], out_groups=[]): Returns ------- - openmc.mgxs.PhotonProductionMatrixXS - Sliced photon production matrix + openmc.mgxs.PhotonTransferMatrixXS + Sliced photon transfer matrix """ production = self._with_energyout_filter() diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 7d3c3d4151b..51c5d4ad92b 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1587,7 +1587,7 @@ def set_scatter_matrix_mgxs(self, scatter, temperature=ROOM_TEMPERATURE_KELVIN, """ check_type('scatter', scatter, (openmc.mgxs.ScatterMatrixXS, - openmc.mgxs.PhotonProductionMatrixXS)) + openmc.mgxs.PhotonTransferMatrixXS)) check_value('energy_groups', scatter.energy_groups, [self.energy_groups]) check_value('domain_type', scatter.domain_type, diff --git a/tests/unit_tests/test_photon_mgxs.py b/tests/unit_tests/test_photon_mgxs.py index a1f1e3eb0cc..b3371a02f5f 100644 --- a/tests/unit_tests/test_photon_mgxs.py +++ b/tests/unit_tests/test_photon_mgxs.py @@ -14,7 +14,7 @@ def _add_tally_data(mgxs): def test_photon_production_matrix_tallies(): material = openmc.Material() groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0]) - production = openmc.mgxs.PhotonProductionMatrixXS( + production = openmc.mgxs.PhotonTransferMatrixXS( material, 'material', groups) assert production.particle_type == openmc.ParticleType.PHOTON @@ -80,7 +80,7 @@ def test_all_mgxs_types_respects_particle_type(): def test_photon_production_matrix_group_transformations(): material = openmc.Material() groups = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0, 1000.0]) - production = openmc.mgxs.PhotonProductionMatrixXS( + production = openmc.mgxs.PhotonTransferMatrixXS( material, 'material', groups) _add_tally_data(production) From e28ce77b6d434141e2042d65b4eeceafdd620c67 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Fri, 18 Sep 2026 11:16:14 -0400 Subject: [PATCH 08/10] Enable multigroup photon transport in MC mode --- include/openmc/mgxs_interface.h | 1 + src/material.cpp | 4 ++-- src/mgxs_interface.cpp | 18 +++++++++++++++--- src/settings.cpp | 5 ----- src/source.cpp | 14 ++++++++++++-- src/tallies/tally.cpp | 6 +++++- 6 files changed, 35 insertions(+), 13 deletions(-) diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index 117ac503d49..dade0843d5e 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -63,6 +63,7 @@ class MgxsInterface { vector> nuc_temps_; // all available temperatures vector default_inverse_velocity_; // approximate default inverse-velocity data + ParticleType particle_type_ {ParticleType::neutron()}; }; namespace data { diff --git a/src/material.cpp b/src/material.cpp index 52cc0ec87bd..cba7faa5577 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -215,7 +215,7 @@ Material::Material(pugi::xml_node node) auto n = names.size(); nuclide_.reserve(n); atom_density_ = tensor::Tensor({n}); - if (settings::photon_transport) + if (settings::photon_transport && settings::run_CE) element_.reserve(n); for (int i = 0; i < n; ++i) { @@ -244,7 +244,7 @@ Material::Material(pugi::xml_node node) // If the corresponding element hasn't been encountered yet and photon // transport will be used, we need to add its symbol to the element_dict - if (settings::photon_transport) { + if (settings::photon_transport && settings::run_CE) { std::string element = to_element(name); // Make sure photon cross section data is available diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 865c565808d..b126644eeca 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -209,6 +209,14 @@ void MgxsInterface::read_header(const std::string& path_cross_sections) // Open file for reading hid_t file_id = file_open(cross_sections_path_, 'r', true); + std::string p_type_str; + if (attribute_exists(file_id, "particle_type")) { + read_attribute(file_id, "particle_type", p_type_str); + if (p_type_str == "photon") { + particle_type_ = ParticleType::photon(); + } + } + ensure_exists(file_id, "energy_groups", true); read_attribute(file_id, "energy_groups", num_energy_groups_); @@ -239,6 +247,10 @@ void MgxsInterface::read_header(const std::string& path_cross_sections) // Calculate approximate default inverse velocity data for (int i = 0; i < energy_bins_.size() - 1; ++i) { + if (particle_type_.is_photon()) { + default_inverse_velocity_.push_back(1.0 / C_LIGHT); + continue; + } double e_min = std::max(energy_bins_[i + 1], 1e-5); double e_max = energy_bins_[i]; double alpha = 1.0 / (C_LIGHT * std::log(e_max / e_min)); @@ -257,9 +269,9 @@ void MgxsInterface::read_header(const std::string& path_cross_sections) void put_mgxs_header_data_to_globals() { // Get the minimum and maximum energies - int neutron = ParticleType::neutron().transport_index(); - data::energy_min[neutron] = data::mg.energy_bins_.back(); - data::energy_max[neutron] = data::mg.energy_bins_.front(); + int particle = data::mg.particle_type_.transport_index(); + data::energy_min[particle] = data::mg.energy_bins_.back(); + data::energy_max[particle] = data::mg.energy_bins_.front(); // Save available XS names to library list, so that when // materials are read, the specified mgxs can be confirmed diff --git a/src/settings.cpp b/src/settings.cpp index 15a7b9c27c9..0d967fcccad 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -611,11 +611,6 @@ void read_settings_xml(pugi::xml_node root) // Check for photon transport if (check_for_node(root, "photon_transport")) { photon_transport = get_node_value_bool(root, "photon_transport"); - - if (!run_CE && photon_transport) { - fatal_error("Photon transport is not currently supported in " - "multigroup mode"); - } } // Check for atomic relaxation diff --git a/src/source.cpp b/src/source.cpp index f951722cc18..a085688c822 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -1224,8 +1224,18 @@ SourceSite sample_external_source(uint64_t* seed) // If running in MG, convert site.E to group if (!settings::run_CE) { - site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(), - data::mg.rev_energy_bins_.end(), site.E); + if (site.particle != data::mg.particle_type_) { + fatal_error(fmt::format( + "Source particle type '{}' does not match the '{}' MGXS library.", + site.particle.str(), data::mg.particle_type_.str())); + } + const auto& bins = data::mg.rev_energy_bins_; + if (site.E < bins.front() || site.E > bins.back()) { + fatal_error(fmt::format( + "Source energy {} eV is outside the MGXS group structure ({}-{} eV).", + site.E, bins.front(), bins.back())); + } + site.E = lower_bound_index(bins.begin(), bins.end(), site.E); site.E = data::mg.num_energy_groups_ - site.E - 1.; } diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 5ffa1b074c6..59148f411c4 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -705,11 +705,15 @@ void Tally::set_scores(const vector& scores) // Make sure all scores are compatible with multigroup mode. if (!settings::run_CE) { - for (auto sc : scores_) + for (auto sc : scores_) { if (sc > 0) fatal_error("Cannot tally " + reaction_name(sc) + " reaction rate " "in multi-group mode"); + if (sc == SCORE_PULSE_HEIGHT) + fatal_error( + "Pulse-height tallies are not supported in multi-group mode."); + } } // Make sure mesh surface tallies contain only current score. From 315debd0e1e4366c149002f6942bb0173c8b2fe5 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Fri, 18 Sep 2026 11:16:48 -0400 Subject: [PATCH 09/10] Add multigroup photon regression test --- tests/regression_tests/mg_photon/__init__.py | 0 .../mg_photon/inputs_true.dat | 46 +++++ .../mg_photon/results_true.dat | 37 ++++ tests/regression_tests/mg_photon/test.py | 194 ++++++++++++++++++ 4 files changed, 277 insertions(+) create mode 100644 tests/regression_tests/mg_photon/__init__.py create mode 100644 tests/regression_tests/mg_photon/inputs_true.dat create mode 100644 tests/regression_tests/mg_photon/results_true.dat create mode 100644 tests/regression_tests/mg_photon/test.py diff --git a/tests/regression_tests/mg_photon/__init__.py b/tests/regression_tests/mg_photon/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/mg_photon/inputs_true.dat b/tests/regression_tests/mg_photon/inputs_true.dat new file mode 100644 index 00000000000..3c06e4c7c6c --- /dev/null +++ b/tests/regression_tests/mg_photon/inputs_true.dat @@ -0,0 +1,46 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + fixed source + 1000 + 20 + + + -5.0 -5.0 -5.0 5.0 5.0 5.0 + + + 1000000.0 1.0 + + + multi-group + true + + + + 1000.0 100000.0 500000.0 2000000.0 + + + photon + + + 1 2 + flux total absorption scatter nu-scatter inverse-velocity + + + diff --git a/tests/regression_tests/mg_photon/results_true.dat b/tests/regression_tests/mg_photon/results_true.dat new file mode 100644 index 00000000000..7776b0f513b --- /dev/null +++ b/tests/regression_tests/mg_photon/results_true.dat @@ -0,0 +1,37 @@ +tally 1: +1.354666E+01 +9.208150E+00 +1.625599E+01 +1.325974E+01 +1.354666E+01 +9.208150E+00 +2.709331E+00 +3.683260E-01 +2.709331E+00 +3.683260E-01 +4.518678E-10 +1.024545E-20 +4.280243E+01 +9.222700E+01 +1.712097E+01 +1.475632E+01 +6.420364E+00 +2.075108E+00 +1.070061E+01 +5.764188E+00 +1.070061E+01 +5.764188E+00 +1.427735E-09 +1.026164E-19 +1.329067E+02 +8.840893E+02 +3.322667E+01 +5.525558E+01 +6.645335E+00 +2.210223E+00 +2.658134E+01 +3.536357E+01 +3.322667E+01 +5.525558E+01 +4.433290E-09 +9.836821E-19 diff --git a/tests/regression_tests/mg_photon/test.py b/tests/regression_tests/mg_photon/test.py new file mode 100644 index 00000000000..db332463071 --- /dev/null +++ b/tests/regression_tests/mg_photon/test.py @@ -0,0 +1,194 @@ +"""Multigroup photon transport with the Monte Carlo solver. + +A reflective cube with a spatially uniform, isotropic source behaves as an +infinite medium, so the group track lengths (volume-integrated fluxes) per +source photon follow from a group balance and can be checked exactly. The +library uses P1 scattering and a top group that produces more photons than it +scatters (a row-constant multiplicity), which exercises angular sampling and +the weight multiplicity used to represent secondary photon production. + +""" +from pathlib import Path + +import numpy as np +import openmc +import pytest + +from openmc.utility_funcs import change_directory +from tests.regression_tests import config +from tests.testing_harness import PyAPITestHarness + +GROUP_EDGES = [1.0e3, 1.0e5, 5.0e5, 2.0e6] + +# Macroscopic cross sections in 1/cm; index 0 is the highest-energy group +TOTAL = np.array([0.25, 0.40, 1.20]) +ABSORPTION = np.array([0.05, 0.15, 1.00]) + +# P0 production (nu-scatter) matrix indexed [g_in, g_out]. The top group +# produces 0.25/cm but only scatters TOTAL - ABSORPTION = 0.20/cm; the extra +# photons stand in for fluorescence and annihilation photons. +PRODUCTION = np.array([ + [0.10, 0.08, 0.07], + [0.00, 0.15, 0.10], + [0.00, 0.00, 0.20], +]) +MU_BAR = 0.2 # average scattering cosine, P1/P0 + +SPEED_OF_LIGHT = 2.99792458e10 # cm/s + + +def _multiplicity(): + """Row-constant multiplicity so each collision yields PRODUCTION/TOTAL""" + ratio = PRODUCTION.sum(axis=1) / (TOTAL - ABSORPTION) + return np.repeat(ratio[:, np.newaxis], len(TOTAL), axis=1) + + +def _exact_track_lengths(): + """Group track length per source photon in an infinite medium. + + Each group balances collisions against source and in-production: + TOTAL[g] * L[g] = S[g] + sum_g' PRODUCTION[g', g] * L[g'] + + """ + source = np.array([1.0, 0.0, 0.0]) # 1 MeV photons are born in group 1 + return np.linalg.solve(np.diag(TOTAL) - PRODUCTION.T, source) + + +def _make_model(): + groups = openmc.mgxs.EnergyGroups(group_edges=GROUP_EDGES) + photon = openmc.XSdata('photon', groups) + photon.order = 1 + photon.set_total(TOTAL) + photon.set_absorption(ABSORPTION) + photon.set_scatter_matrix( + np.stack([PRODUCTION, MU_BAR * PRODUCTION], axis=-1)) + photon.set_multiplicity_matrix(_multiplicity()) + # No inverse velocity is given, so photons must default to 1/c + + library = openmc.MGXSLibrary(groups, particle_type='photon') + library.add_xsdata(photon) + library.export_to_hdf5('mgxs.h5') + + material = openmc.Material() + material.set_density('macro', 1.0) + material.add_macroscopic('photon') + + box = openmc.model.RectangularParallelepiped( + -5.0, 5.0, -5.0, 5.0, -5.0, 5.0, boundary_type='reflective') + cell = openmc.Cell(fill=material, region=-box) + + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.materials = openmc.Materials([material]) + model.materials.cross_sections = 'mgxs.h5' + + tally = openmc.Tally(name='photon') + tally.filters = [ + openmc.EnergyFilter(GROUP_EDGES), + openmc.ParticleFilter('photon'), + ] + tally.scores = ['flux', 'total', 'absorption', 'scatter', 'nu-scatter', + 'inverse-velocity'] + model.tallies.append(tally) + + model.settings.energy_mode = 'multi-group' + model.settings.run_mode = 'fixed source' + model.settings.photon_transport = True + model.settings.batches = 20 + model.settings.particles = 1000 + model.settings.source = openmc.IndependentSource( + particle='photon', + space=openmc.stats.Box((-5.0, -5.0, -5.0), (5.0, 5.0, 5.0)), + energy=openmc.stats.Discrete([1.0e6], [1.0])) + + return model + + +@pytest.fixture +def model(): + return _make_model() + + +class MGPhotonTestHarness(PyAPITestHarness): + """Check exact infinite-medium answers before the regression comparison""" + + def _get_results(self, hash_output=False): + exact = _exact_track_lengths() + expected = { + 'flux': exact, + 'total': TOTAL * exact, + 'absorption': ABSORPTION * exact, + 'scatter': (TOTAL - ABSORPTION) * exact, + 'nu-scatter': PRODUCTION.sum(axis=1) * exact, + } + + with openmc.StatePoint(self._sp_name) as statepoint: + tally = statepoint.get_tally(name='photon') + + # Energy filter bins run from low to high energy, so reverse them + # to match the group ordering of the library + for score, reference in expected.items(): + mean = tally.get_values(scores=[score]).ravel()[::-1] + std_dev = tally.get_values( + scores=[score], value='std_dev').ravel()[::-1] + deviation = np.abs(mean - reference) / std_dev + if np.any(deviation > 5.0): + raise AssertionError( + f"Multigroup photon '{score}' disagrees with the " + f"infinite-medium solution: {mean} vs. {reference} " + f"({deviation} standard deviations)") + + # Photons travel at the speed of light, so the inverse-velocity + # score must equal flux / c in every group + flux = tally.get_values(scores=['flux']).ravel() + inverse_velocity = tally.get_values( + scores=['inverse-velocity']).ravel() + if not np.allclose(SPEED_OF_LIGHT * inverse_velocity, flux, + rtol=1e-10): + raise AssertionError( + 'Multigroup photons are not moving at the speed of light: ' + f'flux / inverse-velocity = {flux / inverse_velocity}') + + return super()._get_results(hash_output) + + def _cleanup(self): + super()._cleanup() + Path('mgxs.h5').unlink(missing_ok=True) + + +def test_mg_photon(model): + harness = MGPhotonTestHarness('statepoint.20.h5', model) + harness.main() + + +@pytest.mark.parametrize( + ('invalid_input', 'error'), + [ + ('energy above library', 'Source energy above range'), + ('energy below library', 'outside the MGXS group structure'), + ('source mismatch', 'does not match'), + ('pulse-height tally', 'Pulse-height tallies are not supported'), + ], +) +def test_mg_photon_validation(tmp_path, invalid_input, error): + with change_directory(tmp_path): + model = _make_model() + source = model.settings.source[0] + if invalid_input == 'energy above library': + source.energy = openmc.stats.Discrete([3.0e6], [1.0]) + elif invalid_input == 'energy below library': + source.energy = openmc.stats.Discrete([500.0], [1.0]) + elif invalid_input == 'source mismatch': + source.particle = 'neutron' + else: + cells = list(model.geometry.get_all_cells().values()) + tally = openmc.Tally() + tally.filters = [ + openmc.CellFilter(cells), + openmc.EnergyFilter(GROUP_EDGES), + ] + tally.scores = ['pulse-height'] + model.tallies.append(tally) + + with pytest.raises(RuntimeError, match=error): + model.run(openmc_exec=config['exe']) From cec81d37572ec8dbd422f765739a73c1b4e0c032 Mon Sep 17 00:00:00 2001 From: Grego01 Date: Wed, 23 Sep 2026 16:17:27 -0400 Subject: [PATCH 10/10] Fix photon MG scoring in tally_scoring and compute multiplicity in xsdata --- src/tallies/tally_scoring.cpp | 16 ++++++++++++++-- src/xsdata.cpp | 20 ++++++++++++++++++++ 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 241cbf0084c..f3eaf911c64 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1772,9 +1772,21 @@ void score_general_mg(Particle& p, int i_tally, int start_index, if (p.event() != TallyEvent::SCATTER) continue; // For scattering production, we need to use the pre-collision weight - // times the multiplicity as the estimate for the number of neutrons - // exiting a reaction with neutrons in the exit channel + // times the multiplicity as the estimate for the number of particles + // exiting a reaction with particles in the exit channel score = (p.wgt_last() - wgt_absorb) * flux; + // Apply the multiplicity of the sampled g_last -> g transfer in the + // material data, which is the factor scatter() applied to the weight. + // Use the same temperature and angle indices that scatter() used. + int t = p.mg_xs_cache().t; + int a = p.mg_xs_cache().a; + double scatt = macro_xs.get_xs( + MgxsType::SCATTER, p.g_last(), &p.g(), nullptr, nullptr, t, a); + if (scatt > 0.0) { + score *= macro_xs.get_xs(MgxsType::NU_SCATTER, p.g_last(), &p.g(), + nullptr, nullptr, t, a) / + scatt; + } // Since we transport based on material data, the angle selected // was not selected from the f(mu) for the nuclide. Therefore // adjust the score by the actual probability for that nuclide. diff --git a/src/xsdata.cpp b/src/xsdata.cpp index 33d063a7b9a..5e8d605f24f 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -119,6 +119,26 @@ void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, for (size_t i = 0; i < total.size(); i++) if (total.data()[i] == 0.0) total.data()[i] = 1.e-10; + + // Photon libraries fold secondary photons into the scatter matrix and store + // no multiplicity. Derive the group-wise one the MC collision game needs. + if (data::mg.particle_type_.is_photon() && + !object_exists(xsdata_grp, "scatter_data/multiplicity_matrix")) { + for (size_t a = 0; a < n_ang; a++) { + for (size_t g = 0; g < energy_groups; g++) { + double production = scatter[a]->scattxs[g]; + if (production <= 0.0) + continue; + double removal = total(a, g) - absorption(a, g); + if (removal <= 0.0) + fatal_error(fmt::format("Photon group {} produces photons but has " + "no scattering to carry them.", + g + 1)); + std::fill(scatter[a]->mult[g].begin(), scatter[a]->mult[g].end(), + production / removal); + } + } + } } //==============================================================================