Fixed a bug in charged particle energy deposition. (#3416)

Co-authored-by: Jonathan Shimwell <drshimwell@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

docs/source/methods/energy_deposition.rst

@@ -25,19 +25,36 @@ KERMA (Kinetic Energy Release in Materials) [Mack97]_ coefficients for reaction
 :math:`\times` cross-section (e.g., eV-barn) and can be used much like a reaction
 cross section for the purpose of tallying energy deposition.
 
-KERMA coefficients can be computed using the energy-balance method with
-a nuclear data processing code like NJOY, which performs the following
-iteration over all reactions :math:`r` for all isotopes :math:`i`
-requested
+KERMA coefficients can be computed using the energy-balance method with a
+nuclear data processing code like NJOY, which estimates the KERMA coefficients
+using the following equation:
 
 .. math::
 
-    k_{i, r}(E) = \left(E + Q_{i, r} - \bar{E}_{i, r, n}
+    k_{i, r}(E) = \left(E + Q_{i, r} - \sum\limits_x \bar{E}_{i, r, x}
+    \right)\sigma_{i, r}(E),
+
+where the summation is over each secondary particle type :math:`x`. This
+equation states that the energy deposited is equal to the energy of the incident
+particle plus the reaction :math:`Q` value less the energy of secondary
+particles that are transported away from the reaction site. For neutron
+interactions, the energy-balance KERMA coefficient is
+
+.. math::
+
+    k_{i, r}(E) = \left(E + Q_{i, r} - \sum\limits_x \bar{E}_{i, r, n}
     - \bar{E}_{i, r, \gamma}\right)\sigma_{i, r}(E),
 
-removing the energy of neutral particles (neutrons and photons) that are
-transported away from the reaction site :math:`\bar{E}`, and the reaction
-:math:`Q` value.
+where :math:`\bar{E}_{i, r, n}` is the average energy of secondary neutrons and
+:math:`\bar{E}_{i, r, \gamma}` is the average energy of secondary photons. For
+photon and charged particle interactions, the :math:`Q` value is zero and thus
+the KERMA coefficient is
+
+.. math::
+    :label: energy-balance-photon
+
+    k_{i, r}(E) = \left(E - \sum\limits_x \bar{E}_{i, r, x}
+    \right)\sigma_{i, r}(E).
 
 -------
 Fission
@@ -120,7 +137,7 @@ run with :math:`N918` reflecting fission heating computed from NJOY.
 This modified heating data is stored as the MT=901 reaction and will be scored
 if ``heating-local`` is included in :attr:`openmc.Tally.scores`.
 
-Coupled neutron-photon transport
+Coupled Neutron-Photon Transport
 --------------------------------
 
 Here, OpenMC instructs ``heatr`` to assume that energy from photons is not
@@ -138,6 +155,50 @@ Let :math:`N301` represent the total heating number returned from this
 This modified heating data is stored as the MT=301 reaction and will be scored
 if ``heating`` is included in :attr:`openmc.Tally.scores`.
 
+Photons and Charged Particles
+-----------------------------
+
+In OpenMC, energy deposition from photons or charged particles is scored using
+the energy balance method based on Equation :eq:`energy-balance-photon`. Special
+consideration is given to electrons and positrons as described below.
+
++++++++++++++++++
+Charged Particles
++++++++++++++++++
+
+OpenMC tracks photons interaction by interaction so the energy deposited in each
+collision is easily attributed back to the nuclide and reaction for which the
+photon interacted with. Charged particles (electrons and photons) aren't tracked
+in the same way. For charged particles, OpenMC assumes that all their energy
+(less the energy of bremsstrahlung radiation) is deposited in the material in
+which they were born. In this way it is harder to trace how much energy should
+be attributed in each nuclide.
+
+According to the CSDA approximation (see :ref:`ttb`) the energy deposited by a
+charged particle with kinetic energy :math:`T` in the :math:`i`-th element can
+be calculated as:
+
+.. math::
+
+    E_{i} = \int_{0}^{R(T)} w_{i}S_{\text{col,i}} dx
+
+where :math:`R(T)` is the CSDA range of the charged particle,
+:math:`S_{\text{col},i}` is the collision stopping power of the charged particle
+in the :math:`i`-th element and :math:`w_i` is the mass fraction of the
+:math:`i`-th element. According to the Bethe formula the collision stopping
+power of the :math:`i`-th element is proportional to :math:`Z_i/A_i`, so the
+fractional collision stopping power from the :math:`i`-th element is:
+
+.. math::
+
+    \frac{w_{i}S_{\text{col},i}(T)}{S_{\text{col}}(T)} =
+    \frac{\frac{w_{i}Z_{i}}{A_{i}}}{\sum_{i}\frac{w_{i}Z_{i}}{A_{i}}} =
+    \frac{\gamma_i Z_{i}}{\sum_{i}\gamma_i Z_{i}}.
+
+where :math:`\gamma_i` is the atomic fraction of the :math:`i`-th element.
+Therefore, the energy deposited by charged particles should be attributed to
+a given element according to its fractional charge density.
+
 ----------
 References
 ----------

include/openmc/material.h

@@ -107,6 +107,10 @@ public:
   //! \return Density in [g/cm^3]
   double density_gpcc() const { return density_gpcc_; }
 
+  //! Get charge density in [e/b-cm]
+  //! \return Charge density in [e/b-cm]
+  double charge_density() const { return charge_density_; };
+
   //! Get name
   //! \return Material name
   const std::string& name() const { return name_; }
@@ -177,6 +181,7 @@ public:
   xt::xtensor<double, 1> atom_density_; //!< Nuclide atom density in [atom/b-cm]
   double density_;                      //!< Total atom density in [atom/b-cm]
   double density_gpcc_;                 //!< Total atom density in [g/cm^3]
+  double charge_density_;               //!< Total charge density in [e/b-cm]
   double volume_ {-1.0};                //!< Volume in [cm^3]
   vector<bool> p0_; //!< Indicate which nuclides are to be treated with
                     //!< iso-in-lab scattering

src/material.cpp

@@ -454,12 +454,15 @@ void Material::normalize_density()
   // Calculate nuclide atom densities
   atom_density_ *= density_;
 
-  // Calculate density in g/cm^3.
+  // Calculate density in [g/cm^3] and charge density in [e/b-cm]
   density_gpcc_ = 0.0;
+  charge_density_ = 0.0;
   for (int i = 0; i < nuclide_.size(); ++i) {
     int i_nuc = nuclide_[i];
     double awr = settings::run_CE ? data::nuclides[i_nuc]->awr_ : 1.0;
+    int z = settings::run_CE ? data::nuclides[i_nuc]->Z_ : 0.0;
     density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO;
+    charge_density_ += atom_density_(i) * z;
   }
 }
 
@@ -982,12 +985,15 @@ void Material::set_density(double density, const std::string& units)
     // Recalculate nuclide atom densities based on given density
     atom_density_ *= density;
 
-    // Calculate density in g/cm^3.
+    // Calculate density in g/cm^3 and charge density in [e/b-cm]
     density_gpcc_ = 0.0;
+    charge_density_ = 0.0;
     for (int i = 0; i < nuclide_.size(); ++i) {
       int i_nuc = nuclide_[i];
       double awr = data::nuclides[i_nuc]->awr_;
+      int z = settings::run_CE ? data::nuclides[i_nuc]->Z_ : 0.0;
       density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO;
+      charge_density_ += atom_density_(i) * z;
     }
   } else if (units == "g/cm3" || units == "g/cc") {
     // Determine factor by which to change densities
@@ -998,6 +1004,7 @@ void Material::set_density(double density, const std::string& units)
     density_gpcc_ = density;
     density_ *= f;
     atom_density_ *= f;
+    charge_density_ *= f;
   } else {
     throw std::invalid_argument {
       "Invalid units '" + std::string(units.data()) + "' specified."};

src/tallies/tally_scoring.cpp

@@ -325,6 +325,39 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux,
   return score;
 }
 
+//! Helper function to obtain particle heating [eV]
+
+double score_particle_heating(const Particle& p, const Tally& tally,
+  double flux, int rxn_bin, int i_nuclide, double atom_density)
+{
+  if (p.type() == ParticleType::neutron)
+    return score_neutron_heating(
+      p, tally, flux, rxn_bin, i_nuclide, atom_density);
+  if (i_nuclide == -1 || i_nuclide == p.event_nuclide() ||
+      p.event_nuclide() == -1) {
+    // Get the pre-collision energy of the particle.
+    auto E = p.E_last();
+    // The energy deposited is the difference between the pre-collision
+    // and post-collision energy...
+    double score = E - p.E();
+    // ...less the energy of any secondary particles since they will be
+    // transported individually later
+    score -= p.bank_second_E();
+    score *= p.wgt_last();
+
+    // if no event_nuclide (charged particle) scale energy deposition by
+    // fractional charge density
+    if (i_nuclide != -1 && p.event_nuclide() == -1) {
+      const auto& mat {model::materials[p.material()]};
+      int z = data::nuclides[i_nuclide]->Z_;
+      auto i = mat->mat_nuclide_index_[i_nuclide];
+      score *= (z * mat->atom_density_[i] / mat->charge_density());
+    }
+    return score;
+  }
+  return 0.0;
+}
+
 //! Helper function for nu-fission tallies with energyout filters.
 //
 //! In this case, we may need to score to multiple bins if there were multiple
@@ -1007,23 +1040,8 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
       break;
 
     case HEATING:
-      if (p.type() == Type::neutron) {
-        score = score_neutron_heating(
-          p, tally, flux, HEATING, i_nuclide, atom_density);
-      } else {
-        if (i_nuclide == -1 || i_nuclide == p.event_nuclide()) {
-          // The energy deposited is the difference between the pre-collision
-          // and post-collision energy...
-          score = E - p.E();
-          // ...less the energy of any secondary particles since they will be
-          // transported individually later
-          score -= p.bank_second_E();
-
-          score *= p.wgt_last();
-        } else {
-          score = 0.0;
-        }
-      }
+      score = score_particle_heating(
+        p, tally, flux, HEATING, i_nuclide, atom_density);
       break;
 
     default:
@@ -1539,19 +1557,8 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
       break;
 
     case HEATING:
-      if (p.type() == Type::neutron) {
-        score = score_neutron_heating(
-          p, tally, flux, HEATING, i_nuclide, atom_density);
-      } else {
-        // The energy deposited is the difference between the pre-collision and
-        // post-collision energy...
-        score = E - p.E();
-        // ...less the energy of any secondary particles since they will be
-        // transported individually later
-        score -= p.bank_second_E();
-
-        score *= p.wgt_last();
-      }
+      score = score_particle_heating(
+        p, tally, flux, HEATING, i_nuclide, atom_density);
       break;
 
     default:

tests/regression_tests/photon_production/results_true.dat

@@ -67,8 +67,8 @@ tally 3:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-0.000000E+00
-0.000000E+00
+1.774484E+05
+3.148794E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -79,8 +79,8 @@ tally 3:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-0.000000E+00
-0.000000E+00
+7.692488E+03
+5.917437E+07
 0.000000E+00
 0.000000E+00
 0.000000E+00

tests/unit_tests/test_data_multipole.py

@@ -52,7 +52,7 @@ def test_from_endf():
     pytest.importorskip('vectfit')
     endf_data = os.environ['OPENMC_ENDF_DATA']
     endf_file = os.path.join(endf_data, 'neutrons', 'n-001_H_001.endf')
-    return openmc.data.WindowedMultipole.from_endf(
+    assert openmc.data.WindowedMultipole.from_endf(
             endf_file, log=True, wmp_options={"n_win": 400, "n_cf": 3})
 
 
@@ -60,5 +60,5 @@ def test_from_endf_search():
     pytest.importorskip('vectfit')
     endf_data = os.environ['OPENMC_ENDF_DATA']
     endf_file = os.path.join(endf_data, 'neutrons', 'n-095_Am_244.endf')
-    return openmc.data.WindowedMultipole.from_endf(
+    assert openmc.data.WindowedMultipole.from_endf(
             endf_file, log=True, wmp_options={"search": True, 'rtol':1e-2})

tests/unit_tests/test_nuclide_heating.py

@@ -0,0 +1,39 @@
+import openmc
+from pytest import approx
+
+
+def test_nuclide_heating(run_in_tmpdir):
+    mat = openmc.Material()
+    mat.add_nuclide("Li6", 0.5)
+    mat.add_nuclide("Li7", 0.5)
+    mat.set_density("g/cm3", 1.0)
+
+    sphere = openmc.Sphere(r=20, boundary_type="reflective")
+    inside_sphere = openmc.Cell(fill=mat, region=-sphere)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([inside_sphere])
+
+    model.settings.particles = 1000
+    model.settings.batches = 1
+    model.settings.photon_transport = True
+    model.settings.electron_treatment = "ttb"
+    model.settings.cutoff = {"energy_photon": 1000}
+    model.settings.run_mode = "fixed source"
+    model.settings.source = openmc.IndependentSource(
+        energy=openmc.stats.delta_function(10.0e6),
+        particle="photon"
+    )
+
+    # Create two tallies, one with heating by nuclide and one with total heating
+    tally1 = openmc.Tally()
+    tally1.scores = ["heating"]
+    tally1.nuclides = mat.get_nuclides()
+    tally2 = openmc.Tally()
+    tally2.scores = ["heating"]
+    model.tallies = [tally1, tally2]
+
+    # Run the model
+    model.run(apply_tally_results=True)
+
+    # Make sure the heating results are consistent
+    assert tally1.mean.sum() == approx(tally2.mean.sum())