diff --git a/docs/source/pythonapi/mgxs.rst b/docs/source/pythonapi/mgxs.rst index 4141aa0a09a..10bd54ae280 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.PhotonTransferMatrixXS openmc.mgxs.ReducedAbsorptionXS openmc.mgxs.ScatterXS openmc.mgxs.ScatterMatrixXS 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/openmc/mgxs/library.py b/openmc/mgxs/library.py index faa83c0481f..71929d74309 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.PhotonTransferMatrixXS): 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..3079f8e22ab 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 = PhotonTransferMatrixXS( + 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 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. + 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.PhotonTransferMatrixXS + Photon transfer 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.PhotonTransferMatrixXS + Sliced photon transfer 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..51c5d4ad92b 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.PhotonTransferMatrixXS)) 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/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. 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); + } + } + } } //============================================================================== 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']) 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..b3371a02f5f --- /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.PhotonTransferMatrixXS( + 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.PhotonTransferMatrixXS( + 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)