Random Ray Adjoint Source Logic Improvement (#3325)

src/random_ray/flat_source_domain.cpp

@@ -1049,14 +1049,21 @@ void FlatSourceDomain::flatten_xs()
 
 void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
 {
-  // Set the external source to 1/forward_flux
-  // The forward flux is given in terms of total for the forward simulation
-  // so we must convert it to a "per batch" quantity
+  // 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.
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions_; sr++) {
     for (int g = 0; g < negroups_; g++) {
-      source_regions_.external_source(sr, g) =
-        1.0 / forward_flux[sr * negroups_ + g];
+      double flux = forward_flux[sr * negroups_ + g];
+      if (flux <= 0.0) {
+        source_regions_.external_source(sr, g) = 0.0;
+      } else {
+        source_regions_.external_source(sr, g) = 1.0 / flux;
+      }
     }
   }
 
@@ -1064,12 +1071,12 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
   // iteration)
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions_; sr++) {
+    int material = source_regions_.material(sr);
+    if (material == MATERIAL_VOID) {
+      continue;
+    }
     for (int g = 0; g < negroups_; g++) {
-      int material = source_regions_.material(sr);
-      if (material == MATERIAL_VOID) {
-        continue;
-      }
-      double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g];
+      double sigma_t = sigma_t_[material * negroups_ + g];
       source_regions_.external_source(sr, g) /= sigma_t;
     }
   }

tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat

@@ -191,7 +191,7 @@
   </geometry>
   <settings>
     <run_mode>fixed source</run_mode>
-    <particles>90</particles>
+    <particles>500</particles>
     <batches>10</batches>
     <inactive>5</inactive>
     <source particle="neutron" strength="3.14" type="independent">
@@ -214,6 +214,7 @@
       </source>
       <volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
       <adjoint>True</adjoint>
+      <volume_estimator>naive</volume_estimator>
     </random_ray>
   </settings>
   <tallies>

tests/regression_tests/random_ray_adjoint_fixed_source/results_true.dat

@@ -1,9 +1,9 @@
 tally 1:
--7.235364E+03
-3.367109E+09
+5.790516E+04
+6.740859E+08
 tally 2:
-4.818311E+05
-6.269371E+10
+6.885455E+04
+9.482551E+08
 tally 3:
-1.515641E+06
-4.598791E+11
+1.956327E+05
+7.654469E+09

tests/regression_tests/random_ray_adjoint_fixed_source/test.py

@@ -16,5 +16,7 @@ class MGXSTestHarness(TolerantPyAPITestHarness):
 def test_random_ray_adjoint_fixed_source():
     model = random_ray_three_region_cube()
     model.settings.random_ray['adjoint'] = True
+    model.settings.random_ray['volume_estimator'] = 'naive'
+    model.settings.particles = 500
     harness = MGXSTestHarness('statepoint.10.h5', model)
     harness.main()