Testing point detectors

CMakeLists.txt

@@ -451,6 +451,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_musurface.cpp
   src/tallies/filter_parent_nuclide.cpp
   src/tallies/filter_particle.cpp
+  src/tallies/filter_point.cpp  
   src/tallies/filter_particle_production.cpp
   src/tallies/filter_polar.cpp
   src/tallies/filter_reaction.cpp

include/openmc/constants.h

@@ -296,9 +296,9 @@ enum class MgxsType {
 
 enum class TallyResult { VALUE, SUM, SUM_SQ, SUM_THIRD, SUM_FOURTH };
 
-enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
+enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT, POINT };
 
-enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
+enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION, NEXT_EVENT };
 
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 

include/openmc/particle.h

@@ -37,6 +37,8 @@ class Particle : public ParticleData {
   //==========================================================================
   // Methods
 
+  double mass() const;
+
   double speed() const;
 
   //! create a secondary particle

include/openmc/particle_data.h

@@ -484,7 +484,7 @@ class GeometryState {
  *   Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18.
  *   https://doi.org/10.1016/j.anucene.2017.11.032.
  */
-class ParticleData : public GeometryState {
+class ParticleData : virtual public GeometryState {
 private:
   //==========================================================================
   // Data members -- see public: below for descriptions

include/openmc/physics.h

@@ -85,7 +85,7 @@ void sample_fission_neutron(
 
 //! handles all reactions with a single secondary neutron (other than fission),
 //! i.e. level scattering, (n,np), (n,na), etc.
-void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p);
+void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle& p);
 
 void sample_secondary_photons(Particle& p, int i_nuclide);
 

include/openmc/position.h

@@ -204,6 +204,17 @@ inline Position operator/(double a, Position b)
   return b /= a;
 }
 
+inline bool operator<(Position a, Position b)
+{
+  if (a[0] != b[0])
+    return (a[0] < b[0]);
+  if (a[1] != b[1])
+    return (a[1] < b[1]);
+  if (a[2] != b[2])
+    return (a[2] < b[2]);
+  return false;
+}
+
 inline Position Position::reflect(Position n) const
 {
   const double projection = n.dot(*this);

include/openmc/ray.h

@@ -1,14 +1,14 @@
 #ifndef OPENMC_RAY_H
 #define OPENMC_RAY_H
 
-#include "openmc/particle_data.h"
+#include "openmc/particle.h"
 #include "openmc/position.h"
 
 namespace openmc {
 
 // Base class that implements ray tracing logic, not necessarily through
 // defined regions of the geometry but also outside of it.
-class Ray : public GeometryState {
+class Ray : virtual public GeometryState {
 
 public:
   // Initialize from location and direction
@@ -32,6 +32,8 @@ class Ray : public GeometryState {
   // Sets the dist_ variable
   void compute_distance();
 
+  virtual void update_distance();
+
 protected:
   // Records how far the ray has traveled
   double traversal_distance_ {0.0};
@@ -46,5 +48,30 @@ class Ray : public GeometryState {
   unsigned event_counter_ {0};
 };
 
+class ParticleRay : public Ray, public Particle {
+
+public:
+  ParticleRay(
+    Position r, Direction u, ParticleType type_, double time_, double E_)
+    : Ray(r, u)
+  {
+    type() = type_;
+    time() = time_;
+    E() = E_;
+  }
+
+  void on_intersection() override;
+
+  // Sets the dist_ variable
+  void update_distance() override;
+
+  const double& traversal_distance() const { return traversal_distance_; }
+  const double& traversal_mfp() const { return traversal_mfp_; }
+
+protected:
+  // Records how much mean free paths the ray traveled
+  double traversal_mfp_ {0.0};
+};
+
 } // namespace openmc
 #endif // OPENMC_RAY_H

include/openmc/simulation.h

@@ -36,12 +36,14 @@ extern "C" double
 extern double log_spacing; //!< lethargy spacing for energy grid searches
 extern "C" int n_lost_particles;   //!< cumulative number of lost particles
 extern "C" bool need_depletion_rx; //!< need to calculate depletion rx?
-extern "C" int restart_batch;      //!< batch at which a restart job resumed
-extern "C" bool satisfy_triggers;  //!< have tally triggers been satisfied?
-extern int ssw_current_file;       //!< current surface source file
-extern "C" int total_gen;          //!< total number of generations simulated
-extern double total_weight;        //!< Total source weight in a batch
-extern int64_t work_per_rank;      //!< number of particles per MPI rank
+extern bool
+  nonvacuum_boundary_present;     //!< Does the geometry contain non-vacuum b.c.
+extern "C" int restart_batch;     //!< batch at which a restart job resumed
+extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
+extern int ssw_current_file;      //!< current surface source file
+extern "C" int total_gen;         //!< total number of generations simulated
+extern double total_weight;       //!< Total source weight in a batch
+extern int64_t work_per_rank;     //!< number of particles per MPI rank
 
 extern const RegularMesh* entropy_mesh;
 extern const RegularMesh* ufs_mesh;

include/openmc/tallies/filter.h

@@ -40,6 +40,7 @@ enum class FilterType {
   PARENT_NUCLIDE,
   PARTICLE,
   PARTICLE_PRODUCTION,
+  POINT,
   POLAR,
   REACTION,
   SPHERICAL_HARMONICS,

include/openmc/tallies/filter_point.h

@@ -0,0 +1,55 @@
+#ifndef OPENMC_TALLIES_FILTER_POINT_H
+#define OPENMC_TALLIES_FILTER_POINT_H
+
+#include "openmc/position.h"
+#include "openmc/span.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins tally by point detectors
+//==============================================================================
+
+class PointFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~PointFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "point"; }
+  FilterType type() const override { return FilterType::POINT; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const vector<std::pair<Position, double>>& detectors() const
+  {
+    return detectors_;
+  }
+
+  void set_detectors(span<std::pair<Position, double>> detectors);
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<std::pair<Position, double>> detectors_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_PARTICLE_H

include/openmc/tallies/tally.h

@@ -3,6 +3,7 @@
 
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
+#include "openmc/position.h"
 #include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/trigger.h"
@@ -11,6 +12,7 @@
 #include "openmc/tensor.h"
 #include "pugixml.hpp"
 
+#include <set>
 #include <string>
 #include <unordered_map>
 
@@ -209,11 +211,13 @@ extern vector<int> active_analog_tallies;
 extern vector<int> active_tracklength_tallies;
 extern vector<int> active_timed_tracklength_tallies;
 extern vector<int> active_collision_tallies;
+extern vector<int> active_point_tallies;
 extern vector<int> active_meshsurf_tallies;
 extern vector<int> active_surface_tallies;
 extern vector<int> active_pulse_height_tallies;
 extern vector<int32_t> pulse_height_cells;
 extern vector<double> time_grid;
+extern std::set<Position> active_point_detectors;
 
 } // namespace model
 

include/openmc/tallies/tally_scoring.h

@@ -1,9 +1,12 @@
 #ifndef OPENMC_TALLIES_TALLY_SCORING_H
 #define OPENMC_TALLIES_TALLY_SCORING_H
 
+#include "openmc/nuclide.h"
 #include "openmc/particle.h"
+#include "openmc/ray.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/thermal.h"
 
 namespace openmc {
 
@@ -81,6 +84,9 @@ void score_analog_tally_ce(Particle& p);
 //! \param p The particle being tracked
 void score_analog_tally_mg(Particle& p);
 
+void score_tracklength_tally_general(
+  Particle& p, double flux, const vector<int>& tallies);
+
 //! Score tallies using a tracklength estimate of the flux.
 //
 //! This is triggered at every event (surface crossing, lattice crossing, or
@@ -122,6 +128,41 @@ void score_surface_tally(
 //! \param tallies A vector of the indices of the tallies to score to
 void score_pulse_height_tally(Particle& p, const vector<int>& tallies);
 
+void score_point_tally(Particle& p, int i_nuclide, const Reaction& rx,
+  int i_product, Direction* v_t);
+
+void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab,
+  const NuclideMicroXS& micro);
+
+void score_point_tally(SourceSite& site, int source_index);
+
+template<typename PDF>
+void score_point_tally_impl(
+  const Position r, const ParticleType type, const double time, PDF pdffunc)
+{
+  for (auto& det : model::active_point_detectors) {
+    auto u = (det - r);
+    double total_distance = u.norm();
+    u /= total_distance;
+    double E;
+    double pdf = pdffunc(u, E);
+    auto p = ParticleRay(r, u, type, time, E);
+    p.Ray::trace();
+    double distance = p.traversal_distance();
+    if (distance < total_distance)
+      continue;
+    double mfp = p.traversal_mfp();
+    double attenuation = std::exp(-mfp);
+
+    // Save the attenuation for point filter handling
+    p.wgt_last() = p.wgt();
+    p.wgt() *= attenuation;
+
+    double flux = p.wgt_last() * pdf;
+    score_tracklength_tally_general(p, flux, model::active_point_tallies);
+  }
+}
+
 } // namespace openmc
 
 #endif // OPENMC_TALLIES_TALLY_SCORING_H

openmc/filter.py

@@ -26,7 +26,7 @@
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell',
     'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre',
     'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle',
-    'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide',
+    'particleproduction', 'point', 'cellinstance', 'collision', 'time', 'parentnuclide',
     'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
 )
 
@@ -793,6 +793,66 @@ def from_xml_element(cls, elem, **kwargs):
         return cls(bins, filter_id=filter_id)
 
 
+class PointFilter(Filter):
+    """Bins tally events based on point detectors.
+
+    Parameters
+    ----------
+    bins : sequence of tuple[tuple[Real, Real, Real], Real]
+        Point detectors positions and exclusion radii.
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : sequence of tuple[tuple[Real, Real, Real], Real]
+        Point detectors positions and exclusion radii.
+    id : int
+        Unique identifier for the filter
+    num_bins : Integral
+        The number of filter bins
+
+    """
+
+    __hash__ = Filter.__hash__
+
+    def __eq__(self, other):
+        if type(self) is not type(other):
+            return False
+        elif len(self.bins) != len(other.bins):
+            return False
+        else:
+            return all(b1==b2 for b1,b2 in zip(self.bins,other.bins))
+
+    @Filter.bins.setter
+    def bins(self, bins):
+        cv.check_type('bins', bins, Sequence, tuple)
+        for i, item in enumerate(bins):
+            cv.check_type(f'bins[{i}]', item, tuple)
+            cv.check_length(f'bins[{i}]', item, 2, 2)
+            cv.check_type(f'bins[{i}][0]', item[0], tuple, Real)
+            cv.check_length(f'bins[{i}][0]', item[0], 3, 3)
+            cv.check_type(f'bins[{i}][1]', item[1], Real)
+        self._bins = bins
+
+    def to_xml_element(self):
+        """Return XML Element representing the Filter.
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing filter data
+
+        """
+        element = ET.Element('filter')
+        element.set('id', str(self.id))
+        element.set('type', self.short_name.lower())
+
+        subelement = ET.SubElement(element, 'bins')
+        subelement.text = ' '.join(str(b) for item in self.bins 
+                                          for b in list(item[0])+[item[1]])
+        return element    
+
 class ParentNuclideFilter(ParticleFilter):
     """Bins tally events based on the parent nuclide
 

openmc/lib/filter.py

@@ -22,7 +22,7 @@
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
     'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
     'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter',
-    'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PolarFilter',
+    'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PointFilter', 'PolarFilter',
     'ReactionFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter',
     'SurfaceFilter', 'TimeFilter', 'UniverseFilter', 'WeightFilter', 'ZernikeFilter',
     'ZernikeRadialFilter', 'filters'
@@ -607,6 +607,8 @@ def bins(self):
             self._index, particle_i.ctypes.data_as(POINTER(c_int32)))
         return [ParticleType(i) for i in particle_i]
 
+class PointFilter(Filter):
+    filter_type = 'point'
 
 class ParticleProductionFilter(Filter):
     filter_type = 'particleproduction'