Automated generation of FW-CADIS source biasing parameters via Random Ray

CMakeLists.txt

@@ -461,6 +461,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_material.cpp
   src/tallies/filter_materialfrom.cpp
   src/tallies/filter_mesh.cpp
+  src/tallies/filter_meshangular.cpp
   src/tallies/filter_meshborn.cpp
   src/tallies/filter_meshmaterial.cpp
   src/tallies/filter_meshsurface.cpp

include/openmc/capi.h

@@ -138,12 +138,15 @@ int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh);
 int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_mesh_filter_get_translation(int32_t index, double translation[3]);
 int openmc_mesh_filter_set_translation(int32_t index, double translation[3]);
+int openmc_mesh_get_bin(int32_t index, const double xyz[3], int32_t* bin);
 int openmc_mesh_get_id(int32_t index, int32_t* id);
 int openmc_mesh_set_id(int32_t index, int32_t id);
 int openmc_mesh_get_n_elements(int32_t index, size_t* n);
 int openmc_mesh_get_volumes(int32_t index, double* volumes);
 int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz,
   int max_mats, int32_t* materials, double* volumes, double* bboxes);
+int openmc_meshangular_filter_get_mesh(int32_t index, int32_t* index_mesh);
+int openmc_meshangular_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
 int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_new_filter(const char* type, int32_t* index);

include/openmc/constants.h

@@ -35,6 +35,7 @@ constexpr array<int, 2> VERSION_VOXEL {2, 0};
 constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
 constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
 constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
+constexpr array<int, 2> VERSION_SOURCE_BIAS {1, 0};
 constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 2};
 
 // ============================================================================

include/openmc/distribution_multi.h

@@ -8,6 +8,7 @@
 #include "openmc/distribution.h"
 #include "openmc/error.h"
 #include "openmc/position.h"
+#include "openmc/mesh.h"
 
 namespace openmc {
 
@@ -135,6 +136,46 @@ class Monodirectional : public UnitSphereDistribution {
   std::pair<Direction, double> sample(uint64_t* seed) const override;
 };
 
+//==============================================================================
+//! MeshAngular distribution
+//==============================================================================
+
+class MeshAngular : public UnitSphereDistribution {
+public:
+  MeshAngular() {};
+  explicit MeshAngular(pugi::xml_node node);
+
+  //! Sample a direction from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return (sampled Direction, sample weight)
+  std::pair<Direction, double> sample(uint64_t* seed) const override;
+
+  //! Sample the mesh for an element and direction from that element
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled element index and direction within that element
+  std::pair<int32_t, Direction> sample_mesh(uint64_t* seed) const;
+
+  //! Sample a mesh element
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled element index
+  int32_t sample_element_index(uint64_t* seed) const;
+
+  //! Ensure that mesh elements are spherical triangles
+  void check_element_types() const;
+
+  // Accessors
+  const Mesh* mesh() const { return model::meshes.at(mesh_idx_).get(); }
+  int32_t n_sources() const { return this->mesh()->n_bins(); }
+
+  double total_strength() { return this->elem_idx_dist_.integral(); }
+
+private:
+  int32_t mesh_idx_ {C_NONE};
+  DiscreteIndex elem_idx_dist_; //!< Distribution of mesh element indices
+  vector<double> weight_;       //!< Importance weights (empty if unbiased)
+  Direction v_ref_ {1.0, 0.0, 0.0}; //!< reference direction
+};
+
 using UPtrAngle = unique_ptr<UnitSphereDistribution>;
 
 } // namespace openmc

include/openmc/mesh.h

@@ -262,7 +262,7 @@ class Mesh {
   //!                          element
   //! \param[inout] materials Array storing material indices
   //! \param[inout] volumes Array storing volumes
-  void material_volumes(int nx, int ny, int nz, int max_materials,
+  virtual void material_volumes(int nx, int ny, int nz, int max_materials,
     int32_t* materials, double* volumes) const;
 
   //! Determine volume and bounding boxes of materials within each mesh element
@@ -275,7 +275,7 @@ class Mesh {
   //! \param[inout] materials Array storing material indices
   //! \param[inout] volumes Array storing volumes
   //! \param[inout] bboxes Array storing bounding boxes (n_elems, table_size, 6)
-  void material_volumes(int nx, int ny, int nz, int max_materials,
+  virtual void material_volumes(int nx, int ny, int nz, int max_materials,
     int32_t* materials, double* volumes, double* bboxes) const;
 
   //! Determine bounding box of mesh
@@ -826,6 +826,122 @@ 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 {{}, {}};
+  }
+
+  void material_volumes(int nx, int ny, int nz, int max_materials,
+    int32_t* materials, double* volumes) const override
+  {
+    fatal_error("material_volumes() is not supported for angular meshes");
+  }
+
+  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 angular meshes");
+  }
+
+  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");
+  }
+
+  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_;
+};
+
+class UnitSphereTriangularMesh : public AngularMesh {
+public:
+  UnitSphereTriangularMesh() = default;
+  explicit UnitSphereTriangularMesh(vector<double> vertices);
+  UnitSphereTriangularMesh(pugi::xml_node node);
+  UnitSphereTriangularMesh(hid_t group);
+
+  //! Direction and Position are the same type underneath--can override
+  Direction sample_element(int32_t bin, uint64_t* seed) const override;
+
+  //! Point-location on an arbitrary spherical triangular mesh is
+  //! difficult, so we leave this unimplemented.
+  //! get_bin() would only be called on an angular mesh in case of
+  //! anisotropic effects, such as angle-dependent weight windows, which a
+  //! particle encounters during transport.
+  int get_bin(Direction u) const override
+  {
+    fatal_error("get_bin() is not supported for UnitSpereTriangularMesh");
+  }
+
+  int n_bins() const override { return static_cast<int>(vertices_.size() / 9); }
+
+  //! Return surface area of a specific triangle
+  double volume(int bin) const override { return areas_[bin]; }
+
+  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<double> vertices_;
+  vector<double> areas_;
+};
+
 #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 lookup_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};