Fixed a bug in charged particle energy deposition.

docs/source/methods/energy_deposition.rst

@@ -81,6 +81,34 @@ NJOY computes the fission KERMA coefficient using this energy-balance method to
     The energy from delayed neutrons and photons and beta particles is intentionally
     left out from the NJOY calculations.
 
+----------------------------------
+Photon and Charged Particles Kerma
+----------------------------------
+
+In a similar way to neutron heating, the heating rate from photons or charged particles is defined as:
+
+.. math::
+
+    H(E) = \phi(E)\sum_i\rho_i\sum_rk_{i, r}(E),
+
+and has units energy per time, typically eV/s. Here, :math:`k_{i, r}` are the
+KERMA (Kinetic Energy Release in Materials) coefficients for reaction
+:math:`r` of isotope :math:`i`. The KERMA coefficients have units of energy
+: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
+
+.. math::
+
+    k_{i, r}(E) = \left(E - \bar{E}_{i, r, \gamma} - \bar{E}_{i, r, e}\right)\sigma_{i, r}(E),
+
+removing the energy of neutral and charged particles (photons and electrons/positrons) that are
+transported away from the reaction site.
+
+.. note::
+In the KERMA coefficients for photons and charged particles the reaction Q-value is zero.
+
 ---------------------
 OpenMC Implementation
 ---------------------
@@ -138,6 +166,39 @@ 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 Energy Deposition
+-----------------------------------------------
+
+In OpenMC, energy deposition from photons or charged particles is scored in the following way:
+After every collision, the kinetic energy of the incident particle is compared before and after the collision and this energy  removing energy of any secondary products particles is scored as deposited. This algorithm is justified by energy balance and the fact that photons and charged particles reaction Q-value is always zero. 
+
++++++++++++++++++++++++++++++++++++
+Charged Particles Energy Deposition
++++++++++++++++++++++++++++++++++++
+
+OpenMC tracks photons interaction by interaction so the energy deposited in each collision is easily traced back to the nuclide and reaction for which the photon interacted with.
+Charged particles aren't tracked in the same way. For charged particles OpenMC assume that all their energy (less 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 deposited in each nuclide.
+
+According to the CSDA approximation the energy deposited by a charged particle with kinetic energy 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:`\frac{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.
+
+So the energy deposited by charged particles should be divided to elements according to the fractional charge density.
+
 ----------
 References
 ----------

include/openmc/material.h

@@ -56,6 +56,8 @@ class Material {
 
   void calculate_xs(Particle& p) const;
 
+  double charge_density();
+
   //! Assign thermal scattering tables to specific nuclides within the material
   //! so the code knows when to apply bound thermal scattering data
   void init_thermal();

src/material.cpp

@@ -933,6 +933,20 @@ void Material::calculate_photon_xs(Particle& p) const
   }
 }
 
+double Material::charge_density()
+{
+  double val = 0.0;
+  // Add contribution from each nuclide in material
+  for (int i = 0; i < nuclide_.size(); ++i) {
+
+    // Get nuclide index
+    int i_nuclide = nuclide_[i];
+    int z = data::nuclides[i_nuclide]->Z_;
+    val += atom_density_(i) * z;
+  }
+  return val;
+}
+
 void Material::set_id(int32_t id)
 {
   assert(id >= 0 || id == C_NONE);

src/tallies/tally_scoring.cpp

@@ -325,6 +325,44 @@ 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 n = mat->nuclide_.size();
+      int z = data::nuclides[i_nuclide]->Z_;
+      for (int i = 0; i < n; ++i) {
+        if (mat->nuclide_[i] == i_nuclide) {
+          score *= (z * mat->atom_density_[i] / mat->charge_density());
+          break;
+        }
+      }
+    }
+    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 +1045,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 +1562,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/nuclide_heating/inputs_true.dat

@@ -0,0 +1,42 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="0.5" name="Cu63"/>
+      <nuclide ao="0.5" name="Cu65"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <surface boundary="reflective" coeffs="0.0 0.0 0.0 20" id="1" type="sphere"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>10000</particles>
+    <batches>1</batches>
+    <source particle="photon" strength="1.0" type="independent">
+      <space type="point">
+        <parameters>0 0 0</parameters>
+      </space>
+      <angle type="isotropic"/>
+      <energy type="discrete">
+        <parameters>10000000.0 1.0</parameters>
+      </energy>
+    </source>
+    <electron_treatment>ttb</electron_treatment>
+    <photon_transport>true</photon_transport>
+    <cutoff>
+      <energy_photon>1000</energy_photon>
+    </cutoff>
+  </settings>
+  <tallies>
+    <tally id="1">
+      <nuclides>Cu63 Cu65</nuclides>
+      <scores>heating</scores>
+    </tally>
+    <tally id="2">
+      <scores>heating</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/nuclide_heating/results_true.dat

@@ -0,0 +1,8 @@
+tally 1:
+5.004050E+06
+2.504051E+13
+4.995950E+06
+2.495952E+13
+tally 2:
+1.000000E+07
+1.000000E+14

tests/regression_tests/nuclide_heating/test.py

@@ -0,0 +1,53 @@
+from math import pi
+
+import numpy as np
+import openmc
+
+from tests.testing_harness import PyAPITestHarness
+
+
+class NuclideHeatingTestHarness(PyAPITestHarness):
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        mat = openmc.Material()
+        mat.add_nuclide("Cu63", 0.5, "ao")
+        mat.add_nuclide("Cu65", 0.5, "ao")
+        mat.set_density("g/cm3", 1.0)
+        self._model.materials = openmc.Materials([mat])
+
+        sphere = openmc.Sphere(r=20, boundary_type="reflective")
+        inside_sphere = openmc.Cell()
+        inside_sphere.region = -sphere
+        inside_sphere.fill = mat
+        self._model.geometry = openmc.Geometry([inside_sphere])
+
+        source = openmc.IndependentSource()
+        source.space = openmc.stats.Point((0, 0, 0))
+        source.angle = openmc.stats.Isotropic()
+        source.energy = openmc.stats.Discrete([10.0e6], [1.0])
+        source.particle = "photon"
+
+        settings = openmc.Settings()
+        settings.particles = 10000
+        settings.batches = 1
+        settings.photon_transport = True
+        settings.electron_treatment = "ttb"
+        settings.cutoff = {"energy_photon": 1000}
+        settings.run_mode = "fixed source"
+        settings.source = source
+        self._model.settings = settings
+
+        tally1 = openmc.Tally()
+        tally1.scores = ["heating"]
+        tally1.nuclides = ["Cu63", "Cu65"]
+
+        tally2 = openmc.Tally()
+        tally2.scores = ["heating"]
+
+        tallies = openmc.Tallies([tally1, tally2])
+        self._model.tallies = tallies
+
+
+def test_nuclide_heating():
+    harness = NuclideHeatingTestHarness("statepoint.1.h5", model=openmc.Model())
+    harness.main()

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