Extend level scattering to support incident photons

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

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.
 

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

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
 };
 
 //===============================================================================

openmc/data/data.py

@@ -274,13 +274,17 @@
 # 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
 
 # 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] = {}
 

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):

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:

src/distribution_energy.cpp

@@ -1,14 +1,17 @@
 #include "openmc/distribution_energy.h"
 
 #include <algorithm> // for max, min, copy, move
+#include <cmath>     // for sqrt, abs
 #include <cstddef>   // for size_t
 #include <iterator>  // 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));
 }
 
 //==============================================================================