Fix a bug in compton scattering shell and angle selection

include/openmc/photon.h

@@ -87,7 +87,7 @@ class PhotonInteraction {
   tensor::Tensor<double> profile_pdf_;
   tensor::Tensor<double> profile_cdf_;
   tensor::Tensor<double> binding_energy_;
-  tensor::Tensor<double> electron_pdf_;
+  tensor::Tensor<double> electron_cdf_;
 
   // Map subshells from Compton profile data obtained from Biggs et al,
   // "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic

src/photon.cpp

@@ -214,8 +214,13 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   rgroup = open_group(group, "compton_profiles");
 
   // Read electron shell PDF and binding energies
-  read_dataset(rgroup, "num_electrons", electron_pdf_);
-  electron_pdf_ /= electron_pdf_.sum();
+  tensor::Tensor<double> electron_pdf;
+  read_dataset(rgroup, "num_electrons", electron_pdf);
+  electron_pdf /= electron_pdf.sum();
+  electron_cdf_ = tensor::Tensor<double>(electron_pdf.shape());
+  electron_cdf_(0) = electron_pdf(0);
+  for (int i = 1; i < electron_pdf.shape()[0]; ++i)
+    electron_cdf_(i) = electron_cdf_(i - 1) + electron_pdf(i);
   read_dataset(rgroup, "binding_energy", binding_energy_);
 
   // Read Compton profiles
@@ -420,10 +425,28 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler,
   double* alpha_out, double* mu, int* i_shell, uint64_t* seed) const
 {
   double form_factor_xmax = 0.0;
+  int last_shell = binding_energy_.shape()[0] - 1;
+  double E_b = binding_energy_(last_shell);
+  double E = alpha * MASS_ELECTRON_EV;
+  double mu_max = 1 - E_b / (alpha * (E - E_b));
+
   while (true) {
     // Sample Klein-Nishina distribution for trial energy and angle
     std::tie(*alpha_out, *mu) = klein_nishina(alpha, seed);
 
+    // If in every angle we cannot eject an electron
+    // Exit with no shell
+    if (mu_max < -1) {
+      *i_shell = -1;
+      return;
+    }
+
+    if (doppler) {
+      // Reject angles that cannot eject the most loosely bound electron
+      if (*mu > mu_max)
+        continue;
+    }
+
     // Note that the parameter used here does not correspond exactly to the
     // momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the
     // parameter as defined by Hubbell, where the actual data comes from
@@ -455,14 +478,25 @@ void PhotonInteraction::compton_doppler(
   double alpha, double mu, double* E_out, int* i_shell, uint64_t* seed) const
 {
   auto n = data::compton_profile_pz.size();
+  double E = alpha * MASS_ELECTRON_EV;
+  int j_shell = 0;
+  for (double E_b : binding_energy_) {
+    if ((E_b - (E - E_b) * alpha * (1.0 - mu)) < 0.0)
+      break;
+    ++j_shell;
+  }
+
+  double offset = 0.0;
+  if (j_shell > 0)
+    offset = electron_cdf_(j_shell - 1);
 
   int shell; // index for shell
   while (true) {
     // Sample electron shell
-    double rn = prn(seed);
-    double c = 0.0;
-    for (shell = 0; shell < electron_pdf_.size(); ++shell) {
-      c += electron_pdf_(shell);
+    double rn = offset + prn(seed) * (1.0 - offset);
+    double c;
+    for (shell = j_shell; shell < electron_cdf_.size(); ++shell) {
+      c = electron_cdf_(shell);
       if (rn < c)
         break;
     }
@@ -471,19 +505,8 @@ void PhotonInteraction::compton_doppler(
     double E_b = binding_energy_(shell);
 
     // Determine p_z,max
-    double E = alpha * MASS_ELECTRON_EV;
-    if (E < E_b) {
-      *E_out = alpha / (1 + alpha * (1 - mu)) * MASS_ELECTRON_EV;
-      break;
-    }
-
     double pz_max = -FINE_STRUCTURE * (E_b - (E - E_b) * alpha * (1.0 - mu)) /
                     std::sqrt(2.0 * E * (E - E_b) * (1.0 - mu) + E_b * E_b);
-    if (pz_max < 0.0) {
-      *E_out = alpha / (1 + alpha * (1 - mu)) * MASS_ELECTRON_EV;
-      break;
-    }
-
     // Determine profile cdf value corresponding to p_z,max
     double c_max;
     if (pz_max > data::compton_profile_pz(n - 1)) {
@@ -537,29 +560,21 @@ void PhotonInteraction::compton_doppler(
     c = E * E * (momentum_sq - 1.0);
 
     double quad = b * b - 4.0 * a * c;
-    if (quad < 0) {
-      *E_out = alpha / (1 + alpha * (1 - mu)) * MASS_ELECTRON_EV;
-      break;
-    }
+    if (quad < 0)
+      continue;
     quad = std::sqrt(quad);
     double E_out1 = -(b + quad) / (2.0 * a);
     double E_out2 = -(b - quad) / (2.0 * a);
 
+    // If no positive solution -- resample
+    if (std::max(E_out1, E_out2) < 0)
+      continue;
+
     // Determine solution to quadratic equation that is positive
-    if (E_out1 > 0.0) {
-      if (E_out2 > 0.0) {
-        // If both are positive, pick one at random
-        *E_out = prn(seed) < 0.5 ? E_out1 : E_out2;
-      } else {
-        *E_out = E_out1;
-      }
+    if ((E_out1 > 0.0) && (E_out2 > 0.0)) {
+      *E_out = prn(seed) < 0.5 ? E_out1 : E_out2;
     } else {
-      if (E_out2 > 0.0) {
-        *E_out = E_out2;
-      } else {
-        // No positive solution -- resample
-        continue;
-      }
+      *E_out = std::max(E_out1, E_out2);
     }
     if (*E_out < E - E_b)
       break;

tests/regression_tests/atomic_relaxation/results_true.dat

@@ -1,9 +1,9 @@
 tally 1:
-1.956204E+00
-3.826732E+00
-7.918768E+04
-6.270688E+09
+1.955720E+00
+3.824843E+00
+7.895341E+04
+6.233641E+09
 0.000000E+00
 0.000000E+00
-9.208123E+05
-8.478953E+11
+9.210466E+05
+8.483268E+11

tests/regression_tests/photon_production/results_true.dat

@@ -1,8 +1,8 @@
 tally 1:
 8.610000E-01
 7.413210E-01
-9.491000E-01
-9.007908E-01
+9.489000E-01
+9.004112E-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.350047E-01
-5.522722E-02
+2.342279E-01
+5.486270E-02
 0.000000E+00
 0.000000E+00
-2.350047E-01
-5.522722E-02
+2.342279E-01
+5.486270E-02
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -53,28 +53,28 @@ tally 3:
 4.104374E+12
 4.296582E-04
 1.846062E-07
-2.297000E-01
-5.276209E-02
-4.196651E+00
-1.761188E+01
+2.294000E-01
+5.262436E-02
+4.488894E+00
+2.015017E+01
 0.000000E+00
 0.000000E+00
-2.297000E-01
-5.276209E-02
-4.196651E+00
-1.761188E+01
+2.294000E-01
+5.262436E-02
+4.488894E+00
+2.015017E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.774484E+05
-3.148794E+10
+1.774550E+05
+3.149027E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.774484E+05
-3.148794E+10
+1.774550E+05
+3.149027E+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.297000E-01
-5.276209E-02
-4.196651E+00
-1.761188E+01
+2.294000E-01
+5.262436E-02
+4.488894E+00
+2.015017E+01
 0.000000E+00
 0.000000E+00
-2.297000E-01
-5.276209E-02
-4.196651E+00
-1.761188E+01
+2.294000E-01
+5.262436E-02
+4.488894E+00
+2.015017E+01
 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.774484E+05
-3.148794E+10
+1.774550E+05
+3.149027E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00

tests/regression_tests/photon_production_fission/results_true.dat

@@ -1,12 +1,12 @@
 k-combined:
-2.297165E+00 1.955494E-02
+2.331233E+00 3.816178E-02
 tally 1:
-2.696393E+00
-2.423937E+00
+2.733254E+00
+2.491835E+00
 0.000000E+00
 0.000000E+00
-2.696393E+00
-2.423937E+00
+2.733254E+00
+2.491835E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -18,52 +18,52 @@ tally 1:
 0.000000E+00
 0.000000E+00
 tally 2:
-2.672029E+00
-2.380251E+00
-4.288244E+08
-6.130569E+16
+2.705904E+00
+2.441602E+00
+4.339210E+08
+6.278920E+16
 0.000000E+00
 0.000000E+00
-2.672029E+00
-2.380251E+00
-4.288244E+08
-6.130569E+16
+2.705904E+00
+2.441602E+00
+4.339210E+08
+6.278920E+16
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.879204E+05
+1.177358E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.879204E+05
+1.177358E+10
 0.000000E+00
 0.000000E+00
 tally 3:
-2.649127E+00
-2.339294E+00
-4.288244E+08
-6.130569E+16
+2.658522E+00
+2.355935E+00
+4.339210E+08
+6.278920E+16
 0.000000E+00
 0.000000E+00
-2.649127E+00
-2.339294E+00
-4.288244E+08
-6.130569E+16
+2.658522E+00
+2.355935E+00
+4.339210E+08
+6.278920E+16
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.879204E+05
+1.177358E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.879204E+05
+1.177358E+10
 0.000000E+00
 0.000000E+00

tests/regression_tests/photon_source/results_true.dat

@@ -1,5 +1,5 @@
 tally 1:
-2.263761E+02
-5.124615E+04
+2.263259E+02
+5.122343E+04
 0.000000E+00
 0.000000E+00

tests/regression_tests/weightwindows/results_true.dat

@@ -1 +1 @@
-a4a3ccca43666e2ca1e71800201b152cca20c387b93d67522c5339807348dcee5cada9acbed3238f37e2e86e76b374b06988742f07d4ea1b413e4e75d0c180b1
+b51d86a9d3a00713fb6d95690adca37d6af93f08884a233caa88b8a21e2a3985c85b6b4d150dd74b19ccc0d4c3fa973d2114309096642b2ea7e38d78ee1d8186

tests/unit_tests/weightwindows/test.py

@@ -246,7 +246,7 @@ def test_photon_heating(run_in_tmpdir):
 
     model.settings.run_mode = 'fixed source'
     model.settings.batches = 5
-    model.settings.particles = 100
+    model.settings.particles = 1000
 
     tally = openmc.Tally()
     tally.scores = ['heating']