Implementation of Shannon Entropy for Random Ray

docs/source/methods/eigenvalue.rst

@@ -55,15 +55,17 @@ in :ref:`fission-bank-algorithms`.
 Source Convergence Issues
 -------------------------
 
+.. _methods-shannon-entropy:
+
 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

docs/source/methods/random_ray.rst

@@ -733,6 +733,30 @@ 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 <methods-shannon-entropy>`, 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:
+
+.. 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
 boundaries) in the simulation geometry. It is acceptable for multiple cells

docs/source/usersguide/random_ray.rst

@@ -46,7 +46,10 @@ 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 <usersguide_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. 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

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,10 @@ 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<float> 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++) {
 
@@ -247,12 +252,39 @@ 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] = 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 / 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) {
+      // 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;
+    }
+  }
+
+  // Adds entropy value to shared entropy vector in openmc namespace.
+  simulation::entropy.push_back(H);
+
   return k_eff_new;
 }
 

src/settings.cpp

@@ -624,29 +624,36 @@ 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) {
+    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")) {
+      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<RegularMesh*>(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<RegularMesh*>(
+        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 <entropy> element "
-      "is deprecated. Please create a mesh using <mesh> and then reference "
-      "it by specifying its ID in an <entropy_mesh> element.");
+    } else if (check_for_node(root, "entropy")) {
+      fatal_error(
+        "Specifying a Shannon entropy mesh via the <entropy> element "
+        "is deprecated. Please create a mesh using <mesh> and then reference "
+        "it by specifying its ID in an <entropy_mesh> element.");
+    }
   }
-
   // Uniform fission source weighting mesh
   if (check_for_node(root, "ufs_mesh")) {
     auto temp = std::stoi(get_node_value(root, "ufs_mesh"));

src/simulation.cpp

@@ -521,7 +521,8 @@ 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

tests/regression_tests/random_ray_entropy/geometry.xml

@@ -0,0 +1,88 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<geometry>
+  <cell id="1" material="1" universe="1"/>
+  <cell fill="2" id="2" name="assembly" region="-2 1 -4 3 -6 5" universe="3"/>
+  <lattice id="2">
+    <pitch>12.5 12.5 12.5</pitch>
+    <dimension>8 8 8</dimension>
+    <lower_left>0.0 0.0 0.0</lower_left>
+    <universes>
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 
+1 1 1 1 1 1 1 1 </universes>
+  </lattice>
+  <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+  <surface boundary="reflective" coeffs="100.0" id="2" type="x-plane"/>
+  <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+  <surface boundary="reflective" coeffs="100.0" id="4" type="y-plane"/>
+  <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+  <surface boundary="reflective" coeffs="100.0" id="6" type="z-plane"/>
+</geometry>

tests/regression_tests/random_ray_entropy/materials.xml

@@ -0,0 +1,8 @@
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <cross_sections>mgxs.h5</cross_sections>
+  <material id="1" name="Core Material">
+    <density units="macro" value="1.0"/>
+    <macroscopic name="CoreMaterial"/>
+  </material>
+</materials>

tests/regression_tests/random_ray_entropy/results_true.dat

@@ -0,0 +1,13 @@
+k-combined:
+1.000000E+00 0.000000E+00
+entropy:
+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

tests/regression_tests/random_ray_entropy/settings.xml

@@ -0,0 +1,17 @@
+<?xml version='1.0' encoding='UTF-8'?>
+<settings>
+  <run_mode>eigenvalue</run_mode>
+  <particles>100</particles>
+  <batches>10</batches>
+  <inactive>5</inactive>
+  <energy_mode>multi-group</energy_mode>
+  <random_ray>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="box">
+        <parameters>0.0 0.0 0.0 100.0 100.0 100.0</parameters>
+      </space>
+    </source>
+    <distance_inactive>40.0</distance_inactive>
+    <distance_active>400.0</distance_active>
+  </random_ray>
+</settings>

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()