From 922e29e29f890389b9d59e5f3d577b48a6639171 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Mon, 3 Jun 2024 10:15:41 -0500 Subject: [PATCH 01/17] Added Shannon entropy functionality to random ray and associated test --- examples/pincell_multigroup/build_xml.py | 2 +- examples/pincell_random_ray/build_xml.py | 8 + src/random_ray/flat_source_domain.cpp | 27 + src/settings.cpp | 41 +- src/simulation.cpp | 2 +- .../random_ray_entropy/__init__.py | 0 .../random_ray_entropy/geometry.xml | 542 ++++++++++++++++++ .../random_ray_entropy/materials.xml | 8 + .../random_ray_entropy/mgxs.h5 | Bin 0 -> 10664 bytes .../random_ray_entropy/results_true.dat | 13 + .../random_ray_entropy/settings.xml | 17 + .../random_ray_entropy/test.py | 33 ++ 12 files changed, 672 insertions(+), 21 deletions(-) create mode 100644 tests/regression_tests/random_ray_entropy/__init__.py create mode 100644 tests/regression_tests/random_ray_entropy/geometry.xml create mode 100644 tests/regression_tests/random_ray_entropy/materials.xml create mode 100644 tests/regression_tests/random_ray_entropy/mgxs.h5 create mode 100644 tests/regression_tests/random_ray_entropy/results_true.dat create mode 100644 tests/regression_tests/random_ray_entropy/settings.xml create mode 100644 tests/regression_tests/random_ray_entropy/test.py diff --git a/examples/pincell_multigroup/build_xml.py b/examples/pincell_multigroup/build_xml.py index 019ae6a113f..ca128a25310 100644 --- a/examples/pincell_multigroup/build_xml.py +++ b/examples/pincell_multigroup/build_xml.py @@ -90,7 +90,7 @@ # Create a region represented as the inside of a rectangular prism pitch = 1.26 -box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective') +box = openmc.model.RectangularPrism(pitch/2, pitch/2, boundary_type='reflective') # Instantiate Cells fuel = openmc.Cell(fill=uo2, region=-fuel_or, name='fuel') diff --git a/examples/pincell_random_ray/build_xml.py b/examples/pincell_random_ray/build_xml.py index b3dd8020a51..9a5fa4f423d 100644 --- a/examples/pincell_random_ray/build_xml.py +++ b/examples/pincell_random_ray/build_xml.py @@ -161,6 +161,14 @@ settings.random_ray['distance_inactive'] = 40.0 settings.random_ray['distance_active'] = 400.0 +entropy_mesh = openmc.RegularMesh() +entropy_mesh.lower_left = (-pitch/2, -pitch/2, -pitch/2) +entropy_mesh.upper_right = (pitch/2, pitch/2, pitch/2) +entropy_mesh.dimension = (8, 8, 8) +settings.entropy_mesh = entropy_mesh + +settings.entropy_on = True + settings.export_to_xml() ############################################################################### diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 5e1194aa92a..2fbe0518e50 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1,6 +1,7 @@ #include "openmc/random_ray/flat_source_domain.h" #include "openmc/cell.h" +#include "openmc/eigenvalue.h" #include "openmc/geometry.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" @@ -225,6 +226,9 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const const int t = 0; const int a = 0; + // Vector for gathering fission source terms for Shannon entropy calculation + vector p(n_source_regions_, 0.0f); + #pragma omp parallel for reduction(+ : fission_rate_old, fission_rate_new) for (int sr = 0; sr < n_source_regions_; sr++) { @@ -249,10 +253,33 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const fission_rate_old += sr_fission_source_old * volume; fission_rate_new += sr_fission_source_new * volume; + + // Assigning the fission source to the vector for entropy calculation. + p[sr] = sr_fission_source_new; } double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); + // Calculating sum for entropy normalization. + double sum = 0.0; + for (float num : p) { + sum += num; + } + + // Normalize to total weight of bank sites. + for (float& num : p) { + num /= sum; + } + + // Sum values to obtain Shannon entropy. + double H = 0.0; + for (int i = 0; i < n_source_regions_; i++) { + H -= p[i] * std::log(p[i]) / std::log(2.0); + } + + // Adds entropy value to shared entropy vector in openmc namespace. + simulation::entropy.push_back(H); + return k_eff_new; } diff --git a/src/settings.cpp b/src/settings.cpp index 6790f2cc0f7..539abad10a9 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -624,29 +624,32 @@ void read_settings_xml(pugi::xml_node root) } // Shannon Entropy mesh - if (check_for_node(root, "entropy_mesh")) { - int temp = std::stoi(get_node_value(root, "entropy_mesh")); - if (model::mesh_map.find(temp) == model::mesh_map.end()) { - fatal_error(fmt::format( - "Mesh {} specified for Shannon entropy does not exist.", temp)); - } + if (solver_type == SolverType::RANDOM_RAY) { + entropy_on = true; + } else if (solver_type == SolverType::MONTE_CARLO) { + if (check_for_node(root, "entropy_mesh")) { + int temp = std::stoi(get_node_value(root, "entropy_mesh")); + if (model::mesh_map.find(temp) == model::mesh_map.end()) { + fatal_error(fmt::format( + "Mesh {} specified for Shannon entropy does not exist.", temp)); + } - auto* m = - dynamic_cast(model::meshes[model::mesh_map.at(temp)].get()); - if (!m) - fatal_error("Only regular meshes can be used as an entropy mesh"); - simulation::entropy_mesh = m; + auto* m = + dynamic_cast(model::meshes[model::mesh_map.at(temp)].get()); + if (!m) + fatal_error("Only regular meshes can be used as an entropy mesh"); + simulation::entropy_mesh = m; - // Turn on Shannon entropy calculation - entropy_on = true; + // Turn on Shannon entropy calculation + entropy_on = true; - } else if (check_for_node(root, "entropy")) { - fatal_error( - "Specifying a Shannon entropy mesh via the element " - "is deprecated. Please create a mesh using and then reference " - "it by specifying its ID in an element."); + } else if (check_for_node(root, "entropy")) { + fatal_error( + "Specifying a Shannon entropy mesh via the element " + "is deprecated. Please create a mesh using and then reference " + "it by specifying its ID in an element."); + } } - // Uniform fission source weighting mesh if (check_for_node(root, "ufs_mesh")) { auto temp = std::stoi(get_node_value(root, "ufs_mesh")); diff --git a/src/simulation.cpp b/src/simulation.cpp index 527d98db4f1..eb18b9c3e5d 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -521,7 +521,7 @@ void finalize_generation() if (settings::run_mode == RunMode::EIGENVALUE) { // Calculate shannon entropy - if (settings::entropy_on) + if (settings::entropy_on && settings::solver_type == SolverType::MONTE_CARLO) shannon_entropy(); // Collect results and statistics diff --git a/tests/regression_tests/random_ray_entropy/__init__.py b/tests/regression_tests/random_ray_entropy/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_entropy/geometry.xml b/tests/regression_tests/random_ray_entropy/geometry.xml new file mode 100644 index 00000000000..9a0e6640193 --- /dev/null +++ b/tests/regression_tests/random_ray_entropy/geometry.xml @@ -0,0 +1,542 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/regression_tests/random_ray_entropy/materials.xml b/tests/regression_tests/random_ray_entropy/materials.xml new file mode 100644 index 00000000000..5a6f93414bb --- /dev/null +++ b/tests/regression_tests/random_ray_entropy/materials.xml @@ -0,0 +1,8 @@ + + + mgxs.h5 + + + + + diff --git a/tests/regression_tests/random_ray_entropy/mgxs.h5 b/tests/regression_tests/random_ray_entropy/mgxs.h5 new file mode 100644 index 0000000000000000000000000000000000000000..6ce80bf713a8d5ff29700ba59e982e18f99b5403 GIT binary patch literal 10664 zcmeHN%}*0S6o1<)ELx=~2P1yeL;JB5M`%OA4UJMuu{yU;m}}jFYub= zpCp>Wc%w?`qzL4JLLk-KJs<`=jH9M_TERFYFm72kyJxt}HjSLnQCudeLWQygKO6SV(mu^{#vDtpJMgnd zH{3=&5sl;2%Eu!cjybM9GwIIQ3|eA|^>K9`h0dQi?=_$Ct!2L!-+%vhQ~3XQ!s+(o zk0k{bjN?Gc@yAL;^%%{l`++9`bj_V4+@MK6%Q@D1c^N!HbJV{c83gv_Q{XH7@ zLUFe!NP%DEcLnpS4N4!}|Be9huotfg$Yp`!8G*|jVO6cwN1#3`^AJf<+$ z;db}k2Yj;uLxT^6kOlNSNdC@?g5Pacu$ja1uHg#02Ep%5%2c{A9m}hQ&q)QJ zb5gOo1FX=yIu|vppSS1FZtb80ix}W%XQ$2Y2?A0_%E1fc90x*l}8~J3RX;`i{jvLfPglKxS;(OKWW1?go4)haasERZ?)r} zg}UG$wGRwK45412!0~}z?V16K?LHMyO*|C8v~|7xQd)UjRl3eE`ovRdPeNFIJeKMT Fk6-_XIk5l$ literal 0 HcmV?d00001 diff --git a/tests/regression_tests/random_ray_entropy/results_true.dat b/tests/regression_tests/random_ray_entropy/results_true.dat new file mode 100644 index 00000000000..1c24e3b1adf --- /dev/null +++ b/tests/regression_tests/random_ray_entropy/results_true.dat @@ -0,0 +1,13 @@ +k-combined: +1.000000E+00 0.000000E+00 +entropy: +9.000000E+00 +9.000000E+00 +9.000000E+00 +9.000000E+00 +9.000000E+00 +9.000000E+00 +9.000000E+00 +9.000000E+00 +9.000000E+00 +9.000000E+00 diff --git a/tests/regression_tests/random_ray_entropy/settings.xml b/tests/regression_tests/random_ray_entropy/settings.xml new file mode 100644 index 00000000000..c80c489acc2 --- /dev/null +++ b/tests/regression_tests/random_ray_entropy/settings.xml @@ -0,0 +1,17 @@ + + + eigenvalue + 200 + 10 + 5 + multi-group + + + + -50.0 -50.0 -1 50.0 50.0 1 + + + 40.0 + 400.0 + + diff --git a/tests/regression_tests/random_ray_entropy/test.py b/tests/regression_tests/random_ray_entropy/test.py new file mode 100644 index 00000000000..a3cba65ad0b --- /dev/null +++ b/tests/regression_tests/random_ray_entropy/test.py @@ -0,0 +1,33 @@ +import glob +import os + +from openmc import StatePoint + +from tests.testing_harness import TestHarness + + +class EntropyTestHarness(TestHarness): + def _get_results(self): + """Digest info in the statepoint and return as a string.""" + # Read the statepoint file. + statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0] + with StatePoint(statepoint) as sp: + # Write out k-combined. + outstr = 'k-combined:\n' + outstr += '{:12.6E} {:12.6E}\n'.format(sp.keff.n, sp.keff.s) + + # Write out entropy data. + outstr += 'entropy:\n' + results = ['{:12.6E}'.format(x) for x in sp.entropy] + outstr += '\n'.join(results) + '\n' + + return outstr + +''' +# This test is adapted from "Monte Carlo power iteration: Entropy and spatial correlations," +M. Nowak et al. The cross sections are defined explicitly so that the value for entropy +is exactly 9 and the eigenvalue is exactly 1. +''' +def test_entropy(): + harness = EntropyTestHarness('statepoint.10.h5') + harness.main() From 6fd0bbee2e64613d3e0bb4bd18d4ab8ca6407707 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Mon, 3 Jun 2024 12:54:20 -0500 Subject: [PATCH 02/17] Added parallelization and removed erroneous edits. --- examples/pincell_multigroup/build_xml.py | 2 +- examples/pincell_random_ray/build_xml.py | 8 ----- src/random_ray/flat_source_domain.cpp | 38 +++++++++++++----------- src/settings.cpp | 5 ++++ 4 files changed, 26 insertions(+), 27 deletions(-) diff --git a/examples/pincell_multigroup/build_xml.py b/examples/pincell_multigroup/build_xml.py index ca128a25310..019ae6a113f 100644 --- a/examples/pincell_multigroup/build_xml.py +++ b/examples/pincell_multigroup/build_xml.py @@ -90,7 +90,7 @@ # Create a region represented as the inside of a rectangular prism pitch = 1.26 -box = openmc.model.RectangularPrism(pitch/2, pitch/2, boundary_type='reflective') +box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective') # Instantiate Cells fuel = openmc.Cell(fill=uo2, region=-fuel_or, name='fuel') diff --git a/examples/pincell_random_ray/build_xml.py b/examples/pincell_random_ray/build_xml.py index 9a5fa4f423d..b3dd8020a51 100644 --- a/examples/pincell_random_ray/build_xml.py +++ b/examples/pincell_random_ray/build_xml.py @@ -161,14 +161,6 @@ settings.random_ray['distance_inactive'] = 40.0 settings.random_ray['distance_active'] = 400.0 -entropy_mesh = openmc.RegularMesh() -entropy_mesh.lower_left = (-pitch/2, -pitch/2, -pitch/2) -entropy_mesh.upper_right = (pitch/2, pitch/2, pitch/2) -entropy_mesh.dimension = (8, 8, 8) -settings.entropy_mesh = entropy_mesh - -settings.entropy_on = True - settings.export_to_xml() ############################################################################### diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 2fbe0518e50..1d1e45931de 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -228,8 +228,9 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const // Vector for gathering fission source terms for Shannon entropy calculation vector p(n_source_regions_, 0.0f); + double sum = 0.0; -#pragma omp parallel for reduction(+ : fission_rate_old, fission_rate_new) +#pragma omp parallel for reduction(+ : fission_rate_old, fission_rate_new, sum) for (int sr = 0; sr < n_source_regions_; sr++) { // If simulation averaged volume is zero, don't include this cell @@ -256,30 +257,31 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const // Assigning the fission source to the vector for entropy calculation. p[sr] = sr_fission_source_new; + + // Calculating sum for entropy normalization. + sum += sr_fission_source_new; } double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); - // Calculating sum for entropy normalization. - double sum = 0.0; - for (float num : p) { - sum += num; - } - - // Normalize to total weight of bank sites. - for (float& num : p) { - num /= sum; - } + if (settings::entropy_on) { + double H = 0.0; + + // defining an inverse sum for better performance + double inverse_sum = 1.0/sum; + + #pragma omp parallel for reduction(+ : H) + for (int i = 0; i < n_source_regions_; i++) { + // Normalize to total weight of bank sites. + p[i] *= inverse_sum; + // Sum values to obtain Shannon entropy. + H -= p[i] * std::log(p[i]) / std::log(2.0); + } - // Sum values to obtain Shannon entropy. - double H = 0.0; - for (int i = 0; i < n_source_regions_; i++) { - H -= p[i] * std::log(p[i]) / std::log(2.0); + // Adds entropy value to shared entropy vector in openmc namespace. + simulation::entropy.push_back(H); } - // Adds entropy value to shared entropy vector in openmc namespace. - simulation::entropy.push_back(H); - return k_eff_new; } diff --git a/src/settings.cpp b/src/settings.cpp index 539abad10a9..ef49869495f 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -625,6 +625,11 @@ void read_settings_xml(pugi::xml_node root) // Shannon Entropy mesh if (solver_type == SolverType::RANDOM_RAY) { + if (check_for_node(root, "entropy_mesh")) { + fatal_error( + "Random ray uses FSRs to compute the Shannon entropy. " + "No user-defined entropy mesh is supported."); + } entropy_on = true; } else if (solver_type == SolverType::MONTE_CARLO) { if (check_for_node(root, "entropy_mesh")) { From 227ed795d51d250b1f812f0106512db206b5eea0 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Mon, 3 Jun 2024 17:02:56 -0500 Subject: [PATCH 03/17] Small edit to reflect entropy always being calculated in random ray and performance improvement. --- src/random_ray/flat_source_domain.cpp | 31 +++++++++++++-------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 1d1e45931de..f4d161d42c3 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -262,26 +262,25 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const sum += sr_fission_source_new; } - double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); - - if (settings::entropy_on) { - double H = 0.0; + double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); + double H = 0.0; // defining an inverse sum for better performance - double inverse_sum = 1.0/sum; - - #pragma omp parallel for reduction(+ : H) - for (int i = 0; i < n_source_regions_; i++) { - // Normalize to total weight of bank sites. - p[i] *= inverse_sum; - // Sum values to obtain Shannon entropy. - H -= p[i] * std::log(p[i]) / std::log(2.0); - } - - // Adds entropy value to shared entropy vector in openmc namespace. - simulation::entropy.push_back(H); + double inverse_sum = 1.0/sum; + // defining inverse log(2) for better performance + const double inv_log2 = 1.0/std::log(2.0); + +#pragma omp parallel for reduction(+ : H) + for (int i = 0; i < n_source_regions_; i++) { + // Normalize to total weight of bank sites. + p[i] *= inverse_sum; + // Sum values to obtain Shannon entropy. + H -= p[i] * std::log(p[i]) * inv_log2; } + // Adds entropy value to shared entropy vector in openmc namespace. + simulation::entropy.push_back(H); + return k_eff_new; } From d2ba8675fa9bbbe77aef44e932fbc477b7677962 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Wed, 5 Jun 2024 15:27:28 -0500 Subject: [PATCH 04/17] Fixed bug where Shannon entropy calculation included FSRs with no fission source, resulting in nan entropy values. --- src/random_ray/flat_source_domain.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index f4d161d42c3..7ff4edaead0 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -272,10 +272,13 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const #pragma omp parallel for reduction(+ : H) for (int i = 0; i < n_source_regions_; i++) { - // Normalize to total weight of bank sites. - p[i] *= inverse_sum; - // Sum values to obtain Shannon entropy. - H -= p[i] * std::log(p[i]) * inv_log2; + // Only if FSR has fission source + if (p[i] != 0.0f) { + // Normalize to total weight of bank sites. + p[i] *= inverse_sum; + // Sum values to obtain Shannon entropy. + H -= p[i] * std::log(p[i]) * inv_log2; + } } // Adds entropy value to shared entropy vector in openmc namespace. From 9a334aa9aa768df2af5bd1ffafeed89ffd6f058d Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Wed, 5 Jun 2024 20:45:40 -0500 Subject: [PATCH 05/17] Updated geometry in random_ray_entropy test for clarity. --- .../random_ray_entropy/geometry.xml | 624 +++--------------- 1 file changed, 85 insertions(+), 539 deletions(-) diff --git a/tests/regression_tests/random_ray_entropy/geometry.xml b/tests/regression_tests/random_ray_entropy/geometry.xml index 9a0e6640193..4c87bbbfb9a 100644 --- a/tests/regression_tests/random_ray_entropy/geometry.xml +++ b/tests/regression_tests/random_ray_entropy/geometry.xml @@ -1,542 +1,88 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + 12.5 12.5 12.5 + 8 8 8 + 0.0 0.0 0.0 + +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 + +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 + +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 + +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 + +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 + +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 + +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 + +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 +1 1 1 1 1 1 1 1 + + + + + + + From 2b07cda1de44e41ed19499529914df493a013276 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Wed, 5 Jun 2024 21:33:57 -0500 Subject: [PATCH 06/17] Adding documentation on Shannon entropy in TRRM. --- docs/source/methods/random_ray.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index aa436784766..616ae5afaec 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -733,6 +733,14 @@ How are Tallies Handled? Most tallies, filters, and scores that you would expect to work with a multigroup solver like random ray should work. For example, you can define 3D mesh tallies with energy filters and flux, fission, and nu-fission scores, etc. + +Additionally, as :math:`k_{eff}` is updated at each generation, the fission +source at each FSR is used to compute the Shannon entropy. This follows the +:ref:`same procedure for computing Shannon entropy in continuous-energy or +multigroup Monte Carlo simulations <_method-shannon-entropy>`, except that +fission sites at FSRs are considered, rather than fission sites of user-defined +regular meshes. + There are some restrictions though. For starters, it is assumed that all filter mesh boundaries will conform to physical surface boundaries (or lattice boundaries) in the simulation geometry. It is acceptable for multiple cells From 5e66462bf7af57dbd9c3452a4ce4e0fc26dda712 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Thu, 6 Jun 2024 14:54:25 -0500 Subject: [PATCH 07/17] Fixed and updated additions to the documentation on Shannon entropy in TRRM. --- docs/source/methods/eigenvalue.rst | 2 ++ docs/source/methods/random_ray.rst | 2 +- docs/source/usersguide/random_ray.rst | 4 +++- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/source/methods/eigenvalue.rst b/docs/source/methods/eigenvalue.rst index 41bf8654920..99cf3566ea8 100644 --- a/docs/source/methods/eigenvalue.rst +++ b/docs/source/methods/eigenvalue.rst @@ -55,6 +55,8 @@ in :ref:`fission-bank-algorithms`. Source Convergence Issues ------------------------- +.. _methods-shannon-entropy: + Diagnosing Convergence with Shannon Entropy ------------------------------------------- diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index 616ae5afaec..a72f780a96d 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -737,7 +737,7 @@ mesh tallies with energy filters and flux, fission, and nu-fission scores, etc. Additionally, as :math:`k_{eff}` is updated at each generation, the fission source at each FSR is used to compute the Shannon entropy. This follows the :ref:`same procedure for computing Shannon entropy in continuous-energy or -multigroup Monte Carlo simulations <_method-shannon-entropy>`, except that +multigroup Monte Carlo simulations `, except that fission sites at FSRs are considered, rather than fission sites of user-defined regular meshes. diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index dad86c826e0..29de2a4caff 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -46,7 +46,9 @@ magnitude or more. For instance, it may be reasonable to only use 50 inactive batches for a light water reactor simulation with Monte Carlo, but random ray might require 500 or more inactive batches. Similar to Monte Carlo, :ref:`Shannon entropy ` can be used to gauge whether the -combined scattering and fission source has fully developed. +combined scattering and fission source has fully developed. The Shannon entropy +is calculated automatically near the end of each batch and is printed to the +statepoint file. Similar to Monte Carlo, active batches are used in the random ray solver mode to accumulate and converge statistics on unknown quantities (i.e., the random ray From 5f98c38f44f4647d8af103e48cfd39c4cd79af78 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Thu, 6 Jun 2024 15:02:07 -0500 Subject: [PATCH 08/17] clang format updates --- src/random_ray/flat_source_domain.cpp | 8 ++++---- src/settings.cpp | 9 ++++----- src/simulation.cpp | 3 ++- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 7ff4edaead0..994e7be3eb2 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -257,18 +257,18 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const // Assigning the fission source to the vector for entropy calculation. p[sr] = sr_fission_source_new; - + // Calculating sum for entropy normalization. sum += sr_fission_source_new; } - double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); + double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); double H = 0.0; // defining an inverse sum for better performance - double inverse_sum = 1.0/sum; + double inverse_sum = 1.0 / sum; // defining inverse log(2) for better performance - const double inv_log2 = 1.0/std::log(2.0); + const double inv_log2 = 1.0 / std::log(2.0); #pragma omp parallel for reduction(+ : H) for (int i = 0; i < n_source_regions_; i++) { diff --git a/src/settings.cpp b/src/settings.cpp index ef49869495f..be0edde3dee 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -626,9 +626,8 @@ void read_settings_xml(pugi::xml_node root) // Shannon Entropy mesh if (solver_type == SolverType::RANDOM_RAY) { if (check_for_node(root, "entropy_mesh")) { - fatal_error( - "Random ray uses FSRs to compute the Shannon entropy. " - "No user-defined entropy mesh is supported."); + fatal_error("Random ray uses FSRs to compute the Shannon entropy. " + "No user-defined entropy mesh is supported."); } entropy_on = true; } else if (solver_type == SolverType::MONTE_CARLO) { @@ -639,8 +638,8 @@ void read_settings_xml(pugi::xml_node root) "Mesh {} specified for Shannon entropy does not exist.", temp)); } - auto* m = - dynamic_cast(model::meshes[model::mesh_map.at(temp)].get()); + auto* m = dynamic_cast( + model::meshes[model::mesh_map.at(temp)].get()); if (!m) fatal_error("Only regular meshes can be used as an entropy mesh"); simulation::entropy_mesh = m; diff --git a/src/simulation.cpp b/src/simulation.cpp index eb18b9c3e5d..3bf60fd10d5 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -521,7 +521,8 @@ void finalize_generation() if (settings::run_mode == RunMode::EIGENVALUE) { // Calculate shannon entropy - if (settings::entropy_on && settings::solver_type == SolverType::MONTE_CARLO) + if (settings::entropy_on && + settings::solver_type == SolverType::MONTE_CARLO) shannon_entropy(); // Collect results and statistics From 663242710843bc557e38878faf5cab549403d587 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Thu, 6 Jun 2024 15:19:37 -0500 Subject: [PATCH 09/17] Rewrap on docs. --- docs/source/methods/eigenvalue.rst | 8 ++++---- docs/source/methods/random_ray.rst | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/source/methods/eigenvalue.rst b/docs/source/methods/eigenvalue.rst index 99cf3566ea8..45099cb1ecf 100644 --- a/docs/source/methods/eigenvalue.rst +++ b/docs/source/methods/eigenvalue.rst @@ -62,10 +62,10 @@ Diagnosing Convergence with Shannon Entropy As discussed earlier, it is necessary to converge both :math:`k_{eff}` and the source distribution before any tallies can begin. Moreover, the convergence rate -of the source distribution is in general slower than that of -:math:`k_{eff}`. One should thus examine not only the convergence of -:math:`k_{eff}` but also the convergence of the source distribution in order to -make decisions on when to start active batches. +of the source distribution is in general slower than that of :math:`k_{eff}`. +One should thus examine not only the convergence of :math:`k_{eff}` but also the +convergence of the source distribution in order to make decisions on when to +start active batches. However, the representation of the source distribution makes it a bit more difficult to analyze its convergence. Since :math:`k_{eff}` is a scalar diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index a72f780a96d..9eeec27a6f7 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -736,9 +736,9 @@ mesh tallies with energy filters and flux, fission, and nu-fission scores, etc. Additionally, as :math:`k_{eff}` is updated at each generation, the fission source at each FSR is used to compute the Shannon entropy. This follows the -:ref:`same procedure for computing Shannon entropy in continuous-energy or -multigroup Monte Carlo simulations `, except that -fission sites at FSRs are considered, rather than fission sites of user-defined +:ref:`same procedure for computing Shannon entropy in continuous-energy or +multigroup Monte Carlo simulations `, except that +fission sites at FSRs are considered, rather than fission sites of user-defined regular meshes. There are some restrictions though. For starters, it is assumed that all filter From 3c06976610eaa8e9565b19cc9b25ccff0f20a935 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Mon, 10 Jun 2024 11:33:17 -0500 Subject: [PATCH 10/17] Resolved volume issue and updated test. --- src/random_ray/flat_source_domain.cpp | 21 ++++++++++--------- .../random_ray_entropy/results_true.dat | 20 +++++++++--------- .../random_ray_entropy/settings.xml | 4 ++-- 3 files changed, 23 insertions(+), 22 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 994e7be3eb2..24a35d86e2b 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -230,7 +230,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const vector p(n_source_regions_, 0.0f); double sum = 0.0; -#pragma omp parallel for reduction(+ : fission_rate_old, fission_rate_new, sum) +#pragma omp parallel for reduction(+ : fission_rate_old, fission_rate_new) for (int sr = 0; sr < n_source_regions_; sr++) { // If simulation averaged volume is zero, don't include this cell @@ -252,21 +252,22 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const sr_fission_source_new += nu_sigma_f * scalar_flux_new_[idx]; } + // Temporary FSR fission rate for performance. + double temp_sr_fission_rate = sr_fission_source_new * volume; + + // Calculating old and new fission rates. fission_rate_old += sr_fission_source_old * volume; - fission_rate_new += sr_fission_source_new * volume; + fission_rate_new += temp_sr_fission_rate; // Assigning the fission source to the vector for entropy calculation. - p[sr] = sr_fission_source_new; - - // Calculating sum for entropy normalization. - sum += sr_fission_source_new; + p[sr] = temp_sr_fission_rate; } double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); double H = 0.0; // defining an inverse sum for better performance - double inverse_sum = 1.0 / sum; + double inverse_sum = 1 / fission_rate_new; // defining inverse log(2) for better performance const double inv_log2 = 1.0 / std::log(2.0); @@ -274,10 +275,10 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const for (int i = 0; i < n_source_regions_; i++) { // Only if FSR has fission source if (p[i] != 0.0f) { - // Normalize to total weight of bank sites. - p[i] *= inverse_sum; + // Normalize to total weight of bank sites. p_i for better performance + float p_i = p[i] * inverse_sum; // Sum values to obtain Shannon entropy. - H -= p[i] * std::log(p[i]) * inv_log2; + H -= p_i * std::log(p_i) * inv_log2; } } diff --git a/tests/regression_tests/random_ray_entropy/results_true.dat b/tests/regression_tests/random_ray_entropy/results_true.dat index 1c24e3b1adf..25425453871 100644 --- a/tests/regression_tests/random_ray_entropy/results_true.dat +++ b/tests/regression_tests/random_ray_entropy/results_true.dat @@ -1,13 +1,13 @@ k-combined: 1.000000E+00 0.000000E+00 entropy: -9.000000E+00 -9.000000E+00 -9.000000E+00 -9.000000E+00 -9.000000E+00 -9.000000E+00 -9.000000E+00 -9.000000E+00 -9.000000E+00 -9.000000E+00 +8.863421E+00 +8.933584E+00 +8.960553E+00 +8.967921E+00 +8.976016E+00 +8.981856E+00 +8.983670E+00 +8.986584E+00 +8.987732E+00 +8.988186E+00 diff --git a/tests/regression_tests/random_ray_entropy/settings.xml b/tests/regression_tests/random_ray_entropy/settings.xml index c80c489acc2..81deaa7751d 100644 --- a/tests/regression_tests/random_ray_entropy/settings.xml +++ b/tests/regression_tests/random_ray_entropy/settings.xml @@ -1,14 +1,14 @@ eigenvalue - 200 + 100 10 5 multi-group - -50.0 -50.0 -1 50.0 50.0 1 + 0.0 0.0 0.0 100.0 100.0 100.0 40.0 From e37c76d60d5fba756053056903db3d99156c01b2 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Mon, 10 Jun 2024 13:08:27 -0500 Subject: [PATCH 11/17] Changes to docs to reflect volume-weighted approach. --- docs/source/methods/random_ray.rst | 18 +++++++++++++++++- docs/source/usersguide/random_ray.rst | 3 ++- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index 9eeec27a6f7..d8d69f338bb 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -739,7 +739,23 @@ source at each FSR is used to compute the Shannon entropy. This follows the :ref:`same procedure for computing Shannon entropy in continuous-energy or multigroup Monte Carlo simulations `, except that fission sites at FSRs are considered, rather than fission sites of user-defined -regular meshes. +regular meshes. Thus, the volume-weighted fission rate is considered instead, +and the fraction of source sites is adjusted such that: + +.. math:: + :label: fraction-source-random-ray + + S_i = \frac{\text{Source sites in $i$-th FSR $*$ Volume of $i$-th FSR}}{\text{Total number of FSRs}} + +The Shannon entropy is then computed normally as + +.. math:: + :label: shannon-entropy-random-ray + + H = - \sum_{i=1}^N S_i \log_2 S_i + +where :math:`N` is the number of FSRs. FSRs with no fission source sites are skipped +to avoid taking the logarithm of zero in :eq:`shannon-entropy-random-ray`. There are some restrictions though. For starters, it is assumed that all filter mesh boundaries will conform to physical surface boundaries (or lattice diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 29de2a4caff..a87c2d71899 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -48,7 +48,8 @@ might require 500 or more inactive batches. Similar to Monte Carlo, :ref:`Shannon entropy ` can be used to gauge whether the combined scattering and fission source has fully developed. The Shannon entropy is calculated automatically near the end of each batch and is printed to the -statepoint file. +statepoint file. Unlike Monte Carlo, an entropy mesh does not need to be defined, +as the Shannon entropy is calculated over FSRs using a volume-weighted approach. Similar to Monte Carlo, active batches are used in the random ray solver mode to accumulate and converge statistics on unknown quantities (i.e., the random ray From 0fb5cb9a0bf4dd23f46d0ea8121c40131b90d208 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Fri, 14 Jun 2024 11:09:44 -0500 Subject: [PATCH 12/17] Revisons to code and docs based on review. --- docs/source/methods/eigenvalue.rst | 5 +++ docs/source/methods/random_ray.rst | 50 +++++++++++++++++---------- docs/source/usersguide/random_ray.rst | 20 +++++++---- src/eigenvalue.cpp | 2 +- src/random_ray/flat_source_domain.cpp | 24 ++++++------- src/settings.cpp | 2 +- 6 files changed, 63 insertions(+), 40 deletions(-) diff --git a/docs/source/methods/eigenvalue.rst b/docs/source/methods/eigenvalue.rst index 45099cb1ecf..9a633762772 100644 --- a/docs/source/methods/eigenvalue.rst +++ b/docs/source/methods/eigenvalue.rst @@ -110,6 +110,11 @@ at plots of :math:`k_{eff}` and the Shannon entropy. A number of methods have been proposed (see e.g. [Romano]_, [Ueki]_), but each of these is not without problems. +Shannon entropy is calculated differently for the random ray solver, as +described :ref:`in the random ray theory section +`, as the Shannon entropy does not serve as +a diagnostic tool for convergence of the scattering source. + --------------------------- Uniform Fission Site Method --------------------------- diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index d8d69f338bb..4dddcbccbe2 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -411,6 +411,8 @@ which when partially simplified becomes: Note that there are now four (seemingly identical) volume terms in this equation. +.. _methods-volume-dilemma: + ~~~~~~~~~~~~~~ Volume Dilemma ~~~~~~~~~~~~~~ @@ -734,18 +736,35 @@ Most tallies, filters, and scores that you would expect to work with a multigroup solver like random ray should work. For example, you can define 3D mesh tallies with energy filters and flux, fission, and nu-fission scores, etc. -Additionally, as :math:`k_{eff}` is updated at each generation, the fission -source at each FSR is used to compute the Shannon entropy. This follows the -:ref:`same procedure for computing Shannon entropy in continuous-energy or -multigroup Monte Carlo simulations `, except that -fission sites at FSRs are considered, rather than fission sites of user-defined -regular meshes. Thus, the volume-weighted fission rate is considered instead, -and the fraction of source sites is adjusted such that: +There are some restrictions though. For starters, it is assumed that all filter +mesh boundaries will conform to physical surface boundaries (or lattice +boundaries) in the simulation geometry. It is acceptable for multiple cells +(FSRs) to be contained within a filter mesh cell (e.g., pincell-level or +assembly-level tallies should work), but it is currently left as undefined +behavior if a single simulation cell is able to score to multiple filter mesh +cells. In the future, the capability to fully support mesh tallies may be added +to OpenMC, but for now this restriction needs to be respected. + +.. _methods-shannon-entropy-random-ray: + +----------------------------- +Shannon Entropy in Random Ray +----------------------------- + +As :math:`k_{eff}` is updated at each generation, the fission source at each FSR +is used to compute the Shannon entropy. This follows the :ref:`same procedure +for computing Shannon entropy in continuous-energy or multigroup Monte Carlo +simulations `, except that fission sources at FSRs are +considered, rather than fission sites of user-defined regular meshes. Thus, the +volume-weighted fission rate is considered instead, and the fraction of fission +sources is adjusted such that: .. math:: :label: fraction-source-random-ray - S_i = \frac{\text{Source sites in $i$-th FSR $*$ Volume of $i$-th FSR}}{\text{Total number of FSRs}} + S_i = \frac{\text{Fission source in FSR $i \times$ Volume of FSR + $i$}}{\text{Total fission source}} = \frac{Q_{i,g} \times + V_{i}}{\sum_{i=1}^{g} Q_{i,g}} The Shannon entropy is then computed normally as @@ -754,17 +773,10 @@ The Shannon entropy is then computed normally as H = - \sum_{i=1}^N S_i \log_2 S_i -where :math:`N` is the number of FSRs. FSRs with no fission source sites are skipped -to avoid taking the logarithm of zero in :eq:`shannon-entropy-random-ray`. - -There are some restrictions though. For starters, it is assumed that all filter -mesh boundaries will conform to physical surface boundaries (or lattice -boundaries) in the simulation geometry. It is acceptable for multiple cells -(FSRs) to be contained within a filter mesh cell (e.g., pincell-level or -assembly-level tallies should work), but it is currently left as undefined -behavior if a single simulation cell is able to score to multiple filter mesh -cells. In the future, the capability to fully support mesh tallies may be added -to OpenMC, but for now this restriction needs to be respected. +where :math:`N` is the number of FSRs. FSRs with no fission source (or, +occassionally, negative fission source, :ref:`due to the volume estimator +problem `) are skipped to avoid taking an undefined +logarithm in :eq:`shannon-entropy-random-ray`. --------------------------- Fundamental Sources of Bias diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index a87c2d71899..e5bf27149ff 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -44,12 +44,7 @@ The additional burden of converging the scattering source generally results in a higher requirement for the number of inactive batches---often by an order of magnitude or more. For instance, it may be reasonable to only use 50 inactive batches for a light water reactor simulation with Monte Carlo, but random ray -might require 500 or more inactive batches. Similar to Monte Carlo, -:ref:`Shannon entropy ` can be used to gauge whether the -combined scattering and fission source has fully developed. The Shannon entropy -is calculated automatically near the end of each batch and is printed to the -statepoint file. Unlike Monte Carlo, an entropy mesh does not need to be defined, -as the Shannon entropy is calculated over FSRs using a volume-weighted approach. +might require 500 or more inactive batches. Similar to Monte Carlo, active batches are used in the random ray solver mode to accumulate and converge statistics on unknown quantities (i.e., the random ray @@ -63,6 +58,19 @@ solver:: settings.batches = 1200 settings.inactive = 600 +--------------- +Shannon Entropy +--------------- + +Similar to Monte Carlo, :ref:`Shannon entropy +` can be used to gauge whether the fission +source has fully developed. However, unlike Monte Carlo, the Shannon entropy +computed for random ray does not show convergence of the scattering source. The +Shannon entropy is calculated automatically after each batch and is printed to +the statepoint file. Unlike Monte Carlo, an entropy mesh does not need to be +defined, as the Shannon entropy is calculated over FSRs using a volume-weighted +approach. + ------------------------------- Inactive Ray Length (Dead Zone) ------------------------------- diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index d5410094b18..8669d76f947 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -556,7 +556,7 @@ void shannon_entropy() double H = 0.0; for (auto p_i : p) { if (p_i > 0.0) { - H -= p_i * std::log(p_i) / std::log(2.0); + H -= p_i * std::log2(p_i); } } diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 24a35d86e2b..4f411b2554c 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -228,7 +228,6 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const // Vector for gathering fission source terms for Shannon entropy calculation vector p(n_source_regions_, 0.0f); - double sum = 0.0; #pragma omp parallel for reduction(+ : fission_rate_old, fission_rate_new) for (int sr = 0; sr < n_source_regions_; sr++) { @@ -252,15 +251,16 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const sr_fission_source_new += nu_sigma_f * scalar_flux_new_[idx]; } - // Temporary FSR fission rate for performance. - double temp_sr_fission_rate = sr_fission_source_new * volume; + // Compute total fission rates in FSR + sr_fission_source_old *= volume; + sr_fission_source_new *= volume; - // Calculating old and new fission rates. - fission_rate_old += sr_fission_source_old * volume; - fission_rate_new += temp_sr_fission_rate; + // Accumulate totals + fission_rate_old += sr_fission_source_old; + fission_rate_new += sr_fission_source_new; - // Assigning the fission source to the vector for entropy calculation. - p[sr] = temp_sr_fission_rate; + // Store total fission rate in the FSR for Shannon calculation + p[sr] = sr_fission_source_new; } double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); @@ -268,17 +268,15 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const double H = 0.0; // defining an inverse sum for better performance double inverse_sum = 1 / fission_rate_new; - // defining inverse log(2) for better performance - const double inv_log2 = 1.0 / std::log(2.0); #pragma omp parallel for reduction(+ : H) for (int i = 0; i < n_source_regions_; i++) { - // Only if FSR has fission source - if (p[i] != 0.0f) { + // Only if FSR has non-negative and non-zero fission source + if (p[i] > 0.0f) { // Normalize to total weight of bank sites. p_i for better performance float p_i = p[i] * inverse_sum; // Sum values to obtain Shannon entropy. - H -= p_i * std::log(p_i) * inv_log2; + H -= p_i * std::log2(p_i); } } diff --git a/src/settings.cpp b/src/settings.cpp index be0edde3dee..67ce2f14bb6 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -623,7 +623,7 @@ void read_settings_xml(pugi::xml_node root) } } - // Shannon Entropy mesh + // Shannon entropy if (solver_type == SolverType::RANDOM_RAY) { if (check_for_node(root, "entropy_mesh")) { fatal_error("Random ray uses FSRs to compute the Shannon entropy. " From 1aaca5ad4c2b845d852a786bb705323ef947e177 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Tue, 25 Jun 2024 10:44:06 -0500 Subject: [PATCH 13/17] Changing index on a loop to sr convention. --- src/random_ray/flat_source_domain.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index f12a0afc015..6b34e359c44 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -323,11 +323,11 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const double inverse_sum = 1 / fission_rate_new; #pragma omp parallel for reduction(+ : H) - for (int i = 0; i < n_source_regions_; i++) { + for (int sr = 0; sr < n_source_regions_; sr++) { // Only if FSR has non-negative and non-zero fission source - if (p[i] > 0.0f) { + if (p[sr] > 0.0f) { // Normalize to total weight of bank sites. p_i for better performance - float p_i = p[i] * inverse_sum; + float p_i = p[sr] * inverse_sum; // Sum values to obtain Shannon entropy. H -= p_i * std::log2(p_i); } From c56a11af637854085306c5033c0ccd944ad097d5 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Tue, 25 Jun 2024 10:57:49 -0500 Subject: [PATCH 14/17] Edits to documentation. --- docs/source/methods/eigenvalue.rst | 6 ++++-- docs/source/methods/random_ray.rst | 4 ++-- docs/source/usersguide/random_ray.rst | 10 ++++------ 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/docs/source/methods/eigenvalue.rst b/docs/source/methods/eigenvalue.rst index 4629b015b16..8abcc09574b 100644 --- a/docs/source/methods/eigenvalue.rst +++ b/docs/source/methods/eigenvalue.rst @@ -112,8 +112,10 @@ problems. Shannon entropy is calculated differently for the random ray solver, as described :ref:`in the random ray theory section -`, as the Shannon entropy does not serve as -a diagnostic tool for convergence of the scattering source. +`. Additionally, as the Shannon entropy only +serves as a diagnostic tool for convergence of the fission source distribution, +there is currently no diagnostic to determine if the scattering source +distribution in random ray is converged. --------------------------- Uniform Fission Site Method diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index e65fb16c6f6..14fa6956c0c 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -775,8 +775,8 @@ sources is adjusted such that: :label: fraction-source-random-ray S_i = \frac{\text{Fission source in FSR $i \times$ Volume of FSR - $i$}}{\text{Total fission source}} = \frac{Q_{i,g} \times - V_{i}}{\sum_{i=1}^{g} Q_{i,g}} + $i$}}{\text{Total fission source}} = \frac{Q_{i} V_{i}}{\sum_{i=1}^{i} Q_{i} + V_{i}} The Shannon entropy is then computed normally as diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 48ceec8b766..2a038b12257 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -68,12 +68,10 @@ Shannon Entropy Similar to Monte Carlo, :ref:`Shannon entropy ` can be used to gauge whether the fission -source has fully developed. However, unlike Monte Carlo, the Shannon entropy -computed for random ray does not show convergence of the scattering source. The -Shannon entropy is calculated automatically after each batch and is printed to -the statepoint file. Unlike Monte Carlo, an entropy mesh does not need to be -defined, as the Shannon entropy is calculated over FSRs using a volume-weighted -approach. +source has fully developed. The Shannon entropy is calculated automatically +after each batch and is printed to the statepoint file. Unlike Monte Carlo, an +entropy mesh does not need to be defined, as the Shannon entropy is calculated +over FSRs using a volume-weighted approach. ------------------------------- Inactive Ray Length (Dead Zone) From a5d54eb7b90c668c548c1f0f2a43d5a337dc6d74 Mon Sep 17 00:00:00 2001 From: Ethan Krammer Date: Tue, 25 Jun 2024 11:06:01 -0500 Subject: [PATCH 15/17] clang update on comment in flat_source_domain --- src/random_ray/flat_source_domain.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 6b34e359c44..654f500a448 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -556,7 +556,7 @@ void FlatSourceDomain::random_ray_tally() tally.results_(task.filter_idx, task.score_idx, TallyResult::VALUE) += score; } // end tally task loop - } // end energy group loop + } // end energy group loop // For flux tallies, the total volume of the spatial region is needed // for normalizing the flux. We store this volume in a separate tensor. From f45b48545308c287b811e450fc714ca8d1bb84bd Mon Sep 17 00:00:00 2001 From: John Tramm Date: Wed, 26 Jun 2024 14:47:57 +0000 Subject: [PATCH 16/17] removed troublesome comments causing issue between clang-format v15 and v18 --- src/random_ray/flat_source_domain.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 654f500a448..792e97da0c1 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -555,8 +555,8 @@ void FlatSourceDomain::random_ray_tally() #pragma omp atomic tally.results_(task.filter_idx, task.score_idx, TallyResult::VALUE) += score; - } // end tally task loop - } // end energy group loop + } + } // For flux tallies, the total volume of the spatial region is needed // for normalizing the flux. We store this volume in a separate tensor. From 9d62eb29120c27556dba1d173ab9e4a79fa81ef7 Mon Sep 17 00:00:00 2001 From: Ethan Krammer <68677905+ethankrammer@users.noreply.github.com> Date: Wed, 26 Jun 2024 16:42:04 -0500 Subject: [PATCH 17/17] Fix to erroneous equation. --- docs/source/methods/random_ray.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index 14fa6956c0c..d46d3ad4165 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -775,8 +775,8 @@ sources is adjusted such that: :label: fraction-source-random-ray S_i = \frac{\text{Fission source in FSR $i \times$ Volume of FSR - $i$}}{\text{Total fission source}} = \frac{Q_{i} V_{i}}{\sum_{i=1}^{i} Q_{i} - V_{i}} + $i$}}{\text{Total fission source}} = \frac{Q_{i} V_{i}}{\sum_{i=1}^{i=N} + Q_{i} V_{i}} The Shannon entropy is then computed normally as