Implementation of Shannon Entropy for Random Ray (#3030)

Co-authored-by: Ethan Krammer <ethan@DESKTOP-MGFGK9N>
Co-authored-by: John Tramm <john.tramm@gmail.com>

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
@@ -108,6 +110,13 @@ 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
+<methods-shannon-entropy-random-ray>`. 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
 ---------------------------

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
 ~~~~~~~~~~~~~~
@@ -745,6 +747,7 @@ 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.
+
 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
@@ -754,6 +757,39 @@ 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 <methods-shannon-entropy>`, 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{Fission source in FSR $i \times$ Volume of FSR
+    $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
+
+.. 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 (or,
+occassionally, negative fission source, :ref:`due to the volume estimator
+problem <methods-volume-dilemma>`) are skipped to avoid taking an undefined
+logarithm in :eq:`shannon-entropy-random-ray`.
+
 .. _usersguide_fixed_source_methods:
 
 ------------

docs/source/usersguide/random_ray.rst

@@ -62,6 +62,17 @@ solver::
     settings.batches = 1200
     settings.inactive = 600
 
+---------------
+Shannon Entropy
+---------------
+
+Similar to Monte Carlo, :ref:`Shannon entropy
+<methods-shannon-entropy-random-ray>` can be used to gauge whether the fission
+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)
 -------------------------------

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);
       }
     }
 

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/material.h"
 #include "openmc/message_passing.h"
@@ -278,6 +279,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<float> 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++) {
 
@@ -300,12 +304,38 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
       sr_fission_source_new += nu_sigma_f * scalar_flux_new_[idx];
     }
 
-    fission_rate_old += sr_fission_source_old * volume;
-    fission_rate_new += sr_fission_source_new * volume;
+    // Compute total fission rates in FSR
+    sr_fission_source_old *= volume;
+    sr_fission_source_new *= volume;
+
+    // Accumulate totals
+    fission_rate_old += sr_fission_source_old;
+    fission_rate_new += sr_fission_source_new;
+
+    // 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);
 
+  double H = 0.0;
+  // defining an inverse sum for better performance
+  double inverse_sum = 1 / fission_rate_new;
+
+#pragma omp parallel for reduction(+ : H)
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    // Only if FSR has non-negative and non-zero fission source
+    if (p[sr] > 0.0f) {
+      // Normalize to total weight of bank sites. p_i for better performance
+      float p_i = p[sr] * inverse_sum;
+      // Sum values to obtain Shannon entropy.
+      H -= p_i * std::log2(p_i);
+    }
+  }
+
+  // Adds entropy value to shared entropy vector in openmc namespace.
+  simulation::entropy.push_back(H);
+
   return k_eff_new;
 }
 
@@ -525,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.

src/settings.cpp

@@ -627,30 +627,37 @@ 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));
+  // 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. "
+                  "No user-defined entropy mesh is supported.");
     }
-
-    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;
+  } 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));
+      }
 
-  } 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.");
+      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;
+
+    } 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

@@ -531,7 +531,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()