Random Ray Linear Source Gradient Limiter

docs/source/io_formats/settings.rst

@@ -634,6 +634,20 @@ found in the :ref:`random ray user guide <random_ray>`.
 
     *Default*: None
 
+  :source_shape:
+    Specifies the assumed shape of the source distribution within each
+    source region. Options are "flat", "linear", or "linear_xy".
+
+    *Default*: flat
+
+  :source_gradient_limiter:
+    Specifies whether to rescale linear source gradients as needed so that
+    the source shape modeled within each source region remains non-negative
+    over the region's bounding box, as sampled by the rays that have crossed
+    it (bool). Only used when the source shape is "linear" or "linear_xy".
+
+    *Default*: false
+
   :volume_normalized_flux_tallies:
     Specifies whether to normalize flux tallies by volume (bool). The
     default is 'False'. When enabled, flux tallies will be reported in units

docs/source/methods/random_ray.rst

@@ -1034,6 +1034,64 @@ 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.
 
+.. _methods_random_ray_gradient_limiter:
+
+~~~~~~~~~~~~~~~~~~~~~~~~
+Source Gradient Limiting
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+The fitted source gradient :math:`\boldsymbol{\vec{Q}}_{i,g} =
+\mathbf{M}_i^{-1} \boldsymbol{\vec{q}}_{i,g}` amplifies noise in the fitted
+moments along any thin extent of a region, so a poorly sampled region can
+carry a spuriously steep gradient and emit a negative source over part of
+its extent. Rays crossing that part can carry negative angular flux
+downstream, which optically thin media with scattering ratios near one can
+amplify.
+
+When the source gradient limiter is enabled, each group's gradient is
+rescaled so that the modeled source stays non-negative over the region's
+axis-aligned bounding box. The box is accumulated from the endpoints of
+every ray segment that has crossed the region past the ray's inactive
+length. These lie on the region's boundary except where a ray starts or
+ends inside it. The linear term is lowest at a corner of the box, where it
+reaches
+
+.. math::
+    :label: gradient-limiter-bound
+
+    \sum_{d \in \{x, y, z\}} \; \min_{x_d \in \{x^{\min}_{i,d},\,
+    x^{\max}_{i,d}\}} \left(\boldsymbol{\vec{Q}}_{i,g}\right)_d \left(x_d -
+    r_{\mathrm{c},i,d}\right),
+
+where :math:`x^{\min}_{i}` and :math:`x^{\max}_{i}` are the box bounds,
+:math:`\mathbf{r}_{\mathrm{c},i}` is the centroid, and :math:`d` indexes
+their components. Whenever the flat source :math:`Q_{i,g}` plus this
+minimum is negative, the gradient is scaled by the ratio of the flat source
+to the magnitude of the minimum, so that the modeled source reaches zero at
+that corner. Because the linear term integrates to zero over the region,
+the rescaling preserves the region's mean emission, and gradients that pass
+the test are left untouched. A group whose flat source is not positive has
+its gradient zeroed. Once the region's extreme points along each axis have
+been sampled, the box contains the region and the modeled source is
+non-negative throughout it. The bound is exact for axis-aligned box regions
+and conservative for others: a sphere is limited by up to a factor of
+:math:`\sqrt{3}` more than necessary, and a thin region lying diagonally to
+the axes by much more, as its bounding box is far larger than the region.
+
+This is the treatment `MPACT <Choi-2024_>`_ applies in its limited linear
+source approximation, with the same mean-preserving factor. MPACT finds
+the minimum source exactly, over the entrance and exit points of every
+segment crossing the region, which requires the fixed set of tracks that
+deterministic MOC lays down once. Random ray samples new rays every batch,
+so no such segment set exists when the source is built, and the sampled
+bounding box takes its place.
+
+The limiter is off by default because a steep fit can also be physical, as
+in the optically thick regions of deep-penetration problems, where
+limiting discards real shape information and alters the solution at
+depth. It is best reserved for simulations that negative sources
+destabilize.
+
 .. _methods-shannon-entropy-random-ray:
 
 -----------------------------
@@ -1196,6 +1254,7 @@ in random ray particle transport are:
 .. _Tramm-2020: https://doi.org/10.1051/EPJCONF/202124703021
 .. _Cosgrove-2023: https://doi.org/10.1080/00295639.2023.2270618
 .. _Ferrer-2016: https://doi.org/10.13182/NSE15-6
+.. _Choi-2024: https://doi.org/10.1080/00295639.2023.2224234
 .. _Gunow-2018: https://dspace.mit.edu/handle/1721.1/119030
 
 .. only:: html

docs/source/usersguide/random_ray.rst

@@ -979,6 +979,21 @@ in the :attr:`openmc.Settings.random_ray` dictionary to ``'linear'`` as::
 LS enables the use of coarser mesh discretizations and lower ray populations,
 offsetting the increased computation per ray.
 
+In poorly sampled source regions, fitted gradients can become spuriously
+steep, producing negative sources that may destabilize optically thin,
+scattering-dominated problems. If this occurs, a gradient limiter can be
+enabled as::
+
+    settings.random_ray['source_gradient_limiter'] = True
+
+The limiter rescales a region's gradient as needed so that the modeled
+source stays non-negative over the region's bounding box, as sampled by
+the rays that have crossed it, preserving the region's mean emission. The
+limiter is off by default, as limiting also clips physically steep source
+shapes such as those found in optically thick regions of deep-penetration
+problems; see the :ref:`methods documentation
+<methods_random_ray_gradient_limiter>` for details.
+
 While OpenMC has no specific mode for 2D simulations, such simulations can be
 performed implicitly by leaving one of the dimensions of the geometry unbounded
 or by imposing reflective boundary conditions with no variation in between them

include/openmc/random_ray/flat_source_domain.h

@@ -122,6 +122,9 @@ class FlatSourceDomain {
   static bool volume_normalized_flux_tallies_;
   // If the user wants outputs based on the adjoint flux
   static bool adjoint_requested_;
+  // If the user wants linear source gradients rescaled so the modeled source
+  // stays non-negative over each source region
+  static bool source_gradient_limiter_;
   // The solve currently being executed
   static RandomRaySolve solve_;
   static bool fw_cadis_local_;

include/openmc/random_ray/source_region.h

@@ -167,6 +167,10 @@ class SourceRegionHandle {
   Position* centroid_t_;
   MomentMatrix* mom_matrix_;
   MomentMatrix* mom_matrix_t_;
+  // Bounding box of the ray segment endpoints sampled in this region, kept
+  // only when the source gradient limiter is enabled (see SourceRegion).
+  Position* extent_min_;
+  Position* extent_max_;
   // A set of volume tally tasks. This more complicated data structure is
   // convenient for ensuring that volumes are only tallied once per source
   // region, regardless of how many energy groups are used for tallying.
@@ -259,6 +263,27 @@ class SourceRegionHandle {
   MomentMatrix& mom_matrix_t() { return *mom_matrix_t_; }
   const MomentMatrix mom_matrix_t() const { return *mom_matrix_t_; }
 
+  const Position extent_min() const { return *extent_min_; }
+  const Position extent_max() const { return *extent_max_; }
+
+  // Grows the sampled bounding box to include a point
+  void expand_extent(const Position& p)
+  {
+    // Conditional stores: once the box has converged, no write is made
+    if (p.x < extent_min_->x)
+      extent_min_->x = p.x;
+    if (p.y < extent_min_->y)
+      extent_min_->y = p.y;
+    if (p.z < extent_min_->z)
+      extent_min_->z = p.z;
+    if (p.x > extent_max_->x)
+      extent_max_->x = p.x;
+    if (p.y > extent_max_->y)
+      extent_max_->y = p.y;
+    if (p.z > extent_max_->z)
+      extent_max_->z = p.z;
+  }
+
   std::unordered_set<TallyTask, TallyTask::HashFunctor>& volume_task()
   {
     return *volume_task_;
@@ -372,6 +397,14 @@ class SourceRegion {
   MomentMatrix mom_matrix_t_ {0.0, 0.0, 0.0, 0.0, 0.0,
     0.0}; //!< The spatial moment matrix accumulated over all iterations
 
+  // Bounding box of the ray segment endpoints sampled in this region. Segment
+  // endpoints lie on the region boundary, so the box converges to the
+  // region's true extent. It is accumulated only when the source gradient
+  // limiter is enabled, which bounds the linear source over it. The empty
+  // box has its minimum above its maximum.
+  Position extent_min_ {INFTY, INFTY, INFTY};
+  Position extent_max_ {-INFTY, -INFTY, -INFTY};
+
   // A set of volume tally tasks. This more complicated data structure is
   // convenient for ensuring that volumes are only tallied once per source
   // region, regardless of how many energy groups are used for tallying.
@@ -411,10 +444,10 @@ class SourceRegionContainer {
 public:
   //----------------------------------------------------------------------------
   // Constructors
-  SourceRegionContainer(
-    int negroups, bool is_linear, bool is_adaptive, bool is_strict_adaptive)
+  SourceRegionContainer(int negroups, bool is_linear, bool is_adaptive,
+    bool is_strict_adaptive, bool track_extents)
     : negroups_(negroups), is_linear_(is_linear), is_adaptive_(is_adaptive),
-      is_strict_adaptive_(is_strict_adaptive)
+      is_strict_adaptive_(is_strict_adaptive), track_extents_(track_extents)
   {}
   SourceRegionContainer() = default;
 
@@ -498,6 +531,12 @@ class SourceRegionContainer {
     return mom_matrix_t_[sr];
   }
 
+  Position& extent_min(int64_t sr) { return extent_min_[sr]; }
+  const Position extent_min(int64_t sr) const { return extent_min_[sr]; }
+
+  Position& extent_max(int64_t sr) { return extent_max_[sr]; }
+  const Position extent_max(int64_t sr) const { return extent_max_[sr]; }
+
   MomentArray& source_gradients(int64_t sr, int g)
   {
     return source_gradients_[index(sr, g)];
@@ -676,6 +715,9 @@ class SourceRegionContainer {
   bool is_linear_ {false};
   bool is_adaptive_ {false};
   bool is_strict_adaptive_ {false};
+  // Whether the sampled bounding boxes are stored (linear source with the
+  // source gradient limiter enabled)
+  bool track_extents_ {false};
 
   // SoA storage for scalar fields (one item per source region)
   vector<int> material_;
@@ -701,6 +743,8 @@ class SourceRegionContainer {
   vector<Position> centroid_t_;
   vector<MomentMatrix> mom_matrix_;
   vector<MomentMatrix> mom_matrix_t_;
+  vector<Position> extent_min_;
+  vector<Position> extent_max_;
   // A set of volume tally tasks. This more complicated data structure is
   // convenient for ensuring that volumes are only tallied once per source
   // region, regardless of how many energy groups are used for tallying.

openmc/settings.py

@@ -210,6 +210,12 @@ class Settings:
         :source_shape:
             Assumed shape of the source distribution within each source region.
             Options are 'flat' (default), 'linear', or 'linear_xy'.
+        :source_gradient_limiter:
+            Whether to rescale linear source gradients as needed so that the
+            source shape modeled within each source region remains
+            non-negative over the region's bounding box, as sampled by the
+            rays that have crossed it (bool). The default is 'False'. Only
+            used when the source shape is 'linear' or 'linear_xy'.
         :volume_normalized_flux_tallies:
             Whether to normalize flux tallies by volume (bool). The default is
             'False'. When enabled, flux tallies will be reported in units of
@@ -1426,6 +1432,8 @@ def random_ray(self, random_ray: dict):
                                ('flat', 'linear', 'linear_xy'))
             elif key == 'volume_normalized_flux_tallies':
                 cv.check_type('volume normalized flux tallies', value, bool)
+            elif key == 'source_gradient_limiter':
+                cv.check_type('source gradient limiter', value, bool)
             elif key == 'adjoint':
                 cv.check_type('adjoint', value, bool)
             elif key == 'source_region_meshes':
@@ -2519,6 +2527,10 @@ def _random_ray_from_xml_element(self, root, meshes=None):
                     self.random_ray['adjoint'] = (
                         child.text in ('true', '1')
                     )
+                elif child.tag == 'source_gradient_limiter':
+                    self.random_ray['source_gradient_limiter'] = (
+                        child.text in ('true', '1')
+                    )
                 elif child.tag == 'adjoint_source':
                     self.random_ray['adjoint_source'] = []
                     for subelem in child.findall('source'):

src/random_ray/flat_source_domain.cpp

@@ -33,6 +33,7 @@ RandomRayVolumeEstimator FlatSourceDomain::resolved_volume_estimator_ {
   RandomRayVolumeEstimator::AUTO};
 bool FlatSourceDomain::volume_normalized_flux_tallies_ {false};
 bool FlatSourceDomain::adjoint_requested_ {false};
+bool FlatSourceDomain::source_gradient_limiter_ {false};
 RandomRaySolve FlatSourceDomain::solve_ {RandomRaySolve::FORWARD};
 bool FlatSourceDomain::fw_cadis_local_ {false};
 double FlatSourceDomain::diagonal_stabilization_rho_ {1.0};
@@ -61,8 +62,9 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
   bool is_adaptive = is_adaptive_family(resolved_volume_estimator_);
   bool is_strict_adaptive =
     resolved_volume_estimator_ == RandomRayVolumeEstimator::STRICT_ADAPTIVE;
-  source_regions_ = SourceRegionContainer(
-    negroups_, is_linear, is_adaptive, is_strict_adaptive);
+  // The sampled bounding boxes exist only for the source gradient limiter
+  source_regions_ = SourceRegionContainer(negroups_, is_linear, is_adaptive,
+    is_strict_adaptive, is_linear && source_gradient_limiter_);
 
   // Initialize tally volumes
   if (volume_normalized_flux_tallies_) {

src/random_ray/linear_source_domain.cpp

@@ -1,5 +1,7 @@
 #include "openmc/random_ray/linear_source_domain.h"
 
+#include <algorithm>
+
 #include "openmc/cell.h"
 #include "openmc/geometry.h"
 #include "openmc/material.h"
@@ -133,6 +135,39 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
       srh.source_gradients(g) = {0.0, 0.0, 0.0};
     }
   }
+
+  // If enabled by the user, limit the source gradients so the modeled local
+  // source q(r) = q_flat + (r - centroid) . q_gradient stays non-negative
+  // over the region's bounding box as sampled by the ray segment endpoints.
+  // The linear term is lowest at the box corner each gradient component
+  // points away from, so its minimum is the sum, over the three axes, of the
+  // gradient component times the offset from the centroid to that face.
+  // Once the region's extreme points along each axis have been sampled the
+  // box contains the region, and the modeled source is non-negative
+  // throughout it whenever the flat source covers the dip. When it does not,
+  // the gradient is scaled by their ratio, which preserves the region's mean
+  // emission, since the linear term integrates to zero over the region;
+  // gradients that pass are left untouched. A non-positive flat source
+  // leaves no shape to keep, so its cap is zero and its gradient is scaled
+  // away. A region with no sampled box yet carries no gradient to limit.
+  if (source_gradient_limiter_ && material != MATERIAL_VOID &&
+      srh.extent_min().x <= srh.extent_max().x) {
+    // Offsets from the centroid to the box faces. The centroid is the
+    // length-weighted mean of segment midpoints, all of which lie in the
+    // box, so lo <= 0 <= hi and the dip below is non-negative.
+    Position lo = srh.extent_min() - srh.centroid();
+    Position hi = srh.extent_max() - srh.centroid();
+    for (int g = 0; g < negroups_; g++) {
+      MomentArray& gradient = srh.source_gradients(g);
+      double cap = std::max<double>(srh.source(g), 0.0);
+      double dip = std::max(-gradient.x * lo.x, -gradient.x * hi.x) +
+                   std::max(-gradient.y * lo.y, -gradient.y * hi.y) +
+                   std::max(-gradient.z * lo.z, -gradient.z * hi.z);
+      if (dip > cap) {
+        gradient *= cap / dip;
+      }
+    }
+  }
 }
 
 void LinearSourceDomain::normalize_scalar_flux_and_volumes(

src/random_ray/random_ray.cpp

@@ -642,6 +642,14 @@ void RandomRay::attenuate_flux_linear_source(
     moment_matrix_estimate *= distance;
     srh.mom_matrix() += moment_matrix_estimate;
 
+    // With the source gradient limiter enabled, grow the region's sampled
+    // bounding box with this segment's endpoints, which lie on the region
+    // boundary (or inside it, where the ray starts or ends).
+    if (FlatSourceDomain::source_gradient_limiter_) {
+      srh.expand_extent(r);
+      srh.expand_extent(r + distance * u());
+    }
+
     srh.n_hits() += 1;
   }
 

src/random_ray/random_ray_simulation.cpp

@@ -291,6 +291,7 @@ void openmc_finalize_random_ray()
   FlatSourceDomain::resolved_volume_estimator_ = RandomRayVolumeEstimator::AUTO;
   FlatSourceDomain::volume_normalized_flux_tallies_ = false;
   FlatSourceDomain::adjoint_requested_ = false;
+  FlatSourceDomain::source_gradient_limiter_ = false;
   FlatSourceDomain::solve_ = RandomRaySolve::FORWARD;
   FlatSourceDomain::fw_cadis_local_ = false;
   FlatSourceDomain::fw_cadis_local_targets_.clear();
@@ -693,6 +694,10 @@ void RandomRaySimulation::print_results_random_ray(
       fatal_error("Invalid random ray source shape");
     }
     fmt::print(" Source Shape                      = {}\n", shape);
+    if (RandomRay::source_shape_ != RandomRaySourceShape::FLAT) {
+      fmt::print(" Source Gradient Limiter           = {}\n",
+        FlatSourceDomain::source_gradient_limiter_ ? "ON" : "OFF");
+    }
     std::string sample_method;
     switch (RandomRay::sample_method_) {
     case RandomRaySampleMethod::PRNG: