Avoid negative heating values during pair production and bremsstrahlung (#3426)

Co-authored-by: GuySten <guyste@post.bgu.ac.il>

Author: paulromano

docs/source/methods/energy_deposition.rst

@@ -47,15 +47,17 @@ interactions, the energy-balance KERMA coefficient is
 
 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
+photon and charged particle interactions the KERMA coefficient is
 
 .. math::
     :label: energy-balance-photon
 
-    k_{i, r}(E) = \left(E - \sum\limits_x \bar{E}_{i, r, x}
+    k_{i, r}(E) = \left(E + Q_{i, r} - \sum\limits_x \bar{E}_{i, r, x}
     \right)\sigma_{i, r}(E).
 
+where the :math:`Q` value is zero for all interactions except for pair
+production and positron annihilation.
+
 -------
 Fission
 -------

include/openmc/particle_data.h

@@ -573,7 +573,8 @@ class ParticleData : public GeometryState {
   bool& fission() { return fission_; }            // true if implicit fission
   int& event_nuclide() { return event_nuclide_; } // index of collision nuclide
   const int& event_nuclide() const { return event_nuclide_; }
-  int& event_mt() { return event_mt_; }           // MT number of collision
+  int& event_mt() { return event_mt_; } // MT number of collision
+  const int& event_mt() const { return event_mt_; }
   int& delayed_group() { return delayed_group_; } // delayed group
   const int& parent_nuclide() const { return parent_nuclide_; }
   int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide

src/bremsstrahlung.cpp

@@ -112,6 +112,12 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
                std::pow(a * (c - c_l) / (std::exp(w_l) * p_l) + 1.0, 1.0 / a);
 
     if (w > settings::energy_cutoff[photon]) {
+      // If the energy of the secondary photon is larger than the remaining
+      // energy of the primary particle, adjust it to the remaining energy
+      if (*E_lost + w > p.E()) {
+        w = p.E() - *E_lost;
+      }
+
       // Create secondary photon
       p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon);
       *E_lost += w;

src/tallies/tally_scoring.cpp

@@ -325,6 +325,20 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux,
   return score;
 }
 
+//! Helper function to obtain reaction Q value for photons and charged particles
+double get_reaction_q_value(const Particle& p)
+{
+  if (p.type() == ParticleType::photon && p.event_mt() == PAIR_PROD) {
+    // pair production
+    return -2 * MASS_ELECTRON_EV;
+  } else if (p.type() == ParticleType::positron) {
+    // positron annihilation
+    return 2 * MASS_ELECTRON_EV;
+  } else {
+    return 0.0;
+  }
+}
+
 //! Helper function to obtain particle heating [eV]
 
 double score_particle_heating(const Particle& p, const Tally& tally,
@@ -335,14 +349,17 @@ double score_particle_heating(const Particle& p, const Tally& tally,
       p, tally, flux, rxn_bin, i_nuclide, atom_density);
   if (i_nuclide == -1 || i_nuclide == p.event_nuclide() ||
       p.event_nuclide() == -1) {
+    // For pair production and positron annihilation, we need to account for the
+    // reaction Q value
+    double Q = get_reaction_q_value(p);
+
     // 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();
+
+    // The energy deposited is the sum of the incident energy and the reaction
+    // Q-value less the energy of any outgoing particles
+    double score = E + Q - p.E() - p.bank_second_E();
+
     score *= p.wgt_last();
 
     // if no event_nuclide (charged particle) scale energy deposition by

tests/regression_tests/photon_production/results_true.dat

@@ -55,14 +55,14 @@ tally 3:
 1.846062E-07
 2.297000E-01
 5.276209E-02
-7.055982E+03
-4.978688E+07
+4.196651E+00
+1.761188E+01
 0.000000E+00
 0.000000E+00
 2.297000E-01
 5.276209E-02
-7.055982E+03
-4.978688E+07
+4.196651E+00
+1.761188E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -79,14 +79,14 @@ tally 3:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.692488E+03
-5.917437E+07
+1.474427E+04
+2.173936E+08
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.692488E+03
-5.917437E+07
+1.474427E+04
+2.173936E+08
 0.000000E+00
 0.000000E+00
 tally 4:
@@ -104,14 +104,14 @@ tally 4:
 0.000000E+00
 2.297000E-01
 5.276209E-02
-7.055982E+03
-4.978688E+07
+4.196651E+00
+1.761188E+01
 0.000000E+00
 0.000000E+00
 2.297000E-01
 5.276209E-02
-7.055982E+03
-4.978688E+07
+4.196651E+00
+1.761188E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -134,7 +134,7 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.692488E+03
-5.917437E+07
+1.474427E+04
+2.173936E+08
 0.000000E+00
 0.000000E+00

tests/regression_tests/photon_production_fission/results_true.dat

@@ -32,14 +32,14 @@ tally 2:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.479234E+06
-2.052367E+12
+1.675918E+05
+9.365733E+09
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.479234E+06
-2.052367E+12
+1.675918E+05
+9.365733E+09
 0.000000E+00
 0.000000E+00
 tally 3:
@@ -57,13 +57,13 @@ tally 3:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.479234E+06
-2.052367E+12
+1.675918E+05
+9.365733E+09
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.479234E+06
-2.052367E+12
+1.675918E+05
+9.365733E+09
 0.000000E+00
 0.000000E+00

tests/unit_tests/test_photon_heating.py

@@ -0,0 +1,30 @@
+import openmc
+
+
+def test_negative_positron_heating():
+    m = openmc.Material()
+    m.add_element('Li', 1.0)
+    m.set_density('g/cm3', 10.0)
+
+    surf = openmc.Sphere(r=100.0, boundary_type='reflective')
+    cell = openmc.Cell(fill=m, region=-surf)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+    model.settings.run_mode = 'fixed source'
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point(),
+        energy=openmc.stats.Discrete([5.0e6], [1.0]),
+        particle='photon',
+    )
+    model.settings.particles = 7
+    model.settings.batches = 1
+    model.settings.electron_treatment = 'led'
+    model.settings.seed = 513836
+
+    tally = openmc.Tally()
+    tally.filters = [openmc.ParticleFilter(['photon', 'electron', 'positron'])]
+    tally.scores = ['heating']
+    model.tallies = openmc.Tallies([tally])
+    model.run(apply_tally_results=True)
+
+    assert (tally.mean >= 0.0).all(), "Negative heating detected"