Volume Calculation enhancement including refactoring and real-valued scoring implementation

include/openmc/volume_calc.h

@@ -7,13 +7,19 @@
 #include <vector>
 
 #include "openmc/array.h"
+#include "openmc/cell.h"
 #include "openmc/openmp_interface.h"
 #include "openmc/position.h"
 #include "openmc/tallies/trigger.h"
+#include "openmc/timer.h"
 #include "openmc/vector.h"
 
 #include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
+
+#ifdef OPENMC_MPI
+#include <mpi.h>
+#endif
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -28,57 +34,173 @@ class VolumeCalculation {
 
 public:
   // Aliases, types
-  struct Result {
-    array<double, 2> volume;    //!< Mean/standard deviation of volume
+  struct NuclResult {
     vector<int> nuclides;       //!< Index of nuclides
     vector<double> atoms;       //!< Number of atoms for each nuclide
     vector<double> uncertainty; //!< Uncertainty on number of atoms
-    int iterations; //!< Number of iterations needed to obtain the results
-  };                // Results for a single domain
+  }; // Results of nuclides calculation for a single domain
+
+  //! \brief Tally corresponding to a single material (AoS)
+  struct VolTally {
+    double score;               //!< Current batch scores accumulator
+    array<double, 2> score_acc; //!< Scores and squared scores accumulator
+    int32_t index;              //!< Material ID
+
+    inline VolTally(const int i_material = 0, const double contrib = 0.0,
+      const double score_acc_ = 0.0, const double score2_acc_ = 0.0);
+
+    // private:
+    //! \brief Add batch scores means to a tally
+    //! \param[in] batch_size_1 Inversed batch size
+    //! \param[in] vol_tallies  All tallies
+    inline void finalize_batch(const double batch_size_1);
+
+    //! \brief Pass counters data from a given tally to this
+    //! \param[in] vol_tally Data source
+    inline void assign_tally(const VolTally& vol_tally);
+
+    //! \brief Add counters data from a given tally to this
+    //! \param[in] vol_tally Data source
+    inline void append_tally(const VolTally& vol_tally);
+
+    //! \brief Determines given trigger condition satisfaction for this tally
+    //
+    //! \param[in] trigger_type Type of trigger condition
+    //! \param[in] threshold    Value for trigger condition (either volume
+    //! fraction variance or squared rel. err. dependent on the trigger type)
+    //! \param[in] n_samples    Statistics size
+    //! \return                 True if the trigger condition is satisfied
+    inline bool trigger_state(const TriggerMetric trigger_type,
+      const double threshold, const size_t& n_samples) const;
+  };
+
+  //! \brief Online results of calculation specific for each thread
+  struct CalcResults {
+    uint64_t n_samples; //!< Number of samples
+    int iterations;     //!< Number of iterations needed to obtain the results
+    double cost; //!< Product of spent time and number of threads/processes
+    Timer sampling_time; // Timer for measurment of the simulation
+    vector<vector<VolumeCalculation::VolTally>>
+      vol_tallies; //!< Volume tallies for each domain
+    vector<VolumeCalculation::NuclResult>
+      nuc_results; //!< Nuclides of each domain
+
+    CalcResults(const VolumeCalculation& vol_calc);
+
+    //! \brief Reset all counters
+    void reset();
+
+    //! \brief Append another counters to this
+    void append(const CalcResults& other);
+
+#ifdef OPENMC_MPI
+    //! \brief Collects results from all MPI processes to this
+    void collect_MPI();
+#endif
+  };
 
   // Constructors
   VolumeCalculation(pugi::xml_node node);
-
   VolumeCalculation() = default;
 
   // Methods
 
   //! \brief Stochastically determine the volume of a set of domains along with
-  //! the
-  //!   average number densities of nuclides within the domain
+  //!        the average number densities of nuclides within the domain
   //
+  //! \param[in,out] results Results of calculation entity for filling
   //! \return Vector of results for each user-specified domain
-  vector<Result> execute() const;
+  void execute(CalcResults& results) const;
+
+  //! \brief Rejection estimator
+  inline void score_hit(const Particle& p, CalcResults& results) const;
+
+  //! \brief Check whether a material has already been hit for a given domain.
+  //! If not, add new entries to the vectors
+  //
+  //! \param[in] i_material      Index in global materials vector
+  //! \param[in] contrib         Scoring value
+  //! \param[in,out] vol_tallies Vector of tallies corresponding to each
+  //! material
+  void check_hit(const int32_t i_material, const double contrib,
+    vector<VolTally>& vol_tallies) const;
+
+  //! \brief Reduce vector of volumetric tallies from each thread to a single
+  //! copy
+  //
+  //! \param[in] local_results Results specific to each thread
+  //! \param[out] results      Reduced results
+  void reduce_results(
+    const CalcResults& local_results, CalcResults& results) const;
+
+  //! \brief Print volume calculation results
+  //
+  //! \param[in] results   Full volume calculation results
+  void show_results(const CalcResults& results) const;
+
+  //! \brief Prints a statistics parameter
+  //
+  //! \param[in] label   Name of parameter
+  //! \param[in] units   Units of measure
+  //! \param[in] value   Value of parameter
+  void show_vol_stat(
+    const std::string label, const std::string units, const double value) const;
+
+  //! \brief Prints volume result for a domain
+  //
+  //! \param[in] domain_type  Either material or cell, ect.
+  //! \param[in] domain_id    Number of this domain
+  //! \param[in] region_name  Domain description
+  //! \param[in] mean         Center of confidence interval
+  //! \param[in] stddev       Half-width of confidence interval
+  void show_volume(const std::string domain_type, const int domain_id,
+    const std::string region_name, const double mean,
+    const double stddev) const;
 
   //! \brief Write volume calculation results to HDF5 file
   //
   //! \param[in] filename  Path to HDF5 file to write
-  //! \param[in] results   Vector of results for each domain
-  void to_hdf5(
-    const std::string& filename, const vector<Result>& results) const;
+  //! \param[in] results   Results entity
+  void to_hdf5(const std::string& filename, const CalcResults& results) const;
 
   // Tally filter and map types
   enum class TallyDomain { UNIVERSE, MATERIAL, CELL };
 
   // Data members
-  TallyDomain domain_type_; //!< Type of domain (cell, material, etc.)
-  size_t n_samples_;        //!< Number of samples to use
-  double threshold_ {-1.0}; //!< Error threshold for domain volumes
+  TallyDomain domain_type_;      //!< Type of domain (cell, material, etc.)
+  size_t n_samples_;             //!< Number of samples to use
+  double volume_sample_;         //!< Volume of bounding primitive
+  double threshold_ {-1.0};      //!< Error threshold for domain volumes
+  double threshold_cnd_;         //!< Pre-computed value for trigger condition
+  int max_iterations_ {INT_MAX}; //!< Limit of iterations number (necessary
+                                 //!< maximum value of data type by default)
   TriggerMetric trigger_type_ {
     TriggerMetric::not_active}; //!< Trigger metric for the volume calculation
   Position lower_left_;         //!< Lower-left position of bounding box
   Position upper_right_;        //!< Upper-right position of bounding box
   vector<int> domain_ids_;      //!< IDs of domains to find volumes of
 
+#ifdef OPENMC_MPI
+  //! \brief Creates MPI structs for passing data between threads
+  void initialize_MPI_struct() const;
+
+  //! \brief Deletes MPI structs for passing data between threads
+  void delete_MPI_struct() const;
+#endif
+
 private:
-  //! \brief Check whether a material has already been hit for a given domain.
-  //! If not, add new entries to the vectors
+  constexpr static int _INDEX_TOTAL =
+    -999; //!< Index of zero-element tally for entire domain totals should be
+          //!< out of material ID range
+
+  //! \brief Computes estimated mean and std.dev. for a tally
   //
-  //! \param[in] i_material Index in global materials vector
-  //! \param[in,out] indices Vector of material indices
-  //! \param[in,out] hits Number of hits corresponding to each material
-  void check_hit(
-    int i_material, vector<uint64_t>& indices, vector<uint64_t>& hits) const;
+  //! \param[in] n_samples    Statistic's size
+  //! \param[in] coeff_norm   Normalization coefficient to multiply
+  //! \param[in] vol_tally    Tally
+  //! \return                 Array of mean and stddev
+  array<double, 2> get_tally_results(const size_t& n_samples,
+    const double coeff_norm, const VolTally& vol_tally) const;
 };
 
 //==============================================================================
@@ -93,35 +215,6 @@ extern vector<VolumeCalculation> volume_calcs;
 // Non-member functions
 //==============================================================================
 
-//! Reduce vector of indices and hits from each thread to a single copy
-//
-//! \param[in] local_indices Indices specific to each thread
-//! \param[in] local_hits Hit count specific to each thread
-//! \param[out] indices Reduced vector of indices
-//! \param[out] hits Reduced vector of hits
-template<typename T, typename T2>
-void reduce_indices_hits(const vector<T>& local_indices,
-  const vector<T2>& local_hits, vector<T>& indices, vector<T2>& hits)
-{
-  const int n_threads = num_threads();
-
-#pragma omp for ordered schedule(static, 1)
-  for (int i = 0; i < n_threads; ++i) {
-#pragma omp ordered
-    for (int j = 0; j < local_indices.size(); ++j) {
-      // Check if this material has been added to the master list and if
-      // so, accumulate the number of hits
-      auto it = std::find(indices.begin(), indices.end(), local_indices[j]);
-      if (it == indices.end()) {
-        indices.push_back(local_indices[j]);
-        hits.push_back(local_hits[j]);
-      } else {
-        hits[it - indices.begin()] += local_hits[j];
-      }
-    }
-  }
-}
-
 void free_memory_volume();
 
 } // namespace openmc

openmc/volume.py

@@ -66,13 +66,19 @@ class VolumeCalculation:
         Number of iterations over samples (for calculations with a trigger).
 
         .. versionadded:: 0.12
+    max_iterations : int
+        Limit of the maximal allowed iterations number (optional, for
+        calculations with a trigger).
+
+        .. versionadded:: 0.15.x
 
     """
     def __init__(self, domains, samples, lower_left=None, upper_right=None):
         self._atoms = {}
         self._volumes = {}
         self._threshold = None
         self._trigger_type = None
+        self._max_iterations = None
         self._iterations = None
 
         cv.check_type('domains', domains, Iterable,
@@ -187,6 +193,18 @@ def trigger_type(self, trigger_type):
                        ('variance', 'std_dev', 'rel_err'))
         self._trigger_type = trigger_type
 
+    @property
+    def max_iterations(self):
+        return self._max_iterations
+
+    @max_iterations.setter
+    def max_iterations(self, max_iterations):
+        name = 'volume calculation iterations limit'
+        cv.check_type(name, max_iterations, Integral, none_ok=True)
+        if max_iterations is not None:
+            cv.check_greater_than(name, max_iterations, 0)
+        self._max_iterations = max_iterations
+
     @property
     def iterations(self):
         return self._iterations
@@ -230,7 +248,7 @@ def atoms_dataframe(self):
 
         return pd.DataFrame.from_records(items, columns=columns)
 
-    def set_trigger(self, threshold, trigger_type):
+    def set_trigger(self, threshold, trigger_type, max_iterations=None):
         """Set a trigger on the volume calculation
 
         .. versionadded:: 0.12
@@ -241,9 +259,12 @@ def set_trigger(self, threshold, trigger_type):
             Threshold for the maximum standard deviation of volumes
         trigger_type : {'variance', 'std_dev', 'rel_err'}
             Value type used to halt volume calculation
+        max_iterations : int
+            Maximal allowed number of iterations (optional)
         """
         self.trigger_type = trigger_type
         self.threshold = threshold
+        self.max_iterations = max_iterations
 
     @classmethod
     def from_hdf5(cls, filename):
@@ -270,6 +291,7 @@ def from_hdf5(cls, filename):
 
             threshold = f.attrs.get('threshold')
             trigger_type = f.attrs.get('trigger_type')
+            max_iterations = f.attrs.get('max_iterations')
             iterations = f.attrs.get('iterations', 1)
 
             volumes = {}
@@ -304,7 +326,7 @@ def from_hdf5(cls, filename):
         vol = cls(domains, samples, lower_left, upper_right)
 
         if trigger_type is not None:
-            vol.set_trigger(threshold, trigger_type.decode())
+            vol.set_trigger(threshold, trigger_type.decode(), max_iterations)
 
         vol.iterations = iterations
         vol.volumes = volumes
@@ -355,6 +377,8 @@ def to_xml_element(self):
             trigger_elem = ET.SubElement(element, "threshold")
             trigger_elem.set("type", self.trigger_type)
             trigger_elem.set("threshold", str(self.threshold))
+            if self.max_iterations is not None:
+                trigger_elem.set("max_iterations", str(self.max_iterations))
         return element
 
     @classmethod
@@ -398,6 +422,7 @@ def from_xml_element(cls, elem):
         if trigger_elem is not None:
             trigger_type = get_text(trigger_elem, "type")
             threshold = float(get_text(trigger_elem, "threshold"))
-            vol.set_trigger(threshold, trigger_type)
+            max_iterations = Integral(get_text(trigger_elem, "max_iterations"))
+            vol.set_trigger(threshold, trigger_type, max_iterations)
 
         return vol