Random Ray Source Region Refactor

CMakeLists.txt

@@ -382,6 +382,7 @@ list(APPEND libopenmc_SOURCES
   src/random_ray/flat_source_domain.cpp
   src/random_ray/linear_source_domain.cpp
   src/random_ray/moment_matrix.cpp
+  src/random_ray/source_region.cpp
   src/reaction.cpp
   src/reaction_product.cpp
   src/scattdata.cpp

include/openmc/openmp_interface.h

@@ -36,13 +36,15 @@ inline int thread_num()
 
 class OpenMPMutex {
 public:
-  OpenMPMutex()
+  void init()
   {
 #ifdef _OPENMP
     omp_init_lock(&mutex_);
 #endif
   }
 
+  OpenMPMutex() { init(); }
+
   ~OpenMPMutex()
   {
 #ifdef _OPENMP
@@ -62,14 +64,10 @@ class OpenMPMutex {
   // rather, it produces two different mutexes.
 
   // Copy constructor
-  OpenMPMutex(const OpenMPMutex& other) { OpenMPMutex(); }
+  OpenMPMutex(const OpenMPMutex& other) { init(); }
 
   // Copy assignment operator
-  OpenMPMutex& operator=(const OpenMPMutex& other)
-  {
-    OpenMPMutex();
-    return *this;
-  }
+  OpenMPMutex& operator=(const OpenMPMutex& other) { return *this; }
 
   //! Lock the mutex.
   //

include/openmc/random_ray/flat_source_domain.h

@@ -4,83 +4,12 @@
 #include "openmc/constants.h"
 #include "openmc/openmp_interface.h"
 #include "openmc/position.h"
+#include "openmc/random_ray/source_region.h"
 #include "openmc/source.h"
+#include <unordered_set>
 
 namespace openmc {
 
-//----------------------------------------------------------------------------
-// Helper Functions
-
-// The hash_combine function is the standard hash combine function from boost
-// that is typically used for combining multiple hash values into a single hash
-// as is needed for larger objects being stored in a hash map. The function is
-// taken from:
-// https://www.boost.org/doc/libs/1_55_0/doc/html/hash/reference.html#boost.hash_combine
-// which carries the following license:
-//
-// Boost Software License - Version 1.0 - August 17th, 2003
-// Permission is hereby granted, free of charge, to any person or organization
-// obtaining a copy of the software and accompanying documentation covered by
-// this license (the "Software") to use, reproduce, display, distribute,
-// execute, and transmit the Software, and to prepare derivative works of the
-// Software, and to permit third-parties to whom the Software is furnished to
-// do so, all subject to the following:
-// The copyright notices in the Software and this entire statement, including
-// the above license grant, this restriction and the following disclaimer,
-// must be included in all copies of the Software, in whole or in part, and
-// all derivative works of the Software, unless such copies or derivative
-// works are solely in the form of machine-executable object code generated by
-// a source language processor.
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
-// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
-// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
-// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
-// DEALINGS IN THE SOFTWARE.
-inline void hash_combine(size_t& seed, const size_t v)
-{
-  seed ^= (v + 0x9e3779b9 + (seed << 6) + (seed >> 2));
-}
-
-//----------------------------------------------------------------------------
-// Helper Structs
-
-// A mapping object that is used to map between a specific random ray
-// source region and an OpenMC native tally bin that it should score to
-// every iteration.
-struct TallyTask {
-  int tally_idx;
-  int filter_idx;
-  int score_idx;
-  int score_type;
-  TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type)
-    : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx),
-      score_type(score_type)
-  {}
-  TallyTask() = default;
-
-  // Comparison and Hash operators are defined to allow usage of the
-  // TallyTask struct as a key in an unordered_set
-  bool operator==(const TallyTask& other) const
-  {
-    return tally_idx == other.tally_idx && filter_idx == other.filter_idx &&
-           score_idx == other.score_idx && score_type == other.score_type;
-  }
-
-  struct HashFunctor {
-    size_t operator()(const TallyTask& task) const
-    {
-      size_t seed = 0;
-      hash_combine(seed, task.tally_idx);
-      hash_combine(seed, task.filter_idx);
-      hash_combine(seed, task.score_idx);
-      hash_combine(seed, task.score_type);
-      return seed;
-    }
-  };
-};
-
 /*
  * The FlatSourceDomain class encompasses data and methods for storing
  * scalar flux and source region for all flat source regions in a
@@ -108,15 +37,16 @@ class FlatSourceDomain {
   void random_ray_tally();
   virtual void accumulate_iteration_flux();
   void output_to_vtk() const;
-  virtual void all_reduce_replicated_source_regions();
+  void all_reduce_replicated_source_regions();
   void convert_external_sources();
   void count_external_source_regions();
   void set_adjoint_sources(const vector<double>& forward_flux);
-  virtual void flux_swap();
+  void flux_swap();
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
   void flatten_xs();
   void transpose_scattering_matrix();
+  void serialize_final_fluxes(vector<double>& flux);
 
   //----------------------------------------------------------------------------
   // Static Data members
@@ -129,7 +59,6 @@ class FlatSourceDomain {
 
   //----------------------------------------------------------------------------
   // Public Data members
-
   bool mapped_all_tallies_ {false}; // If all source regions have been visited
 
   int64_t n_source_regions_ {0}; // Total number of source regions in the model
@@ -140,22 +69,6 @@ class FlatSourceDomain {
   // in model::cells
   vector<int64_t> source_region_offsets_;
 
-  // 1D arrays representing values for all source regions
-  vector<OpenMPMutex> lock_;
-  vector<double> volume_;
-  vector<double> volume_t_;
-  vector<int> position_recorded_;
-  vector<Position> position_;
-
-  // 2D arrays stored in 1D representing values for all source regions x energy
-  // groups
-  vector<double> scalar_flux_old_;
-  vector<double> scalar_flux_new_;
-  vector<float> source_;
-  vector<float> external_source_;
-  vector<bool> external_source_present_;
-  vector<double> scalar_flux_final_;
-
   // 2D arrays stored in 1D representing values for all materials x energy
   // groups
   int n_materials_;
@@ -168,20 +81,22 @@ class FlatSourceDomain {
   // groups x energy groups
   vector<double> sigma_s_;
 
+  // The abstract container holding all source region-specific data
+  SourceRegionContainer source_regions_;
+
 protected:
   //----------------------------------------------------------------------------
   // Methods
   void apply_external_source_to_source_region(
-    Discrete* discrete, double strength_factor, int64_t source_region);
+    Discrete* discrete, double strength_factor, int64_t sr);
   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);
-  virtual void set_flux_to_flux_plus_source(
-    int64_t idx, double volume, int material, int g);
-  void set_flux_to_source(int64_t idx);
-  virtual void set_flux_to_old_flux(int64_t idx);
+  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);
 
   //----------------------------------------------------------------------------
   // Private data members
@@ -193,20 +108,6 @@ class FlatSourceDomain {
     simulation_volume_; // Total physical volume of the simulation domain, as
                         // defined by the 3D box of the random ray source
 
-  // 2D array representing values for all source elements x tally
-  // tasks
-  vector<vector<TallyTask>> tally_task_;
-
-  // 1D array representing values for all source regions, with each region
-  // containing a set of volume tally tasks. This more complicated data
-  // structure is convenient for ensuring that volumes are only tallied once per
-  // source region, regardless of how many energy groups are used for tallying.
-  vector<std::unordered_set<TallyTask, TallyTask::HashFunctor>> volume_task_;
-
-  // 1D arrays representing values for all source regions
-  vector<int> material_;
-  vector<double> volume_naive_;
-
   // Volumes for each tally and bin/score combination. This intermediate data
   // structure is used when tallying quantities that must be normalized by
   // volume (i.e., flux). The vector is index by tally index, while the inner 2D
@@ -250,15 +151,6 @@ T convert_to_big_endian(T in)
   return out;
 }
 
-template<typename T>
-void parallel_fill(vector<T>& arr, T value)
-{
-#pragma omp parallel for schedule(static)
-  for (int i = 0; i < arr.size(); i++) {
-    arr[i] = value;
-  }
-}
-
 } // namespace openmc
 
 #endif // OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H

include/openmc/random_ray/linear_source_domain.h

@@ -18,48 +18,22 @@ namespace openmc {
 
 class LinearSourceDomain : public FlatSourceDomain {
 public:
-  //----------------------------------------------------------------------------
-  // Constructors
-  LinearSourceDomain();
-
   //----------------------------------------------------------------------------
   // Methods
   void update_neutron_source(double k_eff) override;
-  double compute_k_eff(double k_eff_old) const;
   void normalize_scalar_flux_and_volumes(
     double total_active_distance_per_iteration) override;
 
   void batch_reset() override;
-  void convert_source_regions_to_tallies();
-  void reset_tally_volumes();
-  void random_ray_tally();
   void accumulate_iteration_flux() override;
   void output_to_vtk() const;
-  void all_reduce_replicated_source_regions() override;
-  void convert_external_sources();
-  void count_external_source_regions();
-  void flux_swap() override;
   double evaluate_flux_at_point(Position r, int64_t sr, int g) const override;
 
-  //----------------------------------------------------------------------------
-  // Public Data members
-
-  vector<MomentArray> source_gradients_;
-  vector<MomentArray> flux_moments_old_;
-  vector<MomentArray> flux_moments_new_;
-  vector<MomentArray> flux_moments_t_;
-  vector<Position> centroid_;
-  vector<Position> centroid_iteration_;
-  vector<Position> centroid_t_;
-  vector<MomentMatrix> mom_matrix_;
-  vector<MomentMatrix> mom_matrix_t_;
-
 protected:
   //----------------------------------------------------------------------------
   // Methods
-  void set_flux_to_flux_plus_source(
-    int64_t idx, double volume, int material, int g) override;
-  void set_flux_to_old_flux(int64_t idx) override;
+  void set_flux_to_flux_plus_source(int64_t sr, double volume, int g) override;
+  void set_flux_to_old_flux(int64_t sr, int g) override;
 
 }; // class LinearSourceDomain