Randomized Quasi-Monte Carlo Sampling in The Random Ray Method

docs/source/io_formats/settings.rst

@@ -439,6 +439,12 @@ found in the :ref:`random ray user guide <random_ray>`.
 
     *Default*: None
 
+  :sample_method:
+    Specifies the method for sampling the starting ray distribution. This
+    element can be set to "prng" or "halton".
+
+    *Default*: prng
+
 ----------------------------------
 ``<resonance_scattering>`` Element
 ----------------------------------

docs/source/usersguide/random_ray.rst

@@ -299,6 +299,24 @@ acceptable ray source for a two-dimensional 2x2 lattice would look like:
     provide physical particle fixed sources in addition to the random ray
     source.
 
+--------------------------
+Quasi-Monte Carlo Sampling
+--------------------------
+
+By default OpenMC will use a pseudorandom number generator (PRNG) to sample ray
+starting locations from a uniform distribution in space and angle.
+Alternatively, a randomized Halton sequence may be sampled from, which is a form
+of Randomized Qusi-Monte Carlo (RQMC) sampling. RQMC sampling with random ray
+has been shown to offer reduced variance as compared to regular PRNG sampling,
+as the Halton sequence offers a more uniform distribution of sampled points.
+Randomized Halton sampling can be enabled as::
+
+  settings.random_ray['sample_method'] = 'halton'
+
+Default behavior using OpenMC's native PRNG can be manually specified as::
+
+  settings.random_ray['sample_method'] = 'prng'
+
 .. _subdivision_fsr:
 
 ----------------------------------

include/openmc/constants.h

@@ -344,6 +344,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY };
 
 enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };
 enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
+enum class RandomRaySampleMethod { PRNG, HALTON };
 
 //==============================================================================
 // Geometry Constants

include/openmc/random_dist.h

@@ -16,6 +16,17 @@ namespace openmc {
 
 double uniform_distribution(double a, double b, uint64_t* seed);
 
+//==============================================================================
+//! Sample an integer from uniform distribution [a, b]
+//
+//! \param a Lower bound of uniform distribution
+//! \param b Upper bound of uniform distribtion
+//! \param seed A pointer to the pseudorandom seed
+//! \return Sampled variate
+//==============================================================================
+
+int64_t uniform_int_distribution(int64_t a, int64_t b, uint64_t* seed);
+
 //==============================================================================
 //! Samples an energy from the Maxwell fission distribution based on a direct
 //! sampling scheme.

include/openmc/random_ray/random_ray.h

@@ -31,13 +31,16 @@ class RandomRay : public Particle {
 
   void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain);
   uint64_t transport_history_based_single_ray();
+  SourceSite sample_prng();
+  SourceSite sample_halton();
 
   //----------------------------------------------------------------------------
   // Static data members
   static double distance_inactive_;          // Inactive (dead zone) ray length
   static double distance_active_;            // Active ray length
   static unique_ptr<Source> ray_source_;     // Starting source for ray sampling
   static RandomRaySourceShape source_shape_; // Flag for linear source
+  static RandomRaySampleMethod sample_method_; // Flag for sampling method
 
   //----------------------------------------------------------------------------
   // Public data members

openmc/settings.py

@@ -173,6 +173,9 @@ class Settings:
         :adjoint:
             Whether to run the random ray solver in adjoint mode (bool). The
             default is 'False'.
+        :sample_method:
+            Sampling method for the ray starting location and direction of travel.
+            Options are `prng` (default) or 'halton`.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -1131,6 +1134,9 @@ def random_ray(self, random_ray: dict):
                 cv.check_type('volume normalized flux tallies', random_ray[key], bool)
             elif key == 'adjoint':
                 cv.check_type('adjoint', random_ray[key], bool)
+            elif key == 'sample_method':
+                cv.check_value('sample method', random_ray[key],
+                               ('prng', 'halton'))
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -1948,6 +1954,8 @@ def _random_ray_from_xml_element(self, root):
                     self.random_ray['adjoint'] = (
                         child.text in ('true', '1')
                     )
+                elif child.tag == 'sample_method':
+                    self.random_ray['sample_method'] = child.text
 
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.

src/random_dist.cpp

@@ -12,6 +12,11 @@ double uniform_distribution(double a, double b, uint64_t* seed)
   return a + (b - a) * prn(seed);
 }
 
+int64_t uniform_int_distribution(int64_t a, int64_t b, uint64_t* seed)
+{
+  return a + static_cast<int64_t>(prn(seed) * (b - a + 1));
+}
+
 double maxwell_spectrum(double T, uint64_t* seed)
 {
   // Set the random numbers

src/random_ray/random_ray.cpp

@@ -9,6 +9,9 @@
 #include "openmc/search.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
+
+#include "openmc/distribution_spatial.h"
+#include "openmc/random_dist.h"
 #include "openmc/source.h"
 
 namespace openmc {
@@ -174,6 +177,65 @@ float exponentialG2(float tau)
   return num / den;
 }
 
+// Implementation of the Fisher-Yates shuffle algorithm.
+// Algorithm adapted from:
+//    https://en.cppreference.com/w/cpp/algorithm/random_shuffle#Version_3
+void fisher_yates_shuffle(vector<int>& arr, uint64_t* seed)
+{
+  // Loop over the array from the last element down to the second
+  for (size_t i = arr.size() - 1; i > 0; --i) {
+    // Generate a random index in the range [0, i]
+    size_t j = uniform_int_distribution(0, i, seed);
+    // Swap arr[i] with arr[j]
+    std::swap(arr[i], arr[j]);
+  }
+}
+
+// Function to generate randomized Halton sequence samples
+//
+// Algorithm adapted from:
+//      A. B. Owen. A randomized halton algorithm in r. Arxiv, 6 2017.
+//      URL https://arxiv.org/abs/1706.02808
+vector<vector<float>> rhalton(int N, int dim, uint64_t* seed, int64_t skip = 0)
+{
+  int b;
+  double b2r;
+  vector<double> ans(N);
+  vector<int> ind(N);
+  vector<int> primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
+  vector<vector<float>> halton(N, vector<float>(dim, 0.0));
+
+  std::iota(ind.begin(), ind.end(), skip);
+
+  for (int D = 0; D < dim; ++D) {
+    b = primes[D];
+    b2r = 1.0 / b;
+    vector<int> res(ind);
+    std::fill(ans.begin(), ans.end(), 0.0);
+
+    while ((1.0 - b2r) < 1.0) {
+      vector<int> dig(N);
+      // randomaly permute a sequence from skip to skip+N
+      vector<int> perm(b);
+      std::iota(perm.begin(), perm.end(), 0);
+      fisher_yates_shuffle(perm, seed);
+
+      for (int i = 0; i < N; ++i) {
+        dig[i] = res[i] % b;
+        ans[i] += perm[dig[i]] * b2r;
+        res[i] = (res[i] - dig[i]) / b;
+      }
+      b2r /= b;
+    }
+
+    for (int i = 0; i < N; ++i) {
+      halton[i][D] = ans[i];
+    }
+  }
+
+  return halton;
+}
+
 //==============================================================================
 // RandomRay implementation
 //==============================================================================
@@ -183,6 +245,7 @@ double RandomRay::distance_inactive_;
 double RandomRay::distance_active_;
 unique_ptr<Source> RandomRay::ray_source_;
 RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT};
+RandomRaySampleMethod RandomRay::sample_method_ {RandomRaySampleMethod::PRNG};
 
 RandomRay::RandomRay()
   : angular_flux_(data::mg.num_energy_groups_),
@@ -509,14 +572,19 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
   // set identifier for particle
   id() = simulation::work_index[mpi::rank] + ray_id;
 
-  // set random number seed
-  int64_t particle_seed =
-    (simulation::current_batch - 1) * settings::n_particles + id();
-  init_particle_seeds(particle_seed, seeds());
-  stream() = STREAM_TRACKING;
+  // generate source site using sample method
+  SourceSite site;
+  switch (sample_method_) {
+  case RandomRaySampleMethod::PRNG:
+    site = sample_prng();
+    break;
+  case RandomRaySampleMethod::HALTON:
+    site = sample_halton();
+    break;
+  default:
+    fatal_error("Unknown sample method for random ray transport.");
+  }
 
-  // Sample from ray source distribution
-  SourceSite site {ray_source_->sample(current_seed())};
   site.E = lower_bound_index(
     data::mg.rev_energy_bins_.begin(), data::mg.rev_energy_bins_.end(), site.E);
   site.E = negroups_ - site.E - 1.;
@@ -545,4 +613,53 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
   }
 }
 
+SourceSite RandomRay::sample_prng()
+{
+  // set random number seed
+  int64_t particle_seed =
+    (simulation::current_batch - 1) * settings::n_particles + id();
+  init_particle_seeds(particle_seed, seeds());
+  stream() = STREAM_TRACKING;
+
+  // Sample from ray source distribution
+  SourceSite site {ray_source_->sample(current_seed())};
+
+  return site;
+}
+
+SourceSite RandomRay::sample_halton()
+{
+  SourceSite site;
+
+  // Set random number seed
+  int64_t batch_seed = (simulation::current_batch - 1) * settings::n_particles;
+  int64_t skip = id();
+  init_particle_seeds(batch_seed, seeds());
+  stream() = STREAM_TRACKING;
+
+  // Calculate next samples in LDS
+  vector<vector<float>> samples = rhalton(1, 5, current_seed(), skip = skip);
+
+  // Get spatial box of ray_source_
+  SpatialBox* sb = dynamic_cast<SpatialBox*>(
+    dynamic_cast<IndependentSource*>(RandomRay::ray_source_.get())->space());
+
+  // Sample spatial distribution
+  Position xi {samples[0][0], samples[0][1], samples[0][2]};
+  Position shift {1e-9, 1e-9, 1e-9};
+  site.r = (sb->lower_left() + shift) +
+           xi * ((sb->upper_right() - shift) - (sb->lower_left() + shift));
+
+  // Sample Polar cosine and azimuthal angles
+  float mu = 2.0 * samples[0][3] - 1.0;
+  float azi = 2.0 * PI * samples[0][4];
+  // Convert to Cartesian coordinates
+  float c = std::sqrt(1.0 - mu * mu);
+  site.u.x = mu;
+  site.u.y = std::cos(azi) * c;
+  site.u.z = std::sin(azi) * c;
+
+  return site;
+}
+
 } // namespace openmc

src/settings.cpp

@@ -306,6 +306,17 @@ void get_run_parameters(pugi::xml_node node_base)
       FlatSourceDomain::adjoint_ =
         get_node_value_bool(random_ray_node, "adjoint");
     }
+    if (check_for_node(random_ray_node, "sample_method")) {
+      std::string temp_str =
+        get_node_value(random_ray_node, "sample_method", true, true);
+      if (temp_str == "prng") {
+        RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
+      } else if (temp_str == "halton") {
+        RandomRay::sample_method_ = RandomRaySampleMethod::HALTON;
+      } else {
+        fatal_error("Unrecognized sample method: " + temp_str);
+      }
+    }
   }
 }