crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-21 06:25:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Implement tally filter for filtering by reaction (#3809)

Browse source
This commit is contained in:
Paul Romano 2026-02-21 10:17:12 -06:00 • committed by GitHub
parent 153281a490
commit 139907c955
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
24 changed files with 661 additions and 160 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -412,6 +412,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

1
docs/source/pythonapi/base.rst
View file

@ -149,6 +149,7 @@ Constructing Tallies
openmc.ZernikeRadialFilter
openmc.ParentNuclideFilter
openmc.ParticleFilter
openmc.ReactionFilter
openmc.MeshMaterialVolumes
openmc.Trigger
openmc.TallyDerivative

2
docs/source/pythonapi/capi.rst
View file

@ -82,7 +82,9 @@ Classes
Nuclide
ParentNuclideFilter
ParticleFilter
ParticleProductionFilter
PolarFilter
ReactionFilter
RectilinearMesh
RegularMesh
SpatialLegendreFilter

221
docs/source/usersguide/tallies.rst
View file

@ -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.**

1
include/openmc/constants.h
View file

@ -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,

7
include/openmc/endf.h
View file

@ -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
//==============================================================================

16
include/openmc/reaction.h
View file

@ -76,11 +76,21 @@ public:
//! \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

1
include/openmc/tallies/filter.h
View file

@ -41,6 +41,7 @@ enum class FilterType {
PARTICLE,
PARTICLE_PRODUCTION,
POLAR,
REACTION,
SPHERICAL_HARMONICS,
SPATIAL_LEGENDRE,
SURFACE,

52
include/openmc/tallies/filter_reaction.h Normal file
View file

@ -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<int>& bins() const { return bins_; }
void set_bins(span<const int> bins);
protected:
//----------------------------------------------------------------------------
// Data members
//! MT numbers to match
vector<int> bins_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_REACTION_H

12
openmc/data/reaction.py
View file

@ -21,6 +21,7 @@ from .function import Tabulated1D, Polynomial
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 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 3: "(n,nonelastic)",
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,t{i - 700})' for i in range(700, 749)})
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)

83
openmc/filter.py
View file

@ -13,6 +13,7 @@ import pandas as pd
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 @@ from ._xml import get_elem_list, get_text
_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 @@ class CollisionFilter(Filter):
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

17
openmc/lib/filter.py
View file

@ -22,9 +22,10 @@ __all__ = [
'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 @@ class ParticleFilter(Filter):
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 @@ _FILTER_TYPE_MAP = {
'musurface': MuSurfaceFilter,
'parentnuclide': ParentNuclideFilter,
'particle': ParticleFilter,
'particleproduction': ParticleProductionFilter,
'polar': PolarFilter,
'reaction': ReactionFilter,
'sphericalharmonics': SphericalHarmonicsFilter,
'spatiallegendre': SpatialLegendreFilter,
'surface': SurfaceFilter,

2
src/chain.cpp
View file

@ -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});
}

77
src/endf.cpp
View file

@ -87,6 +87,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<Function1D> read_function(hid_t group, const char* name)
{
hid_t obj_id = open_object(group, name);

78
src/reaction.cpp
View file

@ -309,6 +309,7 @@ std::unordered_map<int, std::string> 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;
int mt = it->second;
return mt;
}
// 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;
}
// 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;
}

2
src/settings.cpp
View file

@ -964,7 +964,7 @@ void read_settings_xml(pugi::xml_node root)
if (check_for_node(node_ct, "reactions")) {
auto temp = get_node_array<std::string>(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);
}

3
src/tallies/filter.cpp
View file

@ -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<ParticleProductionFilter>(id);
} else if (type == "polar") {
return Filter::create<PolarFilter>(id);
} else if (type == "reaction") {
return Filter::create<ReactionFilter>(id);
} else if (type == "surface") {
return Filter::create<SurfaceFilter>(id);
} else if (type == "spatiallegendre") {

79
src/tallies/filter_reaction.cpp Normal file
View file

@ -0,0 +1,79 @@
#include "openmc/tallies/filter_reaction.h"
#include <fmt/core.h>
#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<std::string>(node, "bins");
// Convert reaction names to MT numbers
vector<int> 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<const int> 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<std::string> 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

8
src/tallies/tally.cpp
View file

@ -171,6 +171,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;
}
}
}
@ -569,7 +573,7 @@ void Tally::set_scores(const vector<std::string>& scores)
}
// Determine integer code for score
int score = reaction_type(score_str);
int score = reaction_tally_mt(score_str);
switch (score) {
case SCORE_FLUX:
@ -785,7 +789,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()) {

0
tests/regression_tests/filter_reaction/__init__.py Normal file
View file

27
tests/regression_tests/filter_reaction/inputs_true.dat Normal file
View file

@ -0,0 +1,27 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" depletable="true">
<density value="10.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1"/>
<surface id="1" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 100.0"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>5</batches>
</settings>
<tallies>
<filter id="1" type="reaction">
<bins>(n,elastic) (n,2n) (n,fission) (n,gamma) (n,total)</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>flux</scores>
</tally>
</tallies>
</model>

13
tests/regression_tests/filter_reaction/results_true.dat Normal file
View file

@ -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

31
tests/regression_tests/filter_reaction/test.py Normal file
View file

@ -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()

87
tests/unit_tests/test_filter_reaction.py Normal file
View file

@ -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)'