Implement S2 directional sampling in the random ray solver (#3811)

include/openmc/constants.h

@@ -366,7 +366,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 };
+enum class RandomRaySampleMethod { PRNG, HALTON, S2 };
 
 //==============================================================================
 // Geometry Constants

include/openmc/random_ray/random_ray.h

@@ -41,6 +41,7 @@ public:
   uint64_t transport_history_based_single_ray();
   SourceSite sample_prng();
   SourceSite sample_halton();
+  SourceSite sample_s2();
 
   //----------------------------------------------------------------------------
   // Static data members

openmc/settings.py

@@ -207,7 +207,11 @@ class Settings:
             default is 'False'.
         :sample_method:
             Sampling method for the ray starting location and direction of
-            travel. Options are `prng` (default) or 'halton`.
+            travel. Options are `prng` (default), `halton`, or `s2`. `s2`
+            modifies the `prng` sampling method such that rays are sampled
+            with directions (-1, 0, 0) or (1, 0, 0). This is used for verification
+            against analytic transport benchmarks which are often derivied with
+            a reduced angular domain.
         :source_region_meshes:
             List of tuples where each tuple contains a mesh and a list of
             domains. Each domain is an instance of openmc.Material, openmc.Cell,
@@ -1351,7 +1355,7 @@ class Settings:
                                 'openmc.Material, openmc.Cell, or openmc.Universe.')
             elif key == 'sample_method':
                 cv.check_value('sample method', value,
-                               ('prng', 'halton'))
+                               ('prng', 'halton', 's2'))
             elif key == 'diagonal_stabilization_rho':
                 cv.check_type('diagonal stabilization rho', value, Real)
                 cv.check_greater_than('diagonal stabilization rho',

src/random_ray/random_ray.cpp

@@ -792,6 +792,9 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
   case RandomRaySampleMethod::HALTON:
     site = sample_halton();
     break;
+  case RandomRaySampleMethod::S2:
+    site = sample_s2();
+    break;
   default:
     fatal_error("Unknown sample method for random ray transport.");
   }
@@ -874,4 +877,27 @@ SourceSite RandomRay::sample_halton()
   return site;
 }
 
+SourceSite RandomRay::sample_s2()
+{
+  // 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;
+
+  // Get spatial component of the ray_source_
+  SpatialDistribution* space =
+    dynamic_cast<IndependentSource*>(RandomRay::ray_source_.get())->space();
+
+  SourceSite site;
+
+  // Sample spatial distribution
+  site.r = space->sample(current_seed()).first;
+
+  // Sample either left or right for S2 (flashlight) transport.
+  site.u = {prn(current_seed()) < 0.5 ? -1 : 1, 0.0, 0.0};
+
+  return site;
+}
+
 } // namespace openmc

src/random_ray/random_ray_simulation.cpp

@@ -576,9 +576,18 @@ void RandomRaySimulation::print_results_random_ray(
       fatal_error("Invalid random ray source shape");
     }
     fmt::print(" Source Shape                      = {}\n", shape);
-    std::string sample_method =
-      (RandomRay::sample_method_ == RandomRaySampleMethod::PRNG) ? "PRNG"
-                                                                 : "Halton";
+    std::string sample_method;
+    switch (RandomRay::sample_method_) {
+    case RandomRaySampleMethod::PRNG:
+      sample_method = "PRNG";
+      break;
+    case RandomRaySampleMethod::HALTON:
+      sample_method = "Halton";
+      break;
+    case RandomRaySampleMethod::S2:
+      sample_method = "PRNG S2";
+      break;
+    }
     fmt::print(" Sample Method                     = {}\n", sample_method);
 
     if (domain_->is_transport_stabilization_needed_) {

src/settings.cpp

@@ -326,6 +326,8 @@ void get_run_parameters(pugi::xml_node node_base)
         RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
       } else if (temp_str == "halton") {
         RandomRay::sample_method_ = RandomRaySampleMethod::HALTON;
+      } else if (temp_str == "s2") {
+        RandomRay::sample_method_ = RandomRaySampleMethod::S2;
       } else {
         fatal_error("Unrecognized sample method: " + temp_str);
       }

tests/regression_tests/random_ray_s2/inputs_true.dat

@@ -0,0 +1,69 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="domain">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="domain"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="5 -6 1 -2 3 -4" universe="1"/>
+    <surface id="1" name="minimum y" type="y-plane" boundary="reflective" coeffs="-5.0"/>
+    <surface id="2" name="maximum y" type="y-plane" boundary="reflective" coeffs="5.0"/>
+    <surface id="3" name="minimum z" type="z-plane" boundary="reflective" coeffs="-5.0"/>
+    <surface id="4" name="maximum z" type="z-plane" boundary="reflective" coeffs="5.0"/>
+    <surface id="5" type="x-plane" boundary="vacuum" coeffs="0.0"/>
+    <surface id="6" type="x-plane" boundary="vacuum" coeffs="40.0"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>100</particles>
+    <batches>30</batches>
+    <inactive>10</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>0.0 -5.0 -5.0 40.0 5.0 5.0</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>1000.0 1.0</parameters>
+      </energy>
+      <constraints>
+        <domain_type>material</domain_type>
+        <domain_ids>1</domain_ids>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_inactive>100.0</distance_inactive>
+      <distance_active>400.0</distance_active>
+      <source type="independent" strength="1.0" particle="neutron">
+        <space type="box">
+          <parameters>0.0 -5.0 -5.0 40.0 5.0 5.0</parameters>
+        </space>
+      </source>
+      <source_shape>flat</source_shape>
+      <sample_method>s2</sample_method>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="1" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+    </random_ray>
+    <mesh id="1">
+      <dimension>10 1 1</dimension>
+      <lower_left>0.0 -5.0 -5.0</lower_left>
+      <upper_right>40.0 5.0 5.0</upper_right>
+    </mesh>
+  </settings>
+  <tallies>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <tally id="1" name="LR">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/random_ray_s2/results_true.dat

@@ -0,0 +1,21 @@
+tally 1:
+1.153280E+01
+6.650271E+00
+1.413926E+01
+9.995935E+00
+1.579543E+01
+1.247479E+01
+1.676986E+01
+1.406141E+01
+1.722054E+01
+1.482734E+01
+1.722054E+01
+1.482734E+01
+1.676986E+01
+1.406141E+01
+1.579543E+01
+1.247479E+01
+1.413926E+01
+9.995935E+00
+1.153280E+01
+6.650271E+00

tests/regression_tests/random_ray_s2/test.py

@@ -0,0 +1,82 @@
+import os
+
+import openmc
+import openmc.model
+import numpy as np
+
+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_s2():
+    NUM_SOURCE_REGIONS = 10
+    L = 40.0
+
+    # Make the simple MGXS library and material.
+    groups = openmc.mgxs.EnergyGroups(group_edges=[1e-5, 20.0e6])
+
+    domain_mat_data = openmc.XSdata('domain', groups)
+    domain_mat_data.order = 0
+    domain_mat_data.set_total([0.1])
+    domain_mat_data.set_absorption([0.1])
+    domain_mat_data.set_scatter_matrix(np.rollaxis(np.array([[[0.0]]]), 0, 3))
+    mg_cross_sections_file = openmc.MGXSLibrary(groups)
+    mg_cross_sections_file.add_xsdatas([domain_mat_data])
+    mg_cross_sections_file.export_to_hdf5()
+
+    domain_data = openmc.Macroscopic('domain')
+    domain_mat = openmc.Material(name='domain')
+    domain_mat.set_density('macro', 1.0)
+    domain_mat.add_macroscopic(domain_data)
+
+    container = openmc.model.RectangularPrism(width=10.0, height=10.0, axis='x',
+                                              origin=(0.0, 0.0), boundary_type='reflective')
+    left = openmc.XPlane(x0 = 0.0, boundary_type='vacuum')
+    right = openmc.XPlane(x0 = L, boundary_type='vacuum')
+    cell = [openmc.Cell(region = +left & -right & -container, fill = domain_mat)]
+
+    model = openmc.model.Model()
+    model.geometry = openmc.Geometry(root=openmc.Universe(cells=cell))
+    model.materials = openmc.Materials([domain_mat])
+    model.materials.cross_sections = './mgxs.h5'
+
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (NUM_SOURCE_REGIONS, 1, 1)
+    mesh.lower_left = (0.0, -5.0, -5.0)
+    mesh.upper_right = (L, 5.0, 5.0)
+
+    tally = openmc.Tally(name="LR")
+    tally.filters = [openmc.MeshFilter(mesh)]
+    tally.scores = ['flux']
+    tally.estimator = 'tracklength'
+    model.tallies.append(tally)
+
+    uniform_dist = openmc.stats.Box((0.0, -5.0, -5.0), (L, 5.0, 5.0))
+    model.settings.source = [
+      openmc.IndependentSource(space=uniform_dist,
+                              energy=openmc.stats.Discrete(x = 1e3, p = 1.0),
+                              constraints={'domains' : [domain_mat]})
+    ]
+    model.settings.energy_mode = "multi-group"
+    model.settings.batches = 30
+    model.settings.inactive = 10
+    model.settings.particles = 100
+    model.settings.run_mode = 'fixed source'
+    model.settings.random_ray['distance_inactive'] = 100.0
+    model.settings.random_ray['distance_active'] = 400.0
+    model.settings.random_ray['ray_source'] = openmc.IndependentSource(space=uniform_dist)
+    model.settings.random_ray['source_shape'] = 'flat'
+    model.settings.random_ray['sample_method'] = 's2'
+    model.settings.random_ray['source_region_meshes'] = [(mesh, [model.geometry.root_universe])]
+
+    model.export_to_model_xml()
+
+    harness = MGXSTestHarness('statepoint.30.h5', model)
+    harness.main()