Potential fix for negative heating from Compton scattering with atomic relaxation

src/physics.cpp

@@ -332,10 +332,18 @@ void sample_photon_reaction(Particle& p)
       alpha, true, &alpha_out, &p.mu(), &i_shell, p.current_seed());
 
     // Determine binding energy of shell. The binding energy is 0.0 if
-    // doppler broadening is not used.
+    // doppler broadening is not used. When atomic relaxation data is
+    // present, use the relaxation binding energy so that the electron
+    // energy and the relaxation cascade are self-consistent (same
+    // approach as the photoelectric effect). Using the Compton profile
+    // binding energy here while the cascade uses relaxation binding
+    // energies can produce negative heating scores.
     double e_b;
     if (i_shell == -1) {
       e_b = 0.0;
+    } else if (settings::atomic_relaxation &&
+               element.subshell_map_[i_shell] >= 0) {
+      e_b = element.shells_[element.subshell_map_[i_shell]].binding_energy;
     } else {
       e_b = element.binding_energy_[i_shell];
     }
@@ -357,6 +365,17 @@ void sample_photon_reaction(Particle& p)
     // relaxation data, use the mapping between the data to find the subshell
     if (settings::atomic_relaxation && i_shell >= 0 &&
         element.subshell_map_[i_shell] >= 0) {
+      // When E_electron < 0, the photon transferred less energy than the
+      // binding energy but the relaxation cascade still releases ~e_b of
+      // energy from the atom. The electron is not created, so bank_second_E
+      // would only contain relaxation products, making the heating score
+      // (E_last - E' - bank_second_E) negative. Adding E_electron (which is
+      // negative) offsets this so the energy balance is correct. This is
+      // equivalent to Geant4's PENELOPE approach of setting
+      // localEnergyDeposit = diffEnergy when eKineticEnergy < 0.
+      if (E_electron < 0.0) {
+        p.bank_second_E() += E_electron;
+      }
       element.atomic_relaxation(element.subshell_map_[i_shell], p);
     }
 

tests/regression_tests/electron_heating/results_true.dat

@@ -1,3 +1,3 @@
 tally 1:
-1.000000E+07
-1.000000E+14
+1.000002E+07
+1.000003E+14

tests/regression_tests/photon_production/results_true.dat

@@ -55,14 +55,14 @@ tally 3:
 1.846062E-07
 2.297000E-01
 5.276209E-02
-4.196651E+00
-1.761188E+01
+4.497835E+00
+2.023052E+01
 0.000000E+00
 0.000000E+00
 2.297000E-01
 5.276209E-02
-4.196651E+00
-1.761188E+01
+4.497835E+00
+2.023052E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -104,14 +104,14 @@ tally 4:
 0.000000E+00
 2.297000E-01
 5.276209E-02
-4.196651E+00
-1.761188E+01
+4.497835E+00
+2.023052E+01
 0.000000E+00
 0.000000E+00
 2.297000E-01
 5.276209E-02
-4.196651E+00
-1.761188E+01
+4.497835E+00
+2.023052E+01
 0.000000E+00
 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
-1.711908E+05
-9.769096E+09
+1.920206E+05
+1.229102E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.920206E+05
+1.229102E+10
 0.000000E+00
 0.000000E+00
 tally 3:
@@ -57,13 +57,13 @@ tally 3:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.920206E+05
+1.229102E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.920206E+05
+1.229102E+10
 0.000000E+00
 0.000000E+00

tests/unit_tests/test_photon_heating.py

@@ -1,6 +1,54 @@
+import numpy as np
+
 import openmc
 
 
+def test_compton_relaxation_heating():
+    """Test that Compton scattering with atomic relaxation does not produce
+    negative photon heating scores. When Doppler broadening causes the photon
+    to transfer less energy than the shell binding energy, the electron is not
+    ejected and atomic relaxation should not occur. Previously, relaxation ran
+    unconditionally, creating secondaries whose energy exceeded the photon's
+    energy transfer, producing negative heating.
+    """
+    # Lead has high-Z (82) with many inner shells and large binding energies,
+    # maximizing the number of Compton events where the energy transfer is
+    # less than the shell binding energy.
+    lead = openmc.Material()
+    lead.add_element('Pb', 1.0)
+    lead.set_density('g/cm3', 11.35)
+
+    box = openmc.model.RectangularParallelepiped(
+        -10, 10, -10, 10, -10, 10, boundary_type='reflective')
+    cell = openmc.Cell(region=-box, fill=lead)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point((0, 0, 0)),
+        energy=openmc.stats.delta_function(200e3),
+        particle='photon')
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 500
+
+    mesh = openmc.RegularMesh.from_domain(model.geometry, dimension=(3, 3, 3))
+    tally = openmc.Tally()
+    tally.scores = ['heating']
+    tally.filters = [
+        openmc.ParticleFilter(['photon']),
+        openmc.MeshFilter(mesh)
+    ]
+    model.tallies = [tally]
+
+    sp_file = model.run()
+    with openmc.StatePoint(sp_file) as sp:
+        tally_mean = sp.tallies[tally.id].mean
+
+    assert np.all(tally_mean >= 0), \
+        "Negative photon heating from Compton + atomic relaxation"
+
+
 def test_negative_positron_heating():
     m = openmc.Material()
     m.add_element('Li', 1.0)