diff --git a/CMakeLists.txt b/CMakeLists.txt index e29c0d81c68..1b4227159b9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -429,6 +429,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_particle.cpp src/tallies/filter_particle_production.cpp src/tallies/filter_polar.cpp + src/tallies/filter_reaction.cpp src/tallies/filter_sph_harm.cpp src/tallies/filter_sptl_legendre.cpp src/tallies/filter_surface.cpp diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 7bd03f379ea..b3911c8b3ba 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -149,6 +149,7 @@ Constructing Tallies openmc.ZernikeRadialFilter openmc.ParentNuclideFilter openmc.ParticleFilter + openmc.ReactionFilter openmc.MeshMaterialVolumes openmc.Trigger openmc.TallyDerivative diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 465df14d86a..dab45481dd9 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -82,7 +82,9 @@ Classes Nuclide ParentNuclideFilter ParticleFilter + ParticleProductionFilter PolarFilter + ReactionFilter RectilinearMesh RegularMesh SpatialLegendreFilter diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index e3b4e508bc9..5b5630556a9 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -105,109 +105,124 @@ The following tables show all valid scores: .. table:: **Reaction scores: units are reactions per source particle.** - +----------------------+---------------------------------------------------+ - |Score | Description | - +======================+===================================================+ - |absorption |Total absorption rate. For incident neutrons, this | - | |accounts for all reactions that do not produce | - | |secondary neutrons as well as fission. For incident| - | |photons, this includes photoelectric and pair | - | |production. | - +----------------------+---------------------------------------------------+ - |elastic |Elastic scattering reaction rate. | - +----------------------+---------------------------------------------------+ - |fission |Total fission reaction rate. | - +----------------------+---------------------------------------------------+ - |scatter |Total scattering rate. | - +----------------------+---------------------------------------------------+ - |total |Total reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,2nd) |(n,2nd) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,2n) |(n,2n) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,3n) |(n,3n) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,na) |(n,n\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,n3a) |(n,n3\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,2na) |(n,2n\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,3na) |(n,3n\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,np) |(n,np) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,n2a) |(n,n2\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,2n2a) |(n,2n2\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,nd) |(n,nd) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,nt) |(n,nt) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,n3He) |(n,n\ :sup:`3`\ He) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,nd2a) |(n,nd2\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,nt2a) |(n,nt2\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,4n) |(n,4n) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,2np) |(n,2np) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,3np) |(n,3np) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,n2p) |(n,n2p) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,n*X*) |Level inelastic scattering reaction rate. The *X* | - | |indicates what which inelastic level, e.g., (n,n3) | - | |is third-level inelastic scattering. | - +----------------------+---------------------------------------------------+ - |(n,nc) |Continuum level inelastic scattering reaction rate.| - +----------------------+---------------------------------------------------+ - |(n,gamma) |Radiative capture reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,p) |(n,p) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,d) |(n,d) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,t) |(n,t) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,3He) |(n,\ :sup:`3`\ He) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,a) |(n,\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,2a) |(n,2\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,3a) |(n,3\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,2p) |(n,2p) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,pa) |(n,p\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,t2a) |(n,t2\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,d2a) |(n,d2\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,pd) |(n,pd) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,pt) |(n,pt) reaction rate. | - +----------------------+---------------------------------------------------+ - |(n,da) |(n,d\ :math:`\alpha`\ ) reaction rate. | - +----------------------+---------------------------------------------------+ - |coherent-scatter |Coherent (Rayleigh) scattering reaction rate. | - +----------------------+---------------------------------------------------+ - |incoherent-scatter |Incoherent (Compton) scattering reaction rate. | - +----------------------+---------------------------------------------------+ - |photoelectric |Photoelectric absorption reaction rate. | - +----------------------+---------------------------------------------------+ - |pair-production |Pair production reaction rate. | - +----------------------+---------------------------------------------------+ - |*Arbitrary integer* |An arbitrary integer is interpreted to mean the | - | |reaction rate for a reaction with a given ENDF MT | - | |number. | - +----------------------+---------------------------------------------------+ + +------------------------+-------------------------------------------------+ + |Score |Description | + +========================+=================================================+ + |absorption |Total absorption rate. For incident neutrons, | + | |this accounts for all reactions that do not | + | |produce secondary neutrons as well as fission. | + | |For incident photons, this includes | + | |photoelectric and pair production. | + +------------------------+-------------------------------------------------+ + |elastic |Elastic scattering reaction rate. | + +------------------------+-------------------------------------------------+ + |fission |Total fission reaction rate. | + +------------------------+-------------------------------------------------+ + |scatter |Total scattering rate. | + +------------------------+-------------------------------------------------+ + |total |Total reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,2nd) |(n,2nd) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,2n) |(n,2n) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,3n) |(n,3n) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,na) |(n,n\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,n3a) |(n,n3\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,2na) |(n,2n\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,3na) |(n,3n\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,np) |(n,np) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,n2a) |(n,n2\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,2n2a) |(n,2n2\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,nd) |(n,nd) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,nt) |(n,nt) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,n3He) |(n,n\ :sup:`3`\ He) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,nd2a) |(n,nd2\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,nt2a) |(n,nt2\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,4n) |(n,4n) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,2np) |(n,2np) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,3np) |(n,3np) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,n2p) |(n,n2p) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,npa) |(n,np\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,n*X*) |Level inelastic scattering reaction rate. The | + | |*X* indicates which inelastic level, e.g., | + | |(n,n3) is third-level inelastic scattering. | + +------------------------+-------------------------------------------------+ + |(n,nc) |Continuum level inelastic scattering | + | |reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,gamma) |Radiative capture reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,p) |(n,p) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,d) |(n,d) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,t) |(n,t) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,3He) |(n,\ :sup:`3`\ He) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,a) |(n,\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,2a) |(n,2\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,3a) |(n,3\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,2p) |(n,2p) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,pa) |(n,p\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,t2a) |(n,t2\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,d2a) |(n,d2\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,pd) |(n,pd) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,pt) |(n,pt) reaction rate. | + +------------------------+-------------------------------------------------+ + |(n,da) |(n,d\ :math:`\alpha`\ ) reaction rate. | + +------------------------+-------------------------------------------------+ + |photon-total |Total photo-atomic reaction rate. | + +------------------------+-------------------------------------------------+ + |coherent-scatter |Coherent (Rayleigh) scattering reaction rate. | + +------------------------+-------------------------------------------------+ + |incoherent-scatter |Incoherent (Compton) scattering reaction rate. | + +------------------------+-------------------------------------------------+ + |photoelectric |Photoelectric absorption reaction rate. | + +------------------------+-------------------------------------------------+ + |photoelectric-*S* |Subshell photoelectric absorption rate for the | + | |*S* shell. For example, "photoelectric-N3" is the| + | |rate for the N3 subshell. | + +------------------------+-------------------------------------------------+ + |pair-production |Pair production reaction rate (total). | + +------------------------+-------------------------------------------------+ + |pair-production-electron|Pair production reaction rate in the electron | + | |field. | + +------------------------+-------------------------------------------------+ + |pair-production-nuclear |Pair production reaction rate in the nuclear | + | |field. | + +------------------------+-------------------------------------------------+ + |*Arbitrary integer* |An arbitrary integer is interpreted to mean the | + | |reaction rate for a reaction with a given ENDF | + | |MT number. | + +------------------------+-------------------------------------------------+ .. table:: **Particle production scores: units are particles produced per source particles.** diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 38238598eba..2adb5ef0fbc 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -227,6 +227,7 @@ enum ReactionType { N_XA = 207, HEATING = 301, DAMAGE_ENERGY = 444, + PHOTON_TOTAL = 501, COHERENT = 502, INCOHERENT = 504, PAIR_PROD_ELEC = 515, diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 4a737eb8816..34fee9758b9 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -33,6 +33,13 @@ bool is_disappearance(int MT); //! \return Whether corresponding reaction is an inelastic scattering reaction bool is_inelastic_scatter(int MT); +//! Determine whether an MT number matches a target MT, considering that the +//! target may be a summation reaction. +//! \param[in] event_mt MT number of the actual event +//! \param[in] target_mt MT number to check against +//! \return Whether event_mt is a component of target_mt (or equal to it) +bool mt_matches(int event_mt, int target_mt); + //============================================================================== //! Abstract one-dimensional function //============================================================================== diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 3314d18666f..e95db1665d7 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -76,11 +76,21 @@ class Reaction { //! \return Name of the corresponding reaction std::string reaction_name(int mt); -//! Return reaction type (MT value) given a reaction name +//! Return MT value for given a reaction name (including special tally MT +//! values) // //! \param[in] name Reaction name -//! \return Corresponding reaction type (MT value) -int reaction_type(std::string name); +//! \return Corresponding MT number or special tally score +int reaction_tally_mt(std::string name); + +//! Return ENDF MT number given a reaction name +// +//! Unlike reaction_tally_mt(), this function always returns a positive ENDF MT +//! number and never a special negative tally score. +//! +//! \param[in] name Reaction name +//! \return Corresponding ENDF MT number +int reaction_mt(const std::string& name); } // namespace openmc diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 0fef97e502f..77b0d9f420d 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -41,6 +41,7 @@ enum class FilterType { PARTICLE, PARTICLE_PRODUCTION, POLAR, + REACTION, SPHERICAL_HARMONICS, SPATIAL_LEGENDRE, SURFACE, diff --git a/include/openmc/tallies/filter_reaction.h b/include/openmc/tallies/filter_reaction.h new file mode 100644 index 00000000000..015729eefb1 --- /dev/null +++ b/include/openmc/tallies/filter_reaction.h @@ -0,0 +1,52 @@ +#ifndef OPENMC_TALLIES_FILTER_REACTION_H +#define OPENMC_TALLIES_FILTER_REACTION_H + +#include "openmc/span.h" +#include "openmc/tallies/filter.h" +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +//! Bins tally events based on the reaction type (MT number). +//============================================================================== + +class ReactionFilter : public Filter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~ReactionFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "reaction"; } + FilterType type() const override { return FilterType::REACTION; } + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle& p, TallyEstimator estimator, + FilterMatch& match) const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //---------------------------------------------------------------------------- + // Accessors + + const vector& bins() const { return bins_; } + void set_bins(span bins); + +protected: + //---------------------------------------------------------------------------- + // Data members + + //! MT numbers to match + vector bins_; +}; + +} // namespace openmc + +#endif // OPENMC_TALLIES_FILTER_REACTION_H diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 65b59582cf8..ab1285b7219 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -21,6 +21,7 @@ from .kalbach_mann import KalbachMann from .laboratory import LaboratoryAngleEnergy from .nbody import NBodyPhaseSpace +from .photon import _SUBSHELLS from .product import Product from .uncorrelated import UncorrelatedAngleEnergy @@ -54,6 +55,10 @@ 198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 203: '(n,Xp)', 204: '(n,Xd)', 205: '(n,Xt)', 206: '(n,X3He)', 207: '(n,Xa)', 301: 'heating', 444: 'damage-energy', + 501: 'photon-total', 502: 'coherent-scatter', + 504: 'incoherent-scatter', 515: 'pair-production-electron', + 516: 'pair-production', 517: 'pair-production-nuclear', + 522: 'photoelectric', 649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)', 849: '(n,ac)', 891: '(n,2nc)', 901: 'heating-local'} REACTION_NAME.update({i: f'(n,n{i - 50})' for i in range(51, 91)}) @@ -63,9 +68,16 @@ REACTION_NAME.update({i: f'(n,3He{i - 750})' for i in range(750, 799)}) REACTION_NAME.update({i: f'(n,a{i - 800})' for i in range(800, 849)}) REACTION_NAME.update({i: f'(n,2n{i - 875})' for i in range(875, 891)}) +REACTION_NAME.update( + {534 + i: f'photoelectric-{shell}' for i, shell in enumerate(_SUBSHELLS[1:])} +) REACTION_MT = {name: mt for mt, name in REACTION_NAME.items()} +REACTION_MT['total'] = 1 +REACTION_MT['elastic'] = 2 REACTION_MT['fission'] = 18 +REACTION_MT['absorption'] = 27 +REACTION_MT['capture'] = 102 FISSION_MTS = (18, 19, 20, 21, 38) diff --git a/openmc/filter.py b/openmc/filter.py index 9a2836e57c8..8fa2e7e626c 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -13,6 +13,7 @@ import openmc import openmc.checkvalue as cv from .cell import Cell +from .data.reaction import REACTION_NAME, REACTION_MT from .material import Material from .mixin import IDManagerMixin from .surface import Surface @@ -22,11 +23,11 @@ _FILTER_TYPES = ( 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', - 'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', - 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', - 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', + 'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', + 'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', + 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide', - 'weight', 'meshborn', 'meshsurface', 'meshmaterial', + 'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction', ) _CURRENT_NAMES = ( @@ -1415,6 +1416,80 @@ def check_bins(self, bins): cv.check_greater_than('filter value', x, 0, equality=True) +class ReactionFilter(Filter): + """Bins tally events based on the reaction type (MT number). + + .. versionadded:: 0.15.4 + + Parameters + ---------- + bins : str, int, or iterable thereof + The reaction types to tally. Can be reaction name strings + (e.g., ``'(n,elastic)'``, ``'(n,gamma)'``) or integer MT numbers + (e.g., 2, 102). Integer MT values are automatically converted to their + canonical string representation. + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + bins : numpy.ndarray of str + Reaction name strings + id : int + Unique identifier for the filter + num_bins : int + The number of filter bins + + """ + + def __init__(self, bins, filter_id=None): + self.bins = bins + self.id = filter_id + + @Filter.bins.setter + def bins(self, bins): + if isinstance(bins, (str, Integral)): + bins = [bins] + elif not isinstance(bins, list): + bins = list(bins) + normalized = [] + for b in bins: + if isinstance(b, Integral): + if int(b) not in REACTION_NAME: + raise ValueError(f"No known reaction for MT={b}") + normalized.append(REACTION_NAME[int(b)]) + elif isinstance(b, str): + if b == 'total': + warnings.warn( + "The reaction name 'total' is ambiguous. Use " + "'(n,total)' for neutron total cross section or " + "'photon-total' for photon total. Interpreting as" + "'(n,total)'.") + if b not in REACTION_MT: + raise ValueError(f"Unknown reaction name '{b}'") + normalized.append(REACTION_NAME[REACTION_MT[b]]) + else: + raise TypeError(f"Expected str or int for reaction filter " + f"bin, got {type(b)}") + self._bins = np.array(normalized, dtype=str) + + @classmethod + def from_hdf5(cls, group, **kwargs): + if group['type'][()].decode() != cls.short_name.lower(): + raise ValueError("Expected HDF5 data for filter type '" + + cls.short_name.lower() + "' but got '" + + group['type'][()].decode() + "' instead") + bins = [b.decode() for b in group['bins'][()]] + filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + return cls(bins, filter_id=filter_id) + + @classmethod + def from_xml_element(cls, elem, **kwargs): + filter_id = int(get_text(elem, "id")) + bins = get_elem_list(elem, "bins", str) or [] + return cls(bins, filter_id=filter_id) + + class RealFilter(Filter): """Tally modifier that describes phase-space and other characteristics diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index dc81c3487b9..574a37443ae 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -22,9 +22,10 @@ 'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter', 'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', - 'ParentNuclideFilter', 'ParticleFilter', 'PolarFilter', 'SphericalHarmonicsFilter', - 'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', 'UniverseFilter', - 'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' + 'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PolarFilter', + 'ReactionFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', + 'SurfaceFilter', 'TimeFilter', 'UniverseFilter', 'WeightFilter', 'ZernikeFilter', + 'ZernikeRadialFilter', 'filters' ] # Tally functions @@ -607,10 +608,18 @@ def bins(self): return [ParticleType(i) for i in particle_i] +class ParticleProductionFilter(Filter): + filter_type = 'particleproduction' + + class PolarFilter(Filter): filter_type = 'polar' +class ReactionFilter(Filter): + filter_type = 'reaction' + + class SphericalHarmonicsFilter(Filter): filter_type = 'sphericalharmonics' @@ -711,7 +720,9 @@ class ZernikeRadialFilter(ZernikeFilter): 'musurface': MuSurfaceFilter, 'parentnuclide': ParentNuclideFilter, 'particle': ParticleFilter, + 'particleproduction': ParticleProductionFilter, 'polar': PolarFilter, + 'reaction': ReactionFilter, 'sphericalharmonics': SphericalHarmonicsFilter, 'spatiallegendre': SpatialLegendreFilter, 'surface': SurfaceFilter, diff --git a/src/chain.cpp b/src/chain.cpp index a60ef12cd6c..4214bc47321 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -42,7 +42,7 @@ ChainNuclide::ChainNuclide(pugi::xml_node node) branching_ratio = std::stod(get_node_value(reaction_node, "branching_ratio")); } - int mt = reaction_type(rx_name); + int mt = reaction_mt(rx_name); reaction_products_[mt].push_back({rx_target, branching_ratio}); } diff --git a/src/endf.cpp b/src/endf.cpp index c0c1d2e7e8b..f1800545b65 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -88,6 +88,83 @@ bool is_inelastic_scatter(int mt) } } +bool mt_matches(int event_mt, int target_mt) +{ + // Direct match + if (event_mt == target_mt) + return true; + + // Check if event_mt is a component of target_mt summation reaction + switch (target_mt) { + case TOTAL_XS: + return event_mt == ELASTIC || mt_matches(event_mt, N_NONELASTIC); + + case N_NONELASTIC: { + static constexpr int components[] = {4, 5, 11, 16, 17, 22, 23, 24, 25, 27, + 28, 29, 30, 32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152, 153, 154, 156, + 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, + 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, + 188, 189, 190, 194, 195, 196, 198, 199, 200}; + for (int mt : components) { + if (mt_matches(event_mt, mt)) + return true; + } + return false; + } + + case N_LEVEL: + // Inelastic scattering levels + return event_mt >= 50 && event_mt <= N_NC; + + case N_2N: + // (n,2n) to excited states + return event_mt >= N_2N0 && event_mt <= N_2NC; + + case N_FISSION: + return is_fission(event_mt); + + case 27: + return is_fission(event_mt) || is_disappearance(event_mt); + + case N_DISAPPEAR: { + return is_disappearance(event_mt); + } + + case N_P: + // (n,p) to excited states + return event_mt >= N_P0 && event_mt <= N_PC; + + case N_D: + // (n,d) to excited states + return event_mt >= N_D0 && event_mt <= N_DC; + + case N_T: + // (n,t) to excited states + return event_mt >= N_T0 && event_mt <= N_TC; + + case N_3HE: + // (n,3He) to excited states + return event_mt >= N_3HE0 && event_mt <= N_3HEC; + + case N_A: + // (n,alpha) to excited states + return event_mt >= N_A0 && event_mt <= N_AC; + + case 501: + return event_mt == 502 || event_mt == 504 || mt_matches(event_mt, 516) || + mt_matches(event_mt, 522); + + case PAIR_PROD: + return event_mt == PAIR_PROD_ELEC || event_mt == PAIR_PROD_NUC; + + case PHOTOELECTRIC: + return event_mt >= 534 && event_mt < 573; + + default: + return false; + } +} + unique_ptr read_function(hid_t group, const char* name) { hid_t obj_id = open_object(group, name); diff --git a/src/reaction.cpp b/src/reaction.cpp index 9ac5a1f528a..c02e0cc407e 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -309,6 +309,7 @@ std::unordered_map REACTION_NAME_MAP { {N_XA, "(n,Xa)"}, {HEATING, "heating"}, {DAMAGE_ENERGY, "damage-energy"}, + {PHOTON_TOTAL, "photon-total"}, {COHERENT, "coherent-scatter"}, {INCOHERENT, "incoherent-scatter"}, {PAIR_PROD_ELEC, "pair-production-electron"}, @@ -346,13 +347,24 @@ void initialize_maps() // Create photoelectric subshells for (int mt = 534; mt <= 572; ++mt) { REACTION_NAME_MAP[mt] = - fmt::format("photoelectric, {} subshell", SUBSHELLS[mt - 534]); + fmt::format("photoelectric-{}", SUBSHELLS[mt - 534]); } // Invert name map to create type map for (const auto& kv : REACTION_NAME_MAP) { REACTION_TYPE_MAP[kv.second] = kv.first; } + + // Alternate names + REACTION_TYPE_MAP["elastic"] = ELASTIC; + REACTION_TYPE_MAP["n2n"] = N_2N; + REACTION_TYPE_MAP["n3n"] = N_3N; + REACTION_TYPE_MAP["n4n"] = N_4N; + REACTION_TYPE_MAP["H1-production"] = N_XP; + REACTION_TYPE_MAP["H2-production"] = N_XD; + REACTION_TYPE_MAP["H3-production"] = N_XT; + REACTION_TYPE_MAP["He3-production"] = N_X3HE; + REACTION_TYPE_MAP["He4-production"] = N_XA; } std::string reaction_name(int mt) @@ -370,62 +382,42 @@ std::string reaction_name(int mt) } } -int reaction_type(std::string name) +int reaction_tally_mt(std::string name) { - // Initialize remainder of name map and all of type map + // All "total" scores should map to the special SCORE_TOTAL + if (name == "total" || name == "(n,total)" || name == "photon-total") + return SCORE_TOTAL; + + // All fission scores should map to the special SCORE_FISSION + if (name == "fission" || name == "(n,fission)") + return SCORE_FISSION; + + // Delegate everything else to reaction_mt() + return reaction_mt(name); +} + +int reaction_mt(const std::string& name) +{ + // Initialize maps if needed if (REACTION_TYPE_MAP.empty()) initialize_maps(); - // (n,total) exists in REACTION_TYPE_MAP for MT=1, but we need this to return - // the special SCORE_TOTAL score - if (name == "(n,total)") - return SCORE_TOTAL; - - // Check if type map has an entry for this reaction name + // Look up directly in type map (no score indirection) auto it = REACTION_TYPE_MAP.find(name); if (it != REACTION_TYPE_MAP.end()) { - return it->second; - } - - // Alternate names for several reactions - if (name == "elastic") { - return ELASTIC; - } else if (name == "n2n") { - return N_2N; - } else if (name == "n3n") { - return N_3N; - } else if (name == "n4n") { - return N_4N; - } else if (name == "H1-production") { - return N_XP; - } else if (name == "H2-production") { - return N_XD; - } else if (name == "H3-production") { - return N_XT; - } else if (name == "He3-production") { - return N_X3HE; - } else if (name == "He4-production") { - return N_XA; + int mt = it->second; + return mt; } - // Assume the given string is a reaction MT number. Make sure it's a natural - // number then return. + // Assume the given string is an MT number int MT = 0; try { MT = std::stoi(name); } catch (const std::invalid_argument& ex) { - throw std::invalid_argument( - "Invalid tally score \"" + name + - "\". See the docs " - "for details: " - "https://docs.openmc.org/en/stable/usersguide/tallies.html#scores"); + throw std::invalid_argument("Unknown reaction name \"" + name + "\"."); } if (MT < 1) - throw std::invalid_argument( - "Invalid tally score \"" + name + - "\". See the docs " - "for details: " - "https://docs.openmc.org/en/stable/usersguide/tallies.html#scores"); + throw std::invalid_argument("Unknown reaction name \"" + name + "\"."); return MT; } diff --git a/src/settings.cpp b/src/settings.cpp index 5b472468fc4..96b18c56bfd 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -962,7 +962,7 @@ void read_settings_xml(pugi::xml_node root) if (check_for_node(node_ct, "reactions")) { auto temp = get_node_array(node_ct, "reactions"); for (const auto& b : temp) { - int reaction_int = reaction_type(b); + int reaction_int = reaction_mt(b); if (reaction_int > 0) { collision_track_config.mt_numbers.insert(reaction_int); } diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index b0e9624c9e2..badb9107733 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -33,6 +33,7 @@ #include "openmc/tallies/filter_particle.h" #include "openmc/tallies/filter_particle_production.h" #include "openmc/tallies/filter_polar.h" +#include "openmc/tallies/filter_reaction.h" #include "openmc/tallies/filter_sph_harm.h" #include "openmc/tallies/filter_sptl_legendre.h" #include "openmc/tallies/filter_surface.h" @@ -151,6 +152,8 @@ Filter* Filter::create(const std::string& type, int32_t id) return Filter::create(id); } else if (type == "polar") { return Filter::create(id); + } else if (type == "reaction") { + return Filter::create(id); } else if (type == "surface") { return Filter::create(id); } else if (type == "spatiallegendre") { diff --git a/src/tallies/filter_reaction.cpp b/src/tallies/filter_reaction.cpp new file mode 100644 index 00000000000..8ee9f3ce87e --- /dev/null +++ b/src/tallies/filter_reaction.cpp @@ -0,0 +1,79 @@ +#include "openmc/tallies/filter_reaction.h" + +#include + +#include "openmc/capi.h" +#include "openmc/endf.h" +#include "openmc/error.h" +#include "openmc/reaction.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +//============================================================================== +// ReactionFilter implementation +//============================================================================== + +void ReactionFilter::from_xml(pugi::xml_node node) +{ + // Read bins as reaction name strings + auto bins_str = get_node_array(node, "bins"); + + // Convert reaction names to MT numbers + vector bins_mt; + bins_mt.reserve(bins_str.size()); + for (const auto& name : bins_str) { + bins_mt.push_back(reaction_mt(name)); + } + + this->set_bins(bins_mt); +} + +void ReactionFilter::set_bins(span bins) +{ + // Clear existing bins + bins_.clear(); + bins_.reserve(bins.size()); + + // Copy bins and build lookup map + for (int64_t i = 0; i < bins.size(); ++i) { + bins_.push_back(bins[i]); + } + + n_bins_ = bins_.size(); +} + +void ReactionFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + // Get the event MT number from the particle + int event_mt = p.event_mt(); + + // Check each bin, considering summation rules + for (int64_t i = 0; i < bins_.size(); ++i) { + if (mt_matches(event_mt, bins_[i])) { + match.bins_.push_back(i); + match.weights_.push_back(1.0); + } + } +} + +void ReactionFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + + // Write bins as reaction name strings for human readability + vector names; + names.reserve(bins_.size()); + for (auto mt : bins_) { + names.push_back(reaction_name(mt)); + } + write_dataset(filter_group, "bins", names); +} + +std::string ReactionFilter::text_label(int bin) const +{ + return reaction_name(bins_[bin]); +} + +} // namespace openmc diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 8cdf3a9050c..fba474cde7f 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -173,6 +173,10 @@ Tally::Tally(pugi::xml_node node) estimator_ = TallyEstimator::COLLISION; } else if (filt_type == FilterType::PARTICLE_PRODUCTION) { estimator_ = TallyEstimator::ANALOG; + } else if (filt_type == FilterType::REACTION) { + if (estimator_ == TallyEstimator::TRACKLENGTH) { + estimator_ = TallyEstimator::COLLISION; + } } } @@ -571,7 +575,7 @@ void Tally::set_scores(const vector& scores) } // Determine integer code for score - int score = reaction_type(score_str); + int score = reaction_tally_mt(score_str); switch (score) { case SCORE_FLUX: @@ -787,7 +791,7 @@ void Tally::init_triggers(pugi::xml_node node) } else { int i_score = 0; for (; i_score < this->scores_.size(); ++i_score) { - if (this->scores_[i_score] == reaction_type(score_str)) + if (this->scores_[i_score] == reaction_tally_mt(score_str)) break; } if (i_score == this->scores_.size()) { diff --git a/tests/regression_tests/filter_reaction/__init__.py b/tests/regression_tests/filter_reaction/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/filter_reaction/inputs_true.dat b/tests/regression_tests/filter_reaction/inputs_true.dat new file mode 100644 index 00000000000..3ab9ab3fb6e --- /dev/null +++ b/tests/regression_tests/filter_reaction/inputs_true.dat @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + eigenvalue + 1000 + 5 + + + + (n,elastic) (n,2n) (n,fission) (n,gamma) (n,total) + + + 1 + flux + + + diff --git a/tests/regression_tests/filter_reaction/results_true.dat b/tests/regression_tests/filter_reaction/results_true.dat new file mode 100644 index 00000000000..61ec1c65fa3 --- /dev/null +++ b/tests/regression_tests/filter_reaction/results_true.dat @@ -0,0 +1,13 @@ +k-combined: +2.265297E+00 7.172807E-03 +tally 1: +8.364316E+01 +1.400120E+03 +1.412312E-01 +4.277446E-03 +5.471503E+00 +7.015849E+00 +0.000000E+00 +0.000000E+00 +1.349840E+02 +3.644648E+03 diff --git a/tests/regression_tests/filter_reaction/test.py b/tests/regression_tests/filter_reaction/test.py new file mode 100644 index 00000000000..3c60fd5173e --- /dev/null +++ b/tests/regression_tests/filter_reaction/test.py @@ -0,0 +1,31 @@ +import openmc + +from tests.testing_harness import PyAPITestHarness + + +def test_filter_reaction(): + model = openmc.Model() + + m = openmc.Material() + m.set_density('g/cm3', 10.0) + m.add_nuclide('U235', 1.0) + model.materials.append(m) + + s = openmc.Sphere(r=100.0, boundary_type='vacuum') + c = openmc.Cell(fill=m, region=-s) + model.geometry = openmc.Geometry([c]) + + # Create a tally with reaction filter + tally = openmc.Tally() + tally.filters = [openmc.ReactionFilter( + ['(n,elastic)', '(n,2n)', '(n,fission)', '(n,gamma)', 'total'] + )] + tally.scores = ['flux'] + model.tallies = openmc.Tallies([tally]) + + # Reduce particles for faster testing + model.settings.particles = 1000 + model.settings.batches = 5 + + harness = PyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/unit_tests/test_filter_reaction.py b/tests/unit_tests/test_filter_reaction.py new file mode 100644 index 00000000000..df8fc9d7b41 --- /dev/null +++ b/tests/unit_tests/test_filter_reaction.py @@ -0,0 +1,87 @@ +import pytest +import openmc + + +def test_reaction_filter_construction_with_strings(): + f = openmc.ReactionFilter(['(n,elastic)', '(n,gamma)']) + assert len(f.bins) == 2 + assert f.bins[0] == '(n,elastic)' + assert f.bins[1] == '(n,gamma)' + + +def test_reaction_filter_construction_with_mt(): + f = openmc.ReactionFilter([2, 102]) + assert len(f.bins) == 2 + assert f.bins[0] == '(n,elastic)' + assert f.bins[1] == '(n,gamma)' + + +def test_reaction_filter_mixed(): + f = openmc.ReactionFilter([2, '(n,gamma)']) + assert f.bins[0] == '(n,elastic)' + assert f.bins[1] == '(n,gamma)' + + +def test_reaction_filter_single_bin_string(): + f = openmc.ReactionFilter('(n,elastic)') + assert len(f.bins) == 1 + assert f.bins[0] == '(n,elastic)' + + +def test_reaction_filter_single_bin_mt(): + f = openmc.ReactionFilter(2) + assert len(f.bins) == 1 + assert f.bins[0] == '(n,elastic)' + + +def test_reaction_filter_single_bin_naming(): + f = openmc.ReactionFilter('total') + assert len(f.bins) == 1 + assert f.bins[0] == '(n,total)' + + +def test_reaction_filter_invalid_mt(): + with pytest.raises(ValueError, match="No known reaction"): + openmc.ReactionFilter([999999]) + + +def test_reaction_filter_invalid_string(): + with pytest.raises(ValueError, match="Unknown reaction name"): + openmc.ReactionFilter(['not-a-reaction']) + + +def test_reaction_filter_invalid_type(): + with pytest.raises(TypeError, match="Expected str or int"): + openmc.ReactionFilter([3.14]) + + +def test_reaction_filter_xml_roundtrip(): + f = openmc.ReactionFilter([2, 102], filter_id=42) + elem = f.to_xml_element() + f2 = openmc.ReactionFilter.from_xml_element(elem) + assert f2.id == 42 + assert len(f2.bins) == 2 + assert f2.bins[0] == '(n,elastic)' + assert f2.bins[1] == '(n,gamma)' + + +def test_reaction_filter_num_bins(): + f = openmc.ReactionFilter(['(n,elastic)', '(n,fission)', '(n,gamma)']) + assert f.num_bins == 3 + + +def test_reaction_filter_repr(): + f = openmc.ReactionFilter([2, 102]) + r = repr(f) + assert 'ReactionFilter' in r + + +def test_reaction_filter_short_name(): + assert openmc.ReactionFilter.short_name == 'Reaction' + + +def test_reaction_filter_total_warning(): + """Test that using 'total' emits a warning about ambiguity.""" + with pytest.warns(UserWarning, match="ambiguous"): + f = openmc.ReactionFilter(['total']) + assert f.bins[0] == '(n,total)'