Random Ray Base Source Region Refactor (#3576)

Author: jtramm

include/openmc/constants.h

@@ -68,6 +68,11 @@ constexpr double MIN_HITS_PER_BATCH {1.5};
 // prevent extremely large adjoint source terms from being generated.
 constexpr double ZERO_FLUX_CUTOFF {1e-22};
 
+// The minimum macroscopic cross section value considered non-void for the
+// random ray solver. Materials with any group with a cross section below this
+// value will be converted to pure void.
+constexpr double MINIMUM_MACRO_XS {1e-6};
+
 // ============================================================================
 // MATH AND PHYSICAL CONSTANTS
 

include/openmc/random_ray/flat_source_domain.h

@@ -27,8 +27,9 @@ class FlatSourceDomain {
 
   //----------------------------------------------------------------------------
   // Methods
-  virtual void update_neutron_source(double k_eff);
-  double compute_k_eff(double k_eff_old) const;
+  virtual void update_single_neutron_source(SourceRegionHandle& srh);
+  virtual void update_all_neutron_sources();
+  void compute_k_eff();
   virtual void normalize_scalar_flux_and_volumes(
     double total_active_distance_per_iteration);
 
@@ -41,7 +42,7 @@ class FlatSourceDomain {
   void output_to_vtk() const;
   void convert_external_sources();
   void count_external_source_regions();
-  void set_adjoint_sources(const vector<double>& forward_flux);
+  void set_adjoint_sources();
   void flux_swap();
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
@@ -54,9 +55,8 @@ class FlatSourceDomain {
     bool is_target_void);
   void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx,
     int32_t target_material_id, bool is_target_void);
-  void prepare_base_source_regions();
   SourceRegionHandle get_subdivided_source_region_handle(
-    int64_t sr, int mesh_bin, Position r, double dist, Direction u);
+    SourceRegionKey sr_key, Position r, Direction u);
   void finalize_discovered_source_regions();
   void apply_transport_stabilization();
   int64_t n_source_regions() const
@@ -67,6 +67,10 @@ class FlatSourceDomain {
   {
     return source_regions_.n_source_regions() * negroups_;
   }
+  int64_t lookup_base_source_region_idx(const GeometryState& p) const;
+  SourceRegionKey lookup_source_region_key(const GeometryState& p) const;
+  int64_t lookup_mesh_bin(int64_t sr, Position r) const;
+  int lookup_mesh_idx(int64_t sr) const;
 
   //----------------------------------------------------------------------------
   // Static Data members
@@ -86,6 +90,7 @@ class FlatSourceDomain {
 
   //----------------------------------------------------------------------------
   // Public Data members
+  double k_eff_ {1.0};              // Eigenvalue
   bool mapped_all_tallies_ {false}; // If all source regions have been visited
 
   int64_t n_external_source_regions_ {0}; // Total number of source regions with
@@ -110,14 +115,6 @@ class FlatSourceDomain {
   // The abstract container holding all source region-specific data
   SourceRegionContainer source_regions_;
 
-  // Base source region container. When source region subdivision via mesh
-  // is in use, this container holds the original (non-subdivided) material
-  // filled cell instance source regions. These are useful as they can be
-  // initialized with external source and mesh domain information ahead of time.
-  // Then, dynamically discovered source regions can be initialized by cloning
-  // their base region.
-  SourceRegionContainer base_source_regions_;
-
   // Parallel hash map holding all source regions discovered during
   // a single iteration. This is a threadsafe data structure that is cleaned
   // out after each iteration and stored in the "source_regions_" container.
@@ -134,8 +131,17 @@ class FlatSourceDomain {
   // 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_;
+  std::unordered_map<SourceRegionKey, vector<int>, SourceRegionKey::HashFunctor>
+    external_point_source_map_;
+
+  // Map that relates a base source region index to the external source index.
+  // This map is used to check if there are any volumetric sources within a
+  // subdivided source region at the time it is discovered.
+  std::unordered_map<int64_t, vector<int>> external_volumetric_source_map_;
+
+  // Map that relates a base source region index to a mesh index. This map
+  // is used to check which subdivision mesh is present in a source region.
+  std::unordered_map<int64_t, int> mesh_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
@@ -147,12 +153,11 @@ class FlatSourceDomain {
   //----------------------------------------------------------------------------
   // Methods
   void apply_external_source_to_source_region(
-    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);
-  void apply_external_source_to_cell_and_children(int32_t i_cell,
-    Discrete* discrete, double strength_factor, int32_t target_material_id);
+    int src_idx, SourceRegionHandle& srh);
+  void apply_external_source_to_cell_instances(int32_t i_cell, int src_idx,
+    int target_material_id, const vector<int32_t>& instances);
+  void apply_external_source_to_cell_and_children(
+    int32_t i_cell, int src_idx, int32_t target_material_id);
   virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g);
   void set_flux_to_source(int64_t sr, int g);
   virtual void set_flux_to_old_flux(int64_t sr, int g);

include/openmc/random_ray/linear_source_domain.h

@@ -20,7 +20,7 @@ class LinearSourceDomain : public FlatSourceDomain {
 public:
   //----------------------------------------------------------------------------
   // Methods
-  void update_neutron_source(double k_eff) override;
+  void update_single_neutron_source(SourceRegionHandle& srh) override;
   void normalize_scalar_flux_and_volumes(
     double total_active_distance_per_iteration) override;
 

include/openmc/random_ray/random_ray.h

@@ -48,7 +48,6 @@ class RandomRay : public Particle {
   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 bool mesh_subdivision_enabled_;     // Flag for mesh subdivision
   static RandomRaySampleMethod sample_method_; // Flag for sampling method
 
   //----------------------------------------------------------------------------

include/openmc/random_ray/random_ray_simulation.h

@@ -21,11 +21,7 @@ class RandomRaySimulation {
   // Methods
   void compute_segment_correction_factors();
   void apply_fixed_sources_and_mesh_domains();
-  void prepare_fixed_sources_adjoint(vector<double>& forward_flux,
-    SourceRegionContainer& forward_source_regions,
-    SourceRegionContainer& forward_base_source_regions,
-    std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>&
-      forward_source_region_map);
+  void prepare_fixed_sources_adjoint();
   void simulate();
   void output_simulation_results() const;
   void instability_check(
@@ -45,9 +41,6 @@ class RandomRaySimulation {
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;
 
-  // Random ray eigenvalue
-  double k_eff_ {1.0};
-
   // Tracks the average FSR miss rate for analysis and reporting
   double avg_miss_rate_ {0.0};
 

include/openmc/random_ray/source_region.h

@@ -308,7 +308,6 @@ class SourceRegion {
   //----------------------------------------------------------------------------
   // Constructors
   SourceRegion(int negroups, bool is_linear);
-  SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr);
   SourceRegion() = default;
 
   //----------------------------------------------------------------------------