Add PointCloud spatial distribution

docker-dev-notes.txt

@@ -0,0 +1,8 @@
+# docker run -v $PWD:/root/OpenMC/openmc -it openmc_dagmc_libmesh
+
+# in docker
+# cd ~/OpenMC/build
+# cmake ../openmc
+# make -j 16 install
+# cd ../openmc
+# pip install .

docs/source/pythonapi/stats.rst

@@ -59,6 +59,7 @@ Spatial Distributions
    openmc.stats.Box
    openmc.stats.Point
    openmc.stats.MeshSpatial
+   openmc.stats.PointCloud
 
 .. autosummary::
    :toctree: generated

docs/source/usersguide/settings.rst

@@ -192,7 +192,9 @@ distributions using spherical or cylindrical coordinates, you can use
 :class:`openmc.stats.SphericalIndependent` or
 :class:`openmc.stats.CylindricalIndependent`, respectively. Meshes can also be
 used to represent spatial distributions with :class:`openmc.stats.MeshSpatial`
-by specifying a mesh and source strengths for each mesh element.
+by specifying a mesh and source strengths for each mesh element. It is also
+possible to define a "cloud" of source points each with a different relative
+probability using :class:`openmc.stats.PointCloud`.
 
 The angular distribution can be set equal to a sub-class of
 :class:`openmc.stats.UnitSphere` such as :class:`openmc.stats.Isotropic`,

include/openmc/distribution_spatial.h

@@ -136,6 +136,33 @@ class MeshSpatial : public SpatialDistribution {
                                 //!< mesh element indices
 };
 
+//==============================================================================
+//! Distribution of points
+//==============================================================================
+
+class PointCloud : public SpatialDistribution {
+public:
+  explicit PointCloud(pugi::xml_node node);
+  explicit PointCloud(
+    std::vector<Position> point_cloud, 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;
+
+  //! Sample a point
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled point index
+  int32_t sample_point_index(uint64_t* seed) const;
+
+  double total_strength() { return this->point_idx_dist_.integral(); }
+
+private:
+  std::vector<Position> point_cloud_;
+  DiscreteIndex point_idx_dist_; //!< Distribution of Position indices
+};
+
 //==============================================================================
 //! Uniform distribution of points over a box
 //==============================================================================

include/openmc/xml_interface.h

@@ -50,6 +50,9 @@ xt::xarray<T> get_node_xarray(
   return xt::adapt(v, shape);
 }
 
+std::vector<Position> get_node_position_array(
+  pugi::xml_node node, const char* name, bool lowercase = false);
+
 Position get_node_position(
   pugi::xml_node node, const char* name, bool lowercase = false);
 

openmc/stats/multivariate.py

@@ -279,6 +279,8 @@ def from_xml_element(cls, elem, meshes=None):
             return Point.from_xml_element(elem)
         elif distribution == 'mesh':
             return MeshSpatial.from_xml_element(elem, meshes)
+        elif distribution == 'cloud':
+            return PointCloud.from_xml_element(elem)
 
 
 class CartesianIndependent(Spatial):
@@ -756,6 +758,118 @@ def from_xml_element(cls, elem, meshes):
         return cls(meshes[mesh_id], strengths, volume_normalized)
 
 
+class PointCloud(Spatial):
+    """Spatial distribution from a point cloud.
+
+    This distribution specifies a discrete list of points, 
+    with corresponding relative probabilities.
+
+    .. versionadded:: 0.15.x
+
+    Parameters
+    ----------
+    positions : Iterable of 3-tuples
+        The points in space to be sampled
+    strengths : iterable of float, optional
+        An iterable of values that represents the relative probabilty of each point.
+
+    Attributes
+    ----------
+    positions: numpy.ndarray (3xN)
+        The points in space to be sampled with shape (3,N)
+    strengths : numpy.ndarray or None
+        An array of relative probabilities for each mesh point
+    """
+
+    def __init__(self, positions, strengths=None):
+        self.positions = positions
+        self.strengths = strengths
+
+    @property
+    def positions(self):
+        return self._positions
+
+    @positions.setter
+    def positions(self, given_positions):
+        if given_positions is None:
+            raise ValueError('No positions were provided')
+        cv.check_iterable_type('position list passed in', given_positions, Real, 2, 2)
+
+        if isinstance(given_positions, list):
+            cv.check_length('first position entry', given_positions[0], 3, 3)
+            self._positions = np.asarray(given_positions)
+        elif isinstance(given_positions, np.ndarray):
+            self._positions = given_positions
+        else:
+            raise ValueError('Unable to interpret source positions')
+
+    @property
+    def strengths(self):
+        return self._strengths
+
+    @strengths.setter
+    def strengths(self, given_strengths):
+        if given_strengths is not None:
+            cv.check_type('strengths array passed in', given_strengths, Iterable, Real)
+            self._strengths = np.asarray(given_strengths, dtype=float).flatten()
+        else:
+            self._strengths = None
+
+    @property
+    def num_strength_bins(self):
+        if self.strengths is None:
+            raise ValueError('Strengths are not set')
+        return self.strengths.size
+
+    def to_xml_element(self):
+        """Return XML representation of the spatial distribution
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing spatial distribution data
+
+        """
+        element = ET.Element('space')
+
+        element.set('type', 'cloud')
+
+        subelement = ET.SubElement(element, 'coords')
+        subelement.text = ' '.join(str(e) for e in self.positions.flatten())
+
+        if self.strengths is not None:
+            subelement = ET.SubElement(element, 'strengths')
+            subelement.text = ' '.join(str(e) for e in self.strengths)
+
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate spatial distribution from an XML element
+
+        Parameters
+        ----------
+        elem : lxml.etree._Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.PointCloud
+            Spatial distribution generated from XML element
+
+
+        """
+        coord_data = get_text(elem, 'coords')
+        positions = np.asarray([float(b) for b in coord_data.split()]).reshape((3,-1))
+
+        strengths = get_text(elem, 'strengths')
+        if strengths is not None:
+            strengths = [float(b) for b in strengths.split()]
+
+        return cls(positions, strengths)
+
+
+
 class Box(Spatial):
     """Uniform distribution of coordinates in a rectangular cuboid.
 

src/distribution_spatial.cpp

@@ -26,6 +26,8 @@ unique_ptr<SpatialDistribution> SpatialDistribution::create(pugi::xml_node node)
     return UPtrSpace {new SphericalIndependent(node)};
   } else if (type == "mesh") {
     return UPtrSpace {new MeshSpatial(node)};
+  } else if (type == "cloud") {
+    return UPtrSpace {new PointCloud(node)};
   } else if (type == "box") {
     return UPtrSpace {new SpatialBox(node)};
   } else if (type == "fission") {
@@ -298,6 +300,54 @@ Position MeshSpatial::sample(uint64_t* seed) const
   return this->sample_mesh(seed).second;
 }
 
+//==============================================================================
+// PointCloud implementation
+//==============================================================================
+
+PointCloud::PointCloud(pugi::xml_node node)
+{
+  std::vector<double> coords;
+
+  if (check_for_node(node, "coords")) {
+    coords = get_node_array<double>(node, "coords");
+  } else {
+    fatal_error("No coordinates were provided for the PointCloud "
+                "spatial distribution");
+  }
+  point_cloud_ = get_node_position_array(node, "coords");
+
+  std::vector<double> strengths(point_cloud_.size(), 1.0);
+
+  if (check_for_node(node, "strengths")) {
+    strengths = get_node_array<double>(node, "strengths");
+    if (strengths.size() != point_cloud_.size()) {
+      fatal_error(
+        fmt::format("Number of entries for the strengths array {} does "
+                    "not match the number of spatial points provided {}.",
+          strengths.size(), point_cloud_.size()));
+    }
+  }
+
+  point_idx_dist_.assign(strengths);
+}
+
+PointCloud::PointCloud(
+  std::vector<Position> point_cloud, gsl::span<const double> strengths)
+{
+  point_cloud_.assign(point_cloud.begin(), point_cloud.end());
+  point_idx_dist_.assign(strengths);
+}
+
+int32_t PointCloud::sample_point_index(uint64_t* seed) const
+{
+  return point_idx_dist_.sample(seed);
+}
+
+Position PointCloud::sample(uint64_t* seed) const
+{
+  return point_cloud_[this->sample_point_index(seed)];
+}
+
 //==============================================================================
 // SpatialBox implementation
 //==============================================================================

src/xml_interface.cpp

@@ -4,6 +4,7 @@
 
 #include "openmc/error.h"
 #include "openmc/string_utils.h"
+#include "openmc/vector.h"
 
 namespace openmc {
 
@@ -48,6 +49,24 @@ bool get_node_value_bool(pugi::xml_node node, const char* name)
   return false;
 }
 
+vector<Position> get_node_position_array(
+  pugi::xml_node node, const char* name, bool lowercase)
+{
+  vector<double> coords = get_node_array<double>(node, name, lowercase);
+  if (coords.size() % 3 != 0) {
+    fatal_error(fmt::format(
+      "Incorect number of coordinates in Position array ({}) for \"{}\"",
+      coords.size(), name));
+  }
+  vector<Position> positions;
+  positions.resize(coords.size() / 3);
+  auto it = coords.begin();
+  for (int i = 0; i < positions.size(); i++) {
+    positions[i] = {*it++, *it++, *it++};
+  }
+  return positions;
+}
+
 Position get_node_position(
   pugi::xml_node node, const char* name, bool lowercase)
 {

tests/unit_tests/test_source.py

@@ -49,6 +49,57 @@ def test_spherical_uniform():
 
     assert isinstance(sph_indep_function, openmc.stats.SphericalIndependent)
 
+def test_point_cloud():
+    point_list = [[1,0,0], [0,1,0], [0,0,1]]
+    positions = np.asarray(point_list)
+    strengths = [1,2,3]
+
+    space = openmc.stats.PointCloud(positions, strengths)
+    np.testing.assert_equal(space.positions, positions)
+    np.testing.assert_equal(space.strengths, strengths)
+
+    energy = openmc.stats.Discrete([1.0e6], [1.0])
+    angle = openmc.stats.Isotropic()
+
+    src = openmc.IndependentSource(space=space, angle=angle, energy=energy)
+    assert src.space == space
+    np.testing.assert_equal(src.space.positions, positions)
+    np.testing.assert_equal(src.space.strengths, strengths)
+
+    elem = src.to_xml_element()
+    src = openmc.IndependentSource.from_xml_element(elem)
+    np.testing.assert_equal(src.space.positions, positions)
+    np.testing.assert_equal(src.space.strengths, strengths)
+
+
+def test_point_cloud_strengths(run_in_tmpdir, sphere_box_model):
+    point_list = [[1,0,0], [0,1,0], [0,0,1]]
+    positions = np.asarray(point_list)
+    strengths = [3,2,1]
+
+    space = openmc.stats.PointCloud(positions, strengths)
+
+    energy = openmc.stats.Discrete([1.0e6], [1.0])
+    angle = openmc.stats.Isotropic()
+
+    src = openmc.IndependentSource(space=space, angle=angle, energy=energy)
+    model = sphere_box_model[0]
+    model.settings.run_mode = 'fixed source'
+    model.settings.source = src
+
+    try:
+        model.init_lib()
+        n_samples = 50_000
+        sites = openmc.lib.sample_external_source(n_samples)
+    finally:
+        model.finalize_lib()
+
+    sites_df = sites.to_dataframe()
+    for i, (strength, coord) in enumerate(zip(strengths, ('x', 'y', 'z'))):
+        sampled_strength = len(sites_df[sites_df[coord] == 1.0]) / n_samples
+        expected_strength = pytest.approx(strength/sum(strengths), abs=0.02)
+        assert sampled_strength == expected_strength, f'Strength incorrect for {point_list[i]}'
+
 
 def test_source_file():
     filename = 'source.h5'