From c427b88731739935828b976a2a9518c2ad7c136f Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 08:45:51 -0500 Subject: [PATCH 01/40] Add PointCloud spatial distribution --- docs/source/pythonapi/stats.rst | 1 + docs/source/usersguide/settings.rst | 4 +- include/openmc/distribution_spatial.h | 28 ++++++ openmc/settings.py | 1 + openmc/stats/multivariate.py | 121 ++++++++++++++++++++++++++ src/distribution_spatial.cpp | 69 +++++++++++++++ tests/unit_tests/test_source.py | 27 ++++++ 7 files changed, 250 insertions(+), 1 deletion(-) diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index b72896c1860..c8318ba8620 100644 --- a/docs/source/pythonapi/stats.rst +++ b/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 diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 4111481a9ea..2bcdb1daadf 100644 --- a/docs/source/usersguide/settings.rst +++ b/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 relatively +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`, diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 28568c890d8..6cefaee39a9 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -136,6 +136,34 @@ class MeshSpatial : public SpatialDistribution { //!< mesh element indices }; +//============================================================================== +//! Distribution of points +//============================================================================== + +class PointCloud : public SpatialDistribution { +public: + explicit PointCloud(pugi::xml_node node); + explicit PointCloud(gsl::span x, gsl::span y, + gsl::span z, gsl::span strengths); + + //! Sample a position from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return Sampled position + Position sample(uint64_t* seed) const override; + + // Accessors + int32_t n_sources() const { return this->mesh()->n_bins(); } + + double total_strength() { return this->elem_idx_dist_.integral(); } + +private: + gsl::span x_, y_, z_; + DiscreteIndex point_idx_dist_; //!< Distribution of + //!< mesh element indices +}; + + + //============================================================================== //! Uniform distribution of points over a box //============================================================================== diff --git a/openmc/settings.py b/openmc/settings.py index 96e6368e4f3..c74e6a90e62 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1172,6 +1172,7 @@ def _create_max_order_subelement(self, root): def _create_source_subelement(self, root, mesh_memo=None): for source in self.source: root.append(source.to_xml_element()) + # PPHW todo: handle pointcloud here if isinstance(source, IndependentSource) and isinstance(source.space, MeshSpatial): path = f"./mesh[@id='{source.space.mesh.id}']" if root.find(path) is None: diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 06c65896f49..8fb92cd2ccc 100644 --- a/openmc/stats/multivariate.py +++ b/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,125 @@ 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, + each with different relative probability. + + .. 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 + ---------- + psoitions: numpy.ndarray (3xN) + The points in space to be sampled + 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 list of 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') + + for idx, axis in enumerate(('x','y','z')): + subelement = ET.SubElement(element, axis) + subelement.text = ' '.joing(str(e) for e in self.positions[idx,:]) + + 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 = {} + + for axis in enumerate(('x','y','z')): + coord_data = get_text(elem, axis) + if coord_data is not None: + coord[axis] = [float(b) for b in coord_data.split] + + positions = np.column_stack([coord[axis] for axis in ('x','y','z')]) + + 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. diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 8dbd7d88706..4ef68713c7b 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -26,6 +26,8 @@ unique_ptr 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,73 @@ Position MeshSpatial::sample(uint64_t* seed) const return this->sample_mesh(seed).second; } + +//============================================================================== +// PointCloud implementation +//============================================================================== + +PointCloud::PointCloud(pugi::xml_node node) +{ + + if (check_for_node(node, "x")) { + x_ = get_node_array(node, "x") + } + if (check_for_node(node, "y")) { + y_ = get_node_array(node, "y") + if (y_.size() != x_.size()) { + fatal_error( + fmt::format("Number of entries for the y-coordinate array {} does " + "not match the number of entries for the x-coordinate {}." + y_.size(), x_.size()) + ) + } + } + if (check_for_node(node, "z")) { + z_ = get_node_array(node, "z") + if (z_.size() != x_.size()) { + fatal_error( + fmt::format("Number of entries for the z coordinate array {} does " + "not match the number of entries for the x-coordinate {}." + z_.size(), x_.size()) + ) + } + } + + std::vector strengths(x_.size(), 1.0); + + if (check_for_node(node, "strengths")) { + strengths = get_node_array(node, "strengths") + if (strengths.size() != x_.size()) { + fatal_error( + fmt::format("Number of entries for the strengths array {} does " + "not match the number of spatial points provided {}." + strengths.size(), x_.size()) + ) + } + } + + point_idx_dist_.assign(strengths); +} + +PointCloud::PointCloud(gsl::span x, gsl::span y, + gsl::span z, gsl::span strengths) + : x_(x), y_(y), z_(z) +{ + 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 +{ + int32_t point_idx = this->sample_point_index(seed); + return {x_[point_idx], y_[point_idx], z_[point_idx]}; +} + + //============================================================================== // SpatialBox implementation //============================================================================== diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 9a19f6f24dd..0f3f4933e0a 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -49,6 +49,33 @@ 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) + + space = openmc.stats.PointCloud(point_list, 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) + assert src.space == space + np.testing.assert_equal(src.space.positions, positions) + np.testing.assert_equal(src.space.strengths, strengths) + def test_source_file(): filename = 'source.h5' From 93afabfb8f298f964cd860d3c65030f59a161dd3 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 08:55:38 -0500 Subject: [PATCH 02/40] clang-format --- include/openmc/distribution_spatial.h | 8 ++--- src/distribution_spatial.cpp | 47 ++++++++++++--------------- 2 files changed, 24 insertions(+), 31 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 6cefaee39a9..4d2f97edfb2 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -143,8 +143,8 @@ class MeshSpatial : public SpatialDistribution { class PointCloud : public SpatialDistribution { public: explicit PointCloud(pugi::xml_node node); - explicit PointCloud(gsl::span x, gsl::span y, - gsl::span z, gsl::span strengths); + explicit PointCloud(gsl::span x, gsl::span y, + gsl::span z, gsl::span strengths); //! Sample a position from the distribution //! \param seed Pseudorandom number seed pointer @@ -159,11 +159,9 @@ class PointCloud : public SpatialDistribution { private: gsl::span x_, y_, z_; DiscreteIndex point_idx_dist_; //!< Distribution of - //!< mesh element indices + //!< mesh element indices }; - - //============================================================================== //! Uniform distribution of points over a box //============================================================================== diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 4ef68713c7b..ced021ed442 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -300,7 +300,6 @@ Position MeshSpatial::sample(uint64_t* seed) const return this->sample_mesh(seed).second; } - //============================================================================== // PointCloud implementation //============================================================================== @@ -309,48 +308,45 @@ PointCloud::PointCloud(pugi::xml_node node) { if (check_for_node(node, "x")) { - x_ = get_node_array(node, "x") + x_ = get_node_array(node, "x"); } if (check_for_node(node, "y")) { - y_ = get_node_array(node, "y") + y_ = get_node_array(node, "y"); if (y_.size() != x_.size()) { - fatal_error( - fmt::format("Number of entries for the y-coordinate array {} does " - "not match the number of entries for the x-coordinate {}." - y_.size(), x_.size()) - ) + fatal_error(fmt::format( + "Number of entries for the y-coordinate array {} does " + "not match the number of entries for the x-coordinate {}." y_.size(), + x_.size())); } } if (check_for_node(node, "z")) { - z_ = get_node_array(node, "z") + z_ = get_node_array(node, "z"); if (z_.size() != x_.size()) { - fatal_error( - fmt::format("Number of entries for the z coordinate array {} does " - "not match the number of entries for the x-coordinate {}." - z_.size(), x_.size()) - ) + fatal_error(fmt::format( + "Number of entries for the z coordinate array {} does " + "not match the number of entries for the x-coordinate {}." z_.size(), + x_.size())); } } std::vector strengths(x_.size(), 1.0); if (check_for_node(node, "strengths")) { - strengths = get_node_array(node, "strengths") - if (strengths.size() != x_.size()) { - fatal_error( - fmt::format("Number of entries for the strengths array {} does " - "not match the number of spatial points provided {}." - strengths.size(), x_.size()) - ) + strengths = get_node_array(node, "strengths"); + if (strengths.size() != x_.size()) { + fatal_error(fmt::format( + "Number of entries for the strengths array {} does " + "not match the number of spatial points provided {}." strengths.size(), + x_.size())); } - } + } point_idx_dist_.assign(strengths); } -PointCloud::PointCloud(gsl::span x, gsl::span y, - gsl::span z, gsl::span strengths) - : x_(x), y_(y), z_(z) +PointCloud::PointCloud(gsl::span x, gsl::span y, + gsl::span z, gsl::span strengths) + : x_(x), y_(y), z_(z) { point_idx_dist_.assign(strengths); } @@ -366,7 +362,6 @@ Position PointCloud::sample(uint64_t* seed) const return {x_[point_idx], y_[point_idx], z_[point_idx]}; } - //============================================================================== // SpatialBox implementation //============================================================================== From 954191daca753be1d4f99e104a264987fd03e850 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 08:58:35 -0500 Subject: [PATCH 03/40] PointCloud should be handled like a standard IndependentSource in settings --- openmc/settings.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/settings.py b/openmc/settings.py index c74e6a90e62..96e6368e4f3 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1172,7 +1172,6 @@ def _create_max_order_subelement(self, root): def _create_source_subelement(self, root, mesh_memo=None): for source in self.source: root.append(source.to_xml_element()) - # PPHW todo: handle pointcloud here if isinstance(source, IndependentSource) and isinstance(source.space, MeshSpatial): path = f"./mesh[@id='{source.space.mesh.id}']" if root.find(path) is None: From 2077a35e4d7a22283ab6149bdce03d404f20312e Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 09:16:02 -0500 Subject: [PATCH 04/40] clean up reference to mesh info in PointCloud --- include/openmc/distribution_spatial.h | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 4d2f97edfb2..32e57709ed6 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -151,10 +151,7 @@ class PointCloud : public SpatialDistribution { //! \return Sampled position Position sample(uint64_t* seed) const override; - // Accessors - int32_t n_sources() const { return this->mesh()->n_bins(); } - - double total_strength() { return this->elem_idx_dist_.integral(); } + double total_strength() { return this->point_idx_dist_.integral(); } private: gsl::span x_, y_, z_; From 74532a747be3c1a7c3579578f7a4eeb97483c5a9 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 09:19:31 -0500 Subject: [PATCH 05/40] missing commas --- src/distribution_spatial.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index ced021ed442..0211ccd466d 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -313,19 +313,19 @@ PointCloud::PointCloud(pugi::xml_node node) if (check_for_node(node, "y")) { y_ = get_node_array(node, "y"); if (y_.size() != x_.size()) { - fatal_error(fmt::format( - "Number of entries for the y-coordinate array {} does " - "not match the number of entries for the x-coordinate {}." y_.size(), - x_.size())); + fatal_error( + fmt::format("Number of entries for the y-coordinate array {} does " + "not match the number of entries for the x-coordinate {}.", + y_.size(), x_.size())); } } if (check_for_node(node, "z")) { z_ = get_node_array(node, "z"); if (z_.size() != x_.size()) { - fatal_error(fmt::format( - "Number of entries for the z coordinate array {} does " - "not match the number of entries for the x-coordinate {}." z_.size(), - x_.size())); + fatal_error( + fmt::format("Number of entries for the z coordinate array {} does " + "not match the number of entries for the x-coordinate {}.", + z_.size(), x_.size())); } } @@ -334,10 +334,10 @@ PointCloud::PointCloud(pugi::xml_node node) if (check_for_node(node, "strengths")) { strengths = get_node_array(node, "strengths"); if (strengths.size() != x_.size()) { - fatal_error(fmt::format( - "Number of entries for the strengths array {} does " - "not match the number of spatial points provided {}." strengths.size(), - x_.size())); + fatal_error( + fmt::format("Number of entries for the strengths array {} does " + "not match the number of spatial points provided {}.", + strengths.size(), x_.size())); } } From 0c78f394fe3622d7eb037531b3877828fe4a3e3c Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 09:23:43 -0500 Subject: [PATCH 06/40] add missing function declaration --- include/openmc/distribution_spatial.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 32e57709ed6..c3e2f4a2694 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -151,6 +151,11 @@ class PointCloud : public SpatialDistribution { //! \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: From a50d7af5905cc5eec5623c7ad8266bb5afd39dde Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 10:16:43 -0500 Subject: [PATCH 07/40] dumb typo --- openmc/stats/multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 8fb92cd2ccc..f8bd63b79ba 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -836,7 +836,7 @@ def to_xml_element(self): for idx, axis in enumerate(('x','y','z')): subelement = ET.SubElement(element, axis) - subelement.text = ' '.joing(str(e) for e in self.positions[idx,:]) + subelement.text = ' '.join(str(e) for e in self.positions[idx,:]) if self.strengths is not None: subelement = ET.SubElement(element, 'strengths') From 733337d49f061444db254c8aa68ea0e2b1d66a45 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 11:08:40 -0500 Subject: [PATCH 08/40] remove stale enumeration --- openmc/stats/multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index f8bd63b79ba..395ddec1e36 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -862,7 +862,7 @@ def from_xml_element(cls, elem): """ coord = {} - for axis in enumerate(('x','y','z')): + for axis in ('x','y','z'): coord_data = get_text(elem, axis) if coord_data is not None: coord[axis] = [float(b) for b in coord_data.split] From 3138b92355546eed5f5c6bc8df0821abc3c0d212 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 11:59:36 -0500 Subject: [PATCH 09/40] another dumb typo --- openmc/stats/multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 395ddec1e36..600c05fcf9e 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -865,7 +865,7 @@ def from_xml_element(cls, elem): for axis in ('x','y','z'): coord_data = get_text(elem, axis) if coord_data is not None: - coord[axis] = [float(b) for b in coord_data.split] + coord[axis] = [float(b) for b in coord_data.split()] positions = np.column_stack([coord[axis] for axis in ('x','y','z')]) From 27ce9fc986cf381e9f0ac5125424abd111dfe905 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 5 Oct 2024 12:56:56 -0500 Subject: [PATCH 10/40] testing object equivalence is invalid --- tests/unit_tests/test_source.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 0f3f4933e0a..31809d3e0f5 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -72,7 +72,6 @@ def test_point_cloud(): elem = src.to_xml_element() src = openmc.IndependentSource.from_xml_element(elem) - assert src.space == space np.testing.assert_equal(src.space.positions, positions) np.testing.assert_equal(src.space.strengths, strengths) From bbb2d4be5849ff8a6e27425c8a862609d939bc37 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Fri, 11 Oct 2024 17:16:11 -0500 Subject: [PATCH 11/40] add input stream operator for Position --- include/openmc/position.h | 2 ++ src/position.cpp | 6 ++++++ 2 files changed, 8 insertions(+) diff --git a/include/openmc/position.h b/include/openmc/position.h index 11ea3764792..ffd48ec935b 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -207,6 +207,8 @@ inline bool operator!=(Position a, Position b) } std::ostream& operator<<(std::ostream& os, Position a); +// input operator assumes space-delimited list of 3 doubles +std::istream& operator>>(std::istream& is, Position a); //============================================================================== //! Type representing a vector direction in Cartesian coordinates diff --git a/src/position.cpp b/src/position.cpp index 5b3613b2897..fd98dcf96f5 100644 --- a/src/position.cpp +++ b/src/position.cpp @@ -88,4 +88,10 @@ std::ostream& operator<<(std::ostream& os, Position r) return os; } +std::istream& operator>>(std::istream& is, Position r) +{ + is >> r.x >> r.y >> r.z; + return is; +} + } // namespace openmc From 3f6e5e312c02cd2fdf35ebc094697609785a4bc6 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Fri, 11 Oct 2024 17:17:09 -0500 Subject: [PATCH 12/40] update XML of list of positions to be flattened list of coords --- include/openmc/distribution_spatial.h | 5 ++- openmc/stats/multivariate.py | 19 ++++------- src/distribution_spatial.cpp | 45 ++++++++++----------------- 3 files changed, 25 insertions(+), 44 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index c3e2f4a2694..202091f4e51 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -159,9 +159,8 @@ class PointCloud : public SpatialDistribution { double total_strength() { return this->point_idx_dist_.integral(); } private: - gsl::span x_, y_, z_; - DiscreteIndex point_idx_dist_; //!< Distribution of - //!< mesh element indices + std::vector point_cloud_; + DiscreteIndex point_idx_dist_; //!< Distribution of Position indices }; //============================================================================== diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 600c05fcf9e..3f7acfbf8d7 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -775,8 +775,8 @@ class PointCloud(Spatial): Attributes ---------- - psoitions: numpy.ndarray (3xN) - The points in space to be sampled + 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 """ @@ -834,9 +834,8 @@ def to_xml_element(self): element.set('type', 'cloud') - for idx, axis in enumerate(('x','y','z')): - subelement = ET.SubElement(element, axis) - subelement.text = ' '.join(str(e) for e in self.positions[idx,:]) + 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') @@ -860,14 +859,8 @@ def from_xml_element(cls, elem): """ - coord = {} - - for axis in ('x','y','z'): - coord_data = get_text(elem, axis) - if coord_data is not None: - coord[axis] = [float(b) for b in coord_data.split()] - - positions = np.column_stack([coord[axis] for axis in ('x','y','z')]) + 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: diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 0211ccd466d..43e4cc86927 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -307,46 +307,36 @@ Position MeshSpatial::sample(uint64_t* seed) const PointCloud::PointCloud(pugi::xml_node node) { - if (check_for_node(node, "x")) { - x_ = get_node_array(node, "x"); - } - if (check_for_node(node, "y")) { - y_ = get_node_array(node, "y"); - if (y_.size() != x_.size()) { - fatal_error( - fmt::format("Number of entries for the y-coordinate array {} does " - "not match the number of entries for the x-coordinate {}.", - y_.size(), x_.size())); - } - } - if (check_for_node(node, "z")) { - z_ = get_node_array(node, "z"); - if (z_.size() != x_.size()) { - fatal_error( - fmt::format("Number of entries for the z coordinate array {} does " - "not match the number of entries for the x-coordinate {}.", - z_.size(), x_.size())); - } + std::vector coords; + + if (check_for_node(node, "coords")) { + point_cloud_ = get_node_array(node, "coords"); + } else { + fatal_error("No coordinates were provided for the PointCloud " + "spatial distribution"); + ) } - std::vector strengths(x_.size(), 1.0); + int32_t num_coords = point_cloud_.size(); + + std::vector strengths(num_coords, 1.0); if (check_for_node(node, "strengths")) { strengths = get_node_array(node, "strengths"); - if (strengths.size() != x_.size()) { + if (strengths.size() != num_coords) { fatal_error( fmt::format("Number of entries for the strengths array {} does " "not match the number of spatial points provided {}.", - strengths.size(), x_.size())); + strengths.size(), num_coords)); } } point_idx_dist_.assign(strengths); } -PointCloud::PointCloud(gsl::span x, gsl::span y, - gsl::span z, gsl::span strengths) - : x_(x), y_(y), z_(z) +PointCloud::PointCloud( + gsl::span point_cloud, gsl::span strengths) + : point_cloud_(point_cloud) { point_idx_dist_.assign(strengths); } @@ -358,8 +348,7 @@ int32_t PointCloud::sample_point_index(uint64_t* seed) const Position PointCloud::sample(uint64_t* seed) const { - int32_t point_idx = this->sample_point_index(seed); - return {x_[point_idx], y_[point_idx], z_[point_idx]}; + return point_cloud_[this->sample_point_index(seed)]; } //============================================================================== From 3a5a9e6362352b362370e6ab26bc75d82c00811c Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Fri, 11 Oct 2024 17:25:47 -0500 Subject: [PATCH 13/40] fix constructor arguments --- include/openmc/distribution_spatial.h | 3 +-- src/distribution_spatial.cpp | 5 ++--- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 202091f4e51..90fa20db80b 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -143,8 +143,7 @@ class MeshSpatial : public SpatialDistribution { class PointCloud : public SpatialDistribution { public: explicit PointCloud(pugi::xml_node node); - explicit PointCloud(gsl::span x, gsl::span y, - gsl::span z, gsl::span strengths); + explicit PointCloud(std::vector point_cloud, gsl::span strengths); //! Sample a position from the distribution //! \param seed Pseudorandom number seed pointer diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 43e4cc86927..be7888fe51f 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -314,7 +314,6 @@ PointCloud::PointCloud(pugi::xml_node node) } else { fatal_error("No coordinates were provided for the PointCloud " "spatial distribution"); - ) } int32_t num_coords = point_cloud_.size(); @@ -335,9 +334,9 @@ PointCloud::PointCloud(pugi::xml_node node) } PointCloud::PointCloud( - gsl::span point_cloud, gsl::span strengths) - : point_cloud_(point_cloud) + std::vector point_cloud, gsl::span strengths) { + point_cloud_.assign(point_cloud.begin(), point_cloud.end()); point_idx_dist_.assign(strengths); } From c0166667afbaf4356a170400238f99d74bfcb5a8 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Fri, 11 Oct 2024 17:27:41 -0500 Subject: [PATCH 14/40] linting --- include/openmc/distribution_spatial.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 90fa20db80b..367a4129205 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -143,7 +143,8 @@ class MeshSpatial : public SpatialDistribution { class PointCloud : public SpatialDistribution { public: explicit PointCloud(pugi::xml_node node); - explicit PointCloud(std::vector point_cloud, gsl::span strengths); + explicit PointCloud( + std::vector point_cloud, gsl::span strengths); //! Sample a position from the distribution //! \param seed Pseudorandom number seed pointer From f8933d74839e3cae78501f2f4f30e22d9ed55b1a Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 12 Oct 2024 12:26:50 -0500 Subject: [PATCH 15/40] back away from istream operator for Position --- docker-dev-notes.txt | 8 ++++++++ include/openmc/position.h | 2 -- src/distribution_spatial.cpp | 21 ++++++++++++++++----- src/position.cpp | 6 ------ 4 files changed, 24 insertions(+), 13 deletions(-) create mode 100644 docker-dev-notes.txt diff --git a/docker-dev-notes.txt b/docker-dev-notes.txt new file mode 100644 index 00000000000..1b600bb6f9d --- /dev/null +++ b/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 . diff --git a/include/openmc/position.h b/include/openmc/position.h index ffd48ec935b..11ea3764792 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -207,8 +207,6 @@ inline bool operator!=(Position a, Position b) } std::ostream& operator<<(std::ostream& os, Position a); -// input operator assumes space-delimited list of 3 doubles -std::istream& operator>>(std::istream& is, Position a); //============================================================================== //! Type representing a vector direction in Cartesian coordinates diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index be7888fe51f..2517fbb87ed 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -310,23 +310,34 @@ PointCloud::PointCloud(pugi::xml_node node) std::vector coords; if (check_for_node(node, "coords")) { - point_cloud_ = get_node_array(node, "coords"); + coords = get_node_array(node, "coords"); } else { fatal_error("No coordinates were provided for the PointCloud " "spatial distribution"); } - int32_t num_coords = point_cloud_.size(); + int32_t extra_coords = coords.size() % 3; - std::vector strengths(num_coords, 1.0); + if (extra_coords != 0) { + fatal_error( + fmt::format("Number of entries for coordinates must be multiple " + "of three; found {} extra entires.", + extra_coords)); + } + + int32_t num_positions = coords.size() / 3; + + point_cloud_.resize(num_positions); + + std::vector strengths(num_positions, 1.0); if (check_for_node(node, "strengths")) { strengths = get_node_array(node, "strengths"); - if (strengths.size() != num_coords) { + if (strengths.size() != num_positions) { fatal_error( fmt::format("Number of entries for the strengths array {} does " "not match the number of spatial points provided {}.", - strengths.size(), num_coords)); + strengths.size(), num_positions)); } } diff --git a/src/position.cpp b/src/position.cpp index fd98dcf96f5..5b3613b2897 100644 --- a/src/position.cpp +++ b/src/position.cpp @@ -88,10 +88,4 @@ std::ostream& operator<<(std::ostream& os, Position r) return os; } -std::istream& operator>>(std::istream& is, Position r) -{ - is >> r.x >> r.y >> r.z; - return is; -} - } // namespace openmc From 85663d3e62c1e328fd560d21cf7347e3872f953a Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sat, 12 Oct 2024 12:42:30 -0500 Subject: [PATCH 16/40] I forgot to store the position data --- src/distribution_spatial.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 2517fbb87ed..503a26f9f39 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -328,6 +328,10 @@ PointCloud::PointCloud(pugi::xml_node node) int32_t num_positions = coords.size() / 3; point_cloud_.resize(num_positions); + for (int32_t pos_idx = 0; pos_idx < num_positions; pos_idx++) { + point_cloud_[pos_idx] = { + coords[3 * pos_idx], coords[3 * pos_idx + 1], coords[3 * pos_idx + 2]}; + } std::vector strengths(num_positions, 1.0); From d8a0fd7c77dfd647ed864deb12a77fadc87175f2 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 15 Oct 2024 16:19:56 -0500 Subject: [PATCH 17/40] Create a new XML function for reading a flat set of Positions from an XML node --- include/openmc/xml_interface.h | 3 +++ src/distribution_spatial.cpp | 24 ++++-------------------- src/xml_interface.cpp | 18 ++++++++++++++++++ 3 files changed, 25 insertions(+), 20 deletions(-) diff --git a/include/openmc/xml_interface.h b/include/openmc/xml_interface.h index bd6554c134a..f49613ecde1 100644 --- a/include/openmc/xml_interface.h +++ b/include/openmc/xml_interface.h @@ -50,6 +50,9 @@ xt::xarray get_node_xarray( return xt::adapt(v, shape); } +std::vector 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); diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 503a26f9f39..41cbb92efff 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -315,33 +315,17 @@ PointCloud::PointCloud(pugi::xml_node node) fatal_error("No coordinates were provided for the PointCloud " "spatial distribution"); } + point_cloud_ = get_node_position_array(node, "coords"); - int32_t extra_coords = coords.size() % 3; - - if (extra_coords != 0) { - fatal_error( - fmt::format("Number of entries for coordinates must be multiple " - "of three; found {} extra entires.", - extra_coords)); - } - - int32_t num_positions = coords.size() / 3; - - point_cloud_.resize(num_positions); - for (int32_t pos_idx = 0; pos_idx < num_positions; pos_idx++) { - point_cloud_[pos_idx] = { - coords[3 * pos_idx], coords[3 * pos_idx + 1], coords[3 * pos_idx + 2]}; - } - - std::vector strengths(num_positions, 1.0); + std::vector strengths(point_cloud_.size(), 1.0); if (check_for_node(node, "strengths")) { strengths = get_node_array(node, "strengths"); - if (strengths.size() != num_positions) { + 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(), num_positions)); + strengths.size(), point_cloud_.size())); } } diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index 6ce465bafb9..af3bb5b156d 100644 --- a/src/xml_interface.cpp +++ b/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,23 @@ bool get_node_value_bool(pugi::xml_node node, const char* name) return false; } +vector +get_node_position_array(pugi::xml_node node, const char* name, bool lowercase) +{ + vector coords = get_node_array(node, name, lowercase); + if (coords.size() % 3 != 0) { + fatal_error(fmt::format( + "Incorect number of coordinates in Position array ({}) for \"{}\"", coords.size(), name)); + } + vector 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) { From 2c7eb6cfefb0a4a35abbc945f6f40138f438858d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 15 Oct 2024 16:21:00 -0500 Subject: [PATCH 18/40] Adding check for source strengths based using sampled source --- tests/unit_tests/test_source.py | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 31809d3e0f5..89353a82379 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -58,10 +58,6 @@ def test_point_cloud(): np.testing.assert_equal(space.positions, positions) np.testing.assert_equal(space.strengths, strengths) - space = openmc.stats.PointCloud(point_list, 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() @@ -75,6 +71,34 @@ def test_point_cloud(): np.testing.assert_equal(src.space.positions, positions) np.testing.assert_equal(src.space.strengths, strengths) +def test_point_cloud_strengths(run_in_tmpdir): + 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) + assembly = openmc.examples.pwr_assembly() + assembly.settings.run_mode = 'fixed source' + assembly.settings.source = src + + try: + assembly.init_lib() + n_samples = 10_000 + sites = openmc.lib.sample_external_source(n_samples) + finally: + assembly.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.01) + assert sampled_strength == expected_strength, f'Strength incorrect for {point_list[i]}' + def test_source_file(): filename = 'source.h5' From 7f2b39a8b7769ff4a568582b569b90bc68c9a9c3 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 15 Oct 2024 16:23:26 -0500 Subject: [PATCH 19/40] C++ formatting --- src/xml_interface.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index af3bb5b156d..5432b21a95a 100644 --- a/src/xml_interface.cpp +++ b/src/xml_interface.cpp @@ -49,13 +49,14 @@ bool get_node_value_bool(pugi::xml_node node, const char* name) return false; } -vector -get_node_position_array(pugi::xml_node node, const char* name, bool lowercase) +vector get_node_position_array( + pugi::xml_node node, const char* name, bool lowercase) { vector coords = get_node_array(node, name, lowercase); if (coords.size() % 3 != 0) { fatal_error(fmt::format( - "Incorect number of coordinates in Position array ({}) for \"{}\"", coords.size(), name)); + "Incorect number of coordinates in Position array ({}) for \"{}\"", + coords.size(), name)); } vector positions; positions.resize(coords.size() / 3); From 3641f99da45f3a54f82709bd7bcdfa342137d411 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 15 Oct 2024 16:24:51 -0500 Subject: [PATCH 20/40] PEP8 --- tests/unit_tests/test_source.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 89353a82379..ee66cc11a69 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -71,6 +71,7 @@ def test_point_cloud(): np.testing.assert_equal(src.space.positions, positions) np.testing.assert_equal(src.space.strengths, strengths) + def test_point_cloud_strengths(run_in_tmpdir): point_list = [[1,0,0], [0,1,0], [0,0,1]] positions = np.asarray(point_list) From 5fafc7da9e6e7569c9330d036a28cf0594f582f5 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 16 Oct 2024 10:09:52 -0500 Subject: [PATCH 21/40] Use a simpler model to cut down on initialization time. Increase tolerance and number of particles. --- tests/unit_tests/test_source.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index ee66cc11a69..5de26020a5a 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -72,7 +72,7 @@ def test_point_cloud(): np.testing.assert_equal(src.space.strengths, strengths) -def test_point_cloud_strengths(run_in_tmpdir): +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] @@ -83,21 +83,21 @@ def test_point_cloud_strengths(run_in_tmpdir): angle = openmc.stats.Isotropic() src = openmc.IndependentSource(space=space, angle=angle, energy=energy) - assembly = openmc.examples.pwr_assembly() - assembly.settings.run_mode = 'fixed source' - assembly.settings.source = src + model = sphere_box_model[0] + model.settings.run_mode = 'fixed source' + model.settings.source = src try: - assembly.init_lib() - n_samples = 10_000 + model.init_lib() + n_samples = 50_000 sites = openmc.lib.sample_external_source(n_samples) finally: - assembly.finalize_lib() + 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.01) + expected_strength = pytest.approx(strength/sum(strengths), abs=0.02) assert sampled_strength == expected_strength, f'Strength incorrect for {point_list[i]}' From ed3b1d3d198d247589739fe7b74900925a94047d Mon Sep 17 00:00:00 2001 From: Paul Wilson Date: Thu, 17 Oct 2024 12:11:48 -0500 Subject: [PATCH 22/40] grammar fix Co-authored-by: Patrick Shriwise --- docs/source/usersguide/settings.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 2bcdb1daadf..3fbfa950fbb 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -193,7 +193,7 @@ distributions using spherical or cylindrical coordinates, you can use :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. It is also -possible to define a "cloud" of source points each with a different relatively +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 From f4b590101b26af284701f98392a4b59cc1a7aae7 Mon Sep 17 00:00:00 2001 From: Paul Wilson Date: Thu, 17 Oct 2024 12:12:31 -0500 Subject: [PATCH 23/40] improve comment grammar Co-authored-by: Patrick Shriwise --- openmc/stats/multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 3f7acfbf8d7..ff91b9e5826 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -762,7 +762,7 @@ class PointCloud(Spatial): """Spatial distribution from a point cloud. This distribution specifies a discrete list of points, - each with different relative probability. + with corresponding relative probabilities. .. versionadded:: 0.15.x From 3dd55a74904e82120b50780b4f0265bc3dd48693 Mon Sep 17 00:00:00 2001 From: Paul Wilson Date: Thu, 17 Oct 2024 12:12:50 -0500 Subject: [PATCH 24/40] capitalization standard Co-authored-by: Patrick Shriwise --- openmc/stats/multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index ff91b9e5826..24ae3e3e2dc 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -768,7 +768,7 @@ class PointCloud(Spatial): Parameters ---------- - positions : iterable of 3-tuples + 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. From d9ca528a93fd1074e76bbf5d2095c85c26588123 Mon Sep 17 00:00:00 2001 From: Paul Wilson Date: Thu, 17 Oct 2024 12:13:23 -0500 Subject: [PATCH 25/40] remove extra blank line Co-authored-by: Patrick Shriwise --- src/distribution_spatial.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 41cbb92efff..e02f74aff6d 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -306,7 +306,6 @@ Position MeshSpatial::sample(uint64_t* seed) const PointCloud::PointCloud(pugi::xml_node node) { - std::vector coords; if (check_for_node(node, "coords")) { From ea462a947dc99ce6d9b739bf2c99f9834e867878 Mon Sep 17 00:00:00 2001 From: Paul Wilson Date: Thu, 17 Oct 2024 12:13:49 -0500 Subject: [PATCH 26/40] improve error message Co-authored-by: Patrick Shriwise --- openmc/stats/multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 24ae3e3e2dc..3ef38841f8c 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -801,7 +801,7 @@ def positions(self, given_positions): elif isinstance(given_positions, np.ndarray): self._positions = given_positions else: - raise ValueError('Unable to interpret list of positions') + raise ValueError('Unable to interpret source positions') @property def strengths(self): From 34879661ef97bb25967a8652443a022b6e332a71 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sun, 20 Oct 2024 19:05:12 -0500 Subject: [PATCH 27/40] only get coords once with new XML function --- src/distribution_spatial.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index e02f74aff6d..35b39e64317 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -309,12 +309,11 @@ PointCloud::PointCloud(pugi::xml_node node) std::vector coords; if (check_for_node(node, "coords")) { - coords = get_node_array(node, "coords"); + point_cloud_ = get_node_position_array(node, "coords"); } else { fatal_error("No coordinates were provided for the PointCloud " "spatial distribution"); } - point_cloud_ = get_node_position_array(node, "coords"); std::vector strengths(point_cloud_.size(), 1.0); From cbe20587cd7cf052427bfaa0c830d995446f7e91 Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sun, 20 Oct 2024 19:06:15 -0500 Subject: [PATCH 28/40] remove extra blank line for formatting --- openmc/stats/multivariate.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 3ef38841f8c..a84c57966c2 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -869,7 +869,6 @@ def from_xml_element(cls, elem): return cls(positions, strengths) - class Box(Spatial): """Uniform distribution of coordinates in a rectangular cuboid. From 8c7aa8a6a478e2df08677d05915a45469795877e Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sun, 20 Oct 2024 19:08:19 -0500 Subject: [PATCH 29/40] allow existing variable checking to throw error --- openmc/stats/multivariate.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index a84c57966c2..c6f272c6f23 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -791,8 +791,6 @@ def positions(self): @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): From 3528f3d03147ec7b25fac022196a2978093d665a Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Sun, 20 Oct 2024 19:10:27 -0500 Subject: [PATCH 30/40] remove this file from the PR --- docker-dev-notes.txt | 8 -------- 1 file changed, 8 deletions(-) delete mode 100644 docker-dev-notes.txt diff --git a/docker-dev-notes.txt b/docker-dev-notes.txt deleted file mode 100644 index 1b600bb6f9d..00000000000 --- a/docker-dev-notes.txt +++ /dev/null @@ -1,8 +0,0 @@ -# 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 . From 8ab1dbb6772ba129d304f5a664c0c60ba2630dae Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Mon, 21 Oct 2024 09:50:45 -0500 Subject: [PATCH 31/40] add point cloud to settings docs and add missing mesh spatial info --- docs/source/io_formats/settings.rst | 53 +++++++++++++++++++++++------ 1 file changed, 42 insertions(+), 11 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 2c6c3265649..85718770865 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -579,24 +579,38 @@ attributes/sub-elements: :type: The type of spatial distribution. Valid options are "box", "fission", - "point", "cartesian", "cylindrical", and "spherical". A "box" spatial - distribution has coordinates sampled uniformly in a parallelepiped. A - "fission" spatial distribution samples locations from a "box" + "point", "cartesian", "cylindrical", "spherical", "mesh", and "cloud". + + A "box" spatial distribution has coordinates sampled uniformly in a + parallelepiped. + + A "fission" spatial distribution samples locations from a "box" distribution but only locations in fissionable materials are accepted. + A "point" spatial distribution has coordinates specified by a triplet. + A "cartesian" spatial distribution specifies independent distributions of - x-, y-, and z-coordinates. A "cylindrical" spatial distribution specifies - independent distributions of r-, phi-, and z-coordinates where phi is the - azimuthal angle and the origin for the cylindrical coordinate system is - specified by origin. A "spherical" spatial distribution specifies - independent distributions of r-, cos_theta-, and phi-coordinates where - cos_theta is the cosine of the angle with respect to the z-axis, phi is - the azimuthal angle, and the sphere is centered on the coordinate - (x0,y0,z0). A "mesh" spatial distribution samples source sites from a mesh element + x-, y-, and z-coordinates. + + A "cylindrical" spatial distribution specifies independent distributions + of r-, phi-, and z-coordinates where phi is the azimuthal angle and the + origin for the cylindrical coordinate system is specified by origin. + + A "spherical" spatial distribution specifies independent distributions of + r-, cos_theta-, and phi-coordinates where cos_theta is the cosine of the + angle with respect to the z-axis, phi is the azimuthal angle, and the + sphere is centered on the coordinate (x0,y0,z0). + + A "mesh" spatial distribution samples source sites from a mesh element based on the relative strengths provided in the node. Source locations within an element are sampled isotropically. If no strengths are provided, the space within the mesh is uniformly sampled. + A "cloud" spatial distribution samples source sites from a list of spatial + positions provided in the nodw, based on the relative strengths provided + in the node. If no strengths are provided, the positions are uniformly + sampled. + *Default*: None :parameters: @@ -662,6 +676,23 @@ attributes/sub-elements: For "cylindrical and "spherical" distributions, this element specifies the coordinates for the origin of the coordinate system. + :mesh_id: + For "mesh" spatial distributions, this element specifies which mesh ID to + use for the geometric description of the mesh. + + :coords: + For "cloud" distributions, this element specifies a list of coordinates for + each of the points in the cloud. + + :strengths: + For "mesh" and "cloud" spatial distributions, this element specifies the + relative source strength of each mesh element or each point in the cloud. + + :volume_normalized: + For "mesh" spatial distrubtions, this element specifies whether the vector + of relative strengths should be multiplied by the mesh element volume. This + is most common if the strengths represent a source per unit volume. + :angle: An element specifying the angular distribution of source sites. This element has the following attributes: From b8d7f7f04f7872c127aac0a041ff98a51f4a7b7a Mon Sep 17 00:00:00 2001 From: "Paul P.H. Wilson" Date: Mon, 21 Oct 2024 09:55:41 -0500 Subject: [PATCH 32/40] clarify that volume_normalized is optional and default false --- docs/source/io_formats/settings.rst | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 85718770865..e340e0abe26 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -689,9 +689,10 @@ attributes/sub-elements: relative source strength of each mesh element or each point in the cloud. :volume_normalized: - For "mesh" spatial distrubtions, this element specifies whether the vector - of relative strengths should be multiplied by the mesh element volume. This - is most common if the strengths represent a source per unit volume. + For "mesh" spatial distrubtions, this optional boolean element specifies + whether the vector of relative strengths should be multiplied by the mesh + element volume. This is most common if the strengths represent a source + per unit volume. Default: false :angle: An element specifying the angular distribution of source sites. This element From 43efcf4f3b029bd0e92e843dc0f484dc4b02559d Mon Sep 17 00:00:00 2001 From: Paul Wilson Date: Mon, 21 Oct 2024 09:56:21 -0500 Subject: [PATCH 33/40] Update docs/source/io_formats/settings.rst Co-authored-by: Patrick Shriwise --- docs/source/io_formats/settings.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index e340e0abe26..44cb7ceb2fa 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -607,7 +607,7 @@ attributes/sub-elements: the space within the mesh is uniformly sampled. A "cloud" spatial distribution samples source sites from a list of spatial - positions provided in the nodw, based on the relative strengths provided + positions provided in the node, based on the relative strengths provided in the node. If no strengths are provided, the positions are uniformly sampled. From 25c4f63a4cfa893eb5a7d55b6644455c7691b528 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 21 Oct 2024 10:19:43 -0500 Subject: [PATCH 34/40] Refactor PointCloud constructor a bit --- src/distribution_spatial.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 35b39e64317..c02005159e1 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -306,8 +306,6 @@ Position MeshSpatial::sample(uint64_t* seed) const PointCloud::PointCloud(pugi::xml_node node) { - std::vector coords; - if (check_for_node(node, "coords")) { point_cloud_ = get_node_position_array(node, "coords"); } else { @@ -315,16 +313,18 @@ PointCloud::PointCloud(pugi::xml_node node) "spatial distribution"); } - std::vector strengths(point_cloud_.size(), 1.0); + std::vector strengths; - if (check_for_node(node, "strengths")) { + if (check_for_node(node, "strengths")) strengths = get_node_array(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())); - } + else + strengths.resize(point_cloud_.size(), 1.0); + + 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); From b6c5c9e72479c32123f4a8a7a2ffccaeca8f065d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 21 Oct 2024 10:20:04 -0500 Subject: [PATCH 35/40] Save a live; save a life --- src/xml_interface.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index 5432b21a95a..8b53815c38c 100644 --- a/src/xml_interface.cpp +++ b/src/xml_interface.cpp @@ -58,8 +58,7 @@ vector get_node_position_array( "Incorect number of coordinates in Position array ({}) for \"{}\"", coords.size(), name)); } - vector positions; - positions.resize(coords.size() / 3); + vector positions(coords.size() / 3); auto it = coords.begin(); for (int i = 0; i < positions.size(); i++) { positions[i] = {*it++, *it++, *it++}; From 603131f4397bb1c788a331bbdec8de9cd7c1ddeb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 12 Nov 2024 23:20:36 -0600 Subject: [PATCH 36/40] Updatest test_point_cloud to catch error --- tests/unit_tests/test_source.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 5de26020a5a..7b932e1a2f9 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -50,9 +50,8 @@ 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] + positions = [(1, 0, 2), (0, 1, 0), (0, 0, 3), (4, 9, 2)] + strengths = [1, 2, 3, 4] space = openmc.stats.PointCloud(positions, strengths) np.testing.assert_equal(space.positions, positions) From 8bb6dd79368a3f316f6f072b629f1daff25f9634 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Nov 2024 09:01:10 -0600 Subject: [PATCH 37/40] Simplify PointCloud C++ implementation --- include/openmc/distribution_spatial.h | 7 ------- src/distribution_spatial.cpp | 8 ++------ src/xml_interface.cpp | 7 ++++--- 3 files changed, 6 insertions(+), 16 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 367a4129205..8ff766d1333 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -151,13 +151,6 @@ class PointCloud : public SpatialDistribution { //! \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 point_cloud_; DiscreteIndex point_idx_dist_; //!< Distribution of Position indices diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index c02005159e1..5c193a95d29 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -337,14 +337,10 @@ PointCloud::PointCloud( 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)]; + int32_t index = point_idx_dist_.sample(seed); + return point_cloud_[index]; } //============================================================================== diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index 8b53815c38c..840d3f5b871 100644 --- a/src/xml_interface.cpp +++ b/src/xml_interface.cpp @@ -58,10 +58,11 @@ vector get_node_position_array( "Incorect number of coordinates in Position array ({}) for \"{}\"", coords.size(), name)); } - vector positions(coords.size() / 3); + vector positions; + positions.reserve(coords.size() / 3); auto it = coords.begin(); - for (int i = 0; i < positions.size(); i++) { - positions[i] = {*it++, *it++, *it++}; + for (size_t i = 0; i < coords.size(); i += 3) { + positions.push_back({coords[i], coords[i + 1], coords[i + 2]}); } return positions; } From 8462db7e5aecef05e017d43f20889de049717302 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Nov 2024 09:01:19 -0600 Subject: [PATCH 38/40] Update docs --- docs/source/io_formats/settings.rst | 36 +++++++++++++++-------------- docs/source/usersguide/settings.rst | 4 ++-- 2 files changed, 21 insertions(+), 19 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 44cb7ceb2fa..34b73fc55c5 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -579,28 +579,28 @@ attributes/sub-elements: :type: The type of spatial distribution. Valid options are "box", "fission", - "point", "cartesian", "cylindrical", "spherical", "mesh", and "cloud". - + "point", "cartesian", "cylindrical", "spherical", "mesh", and "cloud". + A "box" spatial distribution has coordinates sampled uniformly in a - parallelepiped. - + parallelepiped. + A "fission" spatial distribution samples locations from a "box" distribution but only locations in fissionable materials are accepted. - + A "point" spatial distribution has coordinates specified by a triplet. - + A "cartesian" spatial distribution specifies independent distributions of - x-, y-, and z-coordinates. - + x-, y-, and z-coordinates. + A "cylindrical" spatial distribution specifies independent distributions of r-, phi-, and z-coordinates where phi is the azimuthal angle and the - origin for the cylindrical coordinate system is specified by origin. - + origin for the cylindrical coordinate system is specified by origin. + A "spherical" spatial distribution specifies independent distributions of r-, cos_theta-, and phi-coordinates where cos_theta is the cosine of the angle with respect to the z-axis, phi is the azimuthal angle, and the - sphere is centered on the coordinate (x0,y0,z0). - + sphere is centered on the coordinate (x0,y0,z0). + A "mesh" spatial distribution samples source sites from a mesh element based on the relative strengths provided in the node. Source locations within an element are sampled isotropically. If no strengths are provided, @@ -608,7 +608,7 @@ attributes/sub-elements: A "cloud" spatial distribution samples source sites from a list of spatial positions provided in the node, based on the relative strengths provided - in the node. If no strengths are provided, the positions are uniformly + in the node. If no strengths are provided, the positions are uniformly sampled. *Default*: None @@ -681,18 +681,20 @@ attributes/sub-elements: use for the geometric description of the mesh. :coords: - For "cloud" distributions, this element specifies a list of coordinates for - each of the points in the cloud. + For "cloud" distributions, this element specifies a list of coordinates + for each of the points in the cloud. :strengths: - For "mesh" and "cloud" spatial distributions, this element specifies the + For "mesh" and "cloud" spatial distributions, this element specifies the relative source strength of each mesh element or each point in the cloud. :volume_normalized: For "mesh" spatial distrubtions, this optional boolean element specifies whether the vector of relative strengths should be multiplied by the mesh element volume. This is most common if the strengths represent a source - per unit volume. Default: false + per unit volume. + + *Default*: false :angle: An element specifying the angular distribution of source sites. This element diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 3fbfa950fbb..349aa34d07a 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -193,8 +193,8 @@ distributions using spherical or cylindrical coordinates, you can use :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. It is also -possible to define a "cloud" of source points each with a different relative -probability using :class:`openmc.stats.PointCloud`. +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`, From 42eb127ce2baf8583aa5b57bfa0adf3f4d47cf49 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Nov 2024 09:01:44 -0600 Subject: [PATCH 39/40] Make fixes and add typing in PointCloud Python class --- openmc/stats/multivariate.py | 61 ++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 30 deletions(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index c6f272c6f23..2b0e8e55378 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -761,60 +761,62 @@ def from_xml_element(cls, elem, meshes): class PointCloud(Spatial): """Spatial distribution from a point cloud. - This distribution specifies a discrete list of points, - with corresponding relative probabilities. + This distribution specifies a discrete list of points, with corresponding + relative probabilities. - .. versionadded:: 0.15.x + .. versionadded:: 0.15.1 Parameters ---------- - positions : Iterable of 3-tuples + 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. + 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) + positions : numpy.ndarray + The points in space to be sampled with shape (N, 3) strengths : numpy.ndarray or None An array of relative probabilities for each mesh point """ - def __init__(self, positions, strengths=None): + def __init__( + self, + positions: Sequence[Sequence[float]], + strengths: Sequence[float] | None = None + ): self.positions = positions self.strengths = strengths @property - def positions(self): + def positions(self) -> np.ndarray: return self._positions @positions.setter - def positions(self, given_positions): - 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') + def positions(self, positions): + positions = np.array(positions, dtype=float) + if positions.ndim != 2: + raise ValueError('positions must be a 2D array') + elif positions.shape[1] != 3: + raise ValueError('positions must have 3 columns') + self._positions = positions @property - def strengths(self): + def strengths(self) -> np.ndarray: 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 + def strengths(self, strengths): + if strengths is not None: + strengths = np.array(strengths, dtype=float) + if strengths.ndim != 1: + raise ValueError('strengths must be a 1D array') + self._strengths = strengths @property - def num_strength_bins(self): + def num_strength_bins(self) -> int: if self.strengths is None: raise ValueError('Strengths are not set') return self.strengths.size @@ -829,7 +831,6 @@ def to_xml_element(self): """ element = ET.Element('space') - element.set('type', 'cloud') subelement = ET.SubElement(element, 'coords') @@ -842,7 +843,7 @@ def to_xml_element(self): return element @classmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem: ET.Element) -> PointCloud: """Generate spatial distribution from an XML element Parameters @@ -858,7 +859,7 @@ def from_xml_element(cls, elem): """ coord_data = get_text(elem, 'coords') - positions = np.asarray([float(b) for b in coord_data.split()]).reshape((3,-1)) + positions = np.array([float(b) for b in coord_data.split()]).reshape((-1, 3)) strengths = get_text(elem, 'strengths') if strengths is not None: From 66c57d16e53a6ede8e2a23dd57aba3c2ca6ccb81 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Nov 2024 09:46:31 -0600 Subject: [PATCH 40/40] Add test for invalid input to PointCloud --- openmc/stats/multivariate.py | 6 +++-- tests/unit_tests/test_source.py | 40 +++++++++++++++++++-------------- 2 files changed, 27 insertions(+), 19 deletions(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 2b0e8e55378..cd474fa9b28 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -800,7 +800,7 @@ def positions(self, positions): if positions.ndim != 2: raise ValueError('positions must be a 2D array') elif positions.shape[1] != 3: - raise ValueError('positions must have 3 columns') + raise ValueError('Each position must have 3 values') self._positions = positions @property @@ -813,6 +813,8 @@ def strengths(self, strengths): strengths = np.array(strengths, dtype=float) if strengths.ndim != 1: raise ValueError('strengths must be a 1D array') + elif strengths.size != self.positions.shape[0]: + raise ValueError('strengths must have the same length as positions') self._strengths = strengths @property @@ -821,7 +823,7 @@ def num_strength_bins(self) -> int: raise ValueError('Strengths are not set') return self.strengths.size - def to_xml_element(self): + def to_xml_element(self) -> ET.Element: """Return XML representation of the spatial distribution Returns diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 7b932e1a2f9..891f78daf0b 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -1,3 +1,4 @@ +from collections import Counter from math import pi import openmc @@ -57,10 +58,7 @@ def test_point_cloud(): 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) + src = openmc.IndependentSource(space=space) assert src.space == space np.testing.assert_equal(src.space.positions, positions) np.testing.assert_equal(src.space.strengths, strengths) @@ -71,20 +69,28 @@ def test_point_cloud(): 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] +def test_point_cloud_invalid(): + with pytest.raises(ValueError, match='2D'): + openmc.stats.PointCloud([1, 0, 2, 0, 1, 0]) - space = openmc.stats.PointCloud(positions, strengths) + with pytest.raises(ValueError, match='3 values'): + openmc.stats.PointCloud([(1, 0, 2, 3), (4, 5, 2, 3)]) - energy = openmc.stats.Discrete([1.0e6], [1.0]) - angle = openmc.stats.Isotropic() + with pytest.raises(ValueError, match='1D'): + openmc.stats.PointCloud([(1, 0, 2), (4, 5, 2)], [(1, 2), (3, 4)]) + + with pytest.raises(ValueError, match='same length'): + openmc.stats.PointCloud([(1, 0, 2), (4, 5, 2)], [1, 2, 4]) + + +def test_point_cloud_strengths(run_in_tmpdir, sphere_box_model): + positions = [(1., 0., 2.), (0., 1., 0.), (0., 0., 3.), (-1., -1., 2.)] + strengths = [1, 2, 3, 4] + space = openmc.stats.PointCloud(positions, strengths) - src = openmc.IndependentSource(space=space, angle=angle, energy=energy) model = sphere_box_model[0] model.settings.run_mode = 'fixed source' - model.settings.source = src + model.settings.source = openmc.IndependentSource(space=space) try: model.init_lib() @@ -93,11 +99,11 @@ def test_point_cloud_strengths(run_in_tmpdir, sphere_box_model): 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 + count = Counter(s.r for s in sites) + for i, (strength, position) in enumerate(zip(strengths, positions)): + sampled_strength = count[position] / n_samples expected_strength = pytest.approx(strength/sum(strengths), abs=0.02) - assert sampled_strength == expected_strength, f'Strength incorrect for {point_list[i]}' + assert sampled_strength == expected_strength, f'Strength incorrect for {positions[i]}' def test_source_file():