diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 8174108e1ae..ade764bb348 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -584,9 +584,9 @@ Level Inelastic :Object type: Group :Attributes: - **type** (*char[]*) -- 'level' - - **threshold** (*double*) -- Energy threshold in the laboratory - system in eV - - **mass_ratio** (*double*) -- :math:`(A/(A + 1))^2` + - **q_value** (*double*) -- Q value in eV + - **mass** (*double*) -- Nucleus mass A relative to neutron rest mass + - **particle** (*char[]*) -- Incident particle name Continuous Tabular ------------------ diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst index 2b797e3dbcf..aa45f51ac39 100644 --- a/docs/source/methods/neutron_physics.rst +++ b/docs/source/methods/neutron_physics.rst @@ -568,7 +568,7 @@ and the incoming energy: .. math:: :label: level-scattering - E' = \left ( \frac{A}{A+1} \right )^2 \left ( E - \frac{A + 1}{A} Q \right ) + E' = \left ( \frac{A}{A+1} \right )^2 \left ( E - \frac{A + 1}{A} |Q| \right ) where :math:`A` is the mass of the target nucleus measured in neutron masses. diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index d2bd3ac7609..5eafb439f07 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -16,9 +16,9 @@ de-excitation of these atoms can result in the emission of electrons and photons. Electrons themselves also can produce photons by means of bremsstrahlung radiation. -------------------- -Photon Interactions -------------------- +------------------------ +Photoatomic Interactions +------------------------ Coherent (Rayleigh) Scattering ------------------------------ @@ -598,6 +598,29 @@ The inverse transform method is used to sample :math:`\mu_{-}` and The azimuthal angles for the electron and positron are sampled independently and uniformly on :math:`[0, 2\pi)`. +------------------------- +Photonuclear Interactions +------------------------- + +Photonuclear reactions occur when a nucleus absorbs a high-energy photon (gamma-ray), leading to a change in the nucleus. +The absorption of the photon excites the nucleus, and the energy provided by the photon can result in the emission of particles +like neutrons, protons, or other light nuclei. Because photonuclear interactions are determined by nuclear rather than atomic properties, +their data libraries must be specific to each isotope. + +Inelastic Level Scattering +-------------------------- + +It can be shown (see Wattenberg_) that in inelastic level scattering, the outgoing +energy of the neutron :math:`E'` can be related to the Q-value of the reaction +and the incoming energy of the photon: + +.. math:: + :label: photonuclear-level-scattering + + E' = \left ( \frac{A-1}{A} \right ) \left ( E - |Q| - \frac{E^2}{2 m_n \left (A-1 \right )} \right ) + +where :math:`A` is the mass of the target nucleus measured in neutron masses, math:`m_n` is the neutron mass in eV. + ------------------- Secondary Processes ------------------- @@ -734,3 +757,5 @@ emitted photon. .. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf .. _Salvat: https://doi.org/10.1787/32da5043-en + +.. _Wattenberg: https://doi.org/10.1103/PhysRev.71.497 diff --git a/include/openmc/distribution_energy.h b/include/openmc/distribution_energy.h index 9b08ed039dc..1ca9b905b35 100644 --- a/include/openmc/distribution_energy.h +++ b/include/openmc/distribution_energy.h @@ -61,8 +61,10 @@ class LevelInelastic : public EnergyDistribution { double sample(double E, uint64_t* seed) const override; private: - double threshold_; //!< Energy threshold in lab, (A + 1)/A * |Q| - double mass_ratio_; //!< (A/(A+1))^2 + //! The scattering law is E_out = a*(E_in-b-c*E_in*E_in) + double a_; //!< a coefficient of the scattering law + double b_; //!< b coefficient of the scattering law + double c_; //!< c coefficient of the scattering law }; //=============================================================================== diff --git a/openmc/data/data.py b/openmc/data/data.py index 5ecadd37be0..e075ebbe145 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -274,6 +274,7 @@ # Unit conversions EV_PER_MEV = 1.0e6 JOULE_PER_EV = 1.602176634e-19 +EV_PER_AMU = 931.49410372e6 # eV/c^2 per amu # Avogadro's constant AVOGADRO = 6.02214076e23 @@ -281,6 +282,9 @@ # Neutron mass in units of amu NEUTRON_MASS = 1.00866491595 +# Neutron mass in units of eV/c^2 +NEUTRON_MASS_EV = NEUTRON_MASS*EV_PER_AMU + # Used in atomic_mass function as a cache _ATOMIC_MASS: dict[str, float] = {} diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py index 069ab1b9beb..ede97e98ae9 100644 --- a/openmc/data/energy_distribution.py +++ b/openmc/data/energy_distribution.py @@ -1,5 +1,6 @@ from abc import ABC, abstractmethod from collections.abc import Iterable +from math import sqrt, isclose from numbers import Integral, Real from warnings import warn @@ -8,7 +9,7 @@ import openmc.checkvalue as cv from openmc.mixin import EqualityMixin from openmc.stats.univariate import Univariate, Tabular, Discrete, Mixture -from .data import EV_PER_MEV +from .data import EV_PER_MEV, NEUTRON_MASS_EV from .endf import get_tab1_record, get_tab2_record from .function import Tabulated1D, INTERPOLATION_SCHEME @@ -897,42 +898,67 @@ class LevelInelastic(EnergyDistribution): Parameters ---------- - threshold : float - Energy threshold in the laboratory system, :math:`(A + 1)/A * |Q|` - mass_ratio : float - :math:`(A/(A + 1))^2` + q_value : float + Q-value of the reaction. + mass : float + mass of the nucleus in units of neutron mass. + particle : {'neutron', 'photon'} + incident particle type, defaults to neutron Attributes ---------- + q_value : float + Q-value of the reaction. + mass : float + mass of the nucleus in units of neutron mass. + particle : {'neutron', 'photon'} + incident particle type. threshold : float - Energy threshold in the laboratory system, :math:`(A + 1)/A * |Q|` - mass_ratio : float - :math:`(A/(A + 1))^2` - + Energy threshold in the laboratory system """ - def __init__(self, threshold, mass_ratio): + def __init__(self, q_value, mass, particle = 'neutron'): super().__init__() - self.threshold = threshold - self.mass_ratio = mass_ratio + self.q_value = q_value + self.mass = mass + self.particle = particle @property - def threshold(self): - return self._threshold + def q_value(self): + return self._q_value - @threshold.setter - def threshold(self, threshold): - cv.check_type('level inelastic threhsold', threshold, Real) - self._threshold = threshold + @q_value.setter + def q_value(self, q_value): + cv.check_type('level inelastic q_value', q_value, Real) + self._q_value = q_value @property - def mass_ratio(self): - return self._mass_ratio - - @mass_ratio.setter - def mass_ratio(self, mass_ratio): - cv.check_type('level inelastic mass ratio', mass_ratio, Real) - self._mass_ratio = mass_ratio + def mass(self): + return self._mass + + @mass.setter + def mass(self, mass): + cv.check_type('level inelastic mass', mass, Real) + self._mass = mass + + @property + def particle(self): + return self._particle + + @particle.setter + def particle(self, particle): + cv.check_value('product particle type', particle, ['neutron', 'photon']) + self._particle = particle + + @property + def threshold(self): + A = self.mass + Q = self.q_value + if particle == 'neutron': + return (A+1.0)/A*abs(Q) + else: + b = NEUTRON_MASS_EV*(A-1) + return sqrt(b)*(sqrt(b)-sqrt(b-2*abs(Q))) def to_hdf5(self, group): """Write distribution to an HDF5 group @@ -945,8 +971,9 @@ def to_hdf5(self, group): """ group.attrs['type'] = np.bytes_('level') - group.attrs['threshold'] = self.threshold - group.attrs['mass_ratio'] = self.mass_ratio + group.attrs['q_value'] = self.q_value + group.attrs['mass'] = self.mass + group.attrs['particle'] = np.bytes_(self.particle) @classmethod def from_hdf5(cls, group): @@ -963,9 +990,18 @@ def from_hdf5(cls, group): Level inelastic scattering distribution """ - threshold = group.attrs['threshold'] - mass_ratio = group.attrs['mass_ratio'] - return cls(threshold, mass_ratio) + #backwards compatible read: + if 'threshold' in group.attrs: + threshold = group.attrs['threshold'] + mass_ratio = group.attrs['mass_ratio'] + mass = 1.0/(1.0/sqrt(mass_ratio)-1.0) + q_value = -threshold * sqrt(mass_ratio) + return cls(q_value, mass) + + q_value = group.attrs['q_value'] + mass = group.attrs['mass'] + particle = group.attrs['particle'].decode() + return cls(q_value, mass, particle) @classmethod def from_ace(cls, ace, idx): @@ -984,9 +1020,15 @@ def from_ace(cls, ace, idx): Level inelastic scattering distribution """ + particle = {'u':'photon', 'c': 'neutron'}[ace.data_type.value] threshold = ace.xss[idx]*EV_PER_MEV mass_ratio = ace.xss[idx + 1] - return cls(threshold, mass_ratio) + mass = 1.0/(1.0/sqrt(mass_ratio)-1.0) if particle == 'neutron' else 1.0-1.0/mass_ratio + if not isclose(ace.atomic_weight_ratio, mass): + warn("Level inelastic distribution mass parameter does not match ace table mass.") + mass = ace.atomic_weight_ratio + q_value = -threshold * sqrt(mass_ratio) if particle == 'neutron' else -threshold + return cls(q_value, mass, particle = particle) class ContinuousTabular(EnergyDistribution): diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 65b59582cf8..4160bf71add 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -1200,12 +1200,7 @@ def from_endf(cls, ev, mt): # since it can be calculated analytically. Here we determine the # necessary parameters to create a LevelInelastic object dist = UncorrelatedAngleEnergy() - - A = ev.target['mass'] - threshold = (A + 1.)/A*abs(rx.q_value) - mass_ratio = (A/(A + 1.))**2 - dist.energy = LevelInelastic(threshold, mass_ratio) - + dist.energy = LevelInelastic(rx.q_value, ev.target['mass']) neutron.distribution.append(dist) if (4, mt) in ev.section: diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index a4a5ce9e1b6..8450826447c 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -1,14 +1,17 @@ #include "openmc/distribution_energy.h" #include // for max, min, copy, move +#include // for sqrt, abs #include // for size_t #include // for back_inserter #include "xtensor/xview.hpp" +#include "openmc/constants.h" #include "openmc/endf.h" #include "openmc/hdf5_interface.h" #include "openmc/math_functions.h" +#include "openmc/particle.h" #include "openmc/random_dist.h" #include "openmc/random_lcg.h" #include "openmc/search.h" @@ -41,13 +44,35 @@ double DiscretePhoton::sample(double E, uint64_t* seed) const LevelInelastic::LevelInelastic(hid_t group) { - read_attribute(group, "threshold", threshold_); - read_attribute(group, "mass_ratio", mass_ratio_); + // for backwards compatibility: + if (attribute_exists(group, "mass_ratio")) { + read_attribute(group, "threshold", b_); + read_attribute(group, "mass_ratio", a_); + c_ = 0.0; + } else { + double A, Q; + std::string temp; + read_attribute(group, "mass", A); + read_attribute(group, "q_value", Q); + read_attribute(group, "particle", temp); + auto particle = str_to_particle_type(temp); + if (particle == ParticleType::neutron) { + a_ = (A / (A + 1.0)) * (A / (A + 1.0)); + b_ = (A + 1.0) / A * std::abs(Q); + c_ = 0.0; + } else if (particle == ParticleType::photon) { + a_ = (A - 1.0) / A; + b_ = std::abs(Q); + c_ = 1.0 / (2.0 * MASS_NEUTRON_EV * (A - 1.0)); + } else { + fatal_error("Unrecognized particle: " + temp); + } + } } double LevelInelastic::sample(double E, uint64_t* seed) const { - return mass_ratio_ * (E - threshold_); + return a_ * (E - b_ - c_ * (E * E)); } //==============================================================================