From fe7c6cda2d28c98495da96a60725df6f7d37e38a Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 12 Mar 2026 19:03:20 +0200 Subject: [PATCH 01/19] reject shell if it is of higher energy relative to incident photon --- src/photon.cpp | 3 +-- tests/unit_tests/weightwindows/test.py | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 6bbdc928f2c..92124533af6 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -473,8 +473,7 @@ void PhotonInteraction::compton_doppler( // 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; + continue; } double pz_max = -FINE_STRUCTURE * (E_b - (E - E_b) * alpha * (1.0 - mu)) / diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index d6e509522fa..9325eccf667 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/tests/unit_tests/weightwindows/test.py @@ -370,6 +370,6 @@ def test_unstructured_mesh_applied_wws(request, run_in_tmpdir, library): model.settings.weight_windows = wws model.settings.weight_windows_on = True model.settings.run_mode = 'fixed source' - model.settings.particles = 100 + model.settings.particles = 1000 model.settings.batches = 2 model.run() From 06fe3366b8ba22a8f6c77c6602490a2500a112bb Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 12 Mar 2026 19:33:28 +0200 Subject: [PATCH 02/19] fix correct test --- tests/unit_tests/weightwindows/test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index 9325eccf667..c7acb48ef55 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 = 500 tally = openmc.Tally() tally.scores = ['heating'] @@ -370,6 +370,6 @@ def test_unstructured_mesh_applied_wws(request, run_in_tmpdir, library): model.settings.weight_windows = wws model.settings.weight_windows_on = True model.settings.run_mode = 'fixed source' - model.settings.particles = 1000 + model.settings.particles = 100 model.settings.batches = 2 model.run() From ff9d6612d8b18b97d64af49c363cc3337070072c Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 12 Mar 2026 19:50:06 +0200 Subject: [PATCH 03/19] update test --- .../photon_production_fission/results_true.dat | 16 ++++++++-------- tests/unit_tests/weightwindows/test.py | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/regression_tests/photon_production_fission/results_true.dat b/tests/regression_tests/photon_production_fission/results_true.dat index af325d4b02a..83b17ad5641 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.703498E+05 +9.673203E+09 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.711908E+05 -9.769096E+09 +1.703498E+05 +9.673203E+09 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.703498E+05 +9.673203E+09 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.711908E+05 -9.769096E+09 +1.703498E+05 +9.673203E+09 0.000000E+00 0.000000E+00 diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index c7acb48ef55..4096cb1ac5a 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 = 500 + model.settings.particles = 650 tally = openmc.Tally() tally.scores = ['heating'] From 7ea17cf1480681fc5b176dff6fc0ee948436b117 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 12 Mar 2026 22:23:23 +0200 Subject: [PATCH 04/19] various fixes --- src/photon.cpp | 50 ++++++++++++++------------ tests/unit_tests/weightwindows/test.py | 2 +- 2 files changed, 29 insertions(+), 23 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 92124533af6..2830ac27a7c 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -420,10 +420,27 @@ 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)); + // If in every angle we cannot eject an electron + // Exit with no shell + if (mu_max < -1) { + *i_shell = -1; + break; + } + while (true) { // Sample Klein-Nishina distribution for trial energy and angle std::tie(*alpha_out, *mu) = klein_nishina(alpha, seed); + 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 @@ -478,11 +495,8 @@ void PhotonInteraction::compton_doppler( 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; - } - + if (pz_max < 0.0) + continue; // Determine profile cdf value corresponding to p_z,max double c_max; if (pz_max > data::compton_profile_pz(n - 1)) { @@ -536,29 +550,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; diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index 4096cb1ac5a..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 = 650 + model.settings.particles = 1000 tally = openmc.Tally() tally.scores = ['heating'] From c1f93e3efbb9742c6678797aefafe29da93905da Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 12 Mar 2026 22:31:56 +0200 Subject: [PATCH 05/19] fix --- src/photon.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 2830ac27a7c..a75ecb46e77 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -424,17 +424,18 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler, double E_b = binding_energy_(last_shell); double E = alpha * MASS_ELECTRON_EV; double mu_max = 1 - E_b / (alpha * (E - E_b)); - // If in every angle we cannot eject an electron - // Exit with no shell - if (mu_max < -1) { - *i_shell = -1; - break; - } 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) From 86b5d6b5973e4a4f52ac364420f2db0306acc59d Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 12 Mar 2026 22:59:37 +0200 Subject: [PATCH 06/19] update regressions --- .../photon_production/results_true.dat | 56 +++++++++--------- .../results_true.dat | 58 +++++++++---------- .../photon_source/results_true.dat | 4 +- .../weightwindows/results_true.dat | 2 +- 4 files changed, 60 insertions(+), 60 deletions(-) diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index c4f5fafa252..27579b6d043 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.490000E-01 +9.006010E-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.349781E-01 +5.521471E-02 0.000000E+00 0.000000E+00 -2.350047E-01 -5.522722E-02 +2.349781E-01 +5.521471E-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.298000E-01 +5.280804E-02 +4.490316E+00 +2.016294E+01 0.000000E+00 0.000000E+00 -2.297000E-01 -5.276209E-02 -4.196651E+00 -1.761188E+01 +2.298000E-01 +5.280804E-02 +4.490316E+00 +2.016294E+01 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774484E+05 -3.148794E+10 +1.774569E+05 +3.149095E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774484E+05 -3.148794E+10 +1.774569E+05 +3.149095E+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.298000E-01 +5.280804E-02 +4.490316E+00 +2.016294E+01 0.000000E+00 0.000000E+00 -2.297000E-01 -5.276209E-02 -4.196651E+00 -1.761188E+01 +2.298000E-01 +5.280804E-02 +4.490316E+00 +2.016294E+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.774569E+05 +3.149095E+10 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 83b17ad5641..285156ff3b4 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.322402E+00 2.188209E-03 tally 1: -2.696393E+00 -2.423937E+00 +2.727353E+00 +2.479491E+00 0.000000E+00 0.000000E+00 -2.696393E+00 -2.423937E+00 +2.727353E+00 +2.479491E+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.706559E+00 +2.441851E+00 +4.342123E+08 +6.284729E+16 0.000000E+00 0.000000E+00 -2.672029E+00 -2.380251E+00 -4.288244E+08 -6.130569E+16 +2.706559E+00 +2.441851E+00 +4.342123E+08 +6.284729E+16 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.703498E+05 -9.673203E+09 +1.868708E+05 +1.164089E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.703498E+05 -9.673203E+09 +1.868708E+05 +1.164089E+10 0.000000E+00 0.000000E+00 tally 3: -2.649127E+00 -2.339294E+00 -4.288244E+08 -6.130569E+16 +2.652132E+00 +2.344612E+00 +4.342123E+08 +6.284729E+16 0.000000E+00 0.000000E+00 -2.649127E+00 -2.339294E+00 -4.288244E+08 -6.130569E+16 +2.652132E+00 +2.344612E+00 +4.342123E+08 +6.284729E+16 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.703498E+05 -9.673203E+09 +1.868708E+05 +1.164089E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.703498E+05 -9.673203E+09 +1.868708E+05 +1.164089E+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..a1811b5a198 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.262355E+02 +5.118249E+04 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/weightwindows/results_true.dat b/tests/regression_tests/weightwindows/results_true.dat index 897089e0597..907103744aa 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 +350897e7bf6f812f46be832aa93f06e759ca17da4b2e8bb2d04cccc9612b4fef15797d8990e542ffcfff6de4d75587d6a046017b7139d5a0066a8774099bbb15 \ No newline at end of file From 8aa5f4c39c5c5ce8b80c5f5395ac6b353dfa2fbe Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 12 Mar 2026 23:12:38 +0200 Subject: [PATCH 07/19] update atomic relaxation regression test --- .../atomic_relaxation/results_true.dat | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/regression_tests/atomic_relaxation/results_true.dat b/tests/regression_tests/atomic_relaxation/results_true.dat index 6f100dac8f8..328dbb44f4b 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.956668E+00 +3.828548E+00 +7.892384E+04 +6.228972E+09 0.000000E+00 0.000000E+00 -9.208123E+05 -8.478953E+11 +9.210762E+05 +8.483813E+11 From 075de4d50d8925697ad598afbd9f373adba7a0c2 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 13 Mar 2026 00:05:06 +0200 Subject: [PATCH 08/19] speed up code by avoiding rejections --- include/openmc/photon.h | 2 +- src/photon.cpp | 34 +++++++++++++++++++++------------- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 93c6dba53f5..a3813733295 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -87,7 +87,7 @@ class PhotonInteraction { tensor::Tensor profile_pdf_; tensor::Tensor profile_cdf_; tensor::Tensor binding_energy_; - tensor::Tensor electron_pdf_; + tensor::Tensor electron_cdf_; // Map subshells from Compton profile data obtained from Biggs et al, // "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic diff --git a/src/photon.cpp b/src/photon.cpp index a75ecb46e77..7f52f70f004 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -214,8 +214,12 @@ 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 electron_pdf; + read_dataset(rgroup, "num_electrons", electron_pdf); + electron_pdf /= electron_pdf.sum(); + 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 @@ -473,14 +477,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 : electron_cdf_) { + 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; } @@ -489,15 +504,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) { - 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) - continue; // Determine profile cdf value corresponding to p_z,max double c_max; if (pz_max > data::compton_profile_pz(n - 1)) { From c5a602973ccd18d794e300ade9633b48fbe5db40 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 13 Mar 2026 00:20:17 +0200 Subject: [PATCH 09/19] Revert "speed up code by avoiding rejections" This reverts commit 075de4d50d8925697ad598afbd9f373adba7a0c2. --- include/openmc/photon.h | 2 +- src/photon.cpp | 34 +++++++++++++--------------------- 2 files changed, 14 insertions(+), 22 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index a3813733295..93c6dba53f5 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -87,7 +87,7 @@ class PhotonInteraction { tensor::Tensor profile_pdf_; tensor::Tensor profile_cdf_; tensor::Tensor binding_energy_; - tensor::Tensor electron_cdf_; + tensor::Tensor electron_pdf_; // Map subshells from Compton profile data obtained from Biggs et al, // "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic diff --git a/src/photon.cpp b/src/photon.cpp index 7f52f70f004..a75ecb46e77 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -214,12 +214,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) rgroup = open_group(group, "compton_profiles"); // Read electron shell PDF and binding energies - tensor::Tensor electron_pdf; - read_dataset(rgroup, "num_electrons", electron_pdf); - electron_pdf /= electron_pdf.sum(); - 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, "num_electrons", electron_pdf_); + electron_pdf_ /= electron_pdf_.sum(); read_dataset(rgroup, "binding_energy", binding_energy_); // Read Compton profiles @@ -477,25 +473,14 @@ 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 : electron_cdf_) { - 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 = offset + prn(seed) * (1.0 - offset); - double c; - for (shell = j_shell; shell < electron_cdf_.size(); ++shell) { - c = electron_cdf_(shell); + double rn = prn(seed); + double c = 0.0; + for (shell = 0; shell < electron_pdf_.size(); ++shell) { + c += electron_pdf_(shell); if (rn < c) break; } @@ -504,8 +489,15 @@ void PhotonInteraction::compton_doppler( double E_b = binding_energy_(shell); // Determine p_z,max + double E = alpha * MASS_ELECTRON_EV; + 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) + continue; // Determine profile cdf value corresponding to p_z,max double c_max; if (pz_max > data::compton_profile_pz(n - 1)) { From 31e880f64f1736a461a9f2059cb406b62fc58925 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 13 Mar 2026 00:50:48 +0200 Subject: [PATCH 10/19] Reapply "speed up code by avoiding rejections" This reverts commit c5a602973ccd18d794e300ade9633b48fbe5db40. --- include/openmc/photon.h | 2 +- src/photon.cpp | 34 +++++++++++++++++++++------------- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 93c6dba53f5..a3813733295 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -87,7 +87,7 @@ class PhotonInteraction { tensor::Tensor profile_pdf_; tensor::Tensor profile_cdf_; tensor::Tensor binding_energy_; - tensor::Tensor electron_pdf_; + tensor::Tensor electron_cdf_; // Map subshells from Compton profile data obtained from Biggs et al, // "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic diff --git a/src/photon.cpp b/src/photon.cpp index a75ecb46e77..7f52f70f004 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -214,8 +214,12 @@ 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 electron_pdf; + read_dataset(rgroup, "num_electrons", electron_pdf); + electron_pdf /= electron_pdf.sum(); + 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 @@ -473,14 +477,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 : electron_cdf_) { + 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; } @@ -489,15 +504,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) { - 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) - continue; // Determine profile cdf value corresponding to p_z,max double c_max; if (pz_max > data::compton_profile_pz(n - 1)) { From d7f63092179c60f248ccd7a406f889a7e71fc3f0 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 13 Mar 2026 01:05:34 +0200 Subject: [PATCH 11/19] avoid some rejections --- src/photon.cpp | 1 + .../atomic_relaxation/results_true.dat | 12 ++-- .../photon_production/results_true.dat | 56 +++++++++--------- .../results_true.dat | 58 +++++++++---------- .../photon_source/results_true.dat | 4 +- .../weightwindows/results_true.dat | 2 +- .../survival_biasing/results_true.dat | 2 +- 7 files changed, 68 insertions(+), 67 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 7f52f70f004..fb58df39586 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -217,6 +217,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) tensor::Tensor electron_pdf; read_dataset(rgroup, "num_electrons", electron_pdf); electron_pdf /= electron_pdf.sum(); + electron_cdf_ = tensor::Tensor(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); diff --git a/tests/regression_tests/atomic_relaxation/results_true.dat b/tests/regression_tests/atomic_relaxation/results_true.dat index 328dbb44f4b..0714bb2e876 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.956668E+00 -3.828548E+00 -7.892384E+04 -6.228972E+09 +1.956914E+00 +3.829513E+00 +7.896189E+04 +6.234980E+09 0.000000E+00 0.000000E+00 -9.210762E+05 -8.483813E+11 +9.210381E+05 +8.483112E+11 diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 27579b6d043..6d5c896e0f8 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.490000E-01 -9.006010E-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.349781E-01 -5.521471E-02 +2.347443E-01 +5.510488E-02 0.000000E+00 0.000000E+00 -2.349781E-01 -5.521471E-02 +2.347443E-01 +5.510488E-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.298000E-01 -5.280804E-02 -4.490316E+00 -2.016294E+01 +2.297000E-01 +5.276209E-02 +4.492308E+00 +2.018083E+01 0.000000E+00 0.000000E+00 -2.298000E-01 -5.280804E-02 -4.490316E+00 -2.016294E+01 +2.297000E-01 +5.276209E-02 +4.492308E+00 +2.018083E+01 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774569E+05 -3.149095E+10 +1.774489E+05 +3.148811E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774569E+05 -3.149095E+10 +1.774489E+05 +3.148811E+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.298000E-01 -5.280804E-02 -4.490316E+00 -2.016294E+01 +2.297000E-01 +5.276209E-02 +4.492308E+00 +2.018083E+01 0.000000E+00 0.000000E+00 -2.298000E-01 -5.280804E-02 -4.490316E+00 -2.016294E+01 +2.297000E-01 +5.276209E-02 +4.492308E+00 +2.018083E+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.774569E+05 -3.149095E+10 +1.774489E+05 +3.148811E+10 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 285156ff3b4..0ffa0b98650 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.322402E+00 2.188209E-03 +2.304258E+00 1.564323E-02 tally 1: -2.727353E+00 -2.479491E+00 +2.704024E+00 +2.437522E+00 0.000000E+00 0.000000E+00 -2.727353E+00 -2.479491E+00 +2.704024E+00 +2.437522E+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.706559E+00 -2.441851E+00 -4.342123E+08 -6.284729E+16 +2.670940E+00 +2.378429E+00 +4.288120E+08 +6.130738E+16 0.000000E+00 0.000000E+00 -2.706559E+00 -2.441851E+00 -4.342123E+08 -6.284729E+16 +2.670940E+00 +2.378429E+00 +4.288120E+08 +6.130738E+16 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.868708E+05 -1.164089E+10 +1.837955E+05 +1.126059E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.868708E+05 -1.164089E+10 +1.837955E+05 +1.126059E+10 0.000000E+00 0.000000E+00 tally 3: -2.652132E+00 -2.344612E+00 -4.342123E+08 -6.284729E+16 +2.653118E+00 +2.346351E+00 +4.288120E+08 +6.130738E+16 0.000000E+00 0.000000E+00 -2.652132E+00 -2.344612E+00 -4.342123E+08 -6.284729E+16 +2.653118E+00 +2.346351E+00 +4.288120E+08 +6.130738E+16 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.868708E+05 -1.164089E+10 +1.837955E+05 +1.126059E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.868708E+05 -1.164089E+10 +1.837955E+05 +1.126059E+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 a1811b5a198..e94a9044ff0 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.262355E+02 -5.118249E+04 +2.260811E+02 +5.111264E+04 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/weightwindows/results_true.dat b/tests/regression_tests/weightwindows/results_true.dat index 907103744aa..558f3953f9e 100644 --- a/tests/regression_tests/weightwindows/results_true.dat +++ b/tests/regression_tests/weightwindows/results_true.dat @@ -1 +1 @@ -350897e7bf6f812f46be832aa93f06e759ca17da4b2e8bb2d04cccc9612b4fef15797d8990e542ffcfff6de4d75587d6a046017b7139d5a0066a8774099bbb15 \ No newline at end of file +a5923974afc7ba2fe7cf1a967676f4cf1434c3e47373dbae9f8d40ec82b0d5bd76c25a39e5c9d0dc9677c7a9625af42951d0570ca39249483681f9fbf30f10e9 \ No newline at end of file diff --git a/tests/regression_tests/weightwindows/survival_biasing/results_true.dat b/tests/regression_tests/weightwindows/survival_biasing/results_true.dat index b457a245c7c..9547ef49855 100644 --- a/tests/regression_tests/weightwindows/survival_biasing/results_true.dat +++ b/tests/regression_tests/weightwindows/survival_biasing/results_true.dat @@ -1 +1 @@ -386e507008ed3c72c6e1a101aafc01cfaaff2c0b10555558d1eab6332441f7b17264b55002d721151bac2f3e7c1a8aac4c5ed243f5270533d171850f985206ed \ No newline at end of file +ef6e1f56dd3c4949b0e7308d46e46aea66d6ab203eef230f9559027dd1be5a713166ee5eca5d0bca0344d54d9fc664d757b69bf3109a6b17d7c9d401a86fdf0f \ No newline at end of file From 4a7cf7e21b9cc1102bbc8cbfc934dc11cecdf902 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 13 Mar 2026 01:22:47 +0200 Subject: [PATCH 12/19] fix typo --- src/photon.cpp | 2 +- .../atomic_relaxation/results_true.dat | 12 ++-- .../photon_production/results_true.dat | 56 +++++++++--------- .../results_true.dat | 58 +++++++++---------- .../photon_source/results_true.dat | 4 +- .../weightwindows/results_true.dat | 2 +- 6 files changed, 67 insertions(+), 67 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index fb58df39586..9bff2dec9e3 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -480,7 +480,7 @@ void PhotonInteraction::compton_doppler( auto n = data::compton_profile_pz.size(); double E = alpha * MASS_ELECTRON_EV; int j_shell = 0; - for (double E_b : electron_cdf_) { + for (double E_b : binding_energy_) { if ((E_b - (E - E_b) * alpha * (1.0 - mu)) < 0.0) break; ++j_shell; diff --git a/tests/regression_tests/atomic_relaxation/results_true.dat b/tests/regression_tests/atomic_relaxation/results_true.dat index 0714bb2e876..4ac3bb63b41 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.956914E+00 -3.829513E+00 -7.896189E+04 -6.234980E+09 +1.955720E+00 +3.824843E+00 +7.895341E+04 +6.233641E+09 0.000000E+00 0.000000E+00 -9.210381E+05 -8.483112E+11 +9.210466E+05 +8.483268E+11 diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 6d5c896e0f8..7ebb2ac9b20 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.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.347443E-01 -5.510488E-02 +2.342279E-01 +5.486270E-02 0.000000E+00 0.000000E+00 -2.347443E-01 -5.510488E-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.492308E+00 -2.018083E+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.492308E+00 -2.018083E+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.774489E+05 -3.148811E+10 +1.774550E+05 +3.149027E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774489E+05 -3.148811E+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.492308E+00 -2.018083E+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.492308E+00 -2.018083E+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.774489E+05 -3.148811E+10 +1.774550E+05 +3.149027E+10 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 0ffa0b98650..f6e0ec25eca 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.304258E+00 1.564323E-02 +2.331233E+00 3.816178E-02 tally 1: -2.704024E+00 -2.437522E+00 +2.733254E+00 +2.491835E+00 0.000000E+00 0.000000E+00 -2.704024E+00 -2.437522E+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.670940E+00 -2.378429E+00 -4.288120E+08 -6.130738E+16 +2.705904E+00 +2.441602E+00 +4.339210E+08 +6.278920E+16 0.000000E+00 0.000000E+00 -2.670940E+00 -2.378429E+00 -4.288120E+08 -6.130738E+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.837955E+05 -1.126059E+10 +1.878973E+05 +1.177074E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.837955E+05 -1.126059E+10 +1.878973E+05 +1.177074E+10 0.000000E+00 0.000000E+00 tally 3: -2.653118E+00 -2.346351E+00 -4.288120E+08 -6.130738E+16 +2.658522E+00 +2.355935E+00 +4.339210E+08 +6.278920E+16 0.000000E+00 0.000000E+00 -2.653118E+00 -2.346351E+00 -4.288120E+08 -6.130738E+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.837955E+05 -1.126059E+10 +1.878973E+05 +1.177074E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.837955E+05 -1.126059E+10 +1.878973E+05 +1.177074E+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 e94a9044ff0..d91fc52400b 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.260811E+02 -5.111264E+04 +2.263259E+02 +5.122343E+04 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/weightwindows/results_true.dat b/tests/regression_tests/weightwindows/results_true.dat index 558f3953f9e..de23fe4f8b2 100644 --- a/tests/regression_tests/weightwindows/results_true.dat +++ b/tests/regression_tests/weightwindows/results_true.dat @@ -1 +1 @@ -a5923974afc7ba2fe7cf1a967676f4cf1434c3e47373dbae9f8d40ec82b0d5bd76c25a39e5c9d0dc9677c7a9625af42951d0570ca39249483681f9fbf30f10e9 \ No newline at end of file +b51d86a9d3a00713fb6d95690adca37d6af93f08884a233caa88b8a21e2a3985c85b6b4d150dd74b19ccc0d4c3fa973d2114309096642b2ea7e38d78ee1d8186 \ No newline at end of file From db77037d0ca2a0f6c9ac59883d9f479a6e17676f Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 13 Mar 2026 14:47:47 +0200 Subject: [PATCH 13/19] update regression tests --- .../photon_production_fission/results_true.dat | 16 ++++++++-------- .../survival_biasing/results_true.dat | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/regression_tests/photon_production_fission/results_true.dat b/tests/regression_tests/photon_production_fission/results_true.dat index f6e0ec25eca..4a84d9191d0 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.878973E+05 -1.177074E+10 +1.879204E+05 +1.177358E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.878973E+05 -1.177074E+10 +1.879204E+05 +1.177358E+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.878973E+05 -1.177074E+10 +1.879204E+05 +1.177358E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.878973E+05 -1.177074E+10 +1.879204E+05 +1.177358E+10 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/weightwindows/survival_biasing/results_true.dat b/tests/regression_tests/weightwindows/survival_biasing/results_true.dat index 9547ef49855..b457a245c7c 100644 --- a/tests/regression_tests/weightwindows/survival_biasing/results_true.dat +++ b/tests/regression_tests/weightwindows/survival_biasing/results_true.dat @@ -1 +1 @@ -ef6e1f56dd3c4949b0e7308d46e46aea66d6ab203eef230f9559027dd1be5a713166ee5eca5d0bca0344d54d9fc664d757b69bf3109a6b17d7c9d401a86fdf0f \ No newline at end of file +386e507008ed3c72c6e1a101aafc01cfaaff2c0b10555558d1eab6332441f7b17264b55002d721151bac2f3e7c1a8aac4c5ed243f5270533d171850f985206ed \ No newline at end of file From d943dca1b057e47be68c318549d398e88c757746 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 1 Apr 2026 16:14:58 +0300 Subject: [PATCH 14/19] fix --- include/openmc/photon.h | 2 +- src/photon.cpp | 85 ++++++++----------- .../atomic_relaxation/results_true.dat | 12 +-- .../photon_production/results_true.dat | 68 +++++++-------- .../results_true.dat | 58 ++++++------- .../photon_source/results_true.dat | 4 +- .../weightwindows/results_true.dat | 2 +- 7 files changed, 108 insertions(+), 123 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index a3813733295..93c6dba53f5 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -87,7 +87,7 @@ class PhotonInteraction { tensor::Tensor profile_pdf_; tensor::Tensor profile_cdf_; tensor::Tensor binding_energy_; - tensor::Tensor electron_cdf_; + tensor::Tensor electron_pdf_; // Map subshells from Compton profile data obtained from Biggs et al, // "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic diff --git a/src/photon.cpp b/src/photon.cpp index 9bff2dec9e3..daa0c6bc087 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -214,13 +214,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) rgroup = open_group(group, "compton_profiles"); // Read electron shell PDF and binding energies - tensor::Tensor electron_pdf; - read_dataset(rgroup, "num_electrons", electron_pdf); - electron_pdf /= electron_pdf.sum(); - electron_cdf_ = tensor::Tensor(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, "num_electrons", electron_pdf_); + electron_pdf_ /= electron_pdf_.sum(); read_dataset(rgroup, "binding_energy", binding_energy_); // Read Compton profiles @@ -428,25 +423,13 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler, 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)); + if (E < E_b) + fatal_error("Cannot eject any electron"); 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 @@ -478,42 +461,37 @@ 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 = offset + prn(seed) * (1.0 - offset); - double c; - for (shell = j_shell; shell < electron_cdf_.size(); ++shell) { - c = electron_cdf_(shell); + double rn = prn(seed); + double c = 0.0; + for (shell = 0; shell < electron_pdf_.size(); ++shell) { + c += electron_pdf_(shell); 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 + // 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); + // 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)) { + // TODO: handle linear extrapolation in lin-log scales 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); @@ -522,10 +500,11 @@ 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); } } @@ -561,22 +540,28 @@ void PhotonInteraction::compton_doppler( double quad = b * b - 4.0 * a * c; if (quad < 0) - continue; + fatal_error("Cannot doppler broaden."); 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) && (E_out2 > 0.0)) { - *E_out = prn(seed) < 0.5 ? E_out1 : E_out2; + 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 { - *E_out = std::max(E_out1, E_out2); + if (E_out2 > 0.0) { + *E_out = E_out2; + } else { + // No positive solution -- resample + continue; + } } - if (*E_out < E - E_b) + 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 4ac3bb63b41..72d2fc5db14 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.955720E+00 -3.824843E+00 -7.895341E+04 -6.233641E+09 +1.967531E+00 +3.871180E+00 +7.966753E+04 +6.346915E+09 0.000000E+00 0.000000E+00 -9.210466E+05 -8.483268E+11 +9.203325E+05 +8.470119E+11 diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 7ebb2ac9b20..b78a9edc25f 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.489000E-01 -9.004112E-01 +9.520000E-01 +9.063040E-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.342279E-01 -5.486270E-02 +2.337965E-01 +5.466081E-02 0.000000E+00 0.000000E+00 -2.342279E-01 -5.486270E-02 +2.337965E-01 +5.466081E-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.294000E-01 -5.262436E-02 -4.488894E+00 -2.015017E+01 +2.253000E-01 +5.076009E-02 +4.053867E+00 +1.643384E+01 0.000000E+00 0.000000E+00 -2.294000E-01 -5.262436E-02 -4.488894E+00 -2.015017E+01 +2.253000E-01 +5.076009E-02 +4.053867E+00 +1.643384E+01 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774550E+05 -3.149027E+10 +1.762738E+05 +3.107245E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.774550E+05 -3.149027E+10 +1.762738E+05 +3.107245E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.474427E+04 -2.173936E+08 +1.438891E+04 +2.070408E+08 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.474427E+04 -2.173936E+08 +1.438891E+04 +2.070408E+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.294000E-01 -5.262436E-02 -4.488894E+00 -2.015017E+01 +2.253000E-01 +5.076009E-02 +4.053867E+00 +1.643384E+01 0.000000E+00 0.000000E+00 -2.294000E-01 -5.262436E-02 -4.488894E+00 -2.015017E+01 +2.253000E-01 +5.076009E-02 +4.053867E+00 +1.643384E+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.774550E+05 -3.149027E+10 +1.762738E+05 +3.107245E+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.438891E+04 +2.070408E+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 4a84d9191d0..d643855c08a 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.331233E+00 3.816178E-02 +2.313118E+00 2.756443E-02 tally 1: -2.733254E+00 -2.491835E+00 +2.710112E+00 +2.448973E+00 0.000000E+00 0.000000E+00 -2.733254E+00 -2.491835E+00 +2.710112E+00 +2.448973E+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.705904E+00 -2.441602E+00 -4.339210E+08 -6.278920E+16 +2.689726E+00 +2.411702E+00 +4.310507E+08 +6.193986E+16 0.000000E+00 0.000000E+00 -2.705904E+00 -2.441602E+00 -4.339210E+08 -6.278920E+16 +2.689726E+00 +2.411702E+00 +4.310507E+08 +6.193986E+16 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.879204E+05 -1.177358E+10 +1.698053E+05 +9.611469E+09 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.879204E+05 -1.177358E+10 +1.698053E+05 +9.611469E+09 0.000000E+00 0.000000E+00 tally 3: -2.658522E+00 -2.355935E+00 -4.339210E+08 -6.278920E+16 +2.658739E+00 +2.356315E+00 +4.310507E+08 +6.193986E+16 0.000000E+00 0.000000E+00 -2.658522E+00 -2.355935E+00 -4.339210E+08 -6.278920E+16 +2.658739E+00 +2.356315E+00 +4.310507E+08 +6.193986E+16 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.879204E+05 -1.177358E+10 +1.698053E+05 +9.611469E+09 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.879204E+05 -1.177358E+10 +1.698053E+05 +9.611469E+09 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 d91fc52400b..6701761d38b 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.263259E+02 -5.122343E+04 +2.275725E+02 +5.178926E+04 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/weightwindows/results_true.dat b/tests/regression_tests/weightwindows/results_true.dat index de23fe4f8b2..626b5e6c070 100644 --- a/tests/regression_tests/weightwindows/results_true.dat +++ b/tests/regression_tests/weightwindows/results_true.dat @@ -1 +1 @@ -b51d86a9d3a00713fb6d95690adca37d6af93f08884a233caa88b8a21e2a3985c85b6b4d150dd74b19ccc0d4c3fa973d2114309096642b2ea7e38d78ee1d8186 \ No newline at end of file +c4d56d4eb3e4d10f9238d34f909f950227a99b262c04ae0afc78f9c2e0a2c07139b88461c4b4c62ed24803c59419183fecc405b6155f553c00be454d632144a2 \ No newline at end of file From 7a7d412bb9fecea16901fbf508918de2436fea06 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 1 Apr 2026 16:28:49 +0300 Subject: [PATCH 15/19] fix2 --- src/photon.cpp | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index daa0c6bc087..7912403067a 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -420,12 +420,6 @@ 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; - if (E < E_b) - fatal_error("Cannot eject any electron"); - while (true) { // Sample Klein-Nishina distribution for trial energy and angle std::tie(*alpha_out, *mu) = klein_nishina(alpha, seed); @@ -487,7 +481,6 @@ void PhotonInteraction::compton_doppler( // Determine profile cdf value corresponding to p_z,max double c_max; if (std::abs(pz_max) > data::compton_profile_pz(n - 1)) { - // TODO: handle linear extrapolation in lin-log scales c_max = profile_cdf_(shell, n - 1); } else { int i = lower_bound_index(data::compton_profile_pz.cbegin(), @@ -538,10 +531,7 @@ void PhotonInteraction::compton_doppler( double b = 2.0 * E * (f - momentum_sq * mu); c = E * E * (momentum_sq - 1.0); - double quad = b * b - 4.0 * a * c; - if (quad < 0) - fatal_error("Cannot doppler broaden."); - quad = std::sqrt(quad); + double quad = std::sqrt(b * b - 4.0 * a * c); double E_out1 = -(b + quad) / (2.0 * a); double E_out2 = -(b - quad) / (2.0 * a); From c0713fba0e0278715c68460dbaef666e72cab1b5 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 1 Apr 2026 19:36:20 +0300 Subject: [PATCH 16/19] fix3 --- src/photon.cpp | 45 ++++++------ .../atomic_relaxation/results_true.dat | 12 ++-- .../photon_production/results_true.dat | 68 +++++++++---------- .../results_true.dat | 58 ++++++++-------- .../photon_source/results_true.dat | 4 +- .../weightwindows/results_true.dat | 2 +- 6 files changed, 94 insertions(+), 95 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 7912403067a..7f603303016 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -503,7 +503,11 @@ void PhotonInteraction::compton_doppler( } // 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); @@ -528,29 +532,24 @@ 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 quad = std::sqrt(b * b - 4.0 * a * c); - 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; - } + double b = 2.0 * (f - momentum_sq * mu); + c = (momentum_sq - 1.0); + double quad = b * b - 4.0 * a * c; + if (quad < 0) { + write_message("Negative quad {}", quad); + continue; } + quad = std::sqrt(quad); + 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 72d2fc5db14..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.967531E+00 -3.871180E+00 -7.966753E+04 -6.346915E+09 +1.920387E+00 +3.687888E+00 +7.854617E+04 +6.169501E+09 0.000000E+00 0.000000E+00 -9.203325E+05 -8.470119E+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 b78a9edc25f..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.520000E-01 -9.063040E-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.337965E-01 -5.466081E-02 +2.362831E-01 +5.582969E-02 0.000000E+00 0.000000E+00 -2.337965E-01 -5.466081E-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.253000E-01 -5.076009E-02 -4.053867E+00 -1.643384E+01 +2.344000E-01 +5.494336E-02 +4.593008E+00 +2.109572E+01 0.000000E+00 0.000000E+00 -2.253000E-01 -5.076009E-02 -4.053867E+00 -1.643384E+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.762738E+05 -3.107245E+10 +1.792949E+05 +3.214667E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.762738E+05 -3.107245E+10 +1.792949E+05 +3.214667E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.438891E+04 -2.070408E+08 +1.452085E+04 +2.108550E+08 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.438891E+04 -2.070408E+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.253000E-01 -5.076009E-02 -4.053867E+00 -1.643384E+01 +2.344000E-01 +5.494336E-02 +4.593008E+00 +2.109572E+01 0.000000E+00 0.000000E+00 -2.253000E-01 -5.076009E-02 -4.053867E+00 -1.643384E+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.762738E+05 -3.107245E+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.438891E+04 -2.070408E+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 d643855c08a..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.313118E+00 2.756443E-02 +2.255088E+00 7.363219E-02 tally 1: -2.710112E+00 -2.448973E+00 +2.640439E+00 +2.329186E+00 0.000000E+00 0.000000E+00 -2.710112E+00 -2.448973E+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.689726E+00 -2.411702E+00 -4.310507E+08 -6.193986E+16 +2.625521E+00 +2.303289E+00 +4.210498E+08 +5.923082E+16 0.000000E+00 0.000000E+00 -2.689726E+00 -2.411702E+00 -4.310507E+08 -6.193986E+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.698053E+05 -9.611469E+09 +1.826906E+05 +1.114925E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.698053E+05 -9.611469E+09 +1.826906E+05 +1.114925E+10 0.000000E+00 0.000000E+00 tally 3: -2.658739E+00 -2.356315E+00 -4.310507E+08 -6.193986E+16 +2.658439E+00 +2.355784E+00 +4.210498E+08 +5.923082E+16 0.000000E+00 0.000000E+00 -2.658739E+00 -2.356315E+00 -4.310507E+08 -6.193986E+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.698053E+05 -9.611469E+09 +1.826906E+05 +1.114925E+10 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.698053E+05 -9.611469E+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 6701761d38b..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.275725E+02 -5.178926E+04 +2.193647E+02 +4.812086E+04 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/weightwindows/results_true.dat b/tests/regression_tests/weightwindows/results_true.dat index 626b5e6c070..a3b10b49b58 100644 --- a/tests/regression_tests/weightwindows/results_true.dat +++ b/tests/regression_tests/weightwindows/results_true.dat @@ -1 +1 @@ -c4d56d4eb3e4d10f9238d34f909f950227a99b262c04ae0afc78f9c2e0a2c07139b88461c4b4c62ed24803c59419183fecc405b6155f553c00be454d632144a2 \ No newline at end of file +bb3b76ee105f20bc9eb763fd04753f6212ab0139197f5be03b927769ac4f45efa77dd22cdf48900062c3eaa8874f463bef11c0dd52321ba8f6da2c12c98247b7 \ No newline at end of file From 9b83c02dfb08c747d2ba317e8e7bfcbf0dcb390f Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 1 Apr 2026 19:42:34 +0300 Subject: [PATCH 17/19] fix4 --- src/photon.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 7f603303016..9b40522c271 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -535,10 +535,8 @@ void PhotonInteraction::compton_doppler( double b = 2.0 * (f - momentum_sq * mu); c = (momentum_sq - 1.0); double quad = b * b - 4.0 * a * c; - if (quad < 0) { - write_message("Negative quad {}", quad); + if (quad < 0) continue; - } quad = std::sqrt(quad); double E_out1 = -(b + quad) / (2.0 * a) * E; double E_out2 = -(b - quad) / (2.0 * a) * E; From 0310503c1af87f388701e61eab7607641dec146f Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 1 Apr 2026 21:47:02 +0300 Subject: [PATCH 18/19] update regression --- .../pulse_height/results_true.dat | 102 +++++++++--------- 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/tests/regression_tests/pulse_height/results_true.dat b/tests/regression_tests/pulse_height/results_true.dat index c57e8ff1c83..6aa7df9923a 100644 --- a/tests/regression_tests/pulse_height/results_true.dat +++ b/tests/regression_tests/pulse_height/results_true.dat @@ -1,18 +1,18 @@ 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 +4.120000E+00 +3.409000E+00 +3.000000E-02 +5.000000E-04 0.000000E+00 0.000000E+00 +1.000000E-02 +1.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 @@ -23,26 +23,36 @@ tally 1: 1.000000E-04 0.000000E+00 0.000000E+00 +0.000000E+00 +0.000000E+00 +4.000000E-02 +1.000000E-03 1.000000E-02 1.000000E-04 -3.000000E-02 -5.000000E-04 -2.000000E-02 -4.000000E-04 2.000000E-02 4.000000E-04 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 +3.000000E-02 +3.000000E-04 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 3.000000E-02 -3.000000E-04 +5.000000E-04 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 1.000000E-02 @@ -53,8 +63,6 @@ tally 1: 0.000000E+00 1.000000E-02 1.000000E-04 -1.000000E-02 -1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -63,112 +71,104 @@ tally 1: 0.000000E+00 1.000000E-02 1.000000E-04 +2.000000E-02 +2.000000E-04 0.000000E+00 0.000000E+00 +4.000000E-02 +6.000000E-04 0.000000E+00 0.000000E+00 +2.000000E-02 +2.000000E-04 0.000000E+00 0.000000E+00 -4.000000E-02 -6.000000E-04 +1.000000E-02 +1.000000E-04 1.000000E-02 1.000000E-04 0.000000E+00 0.000000E+00 +3.000000E-02 +5.000000E-04 2.000000E-02 2.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 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 0.000000E+00 0.000000E+00 +3.000000E-02 +3.000000E-04 1.000000E-02 1.000000E-04 2.000000E-02 4.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 -2.000000E-02 -2.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 +0.000000E+00 +0.000000E+00 1.000000E-02 1.000000E-04 3.000000E-02 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 +0.000000E+00 +0.000000E+00 2.000000E-02 2.000000E-04 1.000000E-02 1.000000E-04 1.000000E-02 1.000000E-04 -3.000000E-02 -5.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 -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 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 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 +2.100000E-01 +1.030000E-02 From 9165cc7176cc2b691583189408d71be39453633e Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 1 Apr 2026 23:37:24 +0300 Subject: [PATCH 19/19] update regression test --- .../pulse_height/inputs_true.dat | 4 +- .../pulse_height/results_true.dat | 348 +++++++++--------- tests/regression_tests/pulse_height/test.py | 6 +- 3 files changed, 179 insertions(+), 179 deletions(-) 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 6aa7df9923a..84fcfc0f2de 100644 --- a/tests/regression_tests/pulse_height/results_true.dat +++ b/tests/regression_tests/pulse_height/results_true.dat @@ -1,174 +1,186 @@ tally 1: -4.120000E+00 -3.409000E+00 -3.000000E-02 -5.000000E-04 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.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 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -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-02 -1.000000E-04 -2.000000E-02 -4.000000E-04 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -3.000000E-02 -3.000000E-04 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -3.000000E-02 -5.000000E-04 -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 -1.000000E-02 -1.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 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -2.000000E-02 -2.000000E-04 -0.000000E+00 -0.000000E+00 -4.000000E-02 +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 -0.000000E+00 -0.000000E+00 -2.000000E-02 -2.000000E-04 -0.000000E+00 -0.000000E+00 -1.000000E-02 -1.000000E-04 -1.000000E-02 -1.000000E-04 -0.000000E+00 -0.000000E+00 -3.000000E-02 +3.600000E-07 5.000000E-04 -2.000000E-02 -2.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 -3.000000E-02 +2.500000E-07 +8.000000E-04 +6.400000E-07 3.000000E-04 -1.000000E-02 -1.000000E-04 -2.000000E-02 -4.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 -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 -3.000000E-02 -5.000000E-04 -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 -0.000000E+00 -0.000000E+00 -2.000000E-02 +9.000000E-08 +3.000000E-04 +9.000000E-08 2.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 -0.000000E+00 +4.000000E-08 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 -1.000000E-02 1.000000E-04 +1.000000E-08 0.000000E+00 0.000000E+00 0.000000E+00 @@ -185,17 +197,5 @@ tally 1: 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 -2.100000E-01 -1.030000E-02 +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()