Mesh Source Class (#2759)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
Co-authored-by: Jonathan Shimwell <drshimwell@gmail.com>

.github/workflows/format-check.yml

@@ -1,6 +1,13 @@
 name: C++ Format Check
 
-on: pull_request
+on:
+  # allow workflow to be run manually
+  workflow_dispatch:
+
+  pull_request:
+    branches:
+      - develop
+      - master
 
 jobs:
   cpp-linter:

docs/source/io_formats/settings.rst

@@ -473,6 +473,8 @@ pseudo-random number generator.
 
   *Default*: 1
 
+.. _source_element:
+
 --------------------
 ``<source>`` Element
 --------------------
@@ -491,7 +493,8 @@ attributes/sub-elements:
     *Default*: 1.0
 
   :type:
-    Indicator of source type. One of ``independent``, ``file``, or ``compiled``.
+    Indicator of source type. One of ``independent``, ``file``, ``compiled``, or ``mesh``.
+    The type of the source will be determined by this attribute if it is present.
 
   :particle:
     The source particle type, either ``neutron`` or ``photon``.
@@ -664,6 +667,14 @@ attributes/sub-elements:
 
     *Default*: false
 
+  :mesh:
+    For mesh sources, this indicates the ID of the corresponding mesh.
+
+  :source:
+    For mesh sources, this sub-element specifies the source for an individual
+    mesh element and follows the format for :ref:`source_element`. The number of
+    ``<source>`` sub-elements should correspond to the number of mesh elements.
+
 .. _univariate:
 
 Univariate Probability Distributions

docs/source/pythonapi/base.rst

@@ -25,6 +25,7 @@ Simulation Settings
    openmc.IndependentSource
    openmc.FileSource
    openmc.CompiledSource
+   openmc.MeshSource
    openmc.SourceParticle
    openmc.VolumeCalculation
    openmc.Settings

include/openmc/distribution.h

@@ -7,6 +7,7 @@
 #include <cstddef> // for size_t
 
 #include "pugixml.hpp"
+#include <gsl/gsl-lite.hpp>
 
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
@@ -43,9 +44,9 @@ class DiscreteIndex {
 public:
   DiscreteIndex() {};
   DiscreteIndex(pugi::xml_node node);
-  DiscreteIndex(const double* p, int n);
+  DiscreteIndex(gsl::span<const double> p);
 
-  void assign(const double* p, int n);
+  void assign(gsl::span<const double> p);
 
   //! Sample a value from the distribution
   //! \param seed Pseudorandom number seed pointer
@@ -77,7 +78,7 @@ class DiscreteIndex {
 class Discrete : public Distribution {
 public:
   explicit Discrete(pugi::xml_node node);
-  Discrete(const double* x, const double* p, int n);
+  Discrete(const double* x, const double* p, size_t n);
 
   //! Sample a value from the distribution
   //! \param seed Pseudorandom number seed pointer

include/openmc/distribution_multi.h

@@ -22,6 +22,8 @@ class UnitSphereDistribution {
   explicit UnitSphereDistribution(pugi::xml_node node);
   virtual ~UnitSphereDistribution() = default;
 
+  static unique_ptr<UnitSphereDistribution> create(pugi::xml_node node);
+
   //! Sample a direction from the distribution
   //! \param seed Pseudorandom number seed pointer
   //! \return Direction sampled

include/openmc/distribution_spatial.h

@@ -19,6 +19,8 @@ class SpatialDistribution {
 
   //! Sample a position from the distribution
   virtual Position sample(uint64_t* seed) const = 0;
+
+  static unique_ptr<SpatialDistribution> create(pugi::xml_node node);
 };
 
 //==============================================================================
@@ -102,16 +104,27 @@ class SphericalIndependent : public SpatialDistribution {
 class MeshSpatial : public SpatialDistribution {
 public:
   explicit MeshSpatial(pugi::xml_node node);
+  explicit MeshSpatial(int32_t mesh_id, gsl::span<const double> strengths);
 
   //! Sample a position from the distribution
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled position
   Position sample(uint64_t* seed) const override;
 
-  const Mesh* mesh() const { return model::meshes.at(mesh_idx_).get(); }
+  //! Sample the mesh for an element and position within that element
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled element index and position within that element
+  std::pair<int32_t, Position> sample_mesh(uint64_t* seed) const;
 
+  //! For unstructured meshes, ensure that elements are all linear tetrahedra
+  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

include/openmc/mesh.h

@@ -10,6 +10,7 @@
 #include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
 
+#include "openmc/error.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/particle.h"
 #include "openmc/position.h"
@@ -81,12 +82,12 @@ class Mesh {
   //! Return a position in the local coordinates of the mesh
   virtual Position local_coords(const Position& r) const { return r; };
 
-  //! Sample a mesh volume using a certain seed
+  //! Sample a position within a mesh element
   //
-  //! \param[in] seed Seed to use for random sampling
-  //! \param[in] bin Bin value of the tet sampled
-  //! \return sampled position within tet
-  virtual Position sample(uint64_t* seed, int32_t bin) const = 0;
+  //! \param[in] bin Bin value of the mesh element sampled
+  //! \param[inout] seed Seed to use for random sampling
+  //! \return sampled position within mesh element
+  virtual Position sample_element(int32_t bin, uint64_t* seed) const = 0;
 
   //! Determine which bins were crossed by a particle
   //
@@ -183,7 +184,12 @@ class StructuredMesh : public Mesh {
     }
   };
 
-  Position sample(uint64_t* seed, int32_t bin) const override;
+  Position sample_element(int32_t bin, uint64_t* seed) const override
+  {
+    return sample_element(get_indices_from_bin(bin), seed);
+  };
+
+  virtual Position sample_element(const MeshIndex& ijk, uint64_t* seed) const;
 
   int get_bin(Position r) const override;
 
@@ -240,6 +246,30 @@ class StructuredMesh : public Mesh {
   //! \param[in] i Direction index
   virtual int get_index_in_direction(double r, int i) const = 0;
 
+  //! Get the coordinate for the mesh grid boundary in the positive direction
+  //!
+  //! \param[in] ijk Array of mesh indices
+  //! \param[in] i Direction index
+  virtual double positive_grid_boundary(const MeshIndex& ijk, int i) const
+  {
+    auto msg =
+      fmt::format("Attempting to call positive_grid_boundary on a {} mesh.",
+        get_mesh_type());
+    fatal_error(msg);
+  };
+
+  //! Get the coordinate for the mesh grid boundary in the negative direction
+  //!
+  //! \param[in] ijk Array of mesh indices
+  //! \param[in] i Direction index
+  virtual double negative_grid_boundary(const MeshIndex& ijk, int i) const
+  {
+    auto msg =
+      fmt::format("Attempting to call negative_grid_boundary on a {} mesh.",
+        get_mesh_type());
+    fatal_error(msg);
+  };
+
   //! Get the closest distance from the coordinate r to the grid surface
   //! in i direction  that bounds mesh cell ijk and that is larger than l
   //! The coordinate r does not have to be inside the mesh cell ijk. In
@@ -322,18 +352,17 @@ class RegularMesh : public StructuredMesh {
 
   void to_hdf5(hid_t group) const override;
 
-  // New methods
   //! Get the coordinate for the mesh grid boundary in the positive direction
   //!
   //! \param[in] ijk Array of mesh indices
   //! \param[in] i Direction index
-  double positive_grid_boundary(const MeshIndex& ijk, int i) const;
+  double positive_grid_boundary(const MeshIndex& ijk, int i) const override;
 
   //! Get the coordinate for the mesh grid boundary in the negative direction
   //!
   //! \param[in] ijk Array of mesh indices
   //! \param[in] i Direction index
-  double negative_grid_boundary(const MeshIndex& ijk, int i) const;
+  double negative_grid_boundary(const MeshIndex& ijk, int i) const override;
 
   //! Count number of bank sites in each mesh bin / energy bin
   //
@@ -373,18 +402,17 @@ class RectilinearMesh : public StructuredMesh {
 
   void to_hdf5(hid_t group) const override;
 
-  // New methods
   //! Get the coordinate for the mesh grid boundary in the positive direction
   //!
   //! \param[in] ijk Array of mesh indices
   //! \param[in] i Direction index
-  double positive_grid_boundary(const MeshIndex& ijk, int i) const;
+  double positive_grid_boundary(const MeshIndex& ijk, int i) const override;
 
   //! Get the coordinate for the mesh grid boundary in the negative direction
   //!
   //! \param[in] ijk Array of mesh indices
   //! \param[in] i Direction index
-  double negative_grid_boundary(const MeshIndex& ijk, int i) const;
+  double negative_grid_boundary(const MeshIndex& ijk, int i) const override;
 
   //! Return the volume for a given mesh index
   double volume(const MeshIndex& ijk) const override;
@@ -410,6 +438,8 @@ class CylindricalMesh : public PeriodicStructuredMesh {
 
   static const std::string mesh_type;
 
+  Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override;
+
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
     const Position& r0, const Direction& u, double l) const override;
 
@@ -420,10 +450,16 @@ class CylindricalMesh : public PeriodicStructuredMesh {
 
   double volume(const MeshIndex& ijk) const override;
 
-  array<vector<double>, 3> grid_;
+  // grid accessors
+  double r(int i) const { return grid_[0][i]; }
+  double phi(int i) const { return grid_[1][i]; }
+  double z(int i) const { return grid_[2][i]; }
 
   int set_grid();
 
+  // Data members
+  array<vector<double>, 3> grid_;
+
 private:
   double find_r_crossing(
     const Position& r, const Direction& u, double l, int shell) const;
@@ -466,6 +502,8 @@ class SphericalMesh : public PeriodicStructuredMesh {
 
   static const std::string mesh_type;
 
+  Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override;
+
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
     const Position& r0, const Direction& u, double l) const override;
 
@@ -474,10 +512,15 @@ class SphericalMesh : public PeriodicStructuredMesh {
 
   void to_hdf5(hid_t group) const override;
 
-  array<vector<double>, 3> grid_;
+  double r(int i) const { return grid_[0][i]; }
+  double theta(int i) const { return grid_[1][i]; }
+  double phi(int i) const { return grid_[2][i]; }
 
   int set_grid();
 
+  // Data members
+  array<vector<double>, 3> grid_;
+
 private:
   double find_r_crossing(
     const Position& r, const Direction& u, double l, int shell) const;
@@ -621,7 +664,7 @@ class MOABMesh : public UnstructuredMesh {
 
   // Overridden Methods
 
-  Position sample(uint64_t* seed, int32_t bin) const override;
+  Position sample_element(int32_t bin, uint64_t* seed) const override;
 
   void bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins, vector<double>& lengths) const override;
@@ -789,7 +832,7 @@ class LibMesh : public UnstructuredMesh {
   void bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins, vector<double>& lengths) const override;
 
-  Position sample(uint64_t* seed, int32_t bin) const override;
+  Position sample_element(int32_t bin, uint64_t* seed) const override;
 
   int get_bin(Position r) const override;
 

include/openmc/source.h

@@ -50,6 +50,8 @@ class Source {
 
   // Methods that can be overridden
   virtual double strength() const { return 1.0; }
+
+  static unique_ptr<Source> create(pugi::xml_node node);
 };
 
 //==============================================================================
@@ -101,12 +103,13 @@ class IndependentSource : public Source {
 class FileSource : public Source {
 public:
   // Constructors
-  explicit FileSource(std::string path);
-  explicit FileSource(const vector<SourceSite>& sites) : sites_ {sites} {}
+  explicit FileSource(pugi::xml_node node);
+  explicit FileSource(const std::string& path);
 
   // Methods
   SourceSite sample(uint64_t* seed) const override;
-
+  void load_sites_from_file(
+    const std::string& path); //!< Load source sites from file
 private:
   vector<SourceSite> sites_; //!< Source sites from a file
 };
@@ -118,7 +121,7 @@ class FileSource : public Source {
 class CompiledSourceWrapper : public Source {
 public:
   // Constructors, destructors
-  CompiledSourceWrapper(std::string path, std::string parameters);
+  CompiledSourceWrapper(pugi::xml_node node);
   ~CompiledSourceWrapper();
 
   // Defer implementation to custom source library
@@ -129,13 +132,44 @@ class CompiledSourceWrapper : public Source {
 
   double strength() const override { return compiled_source_->strength(); }
 
+  void setup(const std::string& path, const std::string& parameters);
+
 private:
   void* shared_library_; //!< library from dlopen
   unique_ptr<Source> compiled_source_;
 };
 
 typedef unique_ptr<Source> create_compiled_source_t(std::string parameters);
 
+//==============================================================================
+//! Mesh-based source with different distributions for each element
+//==============================================================================
+
+class MeshSource : public Source {
+public:
+  // Constructors
+  explicit MeshSource(pugi::xml_node node);
+
+  //! Sample from the external source distribution
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Sampled site
+  SourceSite sample(uint64_t* seed) const override;
+
+  // Properties
+  double strength() const override { return space_->total_strength(); }
+
+  // Accessors
+  const std::unique_ptr<Source>& source(int32_t i) const
+  {
+    return sources_.size() == 1 ? sources_[0] : sources_[i];
+  }
+
+private:
+  // Data members
+  unique_ptr<MeshSpatial> space_;           //!< Mesh spatial
+  vector<std::unique_ptr<Source>> sources_; //!< Source distributions
+};
+
 //==============================================================================
 // Functions
 //==============================================================================