Angular flux tallying capabilities for Random Ray

include/openmc/mesh.h

@@ -826,6 +826,82 @@ class UnstructuredMesh : public Mesh {
   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_;
+};
+
 #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

@@ -2,6 +2,7 @@
 #define OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H
 
 #include "openmc/constants.h"
+#include "openmc/mesh.h"
 #include "openmc/openmp_interface.h"
 #include "openmc/position.h"
 #include "openmc/random_ray/parallel_map.h"
@@ -30,12 +31,13 @@ class FlatSourceDomain {
   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();
   virtual void batch_reset();
   void convert_source_regions_to_tallies(int64_t start_sr_id);
+  void initialize_angular_quadrature();
   void reset_tally_volumes();
   void random_ray_tally();
   virtual void accumulate_iteration_flux();
@@ -72,6 +74,10 @@ class FlatSourceDomain {
   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 get_angular_bin(Direction u) const;
+  Direction angular_quadrature_direction(int a) const;
+  bool tally_angular_flux_applies(
+    int cell_idx, int material, int mesh_idx) const;
 
   //----------------------------------------------------------------------------
   // Static Data members
@@ -173,7 +179,15 @@ class FlatSourceDomain {
 
   //----------------------------------------------------------------------------
   // Private data members
-  int negroups_; // Number of energy groups in simulation
+  int negroups_;    // Number of energy groups in simulation
+  int nangles_ {1}; // Number of bins for any angular flux tallies
+  const UnitSpherePointset* angular_mesh_ {nullptr};
+
+  vector<bool> tally_is_angular_;
+  std::unordered_set<int32_t> angular_target_cells_;
+  std::unordered_set<int32_t> angular_target_materials_;
+  std::unordered_set<int32_t> angular_target_meshes_;
+  bool tally_angular_flux_everywhere_ {false};
 
   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

@@ -54,9 +54,13 @@ struct TallyTask {
   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 @@ struct TallyTask {
   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 @@ struct TallyTask {
       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 @@ class SourceRegionHandle {
   //----------------------------------------------------------------------------
   // 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};
 
@@ -159,6 +166,7 @@ class SourceRegionHandle {
   double* volume_naive_;
   int* position_recorded_;
   int* external_source_present_;
+  int* needs_angular_flux_;
   Position* position_;
   Position* centroid_;
   Position* centroid_iteration_;
@@ -180,6 +188,7 @@ class SourceRegionHandle {
   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_;
@@ -236,6 +245,9 @@ class SourceRegionHandle {
     return *external_source_present_;
   }
 
+  int& needs_angular_flux() { return *needs_angular_flux_; }
+  const int needs_angular_flux() const { return *needs_angular_flux_; }
+
   Position& position() { return *position_; }
   const Position position() const { return *position_; }
 
@@ -279,6 +291,15 @@ class SourceRegionHandle {
   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 +336,8 @@ class SourceRegion {
 public:
   //----------------------------------------------------------------------------
   // Constructors
-  SourceRegion(int negroups, bool is_linear);
+  SourceRegion(int negroups, bool is_linear, int nangles = 1,
+    bool needs_angular_flux = false);
   SourceRegion() = default;
 
   //----------------------------------------------------------------------------
@@ -337,7 +359,9 @@ class SourceRegion {
   double volume_naive_ {0.0}; //!< Volume as integrated from this iteration only
   int position_recorded_ {0}; //!< Has the position been recorded yet?
   int external_source_present_ {
-    0};               //!< Is an external source present in this region?
+    0}; //!< Is an external source present in this region?
+  int needs_angular_flux_ {
+    0};               //!< Is angular flux tallying active in this region?
   int is_small_ {0};  //!< Is it "small", receiving < 1.5 hits per iteration?
   int n_hits_ {0};    //!< Number of total hits (ray crossings)
                       // Mesh that subdivides this source region
@@ -375,6 +399,9 @@ class SourceRegion {
                                      //!< 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 +422,8 @@ class SourceRegionContainer {
 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;
 
@@ -450,6 +477,12 @@ class SourceRegionContainer {
     return external_source_present_[sr];
   }
 
+  int& needs_angular_flux(int64_t sr) { return needs_angular_flux_[sr]; }
+  const int needs_angular_flux(int64_t sr) const
+  {
+    return needs_angular_flux_[sr];
+  }
+
   Position& position(int64_t sr) { return position_[sr]; }
   const Position position(int64_t sr) const { return position_[sr]; }
 
@@ -572,6 +605,20 @@ class SourceRegionContainer {
     return scalar_flux_final_[se];
   }
 
+  double& angular_flux_new(int64_t sr, int g, int a)
+  {
+    return angular_flux_new_[angular_flux_offset_[sr] + g * nangles_ + a];
+  }
+  const double angular_flux_new(int64_t sr, int g, int a) const
+  {
+    return angular_flux_new_[angular_flux_offset_[sr] + g * nangles_ + 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 +673,11 @@ class SourceRegionContainer {
   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 angular_flux_new_.size(); }
   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 +688,7 @@ class SourceRegionContainer {
   // 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)
@@ -656,6 +707,7 @@ class SourceRegionContainer {
   vector<double> volume_naive_;
   vector<int> position_recorded_;
   vector<int> external_source_present_;
+  vector<int> needs_angular_flux_;
   vector<Position> position_;
   vector<Position> centroid_;
   vector<Position> centroid_iteration_;
@@ -671,6 +723,8 @@ class SourceRegionContainer {
   vector<double> scalar_flux_old_;
   vector<double> scalar_flux_new_;
   vector<double> scalar_flux_final_;
+  vector<double> angular_flux_new_;
+  vector<int64_t> angular_flux_offset_;
   vector<float> source_;
   vector<float> external_source_;
 

include/openmc/tallies/filter.h

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