Adaptive Volume Estimators

docs/source/io_formats/settings.rst

@@ -626,7 +626,11 @@ found in the :ref:`random ray user guide <random_ray>`.
 
   :volume_estimator:
     Specifies choice of volume estimator for the random ray solver. Options
-    are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'.
+    are 'naive', 'simulation_averaged', 'hybrid', 'adaptive',
+    'strict_adaptive', or 'auto'. The default is 'auto', which selects
+    'adaptive' for standard solves and 'strict_adaptive' for solves whose
+    results feed variance reduction (weight window generation and adjoint
+    workflows).
 
     *Default*: None
 
@@ -1648,6 +1652,7 @@ and 10. The verbosity levels are defined as follows:
   :5: all of the above + file I/O
   :6: all of the above + timing statistics and initialization messages
   :7: all of the above + :math:`k` by generation
+  :8: all of the above + random ray volume-estimator diagnostics
   :9: all of the above + indicate when each particle starts
   :10: all of the above + event information
 

docs/source/methods/random_ray.rst

@@ -511,18 +511,53 @@ 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.
+Thus, the simulation averaged estimator is generally preferred over the naive
+estimator 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
+a best of both worlds result, with the benefits of the unbiased 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.
+with external source terms via use of the naive estimator. If instability is
+encountered despite high ray densities, then the naive estimator may be
+preferable.
+
+OpenMC also features an "adaptive" volume estimator that generalizes the
+hybrid estimator. It uses the simulation averaged estimator by default and
+automatically demotes individual cells to the naive treatment (the naive
+volume and the previous-flux miss treatment) when they show signs of
+instability. The most important case this adds over the hybrid estimator is
+a cell fed almost entirely by in-scatter from other groups, as in the
+optically thin air regions common to shielding problems, where the reduced
+source dwarfs the flux even though no external source is present.
+
+A cell is demoted when it is hit-starved or when its reduced source is
+negative. During the inactive batches, a cell whose reduced source is much
+larger than its scalar flux is also demoted. Finally, beginning at the end
+of the inactive batches and re-evaluated throughout the active phase, a
+cell is demoted permanently if its flux accumulated over the simulation is
+negative or is dominated by sources that do not derive from its own flux.
+Because these last decisions are made from accumulated statistics and are
+never reversed, the estimator choice does not churn with iteration noise
+during the tallied batches. When a linear source shape is in use, demoted
+cells also revert to a flat source representation.
+
+A "strict adaptive" variant is provided for solves whose results feed
+variance reduction, where even a small number of slightly negative flux
+estimates can degrade the adjoint solve and the quality of generated weight
+windows. It runs the same machinery and additionally repairs any negative
+flux estimate each batch, first by recomputing it with the batch's own
+volume and then, if it is still negative, by falling back on the previous
+iterate. A cell that needs the repair repeatedly is demoted outright, which
+keeps the one-sided repair from biasing its flux upward. The repair
+introduces a small conservative bias overall (several hundred pcm on
+typical eigenvalue problems), so the strict variant is not used for
+standard solves.
+
+By default, OpenMC selects the volume estimator automatically ("auto").
+Weight window generation and adjoint solves receive the strict adaptive
+estimator, and all other solves receive the adaptive estimator.
 
 A table that summarizes the pros and cons, as well as recommendations for
 different use cases, is given in the :ref:`volume

docs/source/usersguide/random_ray.rst

@@ -1080,6 +1080,12 @@ following methods are currently available in OpenMC:
      - Description
      - Pros
      - Cons
+   * - ``auto`` (default)
+     - Automatically selects an appropriate estimator for the type of
+       simulation being performed. Most users do not need to consider this
+       setting further.
+     - * No user input needed
+     - * N/A
    * - ``simulation_averaged``
      - Accumulates total active ray lengths in each FSR over all iterations,
        improving the estimate of the volume in each cell each iteration.
@@ -1099,15 +1105,44 @@ following methods are currently available in OpenMC:
          unstable
      - * Biased estimator
        * Requires more rays or longer active ray length to mitigate bias
-   * - ``hybrid`` (default)
+   * - ``hybrid``
      - 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
+     - * Accuracy of the unbiased 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
+   * - ``adaptive``
+     - Generalizes the hybrid estimator. Uses the simulation averaged
+       estimator by default, but automatically (and permanently) demotes
+       individual cells to the naive treatment when their accumulated
+       statistics indicate the simulation averaged estimator is unstable
+       there (e.g., cells dominated by external or in-scatter sources, and
+       hit-starved cells).
+     - * Accuracy of the simulation averaged estimator in most cells
+       * Stable in cases where the simulation averaged and hybrid estimators
+         are not
+       * No parameters to tune
+     - * Benefits from a longer inactive phase to inform the demotion
+         decisions
+   * - ``strict_adaptive``
+     - As ``adaptive``, but additionally applies a per-batch fixup to any
+       negative flux estimate (recomputing it with the batch's own volume,
+       then falling back on the previous iterate) and demotes chronically
+       affected cells to the naive treatment.
+     - * Suppresses the negative flux estimates other estimators can produce
+         in pathological cases
+       * Improves the quality of generated weight windows
+     - * The one-sided fixup introduces a small conservative bias, so it is
+         not recommended where unbiased results are the priority
+
+By default, the ``volume_estimator`` field is set to ``auto``, which selects
+``strict_adaptive`` for solves whose results feed variance reduction (weight
+window generation and adjoint workflows) and ``adaptive`` for all other
+solves. The end-of-simulation output reports which estimator was selected,
+and explicitly setting any other value overrides the automatic selection.
 
 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
@@ -1117,6 +1152,22 @@ estimator, the following code would be used:
 
     settings.random_ray['volume_estimator'] = 'naive'
 
+The ``auto`` setting is the default, as it gives reliable behavior out of
+the box across problem types. The adaptive estimator is especially valuable
+for fixed source and shielding problems, where optically thin, scattering-
+or streaming-dominated regions (for example, the air- or void-filled regions
+of a shielding model) can destabilize the ``hybrid`` and
+``simulation_averaged`` estimators. It detects and stabilizes the affected
+cells automatically while leaving the rest of the problem on the unbiased
+simulation averaged estimator. Because demotions are decided from each
+cell's accumulated statistics rather than from single-iteration values, the
+estimator choice does not churn with iteration noise and avoids the bias
+that per-iteration selection can introduce. Solves that feed variance
+reduction are routed to ``strict_adaptive`` instead, as even a small number
+of slightly negative flux estimates in the near-void regions of pathological
+problems can otherwise degrade the adjoint solve and the quality of
+generated weight windows.
+
 -----------------
 Adjoint Flux Mode
 -----------------

include/openmc/constants.h

@@ -61,20 +61,6 @@ constexpr double RADIAL_MESH_TOL {1e-10};
 // Maximum number of random samples per history
 constexpr int MAX_SAMPLE {100000};
 
-// Avg. number of hits per batch to be defined as a "small"
-// 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};
-
-// The minimum macroscopic cross section value considered non-void for the
-// random ray solver. Materials with any group with a cross section below this
-// value will be converted to pure void.
-constexpr double MINIMUM_MACRO_XS {1e-6};
-
 // Relative dead band applied to weight window comparisons: particles split
 // only above upper * (1 + tol) and roulette only below lower * (1 - tol).
 // Weight window arithmetic can land a particle's weight exactly back on a
@@ -378,7 +364,14 @@ enum class RunMode {
 
 enum class SolverType { MONTE_CARLO, RANDOM_RAY };
 
-enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };
+enum class RandomRayVolumeEstimator {
+  NAIVE,
+  SIMULATION_AVERAGED,
+  HYBRID,
+  ADAPTIVE,
+  STRICT_ADAPTIVE,
+  AUTO
+};
 enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
 enum class RandomRaySampleMethod { PRNG, HALTON, S2 };
 enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT };

include/openmc/random_ray/flat_source_domain.h

@@ -12,6 +12,50 @@
 
 namespace openmc {
 
+// Avg. number of hits per batch to be defined as a "small" source region.
+constexpr double MIN_HITS_PER_BATCH {1.5};
+
+// Strong-source ratio threshold for the adaptive volume estimator. A source
+// region is treated as having a "strong" inhomogeneous source in any group
+// where the reduced source q/Sigma_t exceeds this multiple of the region's
+// scalar flux, indicating a source sustained by an external or in-scatter
+// contribution rather than by the local flux. Such regions are given the
+// naive volume and previous-flux miss treatment. The value sits well inside
+// the range over which benign problems remain untriggered while pathological
+// cells are still caught.
+constexpr double ADAPTIVE_VOLUME_KAPPA {4.0};
+
+// Chronic-negativity demotion thresholds for the strict adaptive volume
+// estimator. A region whose flux has gone negative (before the fixup) in at
+// least max(MIN_COUNT, RATE * current_batch) batches is demoted to the naive
+// volume and previous-flux miss treatment. Without this channel the
+// non-negativity floor would mask the accumulated-flux sign signal that the
+// adaptive demotion relies on, leaving noisy regions to be clipped every
+// batch and biasing their fluxes upward. Demotion instead moves such regions
+// onto an estimator that does not need clipping.
+constexpr int NEGATIVE_FLUX_DEMOTION_MIN_COUNT {3};
+constexpr double NEGATIVE_FLUX_DEMOTION_RATE {0.005};
+
+// 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};
+
+// The minimum macroscopic cross section value considered non-void for the
+// random ray solver. Materials with any group with a cross section below this
+// value will be converted to pure void.
+constexpr double MINIMUM_MACRO_XS {1e-6};
+
+// True for the members of the adaptive volume estimator family: the
+// adaptive estimator, and the strict adaptive estimator, which runs the
+// same machinery plus a per-batch non-negativity enforcement on the flux
+// iterates.
+inline bool is_adaptive_family(RandomRayVolumeEstimator e)
+{
+  return e == RandomRayVolumeEstimator::ADAPTIVE ||
+         e == RandomRayVolumeEstimator::STRICT_ADAPTIVE;
+}
+
 /*
  * The FlatSourceDomain class encompasses data and methods for storing
  * scalar flux and source region for all flat source regions in a
@@ -39,6 +83,7 @@ class FlatSourceDomain {
   void reset_tally_volumes();
   void random_ray_tally();
   virtual void accumulate_iteration_flux();
+  void demotion_step();
   void output_to_vtk() const;
   void convert_external_sources(bool use_adjoint_sources);
   void count_external_source_regions();
@@ -59,7 +104,6 @@ class FlatSourceDomain {
   SourceRegionHandle get_subdivided_source_region_handle(
     SourceRegionKey sr_key, Position r, Direction u);
   void finalize_discovered_source_regions();
-  void apply_transport_stabilization();
   int64_t n_source_regions() const
   {
     return source_regions_.n_source_regions();
@@ -93,7 +137,14 @@ class FlatSourceDomain {
 
   //----------------------------------------------------------------------------
   // Static data members
+  // The volume estimator as configured ("auto" by default). This is set when
+  // the settings are read and is never modified by the solver.
   static RandomRayVolumeEstimator volume_estimator_;
+  // The concrete estimator the solver runs with, assigned at the start of
+  // every random ray solve. It holds the configured value, or for "auto" the
+  // estimator selected for the type of simulation being performed. All
+  // solver code reads this member rather than volume_estimator_.
+  static RandomRayVolumeEstimator resolved_volume_estimator_;
 
   //----------------------------------------------------------------------------
   // Public Data members
@@ -103,6 +154,27 @@ class FlatSourceDomain {
   int64_t n_external_source_regions_ {0}; // Total number of source regions with
                                           // non-zero external source terms
 
+  // Final-batch snapshot of the naive volume treatment, partitioned by
+  // mutually exclusive cause (the cause counts sum to n_final_naive_), for
+  // end-of-simulation reporting. The two demote-only decisions made from the
+  // running accumulated flux (a strong accumulated feed and a negative
+  // accumulated flux) are counted with first priority, so their counts equal
+  // the decisions settled by the final batch. The per-batch strong-source
+  // test and hit-starved causes count the remainder.
+  int64_t n_final_naive_ {0};
+  int64_t n_final_latch_ {0};   // strong source, from the accumulated feed
+  int64_t n_final_strong_ {0};  // strong source, from the per-batch test
+  int64_t n_final_sign_ {0};    // negative accumulated flux
+  int64_t n_final_small_ {0};   // hit-starved
+  int64_t n_final_chronic_ {0}; // chronic negativity (strict adaptive)
+  // Final-batch counts of the strict adaptive estimator's non-negativity
+  // enforcement: regions whose negative batch flux was recomputed with the
+  // batch volume (rescued), and regions floored at the previous iterate
+  // after the rescue was insufficient or unavailable.
+  int64_t n_final_rescued_ {0};
+  int64_t n_final_floored_ {0};
+  bool final_stats_valid_ {false};
+
   // 1D array representing source region starting offset for each OpenMC Cell
   // in model::cells
   vector<int64_t> source_region_offsets_;
@@ -170,6 +242,20 @@ class FlatSourceDomain {
   virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g);
   void set_flux_to_source(int64_t sr, int g);
   virtual void set_flux_to_old_flux(int64_t sr, int g);
+  double flux_additive_term(int64_t sr, int g) const;
+  double stabilized_flux(int64_t sr, int g, double phi_new) const;
+
+  //! Adaptive-estimator "strong source" test. Returns true if, in any group,
+  //! the region's reduced source q/Sigma_t is negative (with a non-negative
+  //! previous-iteration flux) or exceeds ADAPTIVE_VOLUME_KAPPA times the
+  //! (non-negative) previous-iteration scalar flux. The ratio condition is
+  //! only checked when include_ratio is set, which callers do only during
+  //! the inactive batches. The test is shared by the flat volume switch
+  //! (add_source_to_scalar_flux) and the linear gradient fallback
+  //! (update_single_neutron_source), with the region's per-group
+  //! reduced-source and previous-flux arrays passed directly.
+  bool region_has_strong_source(const float* reduced_source,
+    const double* flux_old, bool include_ratio) const;
 
   //----------------------------------------------------------------------------
   // Private data members