diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 34b73fc55c5..24a7cf22c1f 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -439,6 +439,12 @@ found in the :ref:`random ray user guide `. *Default*: None + :sample_method: + Specifies the method for sampling the starting ray distribution. This + element can be set to "prng" or "halton". + + *Default*: prng + ---------------------------------- ```` Element ---------------------------------- diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 2b9cf67240b..b838daef34d 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/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: ---------------------------------- diff --git a/include/openmc/constants.h b/include/openmc/constants.h index a83a4a07d4b..cfcacbabf00 100644 --- a/include/openmc/constants.h +++ b/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 diff --git a/include/openmc/random_dist.h b/include/openmc/random_dist.h index 0fb186edca0..11e88ab8cce 100644 --- a/include/openmc/random_dist.h +++ b/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. diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index 96c38da7b1c..08c2a8488f2 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -31,6 +31,8 @@ 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 @@ -38,6 +40,7 @@ class RandomRay : public Particle { static double distance_active_; // Active ray length static unique_ptr 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 diff --git a/openmc/settings.py b/openmc/settings.py index 77598b204fc..f2eef17258c 100644 --- a/openmc/settings.py +++ b/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. diff --git a/src/random_dist.cpp b/src/random_dist.cpp index 1aa35a689cf..b05b76f99ec 100644 --- a/src/random_dist.cpp +++ b/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(prn(seed) * (b - a + 1)); +} + double maxwell_spectrum(double T, uint64_t* seed) { // Set the random numbers diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 7a359f35664..28572e62099 100644 --- a/src/random_ray/random_ray.cpp +++ b/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& 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> rhalton(int N, int dim, uint64_t* seed, int64_t skip = 0) +{ + int b; + double b2r; + vector ans(N); + vector ind(N); + vector primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}; + vector> halton(N, vector(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 res(ind); + std::fill(ans.begin(), ans.end(), 0.0); + + while ((1.0 - b2r) < 1.0) { + vector dig(N); + // randomaly permute a sequence from skip to skip+N + vector 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 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> samples = rhalton(1, 5, current_seed(), skip = skip); + + // Get spatial box of ray_source_ + SpatialBox* sb = dynamic_cast( + dynamic_cast(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 diff --git a/src/settings.cpp b/src/settings.cpp index 61eda79967a..60263c84653 100644 --- a/src/settings.cpp +++ b/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); + } + } } } diff --git a/tests/regression_tests/random_ray_halton_samples/__init__.py b/tests/regression_tests/random_ray_halton_samples/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_halton_samples/inputs_true.dat b/tests/regression_tests/random_ray_halton_samples/inputs_true.dat new file mode 100644 index 00000000000..624ab495f41 --- /dev/null +++ b/tests/regression_tests/random_ray_halton_samples/inputs_true.dat @@ -0,0 +1,110 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.126 0.126 + 10 10 + -0.63 -0.63 + +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 + + + 1.26 1.26 + 2 2 + -1.26 -1.26 + +2 2 +2 5 + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 10 + 5 + multi-group + + 100.0 + 20.0 + + + -1.26 -1.26 -1 1.26 1.26 1 + + + True + halton + + + + + 2 2 + -1.26 -1.26 + 1.26 1.26 + + + 1 + + + 1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0 + + + 1 2 + flux fission nu-fission + analog + + + diff --git a/tests/regression_tests/random_ray_halton_samples/results_true.dat b/tests/regression_tests/random_ray_halton_samples/results_true.dat new file mode 100644 index 00000000000..651c884b5dc --- /dev/null +++ b/tests/regression_tests/random_ray_halton_samples/results_true.dat @@ -0,0 +1,171 @@ +k-combined: +8.388050E-01 7.383283E-03 +tally 1: +5.033306E+00 +5.072159E+00 +1.917335E+00 +7.360725E-01 +4.666410E+00 +4.360038E+00 +2.851811E+00 +1.629361E+00 +4.365591E-01 +3.818885E-02 +1.062497E+00 +2.262072E-01 +1.697620E+00 +5.829331E-01 +5.639912E-02 +6.427569E-04 +1.372642E-01 +3.807294E-03 +2.376682E+00 +1.151027E+00 +8.060903E-02 +1.323179E-03 +1.961863E-01 +7.837694E-03 +7.145451E+00 +1.037540E+01 +8.551805E-02 +1.486270E-03 +2.081364E-01 +8.803959E-03 +2.053205E+01 +8.469503E+01 +3.235620E-02 +2.102893E-04 +8.006316E-02 +1.287560E-03 +1.326546E+01 +3.519488E+01 +1.867472E-01 +6.975145E-03 +5.194280E-01 +5.396293E-02 +7.558115E+00 +1.142535E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.386210E+00 +2.294414E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.827274E+00 +6.782303E-01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.702858E+00 +1.489752E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +7.475537E+00 +1.133971E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.828685E+01 +6.719608E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.143600E+01 +2.615734E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.590713E+00 +4.224106E+00 +1.705847E+00 +5.831967E-01 +4.151691E+00 +3.454497E+00 +2.730853E+00 +1.495413E+00 +4.072633E-01 +3.325944E-02 +9.911973E-01 +1.970084E-01 +1.664732E+00 +5.598667E-01 +5.385644E-02 +5.858353E-04 +1.310758E-01 +3.470126E-03 +2.312239E+00 +1.088485E+00 +7.662778E-02 +1.195422E-03 +1.864967E-01 +7.080943E-03 +7.105765E+00 +1.025960E+01 +8.287513E-02 +1.396474E-03 +2.017039E-01 +8.272054E-03 +2.099024E+01 +8.854252E+01 +3.191885E-02 +2.048369E-04 +7.898097E-02 +1.254176E-03 +1.355820E+01 +3.676863E+01 +1.815102E-01 +6.590730E-03 +5.048615E-01 +5.098892E-02 +5.093659E+00 +5.192360E+00 +1.874631E+00 +7.031791E-01 +4.562477E+00 +4.165198E+00 +2.870213E+00 +1.650126E+00 +4.244068E-01 +3.608744E-02 +1.032921E+00 +2.137597E-01 +1.703400E+00 +5.873028E-01 +5.464556E-02 +6.042852E-04 +1.329964E-01 +3.579412E-03 +2.389118E+00 +1.163673E+00 +7.832235E-02 +1.251253E-03 +1.906209E-01 +7.411650E-03 +7.162706E+00 +1.042515E+01 +8.273829E-02 +1.391799E-03 +2.013709E-01 +8.244357E-03 +2.043146E+01 +8.383560E+01 +3.096158E-02 +1.924116E-04 +7.661227E-02 +1.178099E-03 +1.314148E+01 +3.454146E+01 +1.771733E-01 +6.279891E-03 +4.927985E-01 +4.858413E-02 diff --git a/tests/regression_tests/random_ray_halton_samples/test.py b/tests/regression_tests/random_ray_halton_samples/test.py new file mode 100644 index 00000000000..478b6502648 --- /dev/null +++ b/tests/regression_tests/random_ray_halton_samples/test.py @@ -0,0 +1,19 @@ +import os + +from openmc.examples import random_ray_lattice + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + +def test_random_ray_halton_samples(): + model = random_ray_lattice() + model.settings.random_ray['sample_method'] = 'halton' + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main()