From 07f5334616c2ca92ad0a04dfa8c900e080509f70 Mon Sep 17 00:00:00 2001
From: John Tramm <john.tramm@gmail.com>
Date: Wed, 2 Apr 2025 00:40:07 -0500
Subject: [PATCH] Random Ray Point Source Locator (#3360)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 docs/source/methods/random_ray.rst            |   3 +-
 docs/source/usersguide/random_ray.rst         |   9 +-
 .../openmc/random_ray/flat_source_domain.h    |   8 +-
 src/random_ray/flat_source_domain.cpp         | 119 ++++++--
 src/random_ray/random_ray_simulation.cpp      |  23 +-
 .../__init__.py                               |   0
 .../inputs_true.dat                           | 253 ++++++++++++++++++
 .../results_true.dat                          |   9 +
 .../random_ray_point_source_locator/test.py   |  44 +++
 9 files changed, 435 insertions(+), 33 deletions(-)
 create mode 100644 tests/regression_tests/random_ray_point_source_locator/__init__.py
 create mode 100644 tests/regression_tests/random_ray_point_source_locator/inputs_true.dat
 create mode 100644 tests/regression_tests/random_ray_point_source_locator/results_true.dat
 create mode 100644 tests/regression_tests/random_ray_point_source_locator/test.py

diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst
index eb22b0544d..41f95f0a9e 100644
--- a/docs/source/methods/random_ray.rst
+++ b/docs/source/methods/random_ray.rst
@@ -1052,7 +1052,8 @@ random ray and Monte Carlo, however.
   regions. Thus, in the OpenMC implementation of random ray, particle sources
   are restricted to being volumetric and isotropic, although different energy
   spectrums are supported. Fixed sources can be applied to specific materials,
-  cells, or universes.
+  cells, or universes. Point sources are "smeared" to fill the volume of the
+  source region that contains the point source coordinate.
 
 - **Inactive batches:** In Monte Carlo, use of a fixed source implies that all
   batches are active batches, as there is no longer a need to develop a fission
diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst
index ce98211966..138ae910c9 100644
--- a/docs/source/usersguide/random_ray.rst
+++ b/docs/source/usersguide/random_ray.rst
@@ -919,9 +919,12 @@ Monte Carlo solver.
 Currently, all of the following conditions must be met for the particle source
 to be valid in random ray mode:
 
-- One or more domain ids must be specified that indicate which cells, universes,
-  or materials the source applies to. This implicitly limits the source type to
-  being volumetric. This is specified via the ``domains`` constraint placed on the
+- Either a point source must be used, or a domain constraint must be specified
+  that indicates which cells, universes, or materials the source applies to. In
+  either case, this implicitly limits the source type to being volumetric, as
+  even in the point source case the source will be "smeared" throughout the
+  source region that contains the point source coordinate. A source domain is
+  specified via the ``domains`` constraint placed on the
   :class:`openmc.IndependentSource` Python class.
 - The source must be isotropic (default for a source)
 - The source must use a discrete (i.e., multigroup) energy distribution. The
diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h
index 656b0354a7..9fdd558750 100644
--- a/include/openmc/random_ray/flat_source_domain.h
+++ b/include/openmc/random_ray/flat_source_domain.h
@@ -131,6 +131,12 @@ public:
   std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
     source_region_map_;
 
+  // Map that relates a SourceRegionKey to the external source index. This map
+  // is used to check if there are any point sources within a subdivided source
+  // region at the time it is discovered.
+  std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
+    point_source_map_;
+
   // If transport corrected MGXS data is being used, there may be negative
   // in-group scattering cross sections that can result in instability in MOC
   // and random ray if used naively. This flag enables a stabilization
@@ -141,7 +147,7 @@ protected:
   //----------------------------------------------------------------------------
   // Methods
   void apply_external_source_to_source_region(
-    Discrete* discrete, double strength_factor, int64_t sr);
+    Discrete* discrete, double strength_factor, SourceRegionHandle& srh);
   void apply_external_source_to_cell_instances(int32_t i_cell,
     Discrete* discrete, double strength_factor, int target_material_id,
     const vector<int32_t>& instances);
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index ff8e5d2327..8ab30abf20 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@@ -946,17 +946,16 @@ void FlatSourceDomain::output_to_vtk() const
 }
 
 void FlatSourceDomain::apply_external_source_to_source_region(
-  Discrete* discrete, double strength_factor, int64_t sr)
+  Discrete* discrete, double strength_factor, SourceRegionHandle& srh)
 {
-  source_regions_.external_source_present(sr) = 1;
+  srh.external_source_present() = 1;
 
   const auto& discrete_energies = discrete->x();
   const auto& discrete_probs = discrete->prob();
 
   for (int i = 0; i < discrete_energies.size(); i++) {
     int g = data::mg.get_group_index(discrete_energies[i]);
-    source_regions_.external_source(sr, g) +=
-      discrete_probs[i] * strength_factor;
+    srh.external_source(g) += discrete_probs[i] * strength_factor;
   }
 }
 
@@ -980,8 +979,9 @@ void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell,
     if (target_material_id == C_NONE ||
         cell_material_id == target_material_id) {
       int64_t source_region = source_region_offsets_[i_cell] + j;
-      apply_external_source_to_source_region(
-        discrete, strength_factor, source_region);
+      SourceRegionHandle srh =
+        source_regions_.get_source_region_handle(source_region);
+      apply_external_source_to_source_region(discrete, strength_factor, srh);
     }
   }
 }
@@ -1023,34 +1023,88 @@ void FlatSourceDomain::convert_external_sources()
 {
   // Loop over external sources
   for (int es = 0; es < model::external_sources.size(); es++) {
+
+    // Extract source information
     Source* s = model::external_sources[es].get();
     IndependentSource* is = dynamic_cast<IndependentSource*>(s);
     Discrete* energy = dynamic_cast<Discrete*>(is->energy());
     const std::unordered_set<int32_t>& domain_ids = is->domain_ids();
-
     double strength_factor = is->strength();
 
-    if (is->domain_type() == Source::DomainType::MATERIAL) {
-      for (int32_t material_id : domain_ids) {
-        for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) {
-          apply_external_source_to_cell_and_children(
-            i_cell, energy, strength_factor, material_id);
+    // If there is no domain constraint specified, then this must be a point
+    // source. In this case, we need to find the source region that contains the
+    // point source and apply or relate it to the external source.
+    if (is->domain_ids().size() == 0) {
+
+      // Extract the point source coordinate and find the base source region at
+      // that point
+      auto sp = dynamic_cast<SpatialPoint*>(is->space());
+      GeometryState gs;
+      gs.r() = sp->r();
+      gs.r_last() = sp->r();
+      gs.u() = {1.0, 0.0, 0.0};
+      bool found = exhaustive_find_cell(gs);
+      if (!found) {
+        fatal_error(fmt::format("Could not find cell containing external "
+                                "point source at {}",
+          sp->r()));
+      }
+      int i_cell = gs.lowest_coord().cell;
+      int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance();
+
+      if (RandomRay::mesh_subdivision_enabled_) {
+        // If mesh subdivision is enabled, we need to determine which subdivided
+        // mesh bin the point source coordinate is in as well
+        int mesh_idx = source_regions_.mesh(sr);
+        int mesh_bin;
+        if (mesh_idx == C_NONE) {
+          mesh_bin = 0;
+        } else {
+          mesh_bin = model::meshes[mesh_idx]->get_bin(gs.r());
         }
+        // With the source region and mesh bin known, we can use the
+        // accompanying SourceRegionKey as a key into a map that stores the
+        // corresponding external source index for the point source. Notably, we
+        // do not actually apply the external source to any source regions here,
+        // as if mesh subdivision is enabled, they haven't actually been
+        // discovered & initilized yet. When discovered, they will read from the
+        // point_source_map to determine if there are any point source terms
+        // that should be applied.
+        SourceRegionKey key {sr, mesh_bin};
+        point_source_map_[key] = es;
+      } else {
+        // If we are not using mesh subdivision, we can apply the external
+        // source directly to the source region as we do for volumetric domain
+        // constraint sources.
+        SourceRegionHandle srh = source_regions_.get_source_region_handle(sr);
+        apply_external_source_to_source_region(energy, strength_factor, srh);
       }
-    } else if (is->domain_type() == Source::DomainType::CELL) {
-      for (int32_t cell_id : domain_ids) {
-        int32_t i_cell = model::cell_map[cell_id];
-        apply_external_source_to_cell_and_children(
-          i_cell, energy, strength_factor, C_NONE);
-      }
-    } else if (is->domain_type() == Source::DomainType::UNIVERSE) {
-      for (int32_t universe_id : domain_ids) {
-        int32_t i_universe = model::universe_map[universe_id];
-        Universe& universe = *model::universes[i_universe];
-        for (int32_t i_cell : universe.cells_) {
+
+    } else {
+      // If not a point source, then use the volumetric domain constraints to
+      // determine which source regions to apply the external source to.
+      if (is->domain_type() == Source::DomainType::MATERIAL) {
+        for (int32_t material_id : domain_ids) {
+          for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) {
+            apply_external_source_to_cell_and_children(
+              i_cell, energy, strength_factor, material_id);
+          }
+        }
+      } else if (is->domain_type() == Source::DomainType::CELL) {
+        for (int32_t cell_id : domain_ids) {
+          int32_t i_cell = model::cell_map[cell_id];
           apply_external_source_to_cell_and_children(
             i_cell, energy, strength_factor, C_NONE);
         }
+      } else if (is->domain_type() == Source::DomainType::UNIVERSE) {
+        for (int32_t universe_id : domain_ids) {
+          int32_t i_universe = model::universe_map[universe_id];
+          Universe& universe = *model::universes[i_universe];
+          for (int32_t i_cell : universe.cells_) {
+            apply_external_source_to_cell_and_children(
+              i_cell, energy, strength_factor, C_NONE);
+          }
+        }
       }
     }
   } // End loop over external sources
@@ -1399,6 +1453,25 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
     sr_key, {base_source_regions_.get_source_region_handle(sr), sr});
   discovered_source_regions_.unlock(sr_key);
   SourceRegionHandle handle {*sr_ptr};
+
+  // Check if the new source region contains a point source and apply it if so
+  auto it2 = point_source_map_.find(sr_key);
+  if (it2 != point_source_map_.end()) {
+    int es = it2->second;
+    auto s = model::external_sources[es].get();
+    auto is = dynamic_cast<IndependentSource*>(s);
+    auto energy = dynamic_cast<Discrete*>(is->energy());
+    double strength_factor = is->strength();
+    apply_external_source_to_source_region(energy, strength_factor, handle);
+    int material = handle.material();
+    if (material != MATERIAL_VOID) {
+      for (int g = 0; g < negroups_; g++) {
+        double sigma_t = sigma_t_[material * negroups_ + g];
+        handle.external_source(g) /= sigma_t;
+      }
+    }
+  }
+
   return handle;
 }
 
diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp
index 60b00a5be7..b831efcd99 100644
--- a/src/random_ray/random_ray_simulation.cpp
+++ b/src/random_ray/random_ray_simulation.cpp
@@ -281,10 +281,21 @@ void validate_random_ray_inputs()
           "allowed in random ray mode.");
       }
 
-      // Validate that a domain ID was specified
-      if (is->domain_ids().size() == 0) {
-        fatal_error("Fixed sources must be specified by domain "
-                    "id (cell, material, or universe) in random ray mode.");
+      // Validate that a domain ID was specified OR that it is a point source
+      auto sp = dynamic_cast<SpatialPoint*>(is->space());
+      if (is->domain_ids().size() == 0 && !sp) {
+        fatal_error("Fixed sources must be point source or spatially "
+                    "constrained by domain id (cell, material, or universe) in "
+                    "random ray mode.");
+      } else if (is->domain_ids().size() > 0 && sp) {
+        // If both a domain constraint and a non-default point source location
+        // are specified, notify user that domain constraint takes precedence.
+        if (sp->r().x == 0.0 && sp->r().y == 0.0 && sp->r().z == 0.0) {
+          warning("Fixed source has both a domain constraint and a point "
+                  "type spatial distribution. The domain constraint takes "
+                  "precedence in random ray mode -- point source coordinate "
+                  "will be ignored.");
+        }
       }
 
       // Check that a discrete energy distribution was used
@@ -393,12 +404,12 @@ RandomRaySimulation::RandomRaySimulation()
 
 void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
 {
+  domain_->apply_meshes();
   if (settings::run_mode == RunMode::FIXED_SOURCE) {
     // Transfer external source user inputs onto random ray source regions
     domain_->convert_external_sources();
     domain_->count_external_source_regions();
   }
-  domain_->apply_meshes();
 }
 
 void RandomRaySimulation::prepare_fixed_sources_adjoint(
@@ -517,6 +528,8 @@ void RandomRaySimulation::simulate()
     finalize_generation();
     finalize_batch();
   } // End random ray power iteration loop
+
+  domain_->count_external_source_regions();
 }
 
 void RandomRaySimulation::output_simulation_results() const
diff --git a/tests/regression_tests/random_ray_point_source_locator/__init__.py b/tests/regression_tests/random_ray_point_source_locator/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat b/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat
new file mode 100644
index 0000000000..9c7cd4ab5b
--- /dev/null
+++ b/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat
@@ -0,0 +1,253 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="absorber">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="absorber"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="infinite source region" universe="1"/>
+    <cell id="2" material="2" name="infinite void region" universe="2"/>
+    <cell id="3" material="3" name="infinite absorber region" universe="3"/>
+    <cell fill="4" id="4" universe="5"/>
+    <cell fill="5" id="5" name="full domain" region="1 -2 3 -4 5 -6" universe="6"/>
+    <lattice id="4">
+      <pitch>2.5 2.5 2.5</pitch>
+      <dimension>12 12 12</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 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 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 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 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 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 </universes>
+    </lattice>
+    <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="2" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="4" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="6" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>90</particles>
+    <batches>30</batches>
+    <inactive>15</inactive>
+    <source particle="neutron" strength="3.14" type="independent">
+      <space type="point">
+        <parameters>2.5 2.5 2.5</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>100.0 1.0</parameters>
+      </energy>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>500.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="6" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+    </random_ray>
+    <mesh id="1">
+      <dimension>30 30 30</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <upper_right>30.0 30.0 30.0</upper_right>
+    </mesh>
+  </settings>
+  <tallies>
+    <filter id="3" type="material">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="material">
+      <bins>2</bins>
+    </filter>
+    <filter id="1" type="material">
+      <bins>3</bins>
+    </filter>
+    <tally id="3" name="Source Tally">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="2" name="Void Tally">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="1" name="Absorber Tally">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/random_ray_point_source_locator/results_true.dat b/tests/regression_tests/random_ray_point_source_locator/results_true.dat
new file mode 100644
index 0000000000..1785dda574
--- /dev/null
+++ b/tests/regression_tests/random_ray_point_source_locator/results_true.dat
@@ -0,0 +1,9 @@
+tally 1:
+2.633900E+00
+2.948207E+00
+tally 2:
+1.440463E-01
+3.294032E-03
+tally 3:
+9.425207E-03
+1.089748E-05
diff --git a/tests/regression_tests/random_ray_point_source_locator/test.py b/tests/regression_tests/random_ray_point_source_locator/test.py
new file mode 100644
index 0000000000..fd3d8a18fe
--- /dev/null
+++ b/tests/regression_tests/random_ray_point_source_locator/test.py
@@ -0,0 +1,44 @@
+import os
+
+import openmc
+from openmc.examples import random_ray_three_region_cube
+
+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_point_source_locator():
+    model = random_ray_three_region_cube()
+
+    # Overlay subdivided SR mesh to reduce resolution from 2.5cm -> 1cm
+    width = 30.0
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (30, 30, 30)
+    mesh.lower_left = (0.0, 0.0, 0.0)
+    mesh.upper_right = (width, width, width)
+    model.settings.random_ray['source_region_meshes'] = [
+        (mesh, [model.geometry.root_universe]),
+    ]
+
+    # Define a point source
+    strengths = [1.0]
+    midpoints = [100.0]
+    energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
+    spatial_distribution = openmc.stats.Point([2.5, 2.5, 2.5])
+    source = openmc.IndependentSource(
+        energy=energy_distribution, space=spatial_distribution, strength=3.14)
+    model.settings.source = [source]
+
+    # Settings
+    model.settings.inactive = 15
+    model.settings.batches = 30
+
+    harness = MGXSTestHarness('statepoint.30.h5', model)
+    harness.main()