Random Ray Eigenvalue Flux Normalization Change

include/openmc/random_ray/flat_source_domain.h

@@ -170,6 +170,9 @@ class FlatSourceDomain {
     simulation_volume_; // Total physical volume of the simulation domain, as
                         // defined by the 3D box of the random ray source
 
+  double
+    fission_rate_; // The system's fission rate (per cm^3), in eigenvalue mode
+
   // Volumes for each tally and bin/score combination. This intermediate data
   // structure is used when tallying quantities that must be normalized by
   // volume (i.e., flux). The vector is index by tally index, while the inner 2D

src/eigenvalue.cpp

@@ -403,6 +403,16 @@ void calculate_average_keff()
         t_value *
         std::sqrt(
           (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1));
+
+      // In some cases (such as an infinite medium problem), random ray
+      // may estimate k exactly and in an unvarying manner between iterations.
+      // In this case, the floating point roundoff between the division and the
+      // power operations may cause an extremely small negative value to occur
+      // inside the sqrt operation, leading to NaN. If this occurs, we check for
+      // it and set the std dev to zero.
+      if (!std::isfinite(simulation::keff_std)) {
+        simulation::keff_std = 0.0;
+      }
     }
   }
 }

src/random_ray/flat_source_domain.cpp

@@ -124,7 +124,9 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
         double chi = chi_[material * negroups_ + g_out];
 
         scatter_source += sigma_s * scalar_flux;
-        fission_source += nu_sigma_f * scalar_flux * chi;
+        if (settings::create_fission_neutrons) {
+          fission_source += nu_sigma_f * scalar_flux * chi;
+        }
       }
       srh.source(g_out) =
         (scatter_source + fission_source * inverse_k_eff) / sigma_t;
@@ -369,6 +371,7 @@ void FlatSourceDomain::compute_k_eff()
   // Adds entropy value to shared entropy vector in openmc namespace.
   simulation::entropy.push_back(H);
 
+  fission_rate_ = fission_rate_new;
   k_eff_ = k_eff_new;
 }
 
@@ -519,12 +522,33 @@ void FlatSourceDomain::reset_tally_volumes()
 // 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 || adjoint_) {
+  // Eigenvalue mode normalization
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    // Normalize fluxes by total number of fission neutrons produced. This
+    // ensures consistent scaling of the eigenvector such that its magnitude is
+    // comparable to the eigenvector produced by the Monte Carlo solver.
+    // Multiplying by the eigenvalue is unintuitive, but it is necessary.
+    // If the eigenvalue is 1.2, per starting source neutron, you will
+    // generate 1.2 neutrons. Thus if we normalize to generating only ONE
+    // neutron in total for the whole domain, then we don't actually have enough
+    // flux to generate the required 1.2 neutrons. We only know the flux
+    // required to generate 1 neutron (which would have required less than one
+    // starting neutron). Thus, you have to scale the flux up by the eigenvalue
+    // such that 1.2 neutrons are generated, so as to be consistent with the
+    // bookkeeping in MC which is all done per starting source neutron (not per
+    // neutron produced).
+    return k_eff_ / (fission_rate_ * simulation_volume_);
+  }
+
+  // If we are in adjoint mode of a fixed source problem, the external
+  // source is already normalized, such that all resulting fluxes are
+  // also normalized.
+  if (adjoint_) {
     return 1.0;
   }
 
+  // Fixed source mode normalization
+
   // 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;

src/random_ray/linear_source_domain.cpp

@@ -68,9 +68,11 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
 
         // Compute source terms for flat and linear components of the flux
         scatter_flat += sigma_s * flux_flat;
-        fission_flat += nu_sigma_f * flux_flat * chi;
         scatter_linear += sigma_s * flux_linear;
-        fission_linear += nu_sigma_f * flux_linear * chi;
+        if (settings::create_fission_neutrons) {
+          fission_flat += nu_sigma_f * flux_flat * chi;
+          fission_linear += nu_sigma_f * flux_linear * chi;
+        }
       }
 
       // Compute the flat source term

tests/regression_tests/random_ray_adjoint_k_eff/results_true.dat

@@ -1,171 +1,171 @@
 k-combined:
 1.006640E+00 1.812969E-03
 tally 1:
-6.684129E+00
-8.939821E+00
-2.685967E+00
-1.443592E+00
+1.208044E+00
+2.920182E-01
+4.854426E-01
+4.715453E-02
 0.000000E+00
 0.000000E+00
-6.358774E+00
-8.091444E+00
-9.687217E-01
-1.878029E-01
+1.149242E+00
+2.643067E-01
+1.750801E-01
+6.134563E-03
 0.000000E+00
 0.000000E+00
-5.963160E+00
-7.117108E+00
-1.932332E-01
-7.473914E-03
+1.077743E+00
+2.324814E-01
+3.492371E-02
+2.441363E-04
 0.000000E+00
 0.000000E+00
-5.137593E+00
-5.283310E+00
-1.714616E-01
-5.884834E-03
-1.086218E-06
-2.361752E-13
-4.857253E+00
-4.719856E+00
-5.689580E-02
-6.476286E-04
-2.989356E-03
-1.787808E-06
-4.830516E+00
-4.666801E+00
-7.203015E-03
-1.037676E-05
-3.620020E+00
-2.620927E+00
-5.161382E+00
-5.328124E+00
-6.786255E-02
-9.210763E-04
-5.531943E+00
-6.120553E+00
-5.414034E+00
-5.864661E+00
+9.285362E-01
+1.725808E-01
+3.098889E-02
+1.922297E-04
+1.963161E-07
+7.714727E-15
+8.778641E-01
+1.541719E-01
+1.028293E-02
+2.115448E-05
+5.402741E-04
+5.839789E-08
+8.730274E-01
+1.524358E-01
+1.301813E-03
+3.389450E-07
+6.542525E-01
+8.560964E-02
+9.328247E-01
+1.740366E-01
+1.226491E-02
+3.008584E-05
+9.997969E-01
+1.999204E-01
+9.784958E-01
+1.915688E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.632338E+00
-6.347626E+00
+1.017952E+00
+2.073461E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.682608E+00
-6.462382E+00
+1.027039E+00
+2.110955E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.310716E+00
-5.645180E+00
+9.598240E-01
+1.844004E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-4.945409E+00
-4.893171E+00
+8.937969E-01
+1.598332E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-4.842688E+00
-4.690352E+00
+8.752275E-01
+1.532052E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-5.117198E+00
-5.237280E+00
+9.248400E-01
+1.710699E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-6.938711E+00
-9.633345E+00
-2.835258E+00
-1.608212E+00
+1.254054E+00
+3.146708E-01
+5.124223E-01
+5.253093E-02
 0.000000E+00
 0.000000E+00
-6.549505E+00
-8.584036E+00
-1.015138E+00
-2.061993E-01
+1.183712E+00
+2.803961E-01
+1.834683E-01
+6.735381E-03
 0.000000E+00
 0.000000E+00
-6.050651E+00
-7.327711E+00
-1.992816E-01
-7.948424E-03
+1.093555E+00
+2.393604E-01
+3.601678E-02
+2.596341E-04
 0.000000E+00
 0.000000E+00
-5.113981E+00
-5.234801E+00
-1.732323E-01
-6.006619E-03
-1.097435E-06
-2.410627E-13
-4.837033E+00
-4.680541E+00
-5.760042E-02
-6.637112E-04
-3.026377E-03
-1.832205E-06
-4.827049E+00
-4.660105E+00
-7.319913E-03
-1.071647E-05
-3.678770E+00
-2.706730E+00
-5.175337E+00
-5.356957E+00
-6.923046E-02
-9.586177E-04
-5.643451E+00
-6.370016E+00
-6.693323E+00
-8.964322E+00
-2.753307E+00
-1.516683E+00
+9.242694E-01
+1.709967E-01
+3.130894E-02
+1.962084E-04
+1.983437E-07
+7.874400E-15
+8.742100E-01
+1.528879E-01
+1.041026E-02
+2.167970E-05
+5.469644E-04
+5.984780E-08
+8.724009E-01
+1.522171E-01
+1.322938E-03
+3.500386E-07
+6.648691E-01
+8.841161E-02
+9.353464E-01
+1.749781E-01
+1.251209E-02
+3.131171E-05
+1.019947E+00
+2.080664E-01
+1.209708E+00
+2.928214E-01
+4.976146E-01
+4.954258E-02
 0.000000E+00
 0.000000E+00
-6.358384E+00
-8.090233E+00
-9.912008E-01
-1.965868E-01
+1.149174E+00
+2.642694E-01
+1.791431E-01
+6.421530E-03
 0.000000E+00
 0.000000E+00
-5.957484E+00
-7.103246E+00
-1.974033E-01
-7.798286E-03
+1.076718E+00
+2.320295E-01
+3.567737E-02
+2.547314E-04
 0.000000E+00
 0.000000E+00
-5.130744E+00
-5.268844E+00
-1.749233E-01
-6.123348E-03
-1.108148E-06
-2.457474E-13
-4.857340E+00
-4.720019E+00
-5.816659E-02
-6.768049E-04
-3.056125E-03
-1.868351E-06
-4.830629E+00
-4.667018E+00
-7.366289E-03
-1.085264E-05
-3.702077E+00
-2.741125E+00
-5.164864E+00
-5.335279E+00
-6.947917E-02
-9.655086E-04
-5.663725E+00
-6.415806E+00
+9.272977E-01
+1.721077E-01
+3.161443E-02
+2.000179E-04
+2.002789E-07
+8.027290E-15
+8.778794E-01
+1.541769E-01
+1.051257E-02
+2.210729E-05
+5.523400E-04
+6.102818E-08
+8.730477E-01
+1.524429E-01
+1.331320E-03
+3.544869E-07
+6.690816E-01
+8.953516E-02
+9.334544E-01
+1.742707E-01
+1.255707E-02
+3.153710E-05
+1.023613E+00
+2.095641E-01

tests/regression_tests/random_ray_auto_convert/material_wise/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-7.356667E-01 6.637270E-03
+7.356667E-01 6.637266E-03

tests/regression_tests/random_ray_auto_convert/stochastic_slab/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-7.551716E-01 8.117378E-03
+7.551716E-01 8.117379E-03