Random Ray Point Source Locator (#3360)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Author: jtramm

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

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

include/openmc/random_ray/flat_source_domain.h

@@ -131,6 +131,12 @@ class FlatSourceDomain {
   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 @@ class FlatSourceDomain {
   //----------------------------------------------------------------------------
   // 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);

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()));
       }
-    } 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);
+      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::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;
 }
 

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