Angular flux tallying capabilities for Random Ray

include/openmc/mesh.h

@@ -1,3 +1,3 @@
 //! \file mesh.h
 //! \brief Mesh types used for tallies, Shannon entropy, CMFD, etc.
 
@@ -826,6 +826,83 @@
   virtual void initialize() = 0;
 };
 
+// Abstract class for meshes over directions (points on the unit sphere), 
+// as opposed to meshes over volumes.
+class AngularMesh : public Mesh {
+public:
+  AngularMesh() { n_dimension_ = 2; }
+  AngularMesh(pugi::xml_node node) : Mesh(node) { n_dimension_ = 2; }
+  AngularMesh(hid_t group) : Mesh(group) { n_dimension_ = 2; }
+
+  // Angular meshes partition direction space, not a volume.
+  // Therefore, "bin crossing" methods are not used; only the "get_bin" method.
+  void bins_crossed(Position r0, Position r1, const Direction& u,
+    vector<int>& bins, vector<double>& lengths) const override
+  {
+    fatal_error("Angular meshes do not support spatial tracklength tallies.");
+  }
+
+  void surface_bins_crossed(Position r0, Position r1, const Direction& u,
+    vector<int>& bins) const override
+  {
+    fatal_error("Angular meshes do not support surface-crossing tallies.");
+  }
+
+  int n_surface_bins() const override { return 0; }
+
+  std::pair<vector<double>, vector<double>> plot(
+    Position plot_ll, Position plot_ur) const override
+  {
+    return {{}, {}};
+  }
+
+  std::string bin_label(int bin) const override
+  {
+    return fmt::format("Element Index ({})", bin);
+  }
+
+  Position lower_left() const override { return {-1., -1., -1.}; }
+  Position upper_right() const override { return {1., 1., 1.}; }
+};
+
+class UnitSpherePointset : public AngularMesh {
+public:
+  UnitSpherePointset() = default;
+  explicit UnitSpherePointset(vector<Direction> points);
+  UnitSpherePointset(pugi::xml_node node);
+  UnitSpherePointset(hid_t group);
+
+  //! TODO: add sampling from within a spherical Voronoi cell
+  Position sample_element(int32_t bin, uint64_t* seed) const override 
+  {
+    fatal_error(
+      "Sampling over a UnitSpherePointset angular mesh is not supported");
+  }
+
+  int get_bin(Direction u) const override;
+
+  int n_bins() const override { return static_cast<int>(points_.size()); }
+
+  double volume(int bin) const override
+  {
+    fatal_error("Volume calculation over UnitSpherePointset is not supported");
+  }
+
+  void material_volumes(int nx, int ny, int nz, int max_materials,
+    int32_t* materials, double* volumes, double* bboxes) const override
+  {
+    fatal_error("material_volumes() is not supported for UnitSpherePointset")
+  }
+
+  std::string get_mesh_type() const override { return mesh_type; }
+  static const std::string mesh_type;
+
+  void to_hdf5_inner(hid_t group) const override;
+
+  vector<Position> points_;
+  vector<double> data_;
+};
+
 #ifdef OPENMC_DAGMC_ENABLED
 
 class MOABMesh : public UnstructuredMesh {

include/openmc/particle.h

@@ -106,6 +106,10 @@ class Particle : public ParticleData {
   void cross_periodic_bc(
     const Surface& surf, Position new_r, Direction new_u, int new_surface);
 
+  //! Reset angular flux tally bin if in RandomRay mode
+  //! (does nothing for MC particles)
+  virtual void direction_changed() {}
+
   //! mark a particle as lost and create a particle restart file
   //! \param message A warning message to display
   virtual void mark_as_lost(const char* message) override;

include/openmc/random_ray/flat_source_domain.h

@@ -1,3 +1,3 @@
 #ifndef OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H
 #define OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H
 
@@ -30,7 +30,7 @@
   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(
+  virtual void normalize_flux_and_volumes(
     double total_active_distance_per_iteration);
 
   int64_t add_source_to_scalar_flux();
@@ -68,10 +68,16 @@
   {
     return source_regions_.n_source_regions() * negroups_;
   }
+  int64_t n_source_angular_elements() const
+  {
+    return source_regions_.n_source_angular_elements();
+  }
   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;
+  int lookup_angular_bin(Direction u) const;
+  Direction angular_mesh_direction(int a) { return angular_bin_angles_[a]; }
 
   //----------------------------------------------------------------------------
   // Static Data members
@@ -174,6 +180,12 @@
   //----------------------------------------------------------------------------
   // Private data members
   int negroups_; // Number of energy groups in simulation
+  int nangles_;  // Number of bins for any angular flux tallies
+
+  vector<double> angular_bin_angles_; // Directions corresponding to each
+                                         // bin midpoint in the angular flux
+                                         // tallying scheme, flattened to 1D: 
+                                         // [x0, y0, z0, x1, y1, z1,...]
 
   double
     simulation_volume_; // Total physical volume of the simulation domain, as

include/openmc/random_ray/linear_source_domain.h

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

include/openmc/random_ray/random_ray.h

@@ -36,12 +36,14 @@ class RandomRay : public Particle {
     SourceRegionHandle& srh, double distance, bool is_active, Position r);
   void attenuate_flux_linear_source_void(
     SourceRegionHandle& srh, double distance, bool is_active, Position r);
+  void direction_changed() override { angular_bin_ = C_NONE; }
 
   void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain);
   uint64_t transport_history_based_single_ray();
   SourceSite sample_prng();
   SourceSite sample_halton();
   SourceSite sample_s2();
+  int angular_bin(); // accessor for current angular quadrature bin index
 
   //----------------------------------------------------------------------------
   // Static data members
@@ -67,6 +69,7 @@ class RandomRay : public Particle {
 
   int negroups_;
   int ntemperature_;
+  int angular_bin_ {C_NONE};
   FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source
                                        // data needed for ray transport
   double distance_travelled_ {0};

include/openmc/random_ray/source_region.h

@@ -1,3 +1,3 @@
 #ifndef OPENMC_RANDOM_RAY_SOURCE_REGION_H
 #define OPENMC_RANDOM_RAY_SOURCE_REGION_H
 
@@ -54,9 +54,13 @@
   int64_t filter_idx;
   int score_idx;
   int score_type;
-  TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type)
+  // Angular quadrature bin index. Defaults to C_NONE for angle-independent, 
+  // scoring via the source region's scalar flux.
+  int angle_bin {C_NONE};
+  TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type,
+    int angle_bin = C_NONE)
     : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx),
-      score_type(score_type)
+      score_type(score_type), angle_bin(angle_bin)
   {}
   TallyTask() = default;
 
@@ -65,7 +69,8 @@
   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;
+           score_idx == other.score_idx && score_type == other.score_type &&
+           angle_bin == other.angle_bin;
   }
 
   struct HashFunctor {
@@ -76,6 +81,7 @@
       hash_combine(seed, task.filter_idx);
       hash_combine(seed, task.score_idx);
       hash_combine(seed, task.score_type);
+      hash_combine(seed, task.angle_bin);
       return seed;
     }
   };
@@ -141,6 +147,7 @@
   //----------------------------------------------------------------------------
   // Public Data members
   int negroups_;
+  int nangles_ {1}; //!< Number of angular bins for angular flux binning
   bool is_numerical_fp_artifact_ {false};
   bool is_linear_ {false};
 
@@ -180,6 +187,7 @@
   float* source_;
   float* external_source_;
   double* scalar_flux_final_;
+  double* angular_flux_new_; //!< only accumulated for ext. source biasing
 
   MomentArray* source_gradients_;
   MomentArray* flux_moments_old_;
@@ -279,6 +287,15 @@
   double& scalar_flux_final(int g) { return scalar_flux_final_[g]; }
   const double scalar_flux_final(int g) const { return scalar_flux_final_[g]; }
 
+  double& angular_flux_new(int g, int a)
+  {
+    return angular_flux_new_[g * nangles_ + a];
+  }
+  const double angular_flux_new(int g, int a) const
+  {
+    return angular_flux_new_[g * nangles_ + a];
+  }
+
   float& source(int g) { return source_[g]; }
   const float source(int g) const { return source_[g]; }
 
@@ -315,7 +332,7 @@
 public:
   //----------------------------------------------------------------------------
   // Constructors
-  SourceRegion(int negroups, bool is_linear);
+  SourceRegion(int negroups, bool is_linear, int nangles = 1);
   SourceRegion() = default;
 
   //----------------------------------------------------------------------------
@@ -375,6 +392,9 @@
                                      //!< active iterations (used for plotting,
                                      //!< or computing adjoint sources)
 
+  vector<double>
+    angular_flux_new_; //!< The angular flux from the current iteration
+
   vector<MomentArray> source_gradients_; //!< The linear source gradients
   vector<MomentArray>
     flux_moments_old_; //!< The linear flux moments from the previous iteration
@@ -395,8 +415,8 @@
 public:
   //----------------------------------------------------------------------------
   // Constructors
-  SourceRegionContainer(int negroups, bool is_linear)
-    : negroups_(negroups), is_linear_(is_linear)
+  SourceRegionContainer(int negroups, bool is_linear, int nangles = 1)
+    : negroups_(negroups), is_linear_(is_linear), nangles_(nangles)
   {}
   SourceRegionContainer() = default;
 
@@ -572,6 +592,20 @@
     return scalar_flux_final_[se];
   }
 
+  double& angular_flux_new(int64_t sr, int g, int a)
+  {
+    return angular_flux_new_[index_angle(sr, g, a)];
+  }
+  const double angular_flux_new(int64_t sr, int g, int a) const
+  {
+    return angular_flux_new_[index_angle(sr, g, a)];
+  }
+  double& angular_flux_new(int64_t sea) { return angular_flux_new_[sea]; }
+  const double angular_flux_new(int64_t sea) const
+  {
+    return angular_flux_new_[sea];
+  }
+
   float& source(int64_t sr, int g) { return source_[index(sr, g)]; }
   const float source(int64_t sr, int g) const { return source_[index(sr, g)]; }
   float& source(int64_t se) { return source_[se]; }
@@ -626,8 +660,14 @@
   void flux_swap();
   int64_t n_source_regions() const { return n_source_regions_; }
   int64_t n_source_elements() const { return n_source_regions_ * negroups_; }
+  int64_t n_source_angular_elements() const
+  {
+    return n_source_regions_ * negroups_ * nangles_;
+  }
   int& negroups() { return negroups_; }
   const int negroups() const { return negroups_; }
+  int& nangles() { return nangles_; }
+  const int nangles() const { return nangles_; }
   bool& is_linear() { return is_linear_; }
   const bool is_linear() const { return is_linear_; }
   SourceRegionHandle get_source_region_handle(int64_t sr);
@@ -638,6 +678,7 @@
   // Private Data Members
   int64_t n_source_regions_ {0};
   int negroups_ {0};
+  int nangles_ {1};
   bool is_linear_ {false};
 
   // SoA storage for scalar fields (one item per source region)
@@ -671,6 +712,7 @@
   vector<double> scalar_flux_old_;
   vector<double> scalar_flux_new_;
   vector<double> scalar_flux_final_;
+  vector<double> angular_flux_new_;
   vector<float> source_;
   vector<float> external_source_;
 
@@ -691,6 +733,10 @@
 
   // Helper function for indexing
   inline int64_t index(int64_t sr, int g) const { return sr * negroups_ + g; }
+  inline int64_t index_angle(int64_t sr, int g, int a) const
+  {
+    return (sr * negroups_ + g) * nangles_ + a;
+  }
 };
 
 } // namespace openmc

include/openmc/tallies/filter.h

@@ -32,6 +32,7 @@ enum class FilterType {
   MATERIAL,
   MATERIALFROM,
   MESH,
+  MESH_ANGULAR,
   MESHBORN,
   MESH_MATERIAL,
   MESH_SURFACE,

include/openmc/tallies/filter_meshangular.h

@@ -0,0 +1,41 @@
+#ifndef OPENMC_TALLIES_FILTER_MESHSURFACE_H
+#define OPENMC_TALLIES_FILTER_MESHSURFACE_H
+
+#include "openmc/tallies/filter_mesh.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Indexes the direction of particle events to a mesh. 
+//==============================================================================
+
+class MeshAngularFilter : public MeshFilter {
+public:
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "meshangular"; }
+  FilterType type() const override { return FilterType::MESH_ANGULAR; }
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  void set_translation(const Position& translation) const override
+  {
+    fatal_error("Angular mesh filters do not permit translation.");
+  }
+
+  void set_translation(const double translation[3]) const override
+  {
+    fatal_error("Angular mesh filters do not permit translation.");
+  }
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MESHSURFACE_H

openmc/filter.py

@@ -27,7 +27,7 @@
     'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre',
     'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle',
     'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide',
-    'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
+    'weight', 'meshangular', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
 )
 
 def _mesh_current_names(mesh):
@@ -1349,6 +1349,14 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
         # Initialize a Pandas DataFrame from the mesh dictionary
         return pd.concat([df, pd.DataFrame(filter_dict)])
 
+class MeshAngularFilter(Filter):
+    """Bins tally events based on incident particle's direction, using 
+    an angular mesh.
+
+    """
+    def __init__(self):
+        pass
+
 
 class CollisionFilter(Filter):
     """Bins tally events based on the number of collisions.