Alternative Random Ray Volume Estimators (openmc-dev#3060)

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

docs/source/methods/random_ray.rst

@@ -411,7 +411,7 @@ which when partially simplified becomes:
 
 Note that there are now four (seemingly identical) volume terms in this equation.
 
-.. _methods-volume-dilemma:
+.. _methods_random_ray_vol:
 
 ~~~~~~~~~~~~~~
 Volume Dilemma
@@ -440,9 +440,11 @@ features stochastic variables (the sums over random ray lengths and angular
 fluxes) in both the numerator and denominator, making it a stochastic ratio
 estimator, which is inherently biased. In practice, usage of the naive estimator
 does result in a biased, but "consistent"  estimator (i.e., it is biased, but
-the bias tends towards zero as the sample size increases). Experimentally, the
-right answer can be obtained with this estimator, though a very fine ray density
-is required to eliminate the bias.
+the bias tends towards zero as the sample size increases). Empirically, this
+bias tends to effect eigenvalue calculations much more significantly than in
+fixed source simulations. Experimentally, the right answer can be obtained with
+this estimator, though for eigenvalue simulations a very fine ray density is
+required to eliminate the bias.
 
 How might we solve the biased ratio estimator problem? While there is no obvious
 way to alter the numerator term (which arises from the characteristic
@@ -463,17 +465,17 @@ replace the actual tracklength that was accumulated inside that FSR each
 iteration with the expected value.
 
 If we know the analytical volumes, then those can be used to directly compute
-the expected value of the tracklength in each cell. However, as the analytical
-volumes are not typically known in OpenMC due to the usage of user-defined
-constructive solid geometry, we need to source this quantity from elsewhere. An
-obvious choice is to simply accumulate the total tracklength through each FSR
-across all iterations (batches) and to use that sum to compute the expected
-average length per iteration, as:
+the expected value of the tracklength in each cell, :math:`L_{avg}`. However, as
+the analytical volumes are not typically known in OpenMC due to the usage of
+user-defined constructive solid geometry, we need to source this quantity from
+elsewhere. An obvious choice is to simply accumulate the total tracklength
+through each FSR across all iterations (batches) and to use that sum to compute
+the expected average length per iteration, as:
 
 .. math::
-    :label: sim_estimator
+    :label: L_avg
 
-       \sum\limits^{}_{i} \ell_i \approx \frac{\sum\limits^{B}_{b}\sum\limits^{N_i}_{r} \ell_{b,r} }{B}
+    \sum\limits^{}_{i} \ell_i \approx L_{avg} = \frac{\sum\limits^{B}_{b}\sum\limits^{N_i}_{r=1} \ell_{b,r} }{B}
 
 where :math:`b` is a single batch in :math:`B` total batches simulated so far.
 
@@ -486,7 +488,7 @@ averaged" estimator is therefore:
 .. math::
     :label: phi_sim
 
-    \phi_{i,g}^{simulation} = \frac{Q_{i,g} }{\Sigma_{t,i,g}} + \frac{\sum\limits_{r=1}^{N_i} \Delta \psi_{r,g}}{\Sigma_{t,i,g} \frac{\sum\limits^{B}_{b}\sum\limits^{N_i}_{r} \ell_{b,r} }{B}}
+    \phi_{i,g}^{simulation} = \frac{Q_{i,g} }{\Sigma_{t,i,g}} + \frac{\sum\limits_{r=1}^{N_i} \Delta \psi_{r,g}}{\Sigma_{t,i,g} L_{avg}}
 
 In practical terms, the "simulation averaged" estimator is virtually
 indistinguishable numerically from use of the true analytical volume to estimate
@@ -500,17 +502,81 @@ in which case the denominator served as a normalization term for the numerator
 integral in Equation :eq:`integral`. Essentially, we have now used a different
 term for the volume in the numerator as compared to the normalizing volume in
 the denominator. The inevitable mismatch (due to noise) between these two
-quantities results in a significant increase in variance. Notably, the same
-problem occurs if using a tracklength estimate based on the analytical volume,
-as again the numerator integral and the normalizing denominator integral no
-longer match on a per-iteration basis.
-
-In practice, the simulation averaged method does completely remove the bias,
-though at the cost of a notable increase in variance. Empirical testing reveals
-that on most problems, the simulation averaged estimator does win out overall in
-numerical performance, as a much coarser quadrature can be used resulting in
-faster runtimes overall. Thus, OpenMC uses the simulation averaged estimator in
-its random ray mode.
+quantities results in a significant increase in variance, and can even result in
+the generation of negative fluxes. Notably, the same problem occurs if using a
+tracklength estimate based on the analytical volume, as again the numerator
+integral and the normalizing denominator integral no longer match on a
+per-iteration basis.
+
+In practice, the simulation averaged method does completely remove the bias seen
+when using the naive estimator, though at the cost of a notable increase in
+variance. Empirical testing reveals that on most eigenvalue problems, the
+simulation averaged estimator does win out overall in numerical performance, as
+a much coarser quadrature can be used resulting in faster runtimes overall.
+Thus, OpenMC uses the simulation averaged estimator as default in its random ray
+mode for eigenvalue solves.
+
+OpenMC also features a "hybrid" volume estimator that uses the naive estimator
+for all regions containing an external (fixed) source term. For all other
+source regions, the "simulation averaged" estimator is used. This typically achieves
+a best of both worlds result, with the benefits of the low bias simulation averaged
+estimator in most regions, while preventing instability and/or large biases in regions
+with external source terms via use of the naive estimator. In general, it is
+recommended to use the "hybrid" estimator, which is the default method used
+in OpenMC. If instability is encountered despite high ray densities, then
+the naive estimator may be preferable.
+
+A table that summarizes the pros and cons, as well as recommendations for
+different use cases, is given in the :ref:`volume
+estimators<usersguide_vol_estimators>` section of the user guide.
+
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+What Happens When a Source Region is Missed?
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Given the stochastic nature of random ray, when low ray densities are used it is
+common for small source regions to occasionally not be hit by any rays in a
+particular power iteration :math:`n`. This naturally collapses the flux estimate
+in that cell for the iteration from Equation :eq:`phi_naive` to:
+
+.. math::
+    :label: phi_missed_one
+
+    \phi_{i,g,n}^{missed} = \frac{Q_{i,g,n} }{\Sigma_{t,i,g}} 
+
+as the streaming operator has gone to zero. While this is obviously innacurate
+as it ignores transport, for most problems where the region is only occasionally
+missed this estimator does not tend to introduce any significant bias.
+
+However, in cases where the total cross section in the region is very small
+(e.g., a void-like material) and where a strong external fixed source has been
+placed, then this treatment causes major issues. In this pathological case, the
+lack of transport forces the entirety of the fixed source to effectively be
+contained and collided within the cell, which for a low cross section region is
+highly unphysical. The net effect is that a very high estimate of the flux
+(often orders of magnitude higher than is expected) is generated that iteration,
+which cannot be washed out even with hundreds or thousands of iterations. Thus,
+huge biases are often seen in spatial tallies containing void-like regions with
+external sources unless a high enough ray density is used such that all source
+regions are always hit each iteration. This is particularly problematic as
+external sources placed in void-like regions are very common in many types of
+fixed source analysis.
+
+For regions where external sources are present, to eliminate this bias it is
+therefore preferable to simply use the previous iteration's estimate of the flux
+in that cell, as:
+
+.. math::
+    :label: phi_missed_two
+
+    \phi_{i,g,n}^{missed} = \phi_{i,g,n-1} .
+
+When linear sources are present, the flux moments from the previous iteration
+are used in the same manner. While this introduces some small degree of
+correlation to the simulation, for miss rates on the order of a few percent the
+correlations are trivial and the bias is eliminated. Thus, in OpenMC the
+previous iteration's scalar flux estimate is applied to cells that are missed
+where there is an external source term present within the cell.
 
 ~~~~~~~~~~~~~~~
 Power Iteration
@@ -563,15 +629,15 @@ total spatial- and energy-integrated fission rate :math:`F^{n-1}` in iteration
 
 Notably, the volume term :math:`V_i` appears in the eigenvalue update equation.
 The same logic applies to the treatment of this term as was discussed earlier.
-In OpenMC, we use the "simulation averaged" volume derived from summing over all
-ray tracklength contributions to a FSR over all iterations and dividing by the
-total integration tracklength to date. Thus, Equation :eq:`fission_source`
-becomes:
+In OpenMC, we use the "simulation averaged" volume (Equation :eq:`L_avg`)
+derived from summing over all ray tracklength contributions to a FSR over all
+iterations and dividing by the total integration tracklength to date. Thus,
+Equation :eq:`fission_source` becomes:
 
 .. math::
     :label: fission_source_volumed
 
-    F^n = \sum\limits^{M}_{i} \left( \frac{\sum\limits^{B}_{b}\sum\limits^{N_i}_{r} \ell_{b,r} }{B} \sum\limits^{G}_{g} \nu \Sigma_f(i, g) \phi^{n}(g) \right)
+    F^n = \sum\limits^{M}_{i} \left( L_{avg} \sum\limits^{G}_{g} \nu \Sigma_f(i, g) \phi^{n}(g) \right)
 
 and a similar substitution can be made to update Equation
 :eq:`fission_source_prev` . In OpenMC, the most up-to-date version of the volume
@@ -965,7 +1031,7 @@ The Shannon entropy is then computed normally as
 
 where :math:`N` is the number of FSRs. FSRs with no fission source (or,
 occassionally, negative fission source, :ref:`due to the volume estimator
-problem <methods-volume-dilemma>`) are skipped to avoid taking an undefined
+problem <methods_random_ray_vol>`) are skipped to avoid taking an undefined
 logarithm in :eq:`shannon-entropy-random-ray`.
 
 .. _usersguide_fixed_source_methods:

docs/source/usersguide/random_ray.rst

@@ -535,6 +535,64 @@ points of 1.0e-2 and 1.0e1.
     # Add fixed source and ray sampling source to settings file
     settings.source = [neutron_source]
 
+.. _usersguide_vol_estimators:
+
+-----------------------------
+Alternative Volume Estimators
+-----------------------------
+
+As discussed in the random ray theory section on :ref:`volume estimators
+<methods_random_ray_vol>`, there are several possible derivations for the scalar
+flux estimate. These options deal with different ways of treating the
+accumulation over ray lengths crossing each FSR (a quantity directly
+proportional to volume), which can be computed using several methods. The
+following methods are currently available in OpenMC:
+
+.. list-table:: Comparison of Estimators
+   :header-rows: 1
+   :widths: 10 30 30 30
+
+   * - Estimator
+     - Description
+     - Pros
+     - Cons
+   * - ``simulation_averaged``
+     - Accumulates total active ray lengths in each FSR over all iterations,
+       improving the estimate of the volume in each cell each iteration. 
+     - * Virtually unbiased after several iterations
+       * Asymptotically approaches the true analytical volume
+       * Typically most efficient in terms of speed vs. accuracy
+     - * Higher variance
+       * Can lead to negative fluxes and numerical instability in pathological
+         cases
+   * - ``naive``
+     - Treats the volume as composed only of the active ray length through each
+       FSR per iteration, being a biased but numerically consistent ratio
+       estimator.
+     - * Low variance
+       * Unlikely to result in negative fluxes
+       * Recommended in cases where the simulation averaged estimator is
+         unstable
+     - * Biased estimator
+       * Requires more rays or longer active ray length to mitigate bias
+   * - ``hybrid`` (default)
+     - Applies the naive estimator to all cells that contain an external (fixed)
+       source contribution. Applies the simulation averaged estimator to all
+       other cells.
+     - * High accuracy/low bias of the simulation averaged estimator in most
+         cells
+       * Stability of the naive estimator in cells with fixed sources
+     - * Can lead to slightly negative fluxes in cells where the simulation
+         averaged estimator is used
+
+These estimators can be selected by setting the ``volume_estimator`` field in the
+:attr:`openmc.Settings.random_ray` dictionary. For example, to use the naive
+estimator, the following code would be used:
+
+::
+
+    settings.random_ray['volume_estimator'] = 'naive'
+
 ---------------------------------------
 Putting it All Together: Example Inputs
 ---------------------------------------

include/openmc/constants.h

@@ -342,6 +342,7 @@ enum class RunMode {
 
 enum class SolverType { MONTE_CARLO, RANDOM_RAY };
 
+enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };
 enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
 
 //==============================================================================

include/openmc/random_ray/flat_source_domain.h

@@ -1,6 +1,7 @@
 #ifndef OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H
 #define OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H
 
+#include "openmc/constants.h"
 #include "openmc/openmp_interface.h"
 #include "openmc/position.h"
 #include "openmc/source.h"
@@ -99,7 +100,8 @@ class FlatSourceDomain {
   double compute_k_eff(double k_eff_old) const;
   virtual void normalize_scalar_flux_and_volumes(
     double total_active_distance_per_iteration);
-  virtual int64_t add_source_to_scalar_flux();
+
+  int64_t add_source_to_scalar_flux();
   virtual void batch_reset();
   void convert_source_regions_to_tallies();
   void reset_tally_volumes();
@@ -117,6 +119,10 @@ class FlatSourceDomain {
   // Static Data members
   static bool volume_normalized_flux_tallies_;
 
+  //----------------------------------------------------------------------------
+  // Static data members
+  static RandomRayVolumeEstimator volume_estimator_;
+
   //----------------------------------------------------------------------------
   // Public Data members
 
@@ -132,18 +138,18 @@ class FlatSourceDomain {
 
   // 1D arrays representing values for all source regions
   vector<OpenMPMutex> lock_;
-  vector<int> was_hit_;
   vector<double> volume_;
   vector<double> volume_t_;
   vector<int> position_recorded_;
   vector<Position> position_;
 
   // 2D arrays stored in 1D representing values for all source regions x energy
   // groups
-  vector<float> scalar_flux_old_;
-  vector<float> scalar_flux_new_;
+  vector<double> scalar_flux_old_;
+  vector<double> scalar_flux_new_;
   vector<float> source_;
   vector<float> external_source_;
+  vector<bool> external_source_present_;
 
 protected:
   //----------------------------------------------------------------------------
@@ -155,6 +161,10 @@ class FlatSourceDomain {
     const vector<int32_t>& instances);
   void apply_external_source_to_cell_and_children(int32_t i_cell,
     Discrete* discrete, double strength_factor, int32_t target_material_id);
+  virtual void set_flux_to_flux_plus_source(
+    int64_t idx, double volume, int material, int g);
+  void set_flux_to_source(int64_t idx);
+  virtual void set_flux_to_old_flux(int64_t idx);
 
   //----------------------------------------------------------------------------
   // Private data members
@@ -178,6 +188,7 @@ class FlatSourceDomain {
 
   // 1D arrays representing values for all source regions
   vector<int> material_;
+  vector<double> volume_naive_;
 
   // 2D arrays stored in 1D representing values for all source regions x energy
   // groups

include/openmc/random_ray/linear_source_domain.h

@@ -28,7 +28,7 @@ class LinearSourceDomain : public FlatSourceDomain {
   double compute_k_eff(double k_eff_old) const;
   void normalize_scalar_flux_and_volumes(
     double total_active_distance_per_iteration) override;
-  int64_t add_source_to_scalar_flux() override;
+
   void batch_reset() override;
   void convert_source_regions_to_tallies();
   void reset_tally_volumes();
@@ -54,6 +54,13 @@ class LinearSourceDomain : public FlatSourceDomain {
   vector<MomentMatrix> mom_matrix_;
   vector<MomentMatrix> mom_matrix_t_;
 
+protected:
+  //----------------------------------------------------------------------------
+  // Methods
+  void set_flux_to_flux_plus_source(
+    int64_t idx, double volume, int material, int g) override;
+  void set_flux_to_old_flux(int64_t idx) override;
+
 }; // class LinearSourceDomain
 
 } // namespace openmc

include/openmc/random_ray/random_ray.h

@@ -43,6 +43,8 @@ class RandomRay : public Particle {
   // Public data members
   vector<float> angular_flux_;
 
+  bool ray_trace_only_ {false}; // If true, only perform geometry operations
+
 private:
   //----------------------------------------------------------------------------
   // Private data members

include/openmc/random_ray/random_ray_simulation.h

@@ -19,6 +19,7 @@ class RandomRaySimulation {
 
   //----------------------------------------------------------------------------
   // Methods
+  void compute_segment_correction_factors();
   void simulate();
   void reduce_simulation_statistics();
   void output_simulation_results() const;

openmc/settings.py

@@ -155,6 +155,10 @@ class Settings:
         :ray_source:
             Starting ray distribution (must be uniform in space and angle) as
             specified by a :class:`openmc.SourceBase` object.
+        :volume_estimator:
+            Choice of volume estimator for the random ray solver. Options are
+            'naive', 'simulation_averaged', or 'hybrid'.
+            The default is 'hybrid'.
         :source_shape:
             Assumed shape of the source distribution within each source
             region. Options are 'flat' (default), 'linear', or 'linear_xy'.
@@ -1091,6 +1095,10 @@ 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_estimator':
+                cv.check_value('volume estimator', random_ray[key],
+                               ('naive', 'simulation_averaged',
+                                'hybrid'))
             elif key == 'source_shape':
                 cv.check_value('source shape', random_ray[key],
                                ('flat', 'linear', 'linear_xy'))
@@ -1889,6 +1897,8 @@ 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_estimator':
+                    self.random_ray['volume_estimator'] = child.text
                 elif child.tag == 'source_shape':
                     self.random_ray['source_shape'] = child.text
                 elif child.tag == 'volume_normalized_flux_tallies':