Removed the ray tracing track-length estimator part of Volume Calculation from PR #3645

Author: vitmog

docs/source/usersguide/volume.rst

@@ -7,15 +7,11 @@ Stochastic Volume Calculations
 .. currentmodule:: openmc
 
 OpenMC has a capability to stochastically determine volumes of cells, materials,
-and universes by two methods. The first method works by overlaying a bounding 
-box, sampling points from within the box, and seeing what fraction of points were
-found in a desired domain. The benefit of doing this stochastically (as opposed
-to equally-spaced points), is that it is possible to give reliable error estimates
-on each stochastic quantity. The second method uses ray tracing in random
-directions from sampled via the first method points, and seeing what fraction of
-rays inside the box were found in a desired domain. This method is predominantly
-intended for low-fraction volume calculations and geometry testing. By default,
-volume calculations are provided using the first method.
+and universes. The method works by overlaying a bounding box, sampling points
+from within the box, and seeing what fraction of points were found in a desired
+domain. The benefit of doing this stochastically (as opposed to equally-spaced
+points), is that it is possible to give reliable error estimates on each
+stochastic quantity.
 
 To specify that a volume calculation be run, you first need to create an
 instance of :class:`openmc.VolumeCalculation`. The constructor takes a list of
@@ -26,16 +22,10 @@ the specified domains::
   lower_left = (-0.62, -0.62, -50.)
   upper_right = (0.62, 0.62, 50.)
   vol_calc = openmc.VolumeCalculation([fuel, clad, moderator], 1000000,
-                                      lower_left, upper_right)     
-
-To specify the same volume calculation but using ray tracing, you need to pass
-an argument::
-
-     vol_calc = openmc.VolumeCalculation([fuel, clad, moderator], 1000000,
-                                      lower_left, upper_right, 'ray')
+                                      lower_left, upper_right)
 
 For domains contained within regions that have simple definitions, OpenMC can
-sometimes automatically determine a bounding box. In this case, the bounding box
+sometimes automatically determine a bounding box. In this case, the last two
 arguments are not necessary. For example,
 
 ::

include/openmc/position.h

@@ -83,7 +83,6 @@ struct Position {
     return x * other.x + y * other.y + z * other.z;
   }
   inline double norm() const { return std::sqrt(x * x + y * y + z * z); }
-  inline double norm_sq() const { return x * x + y * y + z * z; }
   inline Position cross(Position other) const
   {
     return {y * other.z - z * other.y, z * other.x - x * other.z,

include/openmc/volume_calc.h

@@ -9,7 +9,6 @@
 #include "openmc/array.h"
 #include "openmc/cell.h"
 #include "openmc/openmp_interface.h"
-#include "openmc/plot.h"
 #include "openmc/position.h"
 #include "openmc/tallies/trigger.h"
 #include "openmc/timer.h"
@@ -78,9 +77,6 @@ class VolumeCalculation {
   //! \brief Online results of calculation specific for each thread
   struct CalcResults {
     uint64_t n_samples; //!< Number of samples
-    uint64_t n_rays;    //!< Number of traced rays
-    uint64_t n_segs;    //!< Number of traced ray segments
-    uint64_t n_errors;  //!< Number of tracing errors
     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
@@ -116,6 +112,9 @@ class VolumeCalculation {
   //! \return Vector of results for each user-specified domain
   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
   //
@@ -166,8 +165,6 @@ class VolumeCalculation {
 
   // Tally filter and map types
   enum class TallyDomain { UNIVERSE, MATERIAL, CELL };
-  // Volume estimator type
-  enum class EstMode { REJECTION, RAYTRACE };
 
   // Data members
   TallyDomain domain_type_;      //!< Type of domain (cell, material, etc.)
@@ -182,7 +179,6 @@ class VolumeCalculation {
   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
-  EstMode mode_;                //!< Volume estimator type
 
 #ifdef OPENMC_MPI
   //! \brief Creates MPI structs for passing data between threads
@@ -197,14 +193,6 @@ class VolumeCalculation {
     -999; //!< Index of zero-element tally for entire domain totals should be
           //!< out of material ID range
 
-  //! \brief Computes lengths of the two segments of a bounding box chord
-  //
-  //! \param[in] r Location inside bounding sphere splits the chord
-  //! \param[in] u Direction of chord to compute (normalized)
-  //! \return      Array of backward (first) and forward (second) chord segments
-  std::pair<double, double> get_box_chord(
-    const Position& r, const Direction& u) const;
-
   //! \brief Computes estimated mean and std.dev. for a tally
   //
   //! \param[in] n_samples    Statistic's size
@@ -215,63 +203,6 @@ class VolumeCalculation {
     const double coeff_norm, const VolTally& vol_tally) const;
 };
 
-//==============================================================================
-// Volume estimator ray class
-//==============================================================================
-
-class VolEstRay : public Ray {
-  friend class VolumeCalculation; // for access to the CalcResults struct and
-                                  // EstMode variable, not sure that this is
-                                  // a good practice in C++
-
-public:
-  //! \brief Constructs a ray of given length emitted from a given point for the
-  //! volume track length estimator into given tallies
-  //
-  //! \param[in] r   Coordinates of the ray origin
-  //! \param[in] u   Direction of the ray
-  //! \param[in] dist   Target length of the ray
-  //! \param[in] coeff_mult  Coefficient for scores multiplcation, it should
-  //! contain the reciprocal chord lenght and may contain stat. weight, etc.
-  //! \param[in] vol_calc  Volume calculation parameters
-  //! \param[in,out] results  Tallies for scoring during tracing etc.
-  VolEstRay(Position& r, Direction& u, const double dist,
-    const double coeff_mult, const VolumeCalculation& vol_calc,
-    VolumeCalculation::CalcResults& results)
-    : Ray(r, u), traversal_distance_max_(dist), coeff_mult_(coeff_mult),
-      vol_calc_(vol_calc), results_(results)
-  {}
-
-  //! \brief Track length estimator on intersection events
-  virtual void on_intersection() override;
-
-private:
-  const double traversal_distance_max_; //!< Target ray length
-  double traversal_distance_last_ =
-    0.;                               //!< One segment beyond traversal_distance
-  const double coeff_mult_;           //!< Coefficient for scoring multiplcation
-  const VolumeCalculation& vol_calc_; //!< Parameters of volume calculation
-  VolumeCalculation::CalcResults&
-    results_; //!< Tallies for scoring and other statistics collector
-
-  //! \brief Rejection estimator
-  inline void score_hit();
-
-  //! \brief Contributes to volume tallies
-  //
-  //! \param[in] mode    Type of used estimator
-  //! \param[in] score   Contribution value
-  //! \param[in] id_mat  Material ID required due to the material() and
-  //! material_lst() disagreement (see the on_intersection() body for details)
-  inline void vol_scoring(const VolumeCalculation::EstMode mode,
-    const double score, const int id_mat);
-
-  //! \brief Returns a pointer to either current point location cell or last ray
-  //! segment cell
-  inline std::unique_ptr<Cell>& current_cell(
-    VolumeCalculation::EstMode mode, int level);
-};
-
 //==============================================================================
 // Global variables
 //==============================================================================

openmc/volume.py

@@ -32,10 +32,6 @@ class VolumeCalculation:
         Upper-right coordinates of bounding box used to sample points. If this
         argument is not supplied, an attempt is made to automatically determine
         a bounding box.
-    estimator_type : {'hit', 'ray'}
-        Type of volume estimator (default: 'hit').
-
-        .. versionadded:: 0.15.x
 
     Attributes
     ----------
@@ -77,8 +73,7 @@ class VolumeCalculation:
         .. versionadded:: 0.15.x
 
     """
-    def __init__(self, domains, samples, lower_left=None, upper_right=None,
-                 estimator_type='hit'):
+    def __init__(self, domains, samples, lower_left=None, upper_right=None):
         self._atoms = {}
         self._volumes = {}
         self._threshold = None
@@ -136,8 +131,6 @@ def __init__(self, domains, samples, lower_left=None, upper_right=None,
             raise ValueError('Lower-left and upper-right bounding box '
                              'coordinates must be finite.')
 
-        self.estimator_type = estimator_type
-
     @property
     def ids(self):
         return self._ids
@@ -223,16 +216,6 @@ def iterations(self, iterations):
         cv.check_greater_than(name, iterations, 0)
         self._iterations = iterations
 
-    @property
-    def estimator_type(self):
-        return self._estimator_type
-
-    @estimator_type.setter
-    def estimator_type(self, estimator_type):
-        cv.check_value('volume estimator mode', estimator_type,
-                       ('hit', 'ray'))
-        self._estimator_type = estimator_type
-
     @property
     def domain_type(self):
         return self._domain_type
@@ -283,18 +266,6 @@ def set_trigger(self, threshold, trigger_type, max_iterations=None):
         self.threshold = threshold
         self.max_iterations = max_iterations
 
-    def set_estimator(self, estimator_type):
-        """Define type of volume estimator
-
-        .. versionadded:: 0.15
-
-        Parameters
-        ----------
-        estimator_type : {'hit', 'ray'}
-            Either rejection or ray tracing volume estimator
-        """
-        self.estimator_type = estimator_type
-
     @classmethod
     def from_hdf5(cls, filename):
         """Load stochastic volume calculation results from HDF5 file.
@@ -323,8 +294,6 @@ def from_hdf5(cls, filename):
             max_iterations = f.attrs.get('max_iterations')
             iterations = f.attrs.get('iterations', 1)
 
-            estimator_type = f.attrs.get('estimator_type')
-
             volumes = {}
             atoms = {}
             ids = []
@@ -354,8 +323,7 @@ def from_hdf5(cls, filename):
                 domains = [openmc.Universe(uid) for uid in ids]
 
         # Instantiate the class and assign results
-        vol = cls(domains, samples, lower_left, upper_right,
-                  estimator_type.decode())
+        vol = cls(domains, samples, lower_left, upper_right)
 
         if trigger_type is not None:
             vol.set_trigger(threshold, trigger_type.decode(), max_iterations)
@@ -411,8 +379,6 @@ def to_xml_element(self):
             trigger_elem.set("threshold", str(self.threshold))
             if self.max_iterations is not None:
                 trigger_elem.set("max_iterations", str(self.max_iterations))
-        et_elem = ET.SubElement(element, "estimator_type")
-        et_elem.text = str(self.estimator_type)
         return element
 
     @classmethod
@@ -437,7 +403,6 @@ def from_xml_element(cls, elem):
         samples = int(get_text(elem, "samples"))
         lower_left = tuple(get_elem_list(elem, "lower_left", float))
         upper_right = tuple(get_elem_list(elem, "upper_right", float))
-        estimator_type = get_text(elem, "estimator_type")
 
         # Instantiate some throw-away domains that are used by the constructor
         # to assign IDs
@@ -450,7 +415,7 @@ def from_xml_element(cls, elem):
             elif domain_type == 'universe':
                 domains = [openmc.Universe(uid) for uid in ids]
 
-        vol = cls(domains, samples, lower_left, upper_right, estimator_type)
+        vol = cls(domains, samples, lower_left, upper_right)
 
         # Check for trigger
         trigger_elem = elem.find("threshold")

src/plot.cpp

@@ -1695,10 +1695,6 @@ void Ray::trace()
       coord(lev).r() += boundary().distance() * coord(lev).u();
     }
     surface() = boundary().surface();
-    // Initialize last cells from current cell
-    for (int j = 0; j < n_coord(); ++j) {
-      cell_last(j) = coord(j).cell();
-    }
     n_coord_last() = n_coord();
     n_coord() = boundary().coord_level();
     if (boundary().lattice_translation()[0] != 0 ||