Random Ray Normalization Improvements (#3051)

Co-authored-by: Olek <[email protected]>

docs/source/methods/random_ray.rst

@@ -740,6 +740,8 @@ scalar flux value for the FSR).
         global::volume[fsr] += s;
     }
 
+.. _methods_random_tallies:
+
 ------------------------
 How are Tallies Handled?
 ------------------------
@@ -757,6 +759,18 @@ 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.
 
+Flux tallies are handled slightly differently than in Monte Carlo. By default,
+in MC, flux tallies are reported in units of tracklength (cm), so must be
+manually normalized by volume by the user to produce an estimate of flux in
+units of cm\ :sup:`-2`\. Alternatively, MC flux tallies can be normalized via a
+separated volume calculation process as discussed in the :ref:`Volume
+Calculation Section<usersguide_volume>`. In random ray, as the volumes are
+computed on-the-fly as part of the transport process, the flux tallies can
+easily be reported either in units of flux (cm\ :sup:`-2`\) or tracklength (cm).
+By default, the unnormalized flux values (units of cm) will be reported. If the
+user wishes to received volume normalized flux tallies, then an option for this
+is available, as described in the :ref:`User Guide<usersguide_flux_norm>`.
+
 .. _methods-shannon-entropy-random-ray:
 
 -----------------------------

docs/source/usersguide/random_ray.rst

@@ -330,6 +330,8 @@ of a two-dimensional 2x2 reflective pincell lattice:
 
 In the future, automated subdivision of FSRs via mesh overlay may be supported.
 
+.. _usersguide_flux_norm:
+
 -------
 Tallies
 -------
@@ -367,6 +369,25 @@ Note that there is no difference between the analog, tracklength, and collision
 estimators in random ray mode as individual particles are not being simulated.
 Tracklength-style tally estimation is inherent to the random ray method.
 
+As discussed in the random ray theory section on :ref:`Random Ray
+Tallies<methods_random_tallies>`, by default flux tallies in the random ray mode
+are not normalized by the spatial tally volumes such that flux tallies are in
+units of cm. While the volume information is readily available as a byproduct of
+random ray integration, the flux value is reported in unnormalized units of cm
+so that the user will be able to compare "apples to apples" with the default
+flux tallies from the Monte Carlo solver (also reported by default in units of
+cm). If volume normalized flux tallies (in units of cm\ :sup:`-2`) are desired,
+then the user can set the ``volume_normalized_flux_tallies`` field in the
+:attr:`openmc.Settings.random_ray` dictionary to ``True``. An example is given
+below:
+
+::
+
+    settings.random_ray['volume_normalized_flux_tallies'] = True
+
+Note that MC mode flux tallies can also be normalized by volume, as discussed in
+the :ref:`Volume Calculation Section<usersguide_volume>` of the user guide.
+
 --------
 Plotting
 --------

include/openmc/random_ray/flat_source_domain.h

@@ -108,6 +108,11 @@ class FlatSourceDomain {
   void all_reduce_replicated_source_regions();
   void convert_external_sources();
   void count_external_source_regions();
+  double compute_fixed_source_normalization_factor() const;
+
+  //----------------------------------------------------------------------------
+  // Static Data members
+  static bool volume_normalized_flux_tallies_;
 
   //----------------------------------------------------------------------------
   // Public Data members

openmc/settings.py

@@ -157,6 +157,12 @@ class Settings:
         :ray_source:
             Starting ray distribution (must be uniform in space and angle) as
             specified by a :class:`openmc.SourceBase` object.
+        :volume_normalized_flux_tallies:
+            Whether to normalize flux tallies by volume (bool). The default
+            is 'False'. When enabled, flux tallies will be reported in units of
+            cm/cm^3. When disabled, flux tallies will be reported in units
+            of cm (i.e., total distance traveled by neutrons in the spatial
+            tally region).
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -1084,6 +1090,8 @@ def random_ray(self, random_ray: dict):
                                       random_ray[key], 0.0, True)
             elif key == 'ray_source':
                 cv.check_type('random ray source', random_ray[key], SourceBase)
+            elif key == 'volume_normalized_flux_tallies':
+                cv.check_type('volume normalized flux tallies', random_ray[key], bool)
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -1877,6 +1885,10 @@ def _random_ray_from_xml_element(self, root):
                 elif child.tag == 'source':
                     source = SourceBase.from_xml_element(child)
                     self.random_ray['ray_source'] = source
+                elif child.tag == 'volume_normalized_flux_tallies':
+                    self.random_ray['volume_normalized_flux_tallies'] = (
+                        child.text in ('true', '1')
+                    )
 
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.

src/random_ray/flat_source_domain.cpp

@@ -23,6 +23,9 @@ namespace openmc {
 // FlatSourceDomain implementation
 //==============================================================================
 
+// Static Variable Declarations
+bool FlatSourceDomain::volume_normalized_flux_tallies_ {false};
+
 FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
 {
   // Count the number of source regions, compute the cell offset
@@ -81,15 +84,17 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
   }
 
   // Initialize tally volumes
-  tally_volumes_.resize(model::tallies.size());
-  for (int i = 0; i < model::tallies.size(); i++) {
-    //  Get the shape of the 3D result tensor
-    auto shape = model::tallies[i]->results().shape();
-
-    // Create a new 2D tensor with the same size as the first
-    // two dimensions of the 3D tensor
-    tally_volumes_[i] =
-      xt::xtensor<double, 2>::from_shape({shape[0], shape[1]});
+  if (volume_normalized_flux_tallies_) {
+    tally_volumes_.resize(model::tallies.size());
+    for (int i = 0; i < model::tallies.size(); i++) {
+      //  Get the shape of the 3D result tensor
+      auto shape = model::tallies[i]->results().shape();
+
+      // Create a new 2D tensor with the same size as the first
+      // two dimensions of the 3D tensor
+      tally_volumes_[i] =
+        xt::xtensor<double, 2>::from_shape({shape[0], shape[1]});
+    }
   }
 
   // Compute simulation domain volume based on ray source
@@ -462,11 +467,63 @@ void FlatSourceDomain::convert_source_regions_to_tallies()
 // Set the volume accumulators to zero for all tallies
 void FlatSourceDomain::reset_tally_volumes()
 {
+  if (volume_normalized_flux_tallies_) {
 #pragma omp parallel for
-  for (int i = 0; i < tally_volumes_.size(); i++) {
-    auto& tensor = tally_volumes_[i];
-    tensor.fill(0.0); // Set all elements of the tensor to 0.0
+    for (int i = 0; i < tally_volumes_.size(); i++) {
+      auto& tensor = tally_volumes_[i];
+      tensor.fill(0.0); // Set all elements of the tensor to 0.0
+    }
+  }
+}
+
+// In fixed source mode, due to the way that volumetric fixed sources are
+// converted and applied as volumetric sources in one or more source regions,
+// we need to perform an additional normalization step to ensure that the
+// reported scalar fluxes are in units per source neutron. This allows for
+// direct comparison of reported tallies to Monte Carlo flux results.
+// This factor needs to be computed at each iteration, as it is based on the
+// volume estimate of each FSR, which improves over the course of the simulation
+double FlatSourceDomain::compute_fixed_source_normalization_factor() const
+{
+  // If we are not in fixed source mode, then there are no external sources
+  // so no normalization is needed.
+  if (settings::run_mode != RunMode::FIXED_SOURCE) {
+    return 1.0;
+  }
+
+  // Step 1 is to sum over all source regions and energy groups to get the
+  // total external source strength in the simulation.
+  double simulation_external_source_strength = 0.0;
+#pragma omp parallel for reduction(+ : simulation_external_source_strength)
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    int material = material_[sr];
+    double volume = volume_[sr] * simulation_volume_;
+    for (int e = 0; e < negroups_; e++) {
+      // Temperature and angle indices, if using multiple temperature
+      // data sets and/or anisotropic data sets.
+      // TODO: Currently assumes we are only using single temp/single
+      // angle data.
+      const int t = 0;
+      const int a = 0;
+      float sigma_t = data::mg.macro_xs_[material].get_xs(
+        MgxsType::TOTAL, e, nullptr, nullptr, nullptr, t, a);
+      simulation_external_source_strength +=
+        external_source_[sr * negroups_ + e] * sigma_t * volume;
+    }
+  }
+
+  // Step 2 is to determine the total user-specified external source strength
+  double user_external_source_strength = 0.0;
+  for (auto& ext_source : model::external_sources) {
+    user_external_source_strength += ext_source->strength();
   }
+
+  // The correction factor is the ratio of the user-specified external source
+  // strength to the simulation external source strength.
+  double source_normalization_factor =
+    user_external_source_strength / simulation_external_source_strength;
+
+  return source_normalization_factor;
 }
 
 // Tallying in random ray is not done directly during transport, rather,
@@ -492,6 +549,9 @@ void FlatSourceDomain::random_ray_tally()
   const int t = 0;
   const int a = 0;
 
+  double source_normalization_factor =
+    compute_fixed_source_normalization_factor();
+
 // We loop over all source regions and energy groups. For each
 // element, we check if there are any scores needed and apply
 // them.
@@ -510,7 +570,7 @@ void FlatSourceDomain::random_ray_tally()
     double material = material_[sr];
     for (int g = 0; g < negroups_; g++) {
       int idx = sr * negroups_ + g;
-      double flux = scalar_flux_new_[idx];
+      double flux = scalar_flux_new_[idx] * source_normalization_factor;
 
       // Determine numerical score value
       for (auto& task : tally_task_[idx]) {
@@ -561,11 +621,13 @@ void FlatSourceDomain::random_ray_tally()
     // For flux tallies, the total volume of the spatial region is needed
     // for normalizing the flux. We store this volume in a separate tensor.
     // We only contribute to each volume tally bin once per FSR.
-    for (const auto& task : volume_task_[sr]) {
-      if (task.score_type == SCORE_FLUX) {
+    if (volume_normalized_flux_tallies_) {
+      for (const auto& task : volume_task_[sr]) {
+        if (task.score_type == SCORE_FLUX) {
 #pragma omp atomic
-        tally_volumes_[task.tally_idx](task.filter_idx, task.score_idx) +=
-          volume;
+          tally_volumes_[task.tally_idx](task.filter_idx, task.score_idx) +=
+            volume;
+        }
       }
     }
   } // end FSR loop
@@ -575,16 +637,18 @@ void FlatSourceDomain::random_ray_tally()
   // and then scores. For each score, we check the tally data structure to
   // see what index that score corresponds to. If that score is a flux score,
   // then we divide it by volume.
-  for (int i = 0; i < model::tallies.size(); i++) {
-    Tally& tally {*model::tallies[i]};
+  if (volume_normalized_flux_tallies_) {
+    for (int i = 0; i < model::tallies.size(); i++) {
+      Tally& tally {*model::tallies[i]};
 #pragma omp parallel for
-    for (int bin = 0; bin < tally.n_filter_bins(); bin++) {
-      for (int score_idx = 0; score_idx < tally.n_scores(); score_idx++) {
-        auto score_type = tally.scores_[score_idx];
-        if (score_type == SCORE_FLUX) {
-          double vol = tally_volumes_[i](bin, score_idx);
-          if (vol > 0.0) {
-            tally.results_(bin, score_idx, TallyResult::VALUE) /= vol;
+      for (int bin = 0; bin < tally.n_filter_bins(); bin++) {
+        for (int score_idx = 0; score_idx < tally.n_scores(); score_idx++) {
+          auto score_type = tally.scores_[score_idx];
+          if (score_type == SCORE_FLUX) {
+            double vol = tally_volumes_[i](bin, score_idx);
+            if (vol > 0.0) {
+              tally.results_(bin, score_idx, TallyResult::VALUE) /= vol;
+            }
           }
         }
       }
@@ -767,6 +831,9 @@ void FlatSourceDomain::output_to_vtk() const
       }
     }
 
+    double source_normalization_factor =
+      compute_fixed_source_normalization_factor();
+
     // Open file for writing
     std::FILE* plot = std::fopen(filename.c_str(), "wb");
 
@@ -786,7 +853,8 @@ void FlatSourceDomain::output_to_vtk() const
       std::fprintf(plot, "LOOKUP_TABLE default\n");
       for (int fsr : voxel_indices) {
         int64_t source_element = fsr * negroups_ + g;
-        float flux = scalar_flux_final_[source_element];
+        float flux =
+          scalar_flux_final_[source_element] * source_normalization_factor;
         flux /= (settings::n_batches - settings::n_inactive);
         flux = convert_to_big_endian<float>(flux);
         std::fwrite(&flux, sizeof(float), 1, plot);
@@ -819,7 +887,8 @@ void FlatSourceDomain::output_to_vtk() const
       int mat = material_[fsr];
       for (int g = 0; g < negroups_; g++) {
         int64_t source_element = fsr * negroups_ + g;
-        float flux = scalar_flux_final_[source_element];
+        float flux =
+          scalar_flux_final_[source_element] * source_normalization_factor;
         flux /= (settings::n_batches - settings::n_inactive);
         float Sigma_f = data::mg.macro_xs_[mat].get_xs(
           MgxsType::FISSION, g, nullptr, nullptr, nullptr, 0, 0);

src/settings.cpp

@@ -269,6 +269,10 @@ void get_run_parameters(pugi::xml_node node_base)
     } else {
       fatal_error("Specify random ray source in settings XML");
     }
+    if (check_for_node(random_ray_node, "volume_normalized_flux_tallies")) {
+      FlatSourceDomain::volume_normalized_flux_tallies_ =
+        get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies");
+    }
   }
 }
 

tests/regression_tests/random_ray_basic/inputs_true.dat

@@ -85,6 +85,7 @@
           <parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
         </space>
       </source>
+      <volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
     </random_ray>
   </settings>
   <tallies>

tests/regression_tests/random_ray_basic/test.py

@@ -190,6 +190,7 @@ def random_ray_model() -> openmc.Model:
     settings.random_ray['distance_active'] = 100.0
     settings.random_ray['distance_inactive'] = 20.0
     settings.random_ray['ray_source'] = rr_source
+    settings.random_ray['volume_normalized_flux_tallies'] = True
 
     ###############################################################################
     # Define tallies