Naive Flux Linear Source Consistency Fix

docs/source/methods/random_ray.rst

@@ -1034,6 +1034,39 @@ The contents of this section, alongside the equations for the flat source and
 scalar flux, Equations :eq:`source_update` and :eq:`phi_sim` respectively,
 completes the set of equations for LS.
 
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Consistency of the Scalar Flux Estimate
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+One subtlety of the linear source scheme deserves note. Intuitively, the
+issue is a mismatch of statistics: the naive volume treatment of Equation
+:eq:`phi_naive` updates the flux from a single batch's rays, while the
+linear source is anchored to the simulation-averaged centroid, so the two
+halves of the update describe different sets of tracks and the exactness
+the naive treatment promises is quietly broken. Concretely, the transport
+sweep evaluates each region's linear source of Equation :eq:`region_source`
+against the accumulated centroid :math:`\mathbf{r}_{\mathrm{c}}`, but a
+batch's tracks average that source at their own track-length-weighted
+centroid :math:`\mathbf{r}_{\mathrm{c},b}`, so the mean source the batch
+actually integrates is
+
+.. math::
+    :label: batch_sampled_source
+
+    Q_{i,g} + \boldsymbol{\vec{Q}}_{i,g} \cdot \left(\mathbf{r}_{\mathrm{c},b}
+    - \mathbf{r}_{\mathrm{c}}\right)\;.
+
+A flux update that adds back only :math:`Q_{i,g} / \Sigma_{t,i,g}` absorbs
+the difference as gradient-scale noise, which in optically thin scatter-fed
+regions can ignite self-sustaining negative fluxes. OpenMC therefore adds
+back the full batch-sampled source (divided by :math:`\Sigma_{t,i,g}`)
+whenever a region updates with its own batch volume, making that update
+exact for the batch's tracks. Regions updating with the simulation-averaged
+volume keep the original form: each batch's centroid scatters about the
+accumulated centroid it feeds, so the omitted term has no persistent sign
+and its contribution to the accumulated flux shrinks with the number of
+batches, while the original form carries less variance there.
+
 .. _methods-shannon-entropy-random-ray:
 
 -----------------------------

include/openmc/random_ray/flat_source_domain.h

@@ -239,10 +239,11 @@ class FlatSourceDomain {
     int target_material_id, const vector<int32_t>& instances);
   void apply_external_source_to_cell_and_children(
     int32_t i_cell, int src_idx, int32_t target_material_id);
-  virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g);
+  virtual void set_flux_to_flux_plus_source(
+    int64_t sr, double volume, bool batch_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;
+  virtual double flux_additive_term(int64_t sr, int g, bool batch_volume) const;
   double stabilized_flux(int64_t sr, int g, double phi_new) const;
 
   //! Adaptive-estimator "strong source" test. Returns true if, in any group,

include/openmc/random_ray/linear_source_domain.h

@@ -32,8 +32,11 @@ class LinearSourceDomain : public FlatSourceDomain {
 protected:
   //----------------------------------------------------------------------------
   // Methods
-  void set_flux_to_flux_plus_source(int64_t sr, double volume, int g) override;
+  void set_flux_to_flux_plus_source(
+    int64_t sr, double volume, bool batch_volume, int g) override;
   void set_flux_to_old_flux(int64_t sr, int g) override;
+  double flux_additive_term(
+    int64_t sr, int g, bool batch_volume) const override;
 
 }; // class LinearSourceDomain
 

include/openmc/random_ray/source_region.h

@@ -480,6 +480,12 @@ class SourceRegionContainer {
   Position& centroid(int64_t sr) { return centroid_[sr]; }
   const Position centroid(int64_t sr) const { return centroid_[sr]; }
 
+  Position& centroid_offset(int64_t sr) { return centroid_offset_[sr]; }
+  const Position centroid_offset(int64_t sr) const
+  {
+    return centroid_offset_[sr];
+  }
+
   Position& centroid_iteration(int64_t sr) { return centroid_iteration_[sr]; }
   const Position centroid_iteration(int64_t sr) const
   {
@@ -699,6 +705,11 @@ class SourceRegionContainer {
   vector<Position> centroid_;
   vector<Position> centroid_iteration_;
   vector<Position> centroid_t_;
+  // Offset of this batch's track-length-weighted centroid from the
+  // accumulated centroid the transport sweep evaluated the linear source
+  // against, used by the batch-consistent flux update (linear source solver
+  // only)
+  vector<Position> centroid_offset_;
   vector<MomentMatrix> mom_matrix_;
   vector<MomentMatrix> mom_matrix_t_;
   // A set of volume tally tasks. This more complicated data structure is

src/random_ray/flat_source_domain.cpp

@@ -342,13 +342,17 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes(
   }
 }
 
-// The additive term of the flux update for a source region and group. A
-// material region adds its reduced source q/Sigma_t. A void region has no
-// such term and instead adds a bounded contribution from its external
-// source, which is nonzero only in fixed source mode. The same term is used
-// by the strict estimator's rescue, which rescales only the transport part
-// of an update, so the two cannot drift apart.
-double FlatSourceDomain::flux_additive_term(int64_t sr, int g) const
+// The additive term of the flux update for a source region and group,
+// given whether the update divides the transport term by the region's own
+// batch volume (as opposed to its simulation-averaged volume). A material
+// region adds its reduced source q/Sigma_t. A void region has no such term
+// and instead adds a bounded contribution from its external source, which
+// is nonzero only in fixed source mode. The same term is used by the strict
+// estimator's rescue, which rescales only the transport part of an update,
+// so the two cannot drift apart. The linear source solver's term depends on
+// the volume choice (see its override); the flat source term does not.
+double FlatSourceDomain::flux_additive_term(
+  int64_t sr, int g, bool batch_volume) const
 {
   if (source_regions_.material(sr) == MATERIAL_VOID) {
     if (settings::run_mode == RunMode::FIXED_SOURCE) {
@@ -361,7 +365,7 @@ double FlatSourceDomain::flux_additive_term(int64_t sr, int g) const
 }
 
 void FlatSourceDomain::set_flux_to_flux_plus_source(
-  int64_t sr, double volume, int g)
+  int64_t sr, double volume, bool batch_volume, int g)
 {
   int material = source_regions_.material(sr);
   int temp = source_regions_.temperature_idx(sr);
@@ -373,7 +377,8 @@ void FlatSourceDomain::set_flux_to_flux_plus_source(
       source_regions_.density_mult(sr);
     source_regions_.scalar_flux_new(sr, g) /= (sigma_t * volume);
   }
-  source_regions_.scalar_flux_new(sr, g) += flux_additive_term(sr, g);
+  source_regions_.scalar_flux_new(sr, g) +=
+    flux_additive_term(sr, g, batch_volume);
 }
 
 // Applies the "diagonal stabilization" technique developed by Gunow et al.
@@ -645,7 +650,7 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
         // Hit this iteration: the flat source from the previous iteration plus
         // this iteration's transport contribution, normalized by the chosen
         // volume, then stabilized.
-        set_flux_to_flux_plus_source(sr, volume, g);
+        set_flux_to_flux_plus_source(sr, volume, use_naive_volume, g);
         double raw = source_regions_.scalar_flux_new(sr, g);
         double phi = stabilized_flux(sr, g, raw);
         // The strict adaptive estimator applies a per-batch fixup to
@@ -665,9 +670,9 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
         // flat shapes.
         if (is_strict && phi < 0.0) {
           if (volume != volume_iteration) {
-            double additive = flux_additive_term(sr, g);
-            double rescued =
-              (raw - additive) * (volume / volume_iteration) + additive;
+            double rescued = (raw - flux_additive_term(sr, g, false)) *
+                               (volume / volume_iteration) +
+                             flux_additive_term(sr, g, true);
             phi = stabilized_flux(sr, g, rescued);
             region_rescued = true;
           }

src/random_ray/linear_source_domain.cpp

@@ -153,6 +153,22 @@ void LinearSourceDomain::normalize_scalar_flux_and_volumes(
 // update the simulation-averaged cell-wise volume estimates
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions(); sr++) {
+    // Offset of this batch's track centroid from the accumulated centroid
+    // the transport sweep evaluated the linear source against (the value
+    // held here before this batch is folded in). The flux update uses it to
+    // add back the source as actually sampled by this batch's tracks (see
+    // flux_additive_term). Regions with no prior accumulated volume used
+    // per-segment midpoints in place of a centroid, for which the batch
+    // offset is identically zero.
+    if (source_regions_.volume_t(sr) > 0.0 &&
+        source_regions_.volume(sr) > 0.0) {
+      source_regions_.centroid_offset(sr) =
+        source_regions_.centroid_iteration(sr) *
+          (1.0 / source_regions_.volume(sr)) -
+        source_regions_.centroid(sr);
+    } else {
+      source_regions_.centroid_offset(sr) = {0.0, 0.0, 0.0};
+    }
     source_regions_.centroid_t(sr) += source_regions_.centroid_iteration(sr);
     source_regions_.mom_matrix_t(sr) += source_regions_.mom_matrix(sr);
     source_regions_.volume_t(sr) += source_regions_.volume(sr);
@@ -173,15 +189,50 @@ void LinearSourceDomain::normalize_scalar_flux_and_volumes(
   }
 }
 
+// A material region updated with its own batch (naive) volume adds back,
+// on top of its flat source, the gradient part of the source that the
+// batch's rays actually integrated. The sweep evaluates the linear source
+// against the accumulated centroid, but the batch's tracks average it at
+// their own centroid, so the mean emission the tracks sampled is the flat
+// source plus the gradient dotted with the batch centroid offset. Including
+// that term makes the naive-volume update an exact per-batch track
+// identity, in which the flux estimate is the track-length average of the
+// angular flux and so inherits its sign. Omitting it leaves gradient-scale
+// noise in the flux with no flat-source counterpart, which in
+// near-cancellation regions (scattering ratio near one) can exceed the flux
+// itself and ignite self-sustaining negativity.
+//
+// A region updated with the simulation-averaged volume keeps the original
+// form with no added term. Each batch's centroid scatters about the
+// accumulated centroid it feeds, so the omitted term has no persistent
+// sign and its contribution to the accumulated flux shrinks with the
+// number of batches, while the original form samples the residual between
+// the angular flux and the linear source model, whose noise is smaller
+// than the direct form's wherever the model tracks the field. Each volume
+// treatment is thus paired with the update form that is exact or better
+// for it. Void regions have no gradient term and take the flat source
+// treatment.
+double LinearSourceDomain::flux_additive_term(
+  int64_t sr, int g, bool batch_volume) const
+{
+  double term = FlatSourceDomain::flux_additive_term(sr, g, batch_volume);
+  if (batch_volume && source_regions_.material(sr) != MATERIAL_VOID) {
+    term += source_regions_.source_gradients(sr, g).dot(
+      source_regions_.centroid_offset(sr));
+  }
+  return term;
+}
+
 void LinearSourceDomain::set_flux_to_flux_plus_source(
-  int64_t sr, double volume, int g)
+  int64_t sr, double volume, bool batch_volume, int g)
 {
   int material = source_regions_.material(sr);
   if (material == MATERIAL_VOID) {
-    FlatSourceDomain::set_flux_to_flux_plus_source(sr, volume, g);
+    FlatSourceDomain::set_flux_to_flux_plus_source(sr, volume, batch_volume, g);
   } else {
     source_regions_.scalar_flux_new(sr, g) /= volume;
-    source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);
+    source_regions_.scalar_flux_new(sr, g) +=
+      flux_additive_term(sr, g, batch_volume);
   }
   // If a source region is small, then the moments are likely noisy, so we zero
   // them. This is reasonable, given that small regions can get by with a flat

src/random_ray/source_region.cpp

@@ -100,6 +100,7 @@ void SourceRegionContainer::push_back(const SourceRegion& sr)
     centroid_.push_back(sr.centroid_);
     centroid_iteration_.push_back(sr.centroid_iteration_);
     centroid_t_.push_back(sr.centroid_t_);
+    centroid_offset_.push_back({0.0, 0.0, 0.0});
     mom_matrix_.push_back(sr.mom_matrix_);
     mom_matrix_t_.push_back(sr.mom_matrix_t_);
   }
@@ -159,6 +160,7 @@ void SourceRegionContainer::assign(
     centroid_.clear();
     centroid_iteration_.clear();
     centroid_t_.clear();
+    centroid_offset_.clear();
     mom_matrix_.clear();
     mom_matrix_t_.clear();
   }
@@ -262,6 +264,8 @@ void SourceRegionContainer::adjoint_reset()
   std::fill(centroid_iteration_.begin(), centroid_iteration_.end(),
     Position {0.0, 0.0, 0.0});
   std::fill(centroid_t_.begin(), centroid_t_.end(), Position {0.0, 0.0, 0.0});
+  std::fill(
+    centroid_offset_.begin(), centroid_offset_.end(), Position {0.0, 0.0, 0.0});
   std::fill(mom_matrix_.begin(), mom_matrix_.end(),
     MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
   std::fill(mom_matrix_t_.begin(), mom_matrix_t_.end(),

tests/regression_tests/random_ray_fixed_source_linear/linear/results_true.dat

@@ -1,9 +1,9 @@
 tally 1:
-2.339086E+00
-2.747305E-01
+2.347067E+00
+2.759398E-01
 tally 2:
 1.089827E-01
-6.069324E-04
+6.069322E-04
 tally 3:
-7.300831E-03
-2.715940E-06
+7.300833E-03
+2.715941E-06

tests/regression_tests/random_ray_fixed_source_linear/linear_xy/results_true.dat

@@ -1,9 +1,9 @@
 tally 1:
-2.335703E+00
-2.742866E-01
+2.345238E+00
+2.755122E-01
 tally 2:
-1.081884E-01
-5.983316E-04
+1.081885E-01
+5.983320E-04
 tally 3:
-7.295389E-03
-2.711859E-06
+7.295394E-03
+2.711862E-06

tests/regression_tests/random_ray_fixed_source_mesh/linear/results_true.dat

@@ -1,9 +1,9 @@
 tally 1:
-1.780361E+00
-2.137912E-01
+1.768881E+00
+2.095454E-01
 tally 2:
-8.230400E-02
-4.596391E-04
+8.228219E-02
+4.596629E-04
 tally 3:
-5.207797E-03
-1.834531E-06
+5.207806E-03
+1.834565E-06