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..08cf5350d90 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) @@ -159,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: @@ -203,7 +213,11 @@ 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', 'nu-scatter matrix') + else: + self._mgxs_types = all_mgxs_types else: cv.check_iterable_type('mgxs_types', mgxs_types, str) for mgxs_type in mgxs_types: @@ -224,6 +238,20 @@ 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 + if particle_type == openmc.ParticleType.PHOTON: + self._correction = None + @property def domain_type(self): return self._domain_type @@ -513,13 +541,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,8 +1224,18 @@ 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 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, + nuclide=[nuclide], subdomain=subdomain) + 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, @@ -1231,7 +1271,9 @@ def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro', else: using_multiplicity = False - if using_multiplicity: + if self.particle_type == openmc.ParticleType.PHOTON: + pass + elif using_multiplicity: if 'nu-scatter matrix' in self.mgxs_types: nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix') else: @@ -1369,7 +1411,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 +1628,39 @@ def check_library_for_openmc_mgxs(self): error_flag = False + if self.particle_type == openmc.ParticleType.PHOTON: + photon_mgxs_types = { + '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.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', '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.') + error_flag = True + if error_flag: + raise ValueError('Invalid photon MGXS configuration ' + 'encountered.') + return + # 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..e7ba529ecf7 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -177,6 +177,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 +198,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 +266,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 +290,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 +314,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 +481,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 +632,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 +744,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 +778,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 ------- @@ -790,6 +818,13 @@ 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' 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': @@ -835,6 +870,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 +4960,168 @@ 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 = 'nu-scatter' + self._mgxs_type = 'nu-scatter 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 scatter_format(self): + return SCATTER_LEGENDRE + + @property + 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: + 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..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, @@ -2364,6 +2365,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 +2374,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 +2397,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 +2430,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 +2614,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 +2657,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/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/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..130deaa0698 --- /dev/null +++ b/tests/regression_tests/mgxs_photon/inputs_true.dat @@ -0,0 +1,87 @@ + + + + + + + + + + + + + + + + + fixed source + 2000 + 2 + + + 1000000.0 1.0 + + + true + + + + 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..e2d93a7e86c --- /dev/null +++ b/tests/regression_tests/mgxs_photon/test.py @@ -0,0 +1,107 @@ +import hashlib + +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.source = openmc.IndependentSource( + particle='photon', + energy=openmc.stats.delta_function(1.0e6)) + + 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 = 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', 'nu-scatter 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, 'nu-scatter 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(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 new file mode 100644 index 00000000000..a1f1e3eb0cc --- /dev/null +++ b/tests/unit_tests/test_photon_mgxs.py @@ -0,0 +1,102 @@ +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]) + 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_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', 'nu-scatter matrix'], + particle_type='photon') + library.domain_type = 'material' + library.energy_groups = groups + 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' + mg_library.export_to_hdf5(path) + assert openmc.MGXSLibrary.from_hdf5(path).particle_type == \ + openmc.ParticleType.PHOTON + + +def test_all_mgxs_types_respects_particle_type(): + geometry = openmc.Geometry([openmc.Cell()]) + photon_library = openmc.mgxs.Library( + geometry, mgxs_types='all', particle_type='photon') + + 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)