From a5b26de041003d30ee097e83c5540f72f912bb67 Mon Sep 17 00:00:00 2001
From: Olek <45364492+yardasol@users.noreply.github.com>
Date: Mon, 17 Feb 2025 11:02:51 -0600
Subject: [PATCH] Random ray consistency changes (#3298)

---
 openmc/settings.py                            | 20 +++++++-------
 src/random_ray/flat_source_domain.cpp         | 27 +++++--------------
 .../random_ray_adjoint_k_eff/results_true.dat |  2 +-
 3 files changed, 17 insertions(+), 32 deletions(-)

diff --git a/openmc/settings.py b/openmc/settings.py
index 77598b204f..b29393d9ac 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -1110,27 +1110,27 @@ class Settings:
         if not isinstance(random_ray, Mapping):
             raise ValueError(f'Unable to set random_ray from "{random_ray}" '
                              'which is not a dict.')
-        for key in random_ray:
+        for key, value in random_ray.items():
             if key == 'distance_active':
-                cv.check_type('active ray length', random_ray[key], Real)
-                cv.check_greater_than('active ray length', random_ray[key], 0.0)
+                cv.check_type('active ray length', value, Real)
+                cv.check_greater_than('active ray length', value, 0.0)
             elif key == 'distance_inactive':
-                cv.check_type('inactive ray length', random_ray[key], Real)
+                cv.check_type('inactive ray length', value, Real)
                 cv.check_greater_than('inactive ray length',
-                                      random_ray[key], 0.0, True)
+                                      value, 0.0, True)
             elif key == 'ray_source':
-                cv.check_type('random ray source', random_ray[key], SourceBase)
+                cv.check_type('random ray source', value, SourceBase)
             elif key == 'volume_estimator':
-                cv.check_value('volume estimator', random_ray[key],
+                cv.check_value('volume estimator', value,
                                ('naive', 'simulation_averaged',
                                 'hybrid'))
             elif key == 'source_shape':
-                cv.check_value('source shape', random_ray[key],
+                cv.check_value('source shape', value,
                                ('flat', 'linear', 'linear_xy'))
             elif key == 'volume_normalized_flux_tallies':
-                cv.check_type('volume normalized flux tallies', random_ray[key], bool)
+                cv.check_type('volume normalized flux tallies', value, bool)
             elif key == 'adjoint':
-                cv.check_type('adjoint', random_ray[key], bool)
+                cv.check_type('adjoint', value, bool)
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index 39a50ace60..098fd17b08 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@@ -118,7 +118,7 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
 
   double inverse_k_eff = 1.0 / k_eff;
 
-  // Add scattering source
+  // Add scattering + fission source
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions_; sr++) {
     int material = source_regions_.material(sr);
@@ -126,35 +126,20 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
     for (int g_out = 0; g_out < negroups_; g_out++) {
       double sigma_t = sigma_t_[material * negroups_ + g_out];
       double scatter_source = 0.0;
+      double fission_source = 0.0;
 
       for (int g_in = 0; g_in < negroups_; g_in++) {
         double scalar_flux = source_regions_.scalar_flux_old(sr, g_in);
         double sigma_s =
           sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
-        scatter_source += sigma_s * scalar_flux;
-      }
-
-      source_regions_.source(sr, g_out) = scatter_source / sigma_t;
-    }
-  }
-
-  // Add fission source
-#pragma omp parallel for
-  for (int64_t sr = 0; sr < n_source_regions_; sr++) {
-    int material = source_regions_.material(sr);
-
-    for (int g_out = 0; g_out < negroups_; g_out++) {
-      double sigma_t = sigma_t_[material * negroups_ + g_out];
-      double fission_source = 0.0;
-
-      for (int g_in = 0; g_in < negroups_; g_in++) {
-        double scalar_flux = source_regions_.scalar_flux_old(sr, g_in);
         double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
         double chi = chi_[material * negroups_ + g_out];
+
+        scatter_source += sigma_s * scalar_flux;
         fission_source += nu_sigma_f * scalar_flux * chi;
       }
-      source_regions_.source(sr, g_out) +=
-        fission_source * inverse_k_eff / sigma_t;
+      source_regions_.source(sr, g_out) =
+        (scatter_source + fission_source * inverse_k_eff) / sigma_t;
     }
   }
 
diff --git a/tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat b/tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat
index 1690d46e96..657c841b56 100644
--- a/tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat
+++ b/tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat
@@ -1,5 +1,5 @@
 k-combined:
-1.006640E+00 1.812967E-03
+1.006640E+00 1.812969E-03
 tally 1:
 6.684129E+00
 8.939821E+00