Fixing an incorrect computation of CDF of bremsstrahlung photons (#3396)

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

Author: GuySten

src/material.cpp

@@ -777,14 +777,15 @@ void Material::init_bremsstrahlung()
       // Loop over photon energies
       double c = 0.0;
       for (int i = 0; i < j; ++i) {
-        // Integrate the CDF from the PDF using the trapezoidal rule in log-log
-        // space
+        // Integrate the CDF from the PDF using the fact that the PDF is linear
+        // in log-log space
         double w_l = std::log(data::ttb_e_grid(i));
         double w_r = std::log(data::ttb_e_grid(i + 1));
         double x_l = std::log(ttb->pdf(j, i));
         double x_r = std::log(ttb->pdf(j, i + 1));
-
-        c += 0.5 * (w_r - w_l) * (std::exp(w_l + x_l) + std::exp(w_r + x_r));
+        double beta = (x_r - x_l) / (w_r - w_l);
+        double a = beta + 1.0;
+        c += std::exp(w_l + x_l) / a * std::expm1(a * (w_r - w_l));
         ttb->cdf(j, i + 1) = c;
       }
 

tests/regression_tests/photon_production/results_true.dat

@@ -16,12 +16,12 @@ tally 2:
 1.573004E+00
 4.296434E-04
 1.845934E-07
-2.350021E-01
-5.522600E-02
+2.350047E-01
+5.522722E-02
 0.000000E+00
 0.000000E+00
-2.350021E-01
-5.522600E-02
+2.350047E-01
+5.522722E-02
 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.774451E+05
-3.148675E+10
+1.774484E+05
+3.148794E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -85,8 +85,8 @@ tally 3:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-7.691658E+03
-5.916160E+07
+7.692488E+03
+5.917437E+07
 0.000000E+00
 0.000000E+00
 tally 4:
@@ -122,8 +122,8 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.774451E+05
-3.148675E+10
+1.774484E+05
+3.148794E+10
 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.691658E+03
-5.916160E+07
+7.692488E+03
+5.917437E+07
 0.000000E+00
 0.000000E+00

tests/regression_tests/photon_production_fission/results_true.dat

@@ -1,12 +1,12 @@
 k-combined:
-2.270911E+00 4.568134E-02
+2.278476E+00 6.220292E-02
 tally 1:
-2.663476E+00
-2.366726E+00
+2.664071E+00
+2.369250E+00
 0.000000E+00
 0.000000E+00
-2.663476E+00
-2.366726E+00
+2.664071E+00
+2.369250E+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.636709E+00
-2.318646E+00
-4.230052E+08
-5.967879E+16
+2.640755E+00
+2.326201E+00
+4.245217E+08
+6.012128E+16
 0.000000E+00
 0.000000E+00
-2.636709E+00
-2.318646E+00
-4.230052E+08
-5.967879E+16
+2.640755E+00
+2.326201E+00
+4.245217E+08
+6.012128E+16
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.455009E+06
-2.015314E+12
+2.479234E+06
+2.052367E+12
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.455009E+06
-2.015314E+12
+2.479234E+06
+2.052367E+12
 0.000000E+00
 0.000000E+00
 tally 3:
-2.660004E+00
-2.358564E+00
-4.230052E+08
-5.967879E+16
+2.657846E+00
+2.354717E+00
+4.245217E+08
+6.012128E+16
 0.000000E+00
 0.000000E+00
-2.660004E+00
-2.358564E+00
-4.230052E+08
-5.967879E+16
+2.657846E+00
+2.354717E+00
+4.245217E+08
+6.012128E+16
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.455009E+06
-2.015314E+12
+2.479234E+06
+2.052367E+12
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.455009E+06
-2.015314E+12
+2.479234E+06
+2.052367E+12
 0.000000E+00
 0.000000E+00

tests/regression_tests/photon_source/results_true.dat

@@ -1,5 +1,5 @@
 tally 1:
-2.263938E+02
-5.125417E+04
+2.263761E+02
+5.124615E+04
 0.000000E+00
 0.000000E+00

tests/regression_tests/weightwindows/results_true.dat

@@ -1 +1 @@
-952c20fefac374d3b9fd28627fda7d5ae262f8cd7c01a33f0381526204a2287c18e56259a599dc6fdc1b59a2d2866fdfeedea371154c8fa7a69dc5c445b08c2d
+a5880ad9262e8aba90801783891ee74618144101401f06fb46e954e851a3c517ab28d5f0d6e1b2b364844f08d363cba35ab23b63fc81012ea8a6a328755d56c4