From 55c8212954416e6bcf000fc80bfd55d72ce0c465 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Wed, 11 Mar 2026 15:00:54 +0000 Subject: [PATCH] fix negative heating from Compton scattering with atomic relaxation When Doppler broadening causes the Compton-scattered photon to transfer less energy than the shell binding energy, the recoil electron has negative kinetic energy and is not created. However, atomic relaxation still runs and produces secondaries totaling ~e_b. The heating score (E_last - E' - bank_second_E) went negative because bank_second_E contained only relaxation products (~e_b) which exceeded the photon's energy transfer. Fix: when E_electron < 0, add it (a negative value) to bank_second_E before relaxation runs. This offsets the relaxation energy so bank_second_E = E_electron + sum(relaxation) = photon_transfer, giving a heating score of zero. This is equivalent to Geant4 PENELOPE's approach of setting localEnergyDeposit = diffEnergy when eKineticEnergy < 0 (G4PenelopeComptonModel.cc). Also use relaxation binding energy (ENDF/B) instead of Compton profile binding energy (Biggs et al.) for the electron energy calculation when atomic relaxation is enabled, matching the approach already used by the photoelectric effect. This ensures the electron energy and relaxation cascade use a self-consistent binding energy source. Co-Authored-By: Claude Opus 4.6 --- src/physics.cpp | 21 +++++++- .../electron_heating/results_true.dat | 4 +- .../photon_production/results_true.dat | 16 +++---- .../results_true.dat | 16 +++---- tests/unit_tests/test_photon_heating.py | 48 +++++++++++++++++++ 5 files changed, 86 insertions(+), 19 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index 106bd1aa2b2..59259e71d4f 100644 --- a/src/physics.cpp +++ b/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); } diff --git a/tests/regression_tests/electron_heating/results_true.dat b/tests/regression_tests/electron_heating/results_true.dat index 4f54ceaa4d5..268b137b16a 100644 --- a/tests/regression_tests/electron_heating/results_true.dat +++ b/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 diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index c4f5fafa252..249972050ca 100644 --- a/tests/regression_tests/photon_production/results_true.dat +++ b/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 diff --git a/tests/regression_tests/photon_production_fission/results_true.dat b/tests/regression_tests/photon_production_fission/results_true.dat index af325d4b02a..cde059d06da 100644 --- a/tests/regression_tests/photon_production_fission/results_true.dat +++ b/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 diff --git a/tests/unit_tests/test_photon_heating.py b/tests/unit_tests/test_photon_heating.py index 05473f5c6fc..b8900686990 100644 --- a/tests/unit_tests/test_photon_heating.py +++ b/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)