From f31130d36798aef0ca7bed6a73e24a9281c6c51c Mon Sep 17 00:00:00 2001
From: John Tramm <john.tramm@gmail.com>
Date: Wed, 18 Jun 2025 06:34:56 -0500
Subject: [PATCH] Prevent Adjoint Sources from Trending towards Infinity
 (#3449)

---
 include/openmc/constants.h            |  5 ++++
 src/random_ray/flat_source_domain.cpp | 43 +++++++++++++++++++++------
 2 files changed, 39 insertions(+), 9 deletions(-)

diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index 252194528c..ae70560795 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -63,6 +63,11 @@ constexpr int MAX_SAMPLE {100000};
 // source region in the random ray solver
 constexpr double MIN_HITS_PER_BATCH {1.5};
 
+// The minimum flux value to be considered non-zero when computing adjoint
+// sources. Positive values below this cutoff will be treated as zero, so as to
+// prevent extremely large adjoint source terms from being generated.
+constexpr double ZERO_FLUX_CUTOFF {1e-22};
+
 // ============================================================================
 // MATH AND PHYSICAL CONSTANTS
 
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index 9df6513aee..7404a0b29b 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@@ -304,6 +304,13 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
           set_flux_to_source(sr, g);
         }
       }
+      // Halt if NaN implosion is detected
+      if (!std::isfinite(source_regions_.scalar_flux_new(sr, g))) {
+        fatal_error("A source region scalar flux is not finite. "
+                    "This indicates a numerical instability in the "
+                    "simulation. Consider increasing ray density or adjusting "
+                    "the source region mesh.");
+      }
     }
   }
 
@@ -1154,9 +1161,9 @@ void FlatSourceDomain::flatten_xs()
           m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a);
         sigma_t_.push_back(sigma_t);
 
-        double nu_Sigma_f =
+        double nu_sigma_f =
           m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a);
-        nu_sigma_f_.push_back(nu_Sigma_f);
+        nu_sigma_f_.push_back(nu_sigma_f);
 
         double sigma_f =
           m.get_xs(MgxsType::FISSION, g_out, NULL, NULL, NULL, t, a);
@@ -1164,6 +1171,11 @@ void FlatSourceDomain::flatten_xs()
 
         double chi =
           m.get_xs(MgxsType::CHI_PROMPT, g_out, &g_out, NULL, NULL, t, a);
+        if (!std::isfinite(chi)) {
+          // MGXS interface may return NaN in some cases, such as when material
+          // is fissionable but has very small sigma_f.
+          chi = 0.0;
+        }
         chi_.push_back(chi);
 
         for (int g_in = 0; g_in < negroups_; g_in++) {
@@ -1191,17 +1203,30 @@ void FlatSourceDomain::flatten_xs()
 
 void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
 {
-  // Set the external source to 1/forward_flux. If the forward flux is negative
-  // or zero, set the adjoint source to zero, as this is likely a very small
-  // source region that we don't need to bother trying to vector particles
-  // towards. Flux negativity in random ray is not related to the flux being
-  // small in magnitude, but rather due to the source region being physically
-  // small in volume and thus having a noisy flux estimate.
+  // Set the adjoint external source to 1/forward_flux. If the forward flux is
+  // negative, zero, or extremely close to zero, set the adjoint source to zero,
+  // as this is likely a very small source region that we don't need to bother
+  // trying to vector particles towards. In the case of flux "being extremely
+  // close to zero", we define this as being a fixed fraction of the maximum
+  // forward flux, below which we assume the flux would be physically
+  // undetectable.
+
+  // First, find the maximum forward flux value
+  double max_flux = 0.0;
+#pragma omp parallel for reduction(max : max_flux)
+  for (int64_t se = 0; se < n_source_elements(); se++) {
+    double flux = forward_flux[se];
+    if (flux > max_flux) {
+      max_flux = flux;
+    }
+  }
+
+  // Then, compute the adjoint source for each source region
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions(); sr++) {
     for (int g = 0; g < negroups_; g++) {
       double flux = forward_flux[sr * negroups_ + g];
-      if (flux <= 0.0) {
+      if (flux <= ZERO_FLUX_CUTOFF * max_flux) {
         source_regions_.external_source(sr, g) = 0.0;
       } else {
         source_regions_.external_source(sr, g) = 1.0 / flux;