diff --git a/src/photon.cpp b/src/photon.cpp index 6bbdc928f2c..9b40522c271 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -466,31 +466,25 @@ void PhotonInteraction::compton_doppler( if (rn < c) break; } + double E = alpha * MASS_ELECTRON_EV; // Determine binding energy of shell 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; - } + // Resample if photon energy insufficient + if (E < E_b) + continue; 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)) { + if (std::abs(pz_max) > data::compton_profile_pz(n - 1)) { c_max = profile_cdf_(shell, n - 1); } else { int i = lower_bound_index(data::compton_profile_pz.cbegin(), - data::compton_profile_pz.cend(), pz_max); + data::compton_profile_pz.cend(), std::abs(pz_max)); double pz_l = data::compton_profile_pz(i); double pz_r = data::compton_profile_pz(i + 1); double p_l = profile_pdf_(shell, i); @@ -499,16 +493,21 @@ void PhotonInteraction::compton_doppler( if (pz_l == pz_r) { c_max = c_l; } else if (p_l == p_r) { - c_max = c_l + (pz_max - pz_l) * p_l; + c_max = c_l + (std::abs(pz_max) - pz_l) * p_l; } else { double m = (p_l - p_r) / (pz_l - pz_r); - c_max = c_l + (std::pow((m * (pz_max - pz_l) + p_l), 2) - p_l * p_l) / + c_max = c_l + (std::pow((m * (std::abs(pz_max) - pz_l) + p_l), 2) - + p_l * p_l) / (2.0 * m); } } // Sample value on bounded cdf - c = prn(seed) * c_max; + if (pz_max > 0.0) { + c = prn(seed) * c_max; + } else { + c = 1.0 + (c_max - 1.0) * prn(seed); + } // Determine pz corresponding to sampled cdf value tensor::View cdf_shell = profile_cdf_.slice(shell); @@ -533,35 +532,23 @@ void PhotonInteraction::compton_doppler( double momentum_sq = std::pow((pz / FINE_STRUCTURE), 2); double f = 1.0 + alpha * (1.0 - mu); double a = momentum_sq - f * f; - double b = 2.0 * E * (f - momentum_sq * mu); - c = E * E * (momentum_sq - 1.0); - + double b = 2.0 * (f - momentum_sq * mu); + c = (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); - - // 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; - } - } else { - if (E_out2 > 0.0) { - *E_out = E_out2; - } else { - // No positive solution -- resample - continue; - } - } - if (*E_out < E - E_b) + double E_out1 = -(b + quad) / (2.0 * a) * E; + double E_out2 = -(b - quad) / (2.0 * a) * E; + // Determine solution to quadratic equation that is positive and minimal + if ((E_out1 < 0.0) && (E_out2 < 0.0)) + continue; + if ((E_out1 > 0.0) && (E_out2 > 0.0)) + *E_out = std::min(E_out1, E_out2); + else + *E_out = std::max(E_out1, E_out2); + + if (prn(seed) * E <= *E_out) break; } diff --git a/tests/regression_tests/atomic_relaxation/results_true.dat b/tests/regression_tests/atomic_relaxation/results_true.dat index 6f100dac8f8..e99899eedd1 100644 --- a/tests/regression_tests/atomic_relaxation/results_true.dat +++ b/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.920387E+00 +3.687888E+00 +7.854617E+04 +6.169501E+09 0.000000E+00 0.000000E+00 -9.208123E+05 -8.478953E+11 +9.214538E+05 +8.490772E+11 diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index c4f5fafa252..c959f73c311 100644 --- a/tests/regression_tests/photon_production/results_true.dat +++ b/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.493000E-01 +9.011705E-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.362831E-01 +5.582969E-02 0.000000E+00 0.000000E+00 -2.350047E-01 -5.522722E-02 +2.362831E-01 +5.582969E-02 0.000000E+00 0.000000E+00 0.000000E+00 @@ -53,40 +53,40 @@ tally 3: 4.104374E+12 4.296582E-04 1.846062E-07 -2.297000E-01 -5.276209E-02 -4.196651E+00 -1.761188E+01 +2.344000E-01 +5.494336E-02 +4.593008E+00 +2.109572E+01 0.000000E+00 0.000000E+00 -2.297000E-01 -5.276209E-02 -4.196651E+00 -1.761188E+01 +2.344000E-01 +5.494336E-02 +4.593008E+00 +2.109572E+01 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774484E+05 -3.148794E+10 +1.792949E+05 +3.214667E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774484E+05 -3.148794E+10 +1.792949E+05 +3.214667E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.474427E+04 -2.173936E+08 +1.452085E+04 +2.108550E+08 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.474427E+04 -2.173936E+08 +1.452085E+04 +2.108550E+08 0.000000E+00 0.000000E+00 tally 4: @@ -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.344000E-01 +5.494336E-02 +4.593008E+00 +2.109572E+01 0.000000E+00 0.000000E+00 -2.297000E-01 -5.276209E-02 -4.196651E+00 -1.761188E+01 +2.344000E-01 +5.494336E-02 +4.593008E+00 +2.109572E+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.792949E+05 +3.214667E+10 0.000000E+00 0.000000E+00 0.000000E+00 @@ -134,8 +134,8 @@ tally 4: 0.000000E+00 0.000000E+00 0.000000E+00 -1.474427E+04 -2.173936E+08 +1.452085E+04 +2.108550E+08 0.000000E+00 0.000000E+00 tally 5: diff --git a/tests/regression_tests/photon_production_fission/results_true.dat b/tests/regression_tests/photon_production_fission/results_true.dat index af325d4b02a..8c3783ce196 100644 --- a/tests/regression_tests/photon_production_fission/results_true.dat +++ b/tests/regression_tests/photon_production_fission/results_true.dat @@ -1,12 +1,12 @@ k-combined: -2.297165E+00 1.955494E-02 +2.255088E+00 7.363219E-02 tally 1: -2.696393E+00 -2.423937E+00 +2.640439E+00 +2.329186E+00 0.000000E+00 0.000000E+00 -2.696393E+00 -2.423937E+00 +2.640439E+00 +2.329186E+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.625521E+00 +2.303289E+00 +4.210498E+08 +5.923082E+16 0.000000E+00 0.000000E+00 -2.672029E+00 -2.380251E+00 -4.288244E+08 -6.130569E+16 +2.625521E+00 +2.303289E+00 +4.210498E+08 +5.923082E+16 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.711908E+05 -9.769096E+09 +1.826906E+05 +1.114925E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.711908E+05 -9.769096E+09 +1.826906E+05 +1.114925E+10 0.000000E+00 0.000000E+00 tally 3: -2.649127E+00 -2.339294E+00 -4.288244E+08 -6.130569E+16 +2.658439E+00 +2.355784E+00 +4.210498E+08 +5.923082E+16 0.000000E+00 0.000000E+00 -2.649127E+00 -2.339294E+00 -4.288244E+08 -6.130569E+16 +2.658439E+00 +2.355784E+00 +4.210498E+08 +5.923082E+16 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.711908E+05 -9.769096E+09 +1.826906E+05 +1.114925E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.711908E+05 -9.769096E+09 +1.826906E+05 +1.114925E+10 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/photon_source/results_true.dat b/tests/regression_tests/photon_source/results_true.dat index 8d934afc66a..ff23d36d444 100644 --- a/tests/regression_tests/photon_source/results_true.dat +++ b/tests/regression_tests/photon_source/results_true.dat @@ -1,5 +1,5 @@ tally 1: -2.263761E+02 -5.124615E+04 +2.193647E+02 +4.812086E+04 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/pulse_height/inputs_true.dat b/tests/regression_tests/pulse_height/inputs_true.dat index 590928e4355..f38fb8f9688 100644 --- a/tests/regression_tests/pulse_height/inputs_true.dat +++ b/tests/regression_tests/pulse_height/inputs_true.dat @@ -15,8 +15,8 @@ fixed source - 100 - 5 + 10000 + 1 0.0 0.0 0.0 diff --git a/tests/regression_tests/pulse_height/results_true.dat b/tests/regression_tests/pulse_height/results_true.dat index c57e8ff1c83..84fcfc0f2de 100644 --- a/tests/regression_tests/pulse_height/results_true.dat +++ b/tests/regression_tests/pulse_height/results_true.dat @@ -1,186 +1,186 @@ tally 1: -4.140000E+00 -3.443000E+00 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -3.000000E-02 -5.000000E-04 -2.000000E-02 -4.000000E-04 -2.000000E-02 -4.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -3.000000E-02 -3.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -4.000000E-02 +8.122000E-01 +6.596688E-01 +2.300000E-03 +5.290000E-06 +1.400000E-03 +1.960000E-06 +1.200000E-03 +1.440000E-06 +1.400000E-03 +1.960000E-06 +1.600000E-03 +2.560000E-06 +1.500000E-03 +2.250000E-06 +1.000000E-03 +1.000000E-06 +2.100000E-03 +4.410000E-06 +1.400000E-03 +1.960000E-06 +2.300000E-03 +5.290000E-06 +1.900000E-03 +3.610000E-06 +2.200000E-03 +4.840000E-06 +2.100000E-03 +4.410000E-06 +1.100000E-03 +1.210000E-06 +1.400000E-03 +1.960000E-06 +1.900000E-03 +3.610000E-06 +1.200000E-03 +1.440000E-06 +1.700000E-03 +2.890000E-06 +2.600000E-03 +6.760000E-06 +1.600000E-03 +2.560000E-06 +1.500000E-03 +2.250000E-06 +1.800000E-03 +3.240000E-06 +1.500000E-03 +2.250000E-06 +1.900000E-03 +3.610000E-06 +1.300000E-03 +1.690000E-06 +1.800000E-03 +3.240000E-06 +2.000000E-03 +4.000000E-06 +2.200000E-03 +4.840000E-06 +1.400000E-03 +1.960000E-06 +1.900000E-03 +3.610000E-06 +2.300000E-03 +5.290000E-06 +2.100000E-03 +4.410000E-06 +1.100000E-03 +1.210000E-06 +1.000000E-03 +1.000000E-06 +1.600000E-03 +2.560000E-06 +1.600000E-03 +2.560000E-06 +1.900000E-03 +3.610000E-06 +2.100000E-03 +4.410000E-06 +1.500000E-03 +2.250000E-06 +2.400000E-03 +5.760000E-06 +1.600000E-03 +2.560000E-06 +1.800000E-03 +3.240000E-06 +1.800000E-03 +3.240000E-06 +1.900000E-03 +3.610000E-06 +2.000000E-03 +4.000000E-06 +2.200000E-03 +4.840000E-06 +1.600000E-03 +2.560000E-06 +1.200000E-03 +1.440000E-06 +1.500000E-03 +2.250000E-06 +2.100000E-03 +4.410000E-06 +2.300000E-03 +5.290000E-06 +2.600000E-03 +6.760000E-06 +2.300000E-03 +5.290000E-06 +2.200000E-03 +4.840000E-06 +1.300000E-03 +1.690000E-06 +2.100000E-03 +4.410000E-06 +2.100000E-03 +4.410000E-06 +2.800000E-03 +7.840000E-06 +2.700000E-03 +7.290000E-06 +1.300000E-03 +1.690000E-06 +1.900000E-03 +3.610000E-06 +1.700000E-03 +2.890000E-06 +1.200000E-03 +1.440000E-06 +2.700000E-03 +7.290000E-06 +1.900000E-03 +3.610000E-06 +1.700000E-03 +2.890000E-06 +1.100000E-03 +1.210000E-06 +1.700000E-03 +2.890000E-06 +1.900000E-03 +3.610000E-06 +2.300000E-03 +5.290000E-06 +1.400000E-03 +1.960000E-06 +1.400000E-03 +1.960000E-06 +2.800000E-03 +7.840000E-06 +2.700000E-03 +7.290000E-06 +2.100000E-03 +4.410000E-06 +1.900000E-03 +3.610000E-06 +2.100000E-03 +4.410000E-06 +2.500000E-03 +6.250000E-06 +1.700000E-03 +2.890000E-06 +1.200000E-03 +1.440000E-06 +1.000000E-03 +1.000000E-06 6.000000E-04 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -2.000000E-02 -2.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -2.000000E-02 -4.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -2.000000E-02 -2.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -3.000000E-02 +3.600000E-07 5.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -2.000000E-02 -2.000000E-04 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -2.000000E-02 -2.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -2.000000E-02 +2.500000E-07 +8.000000E-04 +6.400000E-07 +3.000000E-04 +9.000000E-08 +3.000000E-04 +9.000000E-08 2.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -3.000000E-02 -5.000000E-04 +4.000000E-08 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -1.000000E-02 1.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -2.000000E-02 -2.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 +1.000000E-08 0.000000E+00 0.000000E+00 0.000000E+00 @@ -197,5 +197,5 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -2.000000E-01 -8.600000E-03 +3.680000E-02 +1.354240E-03 diff --git a/tests/regression_tests/pulse_height/test.py b/tests/regression_tests/pulse_height/test.py index 90d960f6640..81591379e60 100644 --- a/tests/regression_tests/pulse_height/test.py +++ b/tests/regression_tests/pulse_height/test.py @@ -27,8 +27,8 @@ def sphere_model(): # Define settings model.settings.run_mode = 'fixed source' - model.settings.batches = 5 - model.settings.particles = 100 + model.settings.batches = 1 + model.settings.particles = 10000 model.settings.photon_transport = True model.settings.source = openmc.IndependentSource( space=openmc.stats.Point(), @@ -49,5 +49,5 @@ def sphere_model(): def test_pulse_height(sphere_model): - harness = PyAPITestHarness('statepoint.5.h5', sphere_model) + harness = PyAPITestHarness('statepoint.1.h5', sphere_model) harness.main() diff --git a/tests/regression_tests/weightwindows/results_true.dat b/tests/regression_tests/weightwindows/results_true.dat index 897089e0597..a3b10b49b58 100644 --- a/tests/regression_tests/weightwindows/results_true.dat +++ b/tests/regression_tests/weightwindows/results_true.dat @@ -1 +1 @@ -a4a3ccca43666e2ca1e71800201b152cca20c387b93d67522c5339807348dcee5cada9acbed3238f37e2e86e76b374b06988742f07d4ea1b413e4e75d0c180b1 \ No newline at end of file +bb3b76ee105f20bc9eb763fd04753f6212ab0139197f5be03b927769ac4f45efa77dd22cdf48900062c3eaa8874f463bef11c0dd52321ba8f6da2c12c98247b7 \ No newline at end of file diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index d6e509522fa..f3e63dcf867 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/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']