Map Compton subshell data to atomic relaxation data (#3392)

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

Author: amandalund

include/openmc/photon.h

@@ -33,7 +33,6 @@ class ElectronSubshell {
 
   int index_subshell; //!< index in SUBSHELLS
   int threshold;
-  double n_electrons;
   double binding_energy;
   vector<Transition> transitions;
 };
@@ -90,6 +89,11 @@ class PhotonInteraction {
   xt::xtensor<double, 1> binding_energy_;
   xt::xtensor<double, 1> electron_pdf_;
 
+  // Map subshells from Compton profile data obtained from Biggs et al,
+  // "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic
+  // relaxation data
+  xt::xtensor<int, 1> subshell_map_;
+
   // Stopping power data
   double I_; // mean excitation energy
   xt::xtensor<int, 1> n_electrons_;

src/photon.cpp

@@ -165,7 +165,6 @@ PhotonInteraction::PhotonInteraction(hid_t group)
     if (attribute_exists(tgroup, "binding_energy")) {
       has_atomic_relaxation_ = true;
       read_attribute(tgroup, "binding_energy", shell.binding_energy);
-      read_attribute(tgroup, "num_electrons", shell.n_electrons);
     }
 
     // Read subshell cross section
@@ -233,6 +232,28 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   }
   close_group(rgroup);
 
+  // Map Compton subshell data to atomic relaxation data by finding the
+  // subshell with the equivalent binding energy
+  if (has_atomic_relaxation_) {
+    auto is_close = [](double a, double b) {
+      return std::abs(a - b) / a < FP_REL_PRECISION;
+    };
+    subshell_map_ = xt::full_like(binding_energy_, -1);
+    for (int i = 0; i < binding_energy_.size(); ++i) {
+      double E_b = binding_energy_[i];
+      if (i < n_shell && is_close(E_b, shells_[i].binding_energy)) {
+        subshell_map_[i] = i;
+      } else {
+        for (int j = 0; j < n_shell; ++j) {
+          if (is_close(E_b, shells_[j].binding_energy)) {
+            subshell_map_[i] = j;
+            break;
+          }
+        }
+      }
+    }
+  }
+
   // Create Compton profile CDF
   auto n_profile = data::compton_profile_pz.size();
   auto n_shell_compton = profile_pdf_.shape(0);
@@ -423,13 +444,6 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler,
         double E_out;
         this->compton_doppler(alpha, *mu, &E_out, i_shell, seed);
         *alpha_out = E_out / MASS_ELECTRON_EV;
-
-        // It's possible for the Compton profile data to have more shells than
-        // there are in the ENDF data. Make sure the shell index doesn't end up
-        // out of bounds.
-        if (*i_shell >= shells_.size()) {
-          *i_shell = -1;
-        }
       } else {
         *i_shell = -1;
       }

src/physics.cpp

@@ -340,11 +340,11 @@ void sample_photon_reaction(Particle& p)
       p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
     }
 
-    // TODO: Compton subshell data does not match atomic relaxation data
-    // Allow electrons to fill orbital and produce auger electrons
-    // and fluorescent photons
-    if (i_shell >= 0) {
-      element.atomic_relaxation(i_shell, p);
+    // Allow electrons to fill orbital and produce Auger electrons and
+    // fluorescent photons. Since Compton subshell data does not match atomic
+    // relaxation data, use the mapping between the data to find the subshell
+    if (i_shell >= 0 && element.subshell_map_[i_shell] >= 0) {
+      element.atomic_relaxation(element.subshell_map_[i_shell], p);
     }
 
     phi += PI;

tests/regression_tests/photon_production/results_true.dat

@@ -1,8 +1,8 @@
 tally 1:
 8.610000E-01
 7.413210E-01
-9.493000E-01
-9.011705E-01
+9.491000E-01
+9.007908E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -16,12 +16,12 @@ tally 2:
 1.573004E+00
 4.296434E-04
 1.845934E-07
-2.337049E-01
-5.461796E-02
+2.350021E-01
+5.522600E-02
 0.000000E+00
 0.000000E+00
-2.337049E-01
-5.461796E-02
+2.350021E-01
+5.522600E-02
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -53,16 +53,16 @@ tally 3:
 4.104374E+12
 4.296582E-04
 1.846062E-07
-2.286000E-01
-5.225796E-02
-7.054033E+03
-4.975938E+07
+2.297000E-01
+5.276209E-02
+7.055982E+03
+4.978688E+07
 0.000000E+00
 0.000000E+00
-2.286000E-01
-5.225796E-02
-7.054033E+03
-4.975938E+07
+2.297000E-01
+5.276209E-02
+7.055982E+03
+4.978688E+07
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -73,8 +73,8 @@ tally 3:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.764573E+05
-3.113718E+10
+1.774451E+05
+3.148675E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -102,16 +102,16 @@ tally 4:
 4.104374E+12
 0.000000E+00
 0.000000E+00
-2.286000E-01
-5.225796E-02
-7.054033E+03
-4.975938E+07
+2.297000E-01
+5.276209E-02
+7.055982E+03
+4.978688E+07
 0.000000E+00
 0.000000E+00
-2.286000E-01
-5.225796E-02
-7.054033E+03
-4.975938E+07
+2.297000E-01
+5.276209E-02
+7.055982E+03
+4.978688E+07
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -122,8 +122,8 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.764573E+05
-3.113718E+10
+1.774451E+05
+3.148675E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00

tests/regression_tests/weightwindows/results_true.dat

@@ -1 +1 @@
-ebc761815175b25fc95a226174928c226a3ab5dbf3b2a2abf09e079a0b87dee1dee74aa9b9eaec35acd58b8c481d264be7e9b3f052905bcc14ccd76f36b01549
+952c20fefac374d3b9fd28627fda7d5ae262f8cd7c01a33f0381526204a2287c18e56259a599dc6fdc1b59a2d2866fdfeedea371154c8fa7a69dc5c445b08c2d