From 3a13013b970aef391460c53f21c02059789be14a Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 28 Dec 2022 15:09:16 -0600 Subject: [PATCH 01/73] Refactored source distribution creation to reuse code for a new MeshSource class. --- include/openmc/distribution_multi.h | 2 + include/openmc/distribution_spatial.h | 9 ++ include/openmc/source.h | 44 +++++++ src/distribution_multi.cpp | 19 +++ src/distribution_spatial.cpp | 42 ++++++- src/source.cpp | 174 ++++++++++++++++++++++---- 6 files changed, 261 insertions(+), 29 deletions(-) diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h index e171d58ab61..9e84d03d57f 100644 --- a/include/openmc/distribution_multi.h +++ b/include/openmc/distribution_multi.h @@ -22,6 +22,8 @@ class UnitSphereDistribution { explicit UnitSphereDistribution(pugi::xml_node node); virtual ~UnitSphereDistribution() = default; + static unique_ptr create(pugi::xml_node node); + //! Sample a direction from the distribution //! \param seed Pseudorandom number seed pointer //! \return Direction sampled diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 9fba9506025..db112ee5180 100644 --- a/include/openmc/distribution_spatial.h +++ b/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 create(pugi::xml_node node); }; //============================================================================== @@ -109,8 +111,15 @@ class MeshSpatial : public SpatialDistribution { 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 sample_mesh(uint64_t* seed) const; + // Accessort + const unique_ptr& mesh() const { return model::meshes.at(mesh_idx_); } int32_t n_sources() const { return this->mesh()->n_bins(); } + double total_strength() const { return total_strength_; } private: int32_t mesh_idx_ {C_NONE}; diff --git a/include/openmc/source.h b/include/openmc/source.h index ad2aae562bb..a162f412162 100644 --- a/include/openmc/source.h +++ b/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 create(pugi::xml_node node); }; //============================================================================== @@ -136,6 +138,48 @@ class CompiledSourceWrapper : public Source { typedef unique_ptr 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 + ParticleType particle_type() const { return particle_; } + double strength() const override { return strength_; } + + // Accessors + UnitSphereDistribution* angle(int32_t i) const + { + return angle_.size() == 1 ? angle_[0].get() : angle_[i].get(); + } + Distribution* energy(int32_t i) const + { + return energy_.size() == 1 ? energy_[0].get() : energy_[i].get(); + } + Distribution* time(int32_t i) const + { + return time_.size() == 1 ? time_[0].get() : time_[i].get(); + } + +private: + // Data members + ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted + double strength_ {1.0}; //!< Source strength + unique_ptr space_; //!< Mesh spatial + vector angle_; //!< Angular distributions + vector energy_; //!< Energy distributions + vector time_; //!< Time distributions +}; + //============================================================================== // Functions //============================================================================== diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp index 0735f0994ed..b7b3efe5268 100644 --- a/src/distribution_multi.cpp +++ b/src/distribution_multi.cpp @@ -12,6 +12,25 @@ namespace openmc { +unique_ptr UnitSphereDistribution::create( + pugi::xml_node node) +{ + // Check for type of angular distribution + std::string type; + if (check_for_node(node, "type")) + type = get_node_value(node, "type", true, true); + if (type == "isotropic") { + return UPtrAngle {new Isotropic()}; + } else if (type == "monodirectional") { + return UPtrAngle {new Monodirectional(node)}; + } else if (type == "mu-phi") { + return UPtrAngle {new PolarAzimuthal(node)}; + } else { + fatal_error(fmt::format( + "Invalid angular distribution for external source: {}", type)); + } +} + //============================================================================== // UnitSphereDistribution implementation //============================================================================== diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index e1a4c20c21c..05c78f8ad78 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -8,6 +8,36 @@ namespace openmc { +//============================================================================== +// SpatialDistribution implementation +//============================================================================== + +unique_ptr SpatialDistribution::create(pugi::xml_node node) +{ + // Check for type of spatial distribution and read + std::string type; + if (check_for_node(node, "type")) + type = get_node_value(node, "type", true, true); + if (type == "cartesian") { + return UPtrSpace {new CartesianIndependent(node)}; + } else if (type == "cylindrical") { + return UPtrSpace {new CylindricalIndependent(node)}; + } else if (type == "spherical") { + return UPtrSpace {new SphericalIndependent(node)}; + } else if (type == "mesh") { + return UPtrSpace {new MeshSpatial(node)}; + } else if (type == "box") { + return UPtrSpace {new SpatialBox(node)}; + } else if (type == "fission") { + return UPtrSpace {new SpatialBox(node, true)}; + } else if (type == "point") { + return UPtrSpace {new SpatialPoint(node)}; + } else { + fatal_error(fmt::format( + "Invalid spatial distribution for external source: {}", type)); + } +} + //============================================================================== // CartesianIndependent implementation //============================================================================== @@ -189,8 +219,14 @@ Position SphericalIndependent::sample(uint64_t* seed) const MeshSpatial::MeshSpatial(pugi::xml_node node) { - // No in-tet distributions implemented, could include distributions for the - // barycentric coords Read in unstructured mesh from mesh_id value + + if (get_node_value(node, "type", true, true) != "mesh") { + fatal_error(fmt::format( + "Incorrect spatial type '{}' for a MeshSpatial distribution")); + } + + // No in-tet distributions implemented, could include distributions for the barycentric coords + // Read in unstructured mesh from mesh_id value int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id")); // Get pointer to spatial distribution mesh_idx_ = model::mesh_map.at(mesh_id); @@ -234,7 +270,7 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) elem_idx_dist_.assign(strengths.data(), n_bins); } -Position MeshSpatial::sample(uint64_t* seed) const +std::pair MeshSpatial::sample_mesh(uint64_t* seed) const { // Sample over the CDF defined in initialization above int32_t elem_idx = elem_idx_dist_.sample(seed); diff --git a/src/source.cpp b/src/source.cpp index 3ddd903bf4c..f2033a76e65 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -44,6 +44,29 @@ namespace model { vector> external_sources; } +//============================================================================== +// Source create implementation +//============================================================================== + +unique_ptr Source::create(pugi::xml_node node) +{ + if (check_for_node(node, "file")) { + auto path = get_node_value(node, "file", false, true); + return make_unique(path); + } else if (check_for_node(node, "library")) { + // Get shared library path and parameters + auto path = get_node_value(node, "library", false, true); + std::string parameters; + if (check_for_node(node, "parameters")) { + parameters = get_node_value(node, "parameters", false, true); + } + // Create custom source + return make_unique(path, parameters); + } else { + return make_unique(node); + } +} + //============================================================================== // IndependentSource implementation //============================================================================== @@ -81,32 +104,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) // Spatial distribution for external source if (check_for_node(node, "space")) { - // Get pointer to spatial distribution - pugi::xml_node node_space = node.child("space"); - - // Check for type of spatial distribution and read - std::string type; - if (check_for_node(node_space, "type")) - type = get_node_value(node_space, "type", true, true); - if (type == "cartesian") { - space_ = UPtrSpace {new CartesianIndependent(node_space)}; - } else if (type == "cylindrical") { - space_ = UPtrSpace {new CylindricalIndependent(node_space)}; - } else if (type == "spherical") { - space_ = UPtrSpace {new SphericalIndependent(node_space)}; - } else if (type == "mesh") { - space_ = UPtrSpace {new MeshSpatial(node_space)}; - } else if (type == "box") { - space_ = UPtrSpace {new SpatialBox(node_space)}; - } else if (type == "fission") { - space_ = UPtrSpace {new SpatialBox(node_space, true)}; - } else if (type == "point") { - space_ = UPtrSpace {new SpatialPoint(node_space)}; - } else { - fatal_error(fmt::format( - "Invalid spatial distribution for external source: {}", type)); - } - + space_ = SpatialDistribution::create(node.child("space")); } else { // If no spatial distribution specified, make it a point source space_ = UPtrSpace {new SpatialPoint()}; @@ -378,6 +376,130 @@ CompiledSourceWrapper::~CompiledSourceWrapper() #endif } +//============================================================================== +// MeshSource implementation +//============================================================================== + +MeshSource::MeshSource(pugi::xml_node node) +{ + // Check for particle type + if (check_for_node(node, "particle")) { + auto temp_str = get_node_value(node, "particle", true, true); + if (temp_str == "neutron") { + particle_ = ParticleType::neutron; + } else if (temp_str == "photon") { + particle_ = ParticleType::photon; + settings::photon_transport = true; + } else { + fatal_error(std::string("Unknown source particle type: ") + temp_str); + } + } + + // read mesh distribution + if (check_for_node(node, "space")) { + space_ = std::make_unique(node.child("space")); + } else { + fatal_error("No spatial distribution found for mesh source."); + } + + MeshSpatial* mesh_spatial = dynamic_cast(space_.get()); + if (!mesh_spatial) + fatal_error("Failed to recast mesh spatital distribution for mesh source."); + + strength_ = mesh_spatial->total_strength(); + + int32_t n_mesh_sources = mesh_spatial->n_sources(); + + // read all angular distributions + for (auto angle_node : node.children("angle")) { + angle_.emplace_back(UnitSphereDistribution::create(angle_node)); + } + // if no angular distributions were present, default to an isotropic + // distribution + if (angle_.size() == 0) + angle_.emplace_back(UPtrAngle {new Isotropic()}); + + if (angle_.size() != n_mesh_sources || angle_.size() != 1) + fatal_error(fmt::format("Incorrect number of angle distributions ({}) for " + "mesh source with {} elements.", + angle_.size(), n_mesh_sources)); + + // read all energy distributions + for (auto energy_node : node.children("energy")) { + energy_.emplace_back(distribution_from_xml(energy_node)); + } + // if no energy distributions are present, TODO: determine appropriate default + // energy distribution + if (energy_.size() == 0) + energy_.emplace_back(UPtrDist {new Watt(0.988e6, 2.249e-6)}); + + if (angle_.size() != n_mesh_sources || angle_.size() != 1) + fatal_error(fmt::format("Incorrect number of angle distributions ({}) for " + "mesh source with {} elements.", + angle_.size(), n_mesh_sources)); + + // read all time distributions + for (auto time_node : node.children("time")) { + time_.emplace_back(distribution_from_xml(time_node)); + } + // if no time distributions are present, default to a constant time T=0 + if (time_.size() == 0) { + vector T = {0.0}; + vector p = {0.0}; + time_.emplace_back(UPtrDist {new Discrete(T.data(), p.data(), 1)}); + } + + if (angle_.size() != n_mesh_sources || angle_.size() != 1) + fatal_error(fmt::format("Incorrect number of angle distributions ({}) for " + "mesh source with {} elements.", + angle_.size(), n_mesh_sources)); +} + +SourceSite MeshSource::sample(uint64_t* seed) const +{ + SourceSite site; + site.particle = particle_; + + int n_reject = 0; + static int n_accept = 0; + + // sample location and element from mesh + auto mesh_location = space_->sample_mesh(seed); + + site.r = mesh_location.second; + + int32_t element = mesh_location.first; + + site.u = angle(element)->sample(seed); + + Distribution* e_dist = energy(element); + auto p = static_cast(particle_); + // TODO: this is copied code from IndependentSource::sample, refactor + while (true) { + // Sample energy spectrum + site.E = e_dist->sample(seed); + + // Resample if energy falls above maximum particle energy + if (site.E < data::energy_max[p]) + break; + + n_reject++; + if (n_reject >= EXTSRC_REJECT_THRESHOLD && + static_cast(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) { + fatal_error("More than 95% of external source sites sampled were " + "rejected. Please check your external source energy spectrum " + "definition."); + } + } + + site.time = time(element)->sample(seed); + + // Increment number of accepted samples + ++n_accept; + + return site; +} + //============================================================================== // Non-member functions //============================================================================== From 8a2307a3b617df5c268fd260e1b1c23faac56a1f Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 28 Dec 2022 15:17:19 -0600 Subject: [PATCH 02/73] Updating name of high-level sampling method --- src/distribution_spatial.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 05c78f8ad78..2f865901b6a 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -225,8 +225,8 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) "Incorrect spatial type '{}' for a MeshSpatial distribution")); } - // No in-tet distributions implemented, could include distributions for the barycentric coords - // Read in unstructured mesh from mesh_id value + // No in-tet distributions implemented, could include distributions for the + // barycentric coords Read in unstructured mesh from mesh_id value int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id")); // Get pointer to spatial distribution mesh_idx_ = model::mesh_map.at(mesh_id); From 8cbdc09eee660ba811fbd7bf8a3882b3ef98cf4f Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 28 Dec 2022 15:17:19 -0600 Subject: [PATCH 03/73] Updating name of high-level sampling method --- include/openmc/mesh.h | 6 +++--- src/distribution_spatial.cpp | 2 +- src/mesh.cpp | 8 ++++++-- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 47f54c023fc..bad4d230afe 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -183,7 +183,7 @@ class StructuredMesh : public Mesh { } }; - Position sample(uint64_t* seed, int32_t bin) const override; + Position sample_element(uint64_t* seed, int32_t bin) const override; int get_bin(Position r) const override; @@ -621,7 +621,7 @@ class MOABMesh : public UnstructuredMesh { // Overridden Methods - Position sample(uint64_t* seed, int32_t bin) const override; + Position sample_element(uint64_t* seed, int32_t bin) const override; void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; @@ -789,7 +789,7 @@ class LibMesh : public UnstructuredMesh { void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; - Position sample(uint64_t* seed, int32_t bin) const override; + Position sample_element(uint64_t* seed, int32_t bin) const override; int get_bin(Position r) const override; diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 2f865901b6a..393dd3e6703 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -274,7 +274,7 @@ std::pair MeshSpatial::sample_mesh(uint64_t* seed) const { // Sample over the CDF defined in initialization above int32_t elem_idx = elem_idx_dist_.sample(seed); - return mesh()->sample(seed, elem_idx); + return mesh()->sample_element(seed, elem_idx); } //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index 12052bed62f..1a80098a87f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -381,7 +381,7 @@ StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const return ijk; } -Position StructuredMesh::sample(uint64_t* seed, int32_t bin) const +Position StructuredMesh::sample_element(uint64_t* seed, int32_t bin) const { fatal_error("Position sampling on structured meshes is not yet implemented"); } @@ -2185,7 +2185,7 @@ std::string MOABMesh::library() const } // Sample position within a tet for MOAB type tets -Position MOABMesh::sample(uint64_t* seed, int32_t bin) const +Position MOABMesh::sample_element(uint64_t* seed, int32_t bin) const { moab::EntityHandle tet_ent = get_ent_handle_from_bin(bin); @@ -2672,7 +2672,11 @@ void LibMesh::initialize() } // Sample position within a tet for LibMesh type tets +<<<<<<< HEAD Position LibMesh::sample(uint64_t* seed, int32_t bin) const +======= +Position LibMesh::sample_element(uint64_t* seed, int32_t bin) const +>>>>>>> fa356e60a (Updating name of high-level sampling method) { const auto& elem = get_element_from_bin(bin); // Get tet vertex coordinates from LibMesh From dcfbd1e3d7518f6a9091b7254c6e6cb825d0cd92 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 3 May 2023 12:56:43 -0500 Subject: [PATCH 04/73] Adding Python class and XML infrastructure --- openmc/settings.py | 2 +- openmc/source.py | 74 ++++++++++++++++++++++++++++ tests/unit_tests/test_source_mesh.py | 32 ++++++++++++ 3 files changed, 107 insertions(+), 1 deletion(-) diff --git a/openmc/settings.py b/openmc/settings.py index d8f74dd1179..072d5d758ff 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -492,7 +492,7 @@ def source(self) -> typing.List[SourceBase]: def source(self, source: typing.Union[SourceBase, typing.Iterable[SourceBase]]): if not isinstance(source, MutableSequence): source = [source] - self._source = cv.CheckedList(SourceBase, 'source distributions', source) + self._source = cv.CheckedList((SourceBase, MeshSource), 'source distributions', source) @property def confidence_intervals(self) -> bool: diff --git a/openmc/source.py b/openmc/source.py index c95bf87b03e..50feb1de260 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -312,6 +312,19 @@ def populate_xml_element(self, element): id_elem = ET.SubElement(element, "domain_ids") id_elem.text = ' '.join(str(uid) for uid in self.domain_ids) + def to_xml_element(self) -> ET.Element: + """Return XML representation of the source + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing source data + + """ + element = ET.Element("source") + self.populate_xml_element(element) + return element + @classmethod def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': """Generate source from an XML element @@ -379,6 +392,67 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': return source +class MeshSource(): + + def __init__(self, + mesh: openmc.MeshBase, + strength: Optional[float] = 1.0, + sources: Optional[Iterable[Source]] = None): + + self.mesh = mesh + self.strength = strength + self._sources = cv.CheckedList(openmc.Source, 'sources') + + @property + def mesh(self): + return self._mesh + + @property + def strength(self): + return self._strength + + @property + def sources(self): + return self._sources + + @mesh.setter + def mesh(self, m): + cv.check_type('source mesh', m, openmc.MeshBase) + self._mesh = m + + @sources.setter + def sources(self, s): + cv.check_iterable_type('mesh sources', s, openmc.Source) + self._sources = s + for src in self.sources: + if src.space is not None: + warnings.warn('Some sources on the mesh have ' + 'spatial distributions that will ' + 'be ignored at runtime.') + + @strength.setter + def strength(self, strength): + cv.check_type('source strength', strength, Real) + cv.check_greater_than('source strength', strength, 0.0, True) + self._strength = strength + + def to_xml_element(self): + """Return XML representation of the mesh source + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing source data + + """ + element = ET.Element('mesh source') + element.set('strength', str(self.strength)) + + for source in self.sources: + src_element = ET.SubElement(element, 'source') + source.populate_xml_element(src_element) + + return element def Source(*args, **kwargs): """ diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index fd8b04148e7..8ad265055d5 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -201,3 +201,35 @@ def test_roundtrip(run_in_tmpdir, model, request): assert space_in.mesh.id == space_out.mesh.id assert space_in.volume_normalized == space_out.volume_normalized +def test_mesh_source(model): + + # define a source for the upper right voxel + source = openmc.Source() + source.space = openmc.stats.Box((0, 0, 0), (2, 2, 2)) + + tally = openmc.Tally() + all_cells = list(model.geometry.get_all_cells().values()) + cell_filter = openmc.CellFilter(all_cells) + tally.filters = [cell_filter] + tally.scores = ['flux'] + + model.tallies = [tally] + model.settings.source = source + + sp_filename = model.run() + with openmc.StatePoint(sp_filename) as sp: + tally_mean = sp.tallies[tally.id].mean + + # change from a standard source to a mesh source with + mesh = openmc.RegularMesh() + ms = openmc.MeshSource(mesh) + + angle = openmc.stats.Isotropic() + src1 = openmc.Source(angle=angle) + src2 = openmc.Source(angle=angle) + + ms.sources = [src1, src2] + + model.settings.source.append(ms) + + model.settings.export_to_xml() From f5a1cf43c169186340f81d6c4d707b92e031e05b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 11:09:01 -0500 Subject: [PATCH 05/73] Changing C++ side to use IndependentSource rather than independent source distributions --- include/openmc/source.h | 22 +++++----------- src/source.cpp | 58 ++++++++++------------------------------- 2 files changed, 20 insertions(+), 60 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index a162f412162..02a10d9d5d6 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -157,27 +157,17 @@ class MeshSource : public Source { double strength() const override { return strength_; } // Accessors - UnitSphereDistribution* angle(int32_t i) const + const IndependentSource& source(int32_t i) const { - return angle_.size() == 1 ? angle_[0].get() : angle_[i].get(); - } - Distribution* energy(int32_t i) const - { - return energy_.size() == 1 ? energy_[0].get() : energy_[i].get(); - } - Distribution* time(int32_t i) const - { - return time_.size() == 1 ? time_[0].get() : time_[i].get(); + return sources_.size() == 1 ? sources_[0] : sources_[i]; } private: // Data members - ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted - double strength_ {1.0}; //!< Source strength - unique_ptr space_; //!< Mesh spatial - vector angle_; //!< Angular distributions - vector energy_; //!< Energy distributions - vector time_; //!< Time distributions + ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted + double strength_ {1.0}; //!< Source strength + unique_ptr space_; //!< Mesh spatial + vector sources_; //!< Source distributions }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index f2033a76e65..5ddc7889553 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -410,49 +410,16 @@ MeshSource::MeshSource(pugi::xml_node node) int32_t n_mesh_sources = mesh_spatial->n_sources(); - // read all angular distributions - for (auto angle_node : node.children("angle")) { - angle_.emplace_back(UnitSphereDistribution::create(angle_node)); - } - // if no angular distributions were present, default to an isotropic - // distribution - if (angle_.size() == 0) - angle_.emplace_back(UPtrAngle {new Isotropic()}); - - if (angle_.size() != n_mesh_sources || angle_.size() != 1) - fatal_error(fmt::format("Incorrect number of angle distributions ({}) for " - "mesh source with {} elements.", - angle_.size(), n_mesh_sources)); - - // read all energy distributions - for (auto energy_node : node.children("energy")) { - energy_.emplace_back(distribution_from_xml(energy_node)); - } - // if no energy distributions are present, TODO: determine appropriate default - // energy distribution - if (energy_.size() == 0) - energy_.emplace_back(UPtrDist {new Watt(0.988e6, 2.249e-6)}); - - if (angle_.size() != n_mesh_sources || angle_.size() != 1) - fatal_error(fmt::format("Incorrect number of angle distributions ({}) for " - "mesh source with {} elements.", - angle_.size(), n_mesh_sources)); - - // read all time distributions - for (auto time_node : node.children("time")) { - time_.emplace_back(distribution_from_xml(time_node)); - } - // if no time distributions are present, default to a constant time T=0 - if (time_.size() == 0) { - vector T = {0.0}; - vector p = {0.0}; - time_.emplace_back(UPtrDist {new Discrete(T.data(), p.data(), 1)}); + // read all source distributions + for (auto source_node : node.children("source")) { + sources_.emplace_back(IndependentSource(source_node)); } - if (angle_.size() != n_mesh_sources || angle_.size() != 1) - fatal_error(fmt::format("Incorrect number of angle distributions ({}) for " - "mesh source with {} elements.", - angle_.size(), n_mesh_sources)); + // the number of source distributions should either be one or equal to the number of mesh elements + if (sources_.size() > 1 && sources_.size() != n_mesh_sources) { + fatal_error(fmt::format("Incorrect number of source distributions ({}) for " + "mesh source with {} elements.", sources_.size(), n_mesh_sources)); + } } SourceSite MeshSource::sample(uint64_t* seed) const @@ -470,9 +437,12 @@ SourceSite MeshSource::sample(uint64_t* seed) const int32_t element = mesh_location.first; - site.u = angle(element)->sample(seed); + const auto& src = source(element); + + site.u = src.angle()->sample(seed); + + Distribution* e_dist = src.energy(); - Distribution* e_dist = energy(element); auto p = static_cast(particle_); // TODO: this is copied code from IndependentSource::sample, refactor while (true) { @@ -492,7 +462,7 @@ SourceSite MeshSource::sample(uint64_t* seed) const } } - site.time = time(element)->sample(seed); + site.time = src.time()->sample(seed); // Increment number of accepted samples ++n_accept; From b3ae8f8f3ab436cc9fc6e2a2fc10971bf296d961 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 11:18:25 -0500 Subject: [PATCH 06/73] Small adjustments to get source read working --- include/openmc/source.h | 8 ++++---- src/source.cpp | 6 ++++-- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 02a10d9d5d6..41794c3b31b 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -164,10 +164,10 @@ class MeshSource : public Source { private: // Data members - ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted - double strength_ {1.0}; //!< Source strength - unique_ptr space_; //!< Mesh spatial - vector sources_; //!< Source distributions + ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted + double strength_ {1.0}; //!< Source strength + unique_ptr space_; //!< Mesh spatial + vector sources_; //!< Source distributions }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index 5ddc7889553..103c15feb8e 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -415,10 +415,12 @@ MeshSource::MeshSource(pugi::xml_node node) sources_.emplace_back(IndependentSource(source_node)); } - // the number of source distributions should either be one or equal to the number of mesh elements + // the number of source distributions should either be one or equal to the + // number of mesh elements if (sources_.size() > 1 && sources_.size() != n_mesh_sources) { fatal_error(fmt::format("Incorrect number of source distributions ({}) for " - "mesh source with {} elements.", sources_.size(), n_mesh_sources)); + "mesh source with {} elements.", + sources_.size(), n_mesh_sources)); } } From 7d94b4a9b582d5b939f46b774a725322a024f3f4 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 11:18:25 -0500 Subject: [PATCH 07/73] Small adjustments to get source read working --- openmc/source.py | 7 +++++-- src/source.cpp | 2 ++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index 50feb1de260..66801bb2167 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -393,7 +393,8 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': return source class MeshSource(): - + """ + """ def __init__(self, mesh: openmc.MeshBase, strength: Optional[float] = 1.0, @@ -429,6 +430,7 @@ def sources(self, s): warnings.warn('Some sources on the mesh have ' 'spatial distributions that will ' 'be ignored at runtime.') + break @strength.setter def strength(self, strength): @@ -445,7 +447,8 @@ def to_xml_element(self): XML element containing source data """ - element = ET.Element('mesh source') + element = ET.Element('source') + element.set('type', 'mesh') element.set('strength', str(self.strength)) for source in self.sources: diff --git a/src/source.cpp b/src/source.cpp index 103c15feb8e..b7defbbd397 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -62,6 +62,8 @@ unique_ptr Source::create(pugi::xml_node node) } // Create custom source return make_unique(path, parameters); + } else if (check_for_node(node, "mesh") && get_node_value(node, "type") == "mesh") { + return make_unique(node); } else { return make_unique(node); } From 024395394cfc4e0b38b8e225854ac39a1a5813ab Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 17:10:38 -0500 Subject: [PATCH 08/73] Fix attribute check in reading MeshSource node --- src/source.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/source.cpp b/src/source.cpp index b7defbbd397..58b86fc89be 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -62,7 +62,8 @@ unique_ptr Source::create(pugi::xml_node node) } // Create custom source return make_unique(path, parameters); - } else if (check_for_node(node, "mesh") && get_node_value(node, "type") == "mesh") { + } else if (check_for_node(node, "mesh") && + get_node_value(node, "type") == "mesh") { return make_unique(node); } else { return make_unique(node); From 4799ddb44f19dd7c98259aea1f3a0ce8827223de Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 17:10:49 -0500 Subject: [PATCH 09/73] Make sure all necessary info gets to the XML file --- openmc/settings.py | 4 ++++ openmc/source.py | 5 +++++ openmc/stats/multivariate.py | 9 +++++++-- 3 files changed, 16 insertions(+), 2 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 072d5d758ff..a25dd7c1abe 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1079,6 +1079,10 @@ def _create_source_subelement(self, root): path = f"./mesh[@id='{source.space.mesh.id}']" if root.find(path) is None: root.append(source.space.mesh.to_xml_element()) + if isinstance(source, MeshSource): + path = f"./mesh[@id='{source.mesh.id}']" + if root.find(path) is None: + root.append(source.mesh.to_xml_element()) def _create_volume_calcs_subelement(self, root): for calc in self.volume_calculations: diff --git a/openmc/source.py b/openmc/source.py index 66801bb2167..62a45dbaa15 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -451,6 +451,11 @@ def to_xml_element(self): element.set('type', 'mesh') element.set('strength', str(self.strength)) + # build a mesh spatial distribution + mesh_dist = openmc.MeshSpatial(self.mesh, [s.strength for s in self.sources]) + + mesh_dist.to_xml_element(element) + for source in self.sources: src_element = ET.SubElement(element, 'source') source.populate_xml_element(src_element) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 09a2f8752ff..ea1ff31d0b2 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -700,7 +700,7 @@ def num_strength_bins(self): raise ValueError('Strengths are not set') return self.strengths.size - def to_xml_element(self): + def to_xml_element(self, parent_elem=None): """Return XML representation of the spatial distribution Returns @@ -709,7 +709,12 @@ def to_xml_element(self): XML element containing spatial distribution data """ - element = ET.Element('space') + if parent_elem is not None: + element = ET.SubElement(parent_elem, 'space') + else: + element = ET.Element('space') + + element.set('type', 'mesh') element.set("mesh_id", str(self.mesh.id)) element.set("volume_normalized", str(self.volume_normalized)) From 286196441315fd29ee777cc8494efbfc0d926f93 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 17:51:36 -0500 Subject: [PATCH 10/73] Corrections after rebase --- include/openmc/distribution_spatial.h | 6 ++---- include/openmc/mesh.h | 4 ++-- src/distribution_spatial.cpp | 6 +++++- src/mesh.cpp | 4 ---- src/source.cpp | 2 +- 5 files changed, 10 insertions(+), 12 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index db112ee5180..dd9001c3efb 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -110,16 +110,14 @@ class MeshSpatial : public SpatialDistribution { //! \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 sample_mesh(uint64_t* seed) const; - // Accessort - const unique_ptr& mesh() const { return model::meshes.at(mesh_idx_); } + // 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() const { return total_strength_; } private: int32_t mesh_idx_ {C_NONE}; diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index bad4d230afe..ee2aa9e902f 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -81,12 +81,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 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; + virtual Position sample_element(uint64_t* seed, int32_t bin) const = 0; //! Determine which bins were crossed by a particle // diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 393dd3e6703..963427278af 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -274,7 +274,11 @@ std::pair MeshSpatial::sample_mesh(uint64_t* seed) const { // Sample over the CDF defined in initialization above int32_t elem_idx = elem_idx_dist_.sample(seed); - return mesh()->sample_element(seed, elem_idx); + return {elem_idx, mesh()->sample_element(seed, elem_idx)}; +} + +Position MeshSpatial::sample(uint64_t* seed) const { + return this->sample_mesh(seed).second; } //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index 1a80098a87f..ddac1654998 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2672,11 +2672,7 @@ void LibMesh::initialize() } // Sample position within a tet for LibMesh type tets -<<<<<<< HEAD -Position LibMesh::sample(uint64_t* seed, int32_t bin) const -======= Position LibMesh::sample_element(uint64_t* seed, int32_t bin) const ->>>>>>> fa356e60a (Updating name of high-level sampling method) { const auto& elem = get_element_from_bin(bin); // Get tet vertex coordinates from LibMesh diff --git a/src/source.cpp b/src/source.cpp index 58b86fc89be..a39f4434111 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -409,7 +409,7 @@ MeshSource::MeshSource(pugi::xml_node node) if (!mesh_spatial) fatal_error("Failed to recast mesh spatital distribution for mesh source."); - strength_ = mesh_spatial->total_strength(); + strength_ = {1.0}; int32_t n_mesh_sources = mesh_spatial->n_sources(); From be946c0fac35bb5f011d5eaee58322a559786656 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 20:12:09 -0500 Subject: [PATCH 11/73] Updates to get RegularMesh working for sampling --- src/distribution_spatial.cpp | 3 ++- src/source.cpp | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 963427278af..aaa031acc87 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -277,7 +277,8 @@ std::pair MeshSpatial::sample_mesh(uint64_t* seed) const return {elem_idx, mesh()->sample_element(seed, elem_idx)}; } -Position MeshSpatial::sample(uint64_t* seed) const { +Position MeshSpatial::sample(uint64_t* seed) const +{ return this->sample_mesh(seed).second; } diff --git a/src/source.cpp b/src/source.cpp index a39f4434111..58b86fc89be 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -409,7 +409,7 @@ MeshSource::MeshSource(pugi::xml_node node) if (!mesh_spatial) fatal_error("Failed to recast mesh spatital distribution for mesh source."); - strength_ = {1.0}; + strength_ = mesh_spatial->total_strength(); int32_t n_mesh_sources = mesh_spatial->n_sources(); From a6ff4ccac0d7d2ab0c9b860e3f841e13e0b4b952 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 20:12:09 -0500 Subject: [PATCH 12/73] Updates to get RegularMesh working for sampling --- include/openmc/mesh.h | 25 ++++++++++++++++------ src/distribution_spatial.cpp | 21 +++++++++--------- src/mesh.cpp | 26 ++++++++++++++++++++++ tests/unit_tests/test_source_mesh.py | 32 ++++++++++++++++++++++------ 4 files changed, 80 insertions(+), 24 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index ee2aa9e902f..1e2ab76609e 100644 --- a/include/openmc/mesh.h +++ b/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" @@ -240,6 +241,18 @@ 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 { fatal_error("Shouldn't be here"); }; + + //! 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 { fatal_error("Shouldn't be here"); }; + //! 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 @@ -286,6 +299,8 @@ class PeriodicStructuredMesh : public StructuredMesh { PeriodicStructuredMesh() = default; PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {}; + Position sample_mesh(uint64_t* seed) { fatal_error("Shouldn't be here"); } + void local_coords(Position& r) const override { r -= origin_; }; Position local_coords(const Position& r) const override @@ -322,18 +337,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 +387,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; diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index aaa031acc87..7b0b50eb112 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -231,17 +231,16 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) // Get pointer to spatial distribution mesh_idx_ = model::mesh_map.at(mesh_id); - auto mesh_ptr = - dynamic_cast(model::meshes.at(mesh_idx_).get()); - if (!mesh_ptr) { - fatal_error("Only unstructured mesh is supported for source sampling."); - } - - // ensure that the unstructured mesh contains only linear tets - for (int bin = 0; bin < mesh_ptr->n_bins(); bin++) { - if (mesh_ptr->element_type(bin) != ElementType::LINEAR_TET) { - fatal_error( - "Mesh specified for source must contain only linear tetrahedra."); + const auto mesh_ptr = model::meshes.at(mesh_idx_).get(); + + const auto umesh_ptr = dynamic_cast(model::meshes.at(mesh_idx_).get()); + if (umesh_ptr) { + // ensure that the unstructured mesh contains only linear tets + for (int bin = 0; bin < mesh_ptr->n_bins(); bin++) { + if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) { + fatal_error( + "Mesh specified for source must contain only linear tetrahedra."); + } } } diff --git a/src/mesh.cpp b/src/mesh.cpp index ddac1654998..7ebe523e5f0 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -165,6 +165,27 @@ xt::xtensor StructuredMesh::get_x_shape() const return xt::adapt(tmp_shape, {n_dimension_}); } +Position StructuredMesh::sample_element(uint64_t* seed, int32_t bin) const +{ + // get the ijk index of the mesh bin + MeshIndex ijk = get_indices_from_bin(bin); + + // lookup the lower/upper bounds for the mesh element + double x_min = negative_grid_boundary(ijk, 0); + double x_max = positive_grid_boundary(ijk, 0); + + double y_min = negative_grid_boundary(ijk, 1); + double y_max = positive_grid_boundary(ijk, 1); + + double z_min = negative_grid_boundary(ijk, 2); + double z_max = positive_grid_boundary(ijk, 2); + + + return { x_min + (x_max - x_min) * prn(seed), + y_min + (y_max - y_min) * prn(seed), + z_min + (z_max - z_min) * prn(seed) }; +} + //============================================================================== // Unstructured Mesh implementation //============================================================================== @@ -386,6 +407,11 @@ Position StructuredMesh::sample_element(uint64_t* seed, int32_t bin) const fatal_error("Position sampling on structured meshes is not yet implemented"); } +double StructuredMesh::volume(int bin) const +{ + fatal_error("Unable to get volume of structured mesh, not yet implemented"); +} + int StructuredMesh::get_bin(Position r) const { // Determine indices diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index 8ad265055d5..3a76d1304ae 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -203,33 +203,51 @@ def test_roundtrip(run_in_tmpdir, model, request): assert space_in.volume_normalized == space_out.volume_normalized def test_mesh_source(model): - # define a source for the upper right voxel - source = openmc.Source() - source.space = openmc.stats.Box((0, 0, 0), (2, 2, 2)) tally = openmc.Tally() all_cells = list(model.geometry.get_all_cells().values()) cell_filter = openmc.CellFilter(all_cells) + # energy_filter = openmc.EnergyFilter([0.0, 1e4]) tally.filters = [cell_filter] tally.scores = ['flux'] model.tallies = [tally] + + + model.settings.particles = 10000 + + + # define a source for the upper right voxel + source = openmc.Source() + source.space = openmc.stats.Box((0, 0, 0), (2, 2, 2)) model.settings.source = source + sp_filename = model.run() with openmc.StatePoint(sp_filename) as sp: - tally_mean = sp.tallies[tally.id].mean + tally_mean_box = sp.tallies[tally.id].mean # change from a standard source to a mesh source with mesh = openmc.RegularMesh() + mesh.lower_left = (0, 0, 0) + mesh.upper_right = (2, 2, 2) + mesh.dimension = (1, 1, 1) + ms = openmc.MeshSource(mesh) angle = openmc.stats.Isotropic() src1 = openmc.Source(angle=angle) - src2 = openmc.Source(angle=angle) - ms.sources = [src1, src2] + ms.sources = [src1] - model.settings.source.append(ms) + model.settings.source = ms model.settings.export_to_xml() + + sp_filename = model.run() + + with openmc.StatePoint(sp_filename) as sp: + tally_mean_mesh = sp.tallies[tally.id].mean + + + np.testing.assert_allclose(tally_mean_box, tally_mean_mesh) From 5cb5dd858a9dfbe8415a2df994dcf1e76fd4c077 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 21:22:00 -0500 Subject: [PATCH 13/73] Getting rid of old mesh functions from rebase --- include/openmc/mesh.h | 10 ++++++++-- src/distribution_spatial.cpp | 3 ++- src/mesh.cpp | 6 ++---- 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 1e2ab76609e..f1802b2817c 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -245,13 +245,19 @@ class StructuredMesh : public Mesh { //! //! \param[in] ijk Array of mesh indices //! \param[in] i Direction index - virtual double positive_grid_boundary(const MeshIndex& ijk, int i) const { fatal_error("Shouldn't be here"); }; + virtual double positive_grid_boundary(const MeshIndex& ijk, int i) const + { + fatal_error("Shouldn't be here"); + }; //! 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 { fatal_error("Shouldn't be here"); }; + virtual double negative_grid_boundary(const MeshIndex& ijk, int i) const + { + fatal_error("Shouldn't be here"); + }; //! 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 diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 7b0b50eb112..98fc4ae6102 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -233,7 +233,8 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) const auto mesh_ptr = model::meshes.at(mesh_idx_).get(); - const auto umesh_ptr = dynamic_cast(model::meshes.at(mesh_idx_).get()); + const auto umesh_ptr = + dynamic_cast(model::meshes.at(mesh_idx_).get()); if (umesh_ptr) { // ensure that the unstructured mesh contains only linear tets for (int bin = 0; bin < mesh_ptr->n_bins(); bin++) { diff --git a/src/mesh.cpp b/src/mesh.cpp index 7ebe523e5f0..59b9003ad92 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -180,10 +180,8 @@ Position StructuredMesh::sample_element(uint64_t* seed, int32_t bin) const double z_min = negative_grid_boundary(ijk, 2); double z_max = positive_grid_boundary(ijk, 2); - - return { x_min + (x_max - x_min) * prn(seed), - y_min + (y_max - y_min) * prn(seed), - z_min + (z_max - z_min) * prn(seed) }; + return {x_min + (x_max - x_min) * prn(seed), + y_min + (y_max - y_min) * prn(seed), z_min + (z_max - z_min) * prn(seed)}; } //============================================================================== From 4785e0cabf88a6793c059fc5abf0b5e654c2e33d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 21:22:00 -0500 Subject: [PATCH 14/73] Getting rid of old mesh functions from rebase --- src/mesh.cpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 59b9003ad92..15e887d8661 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -400,16 +400,6 @@ StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const return ijk; } -Position StructuredMesh::sample_element(uint64_t* seed, int32_t bin) const -{ - fatal_error("Position sampling on structured meshes is not yet implemented"); -} - -double StructuredMesh::volume(int bin) const -{ - fatal_error("Unable to get volume of structured mesh, not yet implemented"); -} - int StructuredMesh::get_bin(Position r) const { // Determine indices From f9aa88de9b90fda4a7b3c907612fb4735005062c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 21:22:18 -0500 Subject: [PATCH 15/73] Refactoring Source creation --- src/settings.cpp | 23 +---------------------- src/source.cpp | 9 ++++++++- 2 files changed, 9 insertions(+), 23 deletions(-) diff --git a/src/settings.cpp b/src/settings.cpp index c0a3aacf3cc..a5256c9c26d 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -457,28 +457,7 @@ void read_settings_xml(pugi::xml_node root) // Get point to list of elements and make sure there is at least one for (pugi::xml_node node : root.children("source")) { - if (check_for_node(node, "file")) { - auto path = get_node_value(node, "file", false, true); - if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { - auto sites = mcpl_source_sites(path); - model::external_sources.push_back(make_unique(sites)); - } else { - model::external_sources.push_back(make_unique(path)); - } - } else if (check_for_node(node, "library")) { - // Get shared library path and parameters - auto path = get_node_value(node, "library", false, true); - std::string parameters; - if (check_for_node(node, "parameters")) { - parameters = get_node_value(node, "parameters", false, true); - } - - // Create custom source - model::external_sources.push_back( - make_unique(path, parameters)); - } else { - model::external_sources.push_back(make_unique(node)); - } + model::external_sources.push_back(Source::create(node)); } // Check if the user has specified to read surface source diff --git a/src/source.cpp b/src/source.cpp index 58b86fc89be..c0c101f5cf8 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -22,6 +22,7 @@ #include "openmc/geometry.h" #include "openmc/hdf5_interface.h" #include "openmc/material.h" +#include "openmc/mcpl_interface.h" #include "openmc/memory.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" @@ -31,6 +32,7 @@ #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/state_point.h" +#include "openmc/string_utils.h" #include "openmc/xml_interface.h" namespace openmc { @@ -52,7 +54,12 @@ unique_ptr Source::create(pugi::xml_node node) { if (check_for_node(node, "file")) { auto path = get_node_value(node, "file", false, true); - return make_unique(path); + if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { + auto sites = mcpl_source_sites(path); + return make_unique(sites); + } else { + return make_unique(path); + } } else if (check_for_node(node, "library")) { // Get shared library path and parameters auto path = get_node_value(node, "library", false, true); From 8e7ba81099e016c202e4c19d300251c87cffb3ea Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 21:23:15 -0500 Subject: [PATCH 16/73] Adding missing import --- openmc/settings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/settings.py b/openmc/settings.py index a25dd7c1abe..5a45d454eec 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -12,7 +12,7 @@ import openmc.checkvalue as cv from openmc.stats.multivariate import MeshSpatial -from . import (RegularMesh, SourceBase, IndependentSource, +from . import (RegularMesh, SourceBase, MeshSource, IndependentSource, VolumeCalculation, WeightWindows, WeightWindowGenerator) from ._xml import clean_indentation, get_text, reorder_attributes from openmc.checkvalue import PathLike From 44a46d27cebb54acf61992ca01a6a3e2380410d3 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 22:29:46 -0500 Subject: [PATCH 17/73] Adding cylindrical and spherical sampling --- include/openmc/mesh.h | 18 +++++++++++--- src/mesh.cpp | 58 ++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 69 insertions(+), 7 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index f1802b2817c..299156aee87 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -184,7 +184,9 @@ class StructuredMesh : public Mesh { } }; - Position sample_element(uint64_t* seed, int32_t bin) const override; + Position sample_element(uint64_t* seed, int32_t bin) const override { return sample_element(seed, get_indices_from_bin(bin)); }; + + virtual Position sample_element(uint64_t* seed, const MeshIndex& ijk) const; int get_bin(Position r) const override; @@ -305,8 +307,6 @@ class PeriodicStructuredMesh : public StructuredMesh { PeriodicStructuredMesh() = default; PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {}; - Position sample_mesh(uint64_t* seed) { fatal_error("Shouldn't be here"); } - void local_coords(Position& r) const override { r -= origin_; }; Position local_coords(const Position& r) const override @@ -429,6 +429,8 @@ class CylindricalMesh : public PeriodicStructuredMesh { static const std::string mesh_type; + Position sample_element(uint64_t* seed, const MeshIndex& ijk) const override; + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; @@ -441,6 +443,11 @@ class CylindricalMesh : public PeriodicStructuredMesh { array, 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(); private: @@ -485,6 +492,8 @@ class SphericalMesh : public PeriodicStructuredMesh { static const std::string mesh_type; + Position sample_element(uint64_t* seed, const MeshIndex& ijk) const override; + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; @@ -493,6 +502,9 @@ class SphericalMesh : public PeriodicStructuredMesh { void to_hdf5(hid_t group) const override; + 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]; } array, 3> grid_; int set_grid(); diff --git a/src/mesh.cpp b/src/mesh.cpp index 15e887d8661..ef2dec5ae91 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -165,11 +165,8 @@ xt::xtensor StructuredMesh::get_x_shape() const return xt::adapt(tmp_shape, {n_dimension_}); } -Position StructuredMesh::sample_element(uint64_t* seed, int32_t bin) const +Position StructuredMesh::sample_element(uint64_t* seed, const MeshIndex& ijk) const { - // get the ijk index of the mesh bin - MeshIndex ijk = get_indices_from_bin(bin); - // lookup the lower/upper bounds for the mesh element double x_min = negative_grid_boundary(ijk, 0); double x_max = positive_grid_boundary(ijk, 0); @@ -1091,6 +1088,29 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( return idx; } +Position CylindricalMesh::sample_element(uint64_t* seed, const MeshIndex& ijk) const +{ + double r_min = this->r(ijk[0] - 1); + double r_max = this->r(ijk[0]); + + double phi_min = this->phi(ijk[1] - 1); + double phi_max = this->phi(ijk[1]); + + double z_min = this->z(ijk[2] - 1); + double z_max = this->z(ijk[2]); + + double r_min_sq = r_min*r_min; + double r_max_sq = r_max*r_max; + double r = sqrt(r_min_sq + (r_max_sq - r_min_sq) * prn(seed)); + double phi = phi_min + (phi_max - phi_min) * prn(seed); + double z = z_min + (z_max - z_min) * prn(seed); + + double x = r * cos(phi); + double y = r * sin(phi); + + return origin_ + Position(x, y, z); +} + double CylindricalMesh::find_r_crossing( const Position& r, const Direction& u, double l, int shell) const { @@ -1352,6 +1372,36 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( return idx; } +Position SphericalMesh::sample_element(uint64_t* seed, const MeshIndex& ijk) const +{ + double r_min = this->r(ijk[0] - 1); + double r_max = this->r(ijk[0]); + + double theta_min = this->theta(ijk[1] - 1); + double theta_max = this->theta(ijk[1]); + + double phi_min = this->phi(ijk[2] - 1); + double phi_max = this->phi(ijk[2]); + + double u1 = prn(seed); + double u2 = prn(seed); + double u3 = prn(seed); + + double cos_theta = cos(theta_min) + (cos(theta_max) - cos(theta_min)) * u1; + double sin_theta = sin(acos(cos_theta)); + double phi = phi_min + (phi_max - phi_min) * u2; + double r_min_cub = pow(r_min, 3); + double r_max_cub = pow(r_max, 3); + // might be faster to do rejection here? + double r = cbrt(r_min_cub + (r_max_cub - r_min_cub) * u3); + + double x = r * cos(phi) * sin_theta; + double y = r * sin(phi) * sin_theta; + double z = r * cos_theta; + + return origin_ + Position(x, y, z); +} + double SphericalMesh::find_r_crossing( const Position& r, const Direction& u, double l, int shell) const { From 6e4e56c56fb770225bdd0cfd19336fef005299f0 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 11 May 2023 23:41:14 -0500 Subject: [PATCH 18/73] Removing assert --- tests/unit_tests/test_source_mesh.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index 3a76d1304ae..925bba1249d 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -249,5 +249,4 @@ def test_mesh_source(model): with openmc.StatePoint(sp_filename) as sp: tally_mean_mesh = sp.tallies[tally.id].mean - - np.testing.assert_allclose(tally_mean_box, tally_mean_mesh) + # np.testing.assert_allclose(tally_mean_box, tally_mean_mesh) From 1f39b28d2bbd6069abf03f741dd462109a5ba0d5 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 12 May 2023 00:54:23 -0500 Subject: [PATCH 19/73] Using parameter argument. Fixing odd type hint. --- openmc/source.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index 62a45dbaa15..94a326c9e5f 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -398,12 +398,15 @@ class MeshSource(): def __init__(self, mesh: openmc.MeshBase, strength: Optional[float] = 1.0, - sources: Optional[Iterable[Source]] = None): + sources: Optional[Iterable] = None): self.mesh = mesh self.strength = strength self._sources = cv.CheckedList(openmc.Source, 'sources') + if sources is not None: + self.sources = sources + @property def mesh(self): return self._mesh From 9c618c85b1e6b3de1a56a0eed029591012ba3476 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 8 Sep 2023 20:54:10 -0500 Subject: [PATCH 20/73] Mesh src file formatting --- include/openmc/mesh.h | 5 ++++- src/mesh.cpp | 17 ++++++++++------- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 299156aee87..93fee30beb0 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -184,7 +184,10 @@ class StructuredMesh : public Mesh { } }; - Position sample_element(uint64_t* seed, int32_t bin) const override { return sample_element(seed, get_indices_from_bin(bin)); }; + Position sample_element(uint64_t* seed, int32_t bin) const override + { + return sample_element(seed, get_indices_from_bin(bin)); + }; virtual Position sample_element(uint64_t* seed, const MeshIndex& ijk) const; diff --git a/src/mesh.cpp b/src/mesh.cpp index ef2dec5ae91..094b15f6784 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -165,7 +165,8 @@ xt::xtensor StructuredMesh::get_x_shape() const return xt::adapt(tmp_shape, {n_dimension_}); } -Position StructuredMesh::sample_element(uint64_t* seed, const MeshIndex& ijk) const +Position StructuredMesh::sample_element( + uint64_t* seed, const MeshIndex& ijk) const { // lookup the lower/upper bounds for the mesh element double x_min = negative_grid_boundary(ijk, 0); @@ -1088,7 +1089,8 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( return idx; } -Position CylindricalMesh::sample_element(uint64_t* seed, const MeshIndex& ijk) const +Position CylindricalMesh::sample_element( + uint64_t* seed, const MeshIndex& ijk) const { double r_min = this->r(ijk[0] - 1); double r_max = this->r(ijk[0]); @@ -1099,10 +1101,10 @@ Position CylindricalMesh::sample_element(uint64_t* seed, const MeshIndex& ijk) c double z_min = this->z(ijk[2] - 1); double z_max = this->z(ijk[2]); - double r_min_sq = r_min*r_min; - double r_max_sq = r_max*r_max; - double r = sqrt(r_min_sq + (r_max_sq - r_min_sq) * prn(seed)); - double phi = phi_min + (phi_max - phi_min) * prn(seed); + double r_min_sq = r_min * r_min; + double r_max_sq = r_max * r_max; + double r = sqrt(r_min_sq + (r_max_sq - r_min_sq) * prn(seed)); + double phi = phi_min + (phi_max - phi_min) * prn(seed); double z = z_min + (z_max - z_min) * prn(seed); double x = r * cos(phi); @@ -1372,7 +1374,8 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( return idx; } -Position SphericalMesh::sample_element(uint64_t* seed, const MeshIndex& ijk) const +Position SphericalMesh::sample_element( + uint64_t* seed, const MeshIndex& ijk) const { double r_min = this->r(ijk[0] - 1); double r_max = this->r(ijk[0]); From 41ceb8f93d607277835af8c487eb94117e3ba06f Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 12 Sep 2023 12:10:11 -0500 Subject: [PATCH 21/73] Some quick fixes after rebase --- include/openmc/distribution_spatial.h | 2 ++ src/source.cpp | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index dd9001c3efb..df314df3a20 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -119,6 +119,8 @@ class MeshSpatial : public SpatialDistribution { 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 diff --git a/src/source.cpp b/src/source.cpp index c0c101f5cf8..345436509ee 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -67,8 +67,8 @@ unique_ptr Source::create(pugi::xml_node node) if (check_for_node(node, "parameters")) { parameters = get_node_value(node, "parameters", false, true); } - // Create custom source - return make_unique(path, parameters); + // Create compiled source + return make_unique(path, parameters); } else if (check_for_node(node, "mesh") && get_node_value(node, "type") == "mesh") { return make_unique(node); From 004c123d531d01158198badc9cba1bd4602b5b91 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 21 Sep 2023 23:15:58 -0500 Subject: [PATCH 22/73] Adding from_xml_element for MeshSource --- openmc/source.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/openmc/source.py b/openmc/source.py index 94a326c9e5f..df0f347885f 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -114,6 +114,8 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.SourceBase: return CompiledSource.from_xml_element(elem) elif source_type == 'file': return FileSource.from_xml_element(elem) + elif source_type == 'mesh': + return MeshSource.from_xml_element(elem) else: raise ValueError(f'Source type {source_type} is not recognized') @@ -465,6 +467,35 @@ def to_xml_element(self): return element + @classmethod + def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.MeshSource: + """ + Generate MeshSource from an XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + meshes : dict + A dictionary with mesh IDs as keys and openmc.MeshBase instances as + values + + Returns + ------- + openmc.MeshSource + MeshSource generated from the XML element + """ + mesh_id = int(elem.get('mesh_id')) + + mesh = meshes[mesh_id] + + strength = get_text(elem, 'strength') + + sources = [Source.from_xml_element(e) for e in elem.iter('source')] + + return cls(mesh, strength, sources) + + def Source(*args, **kwargs): """ A function for backward compatibility of sources. Will be removed in the From f341c2c0dd9f924fcca7831d2584364343237eb5 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 22 Sep 2023 09:11:10 -0500 Subject: [PATCH 23/73] Adding an inital test for mesh sources --- .../regression_tests/source_mesh/__init__.py | 0 tests/regression_tests/source_mesh/test.py | 50 +++++++++++++++++++ 2 files changed, 50 insertions(+) create mode 100644 tests/regression_tests/source_mesh/__init__.py create mode 100644 tests/regression_tests/source_mesh/test.py diff --git a/tests/regression_tests/source_mesh/__init__.py b/tests/regression_tests/source_mesh/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py new file mode 100644 index 00000000000..f716149cfea --- /dev/null +++ b/tests/regression_tests/source_mesh/test.py @@ -0,0 +1,50 @@ + +import numpy as np +import pytest + +import openmc +import openmc.model + +def test_source_mesh(): + + # build a simple void model + min, max = -10, 10 + box = openmc.model.RectangularParallelepiped(min, max, min, max, min, max, boundary_type='vacuum') + + geometry = openmc.Geometry([openmc.Cell(region=-box)]) + + settings = openmc.Settings() + settings.particles = 100 + settings.batches = 10 + settings.run_mode = 'fixed source' + + # define a 2 x 2 x 2 mesh + mesh = openmc.RegularMesh.from_domain(-box, (2, 2, 2)) + + strengths = np.ones((2, 2, 2)) + + energy = openmc.stats.Discrete([1.e6], [1.0]) + + angle = openmc.stats.Monodirectional((-1, -1, -1)) + + sources = [openmc.Source(energy=energy, angle=angle, strength=0.0) for _ in range(strengths.size)] + + # create only one non-zero strength source in voxel (0, 0, 0) + sources[0].strength = 1.0 + + mesh_source = openmc.MeshSource(mesh, sources=sources) + + settings.source = mesh_source + + model = openmc.Model(geometry=geometry, settings=settings) + + # tally the flux on the mesh + mesh_filter = openmc.MeshFilter(mesh) + tally = openmc.Tally() + tally.filters = [mesh_filter] + tally.scores = ['flux'] + + model.tallies = openmc.Tallies([tally]) + + model.run() + From a8a4be353a27dd4bbae100f85e5a2c38e7b7f956 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 22 Sep 2023 09:12:14 -0500 Subject: [PATCH 24/73] Updating check for source type --- openmc/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index df0f347885f..995bb1656bd 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -428,7 +428,7 @@ def mesh(self, m): @sources.setter def sources(self, s): - cv.check_iterable_type('mesh sources', s, openmc.Source) + cv.check_iterable_type('mesh sources', s, openmc.SourceBase) self._sources = s for src in self.sources: if src.space is not None: From b744406efff5e7e339f633949979cd18a9c31d74 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 22 Sep 2023 13:50:30 -0500 Subject: [PATCH 25/73] Correcting detection of a mesh source --- src/source.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 345436509ee..50c170a1dec 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -69,8 +69,7 @@ unique_ptr Source::create(pugi::xml_node node) } // Create compiled source return make_unique(path, parameters); - } else if (check_for_node(node, "mesh") && - get_node_value(node, "type") == "mesh") { + } else if (get_node_value(node, "type") == "mesh") { return make_unique(node); } else { return make_unique(node); From 71a43edd4bfc8a51d096b10b7321ba2e4a5409d3 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 22 Sep 2023 13:51:09 -0500 Subject: [PATCH 26/73] Adding proper test check and explanation of test --- tests/regression_tests/source_mesh/test.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index f716149cfea..f8695d63d9e 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -6,7 +6,6 @@ import openmc.model def test_source_mesh(): - # build a simple void model min, max = -10, 10 box = openmc.model.RectangularParallelepiped(min, max, min, max, min, max, boundary_type='vacuum') @@ -25,11 +24,17 @@ def test_source_mesh(): energy = openmc.stats.Discrete([1.e6], [1.0]) - angle = openmc.stats.Monodirectional((-1, -1, -1)) + # create sources with only one non-zero strength for the source in the + # mesh voxel occupyting the lowest octant. Direct source particles straight + # out of the problem from there. This demonstrates that + # 1) particles are only being sourced within the intented mesh voxel based on source strength + # 2) particles are respecting the angle distributions assigned to each voxel + vec = np.ones((3,)) + vec /= -np.linalg.norm(vec) + angle = openmc.stats.Monodirectional(vec) sources = [openmc.Source(energy=energy, angle=angle, strength=0.0) for _ in range(strengths.size)] - # create only one non-zero strength source in voxel (0, 0, 0) sources[0].strength = 1.0 mesh_source = openmc.MeshSource(mesh, sources=sources) @@ -46,5 +51,12 @@ def test_source_mesh(): model.tallies = openmc.Tallies([tally]) - model.run() + sp_file = model.run() + + with openmc.StatePoint(sp_file) as sp: + tally_out = sp.get_tally(id=tally.id) + + assert tally_out.mean[0] != 0 + assert np.all(tally_out.mean[1:] == 0) + \ No newline at end of file From 119e5f3fd11cc96cacb7d520eb66562a9e12a0ad Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 22 Sep 2023 21:21:19 -0500 Subject: [PATCH 27/73] Improvement to mesh source test --- src/distribution_spatial.cpp | 2 +- tests/regression_tests/source_mesh/test.py | 63 ++++++++++++++-------- 2 files changed, 42 insertions(+), 23 deletions(-) diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 98fc4ae6102..2fcd17b1008 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -272,7 +272,7 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) std::pair MeshSpatial::sample_mesh(uint64_t* seed) const { - // Sample over the CDF defined in initialization above + // Sample the CDF defined in initialization above int32_t elem_idx = elem_idx_dist_.sample(seed); return {elem_idx, mesh()->sample_element(seed, elem_idx)}; } diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index f8695d63d9e..e64e1636627 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -5,8 +5,11 @@ import openmc import openmc.model + def test_source_mesh(): - # build a simple void model + """ + A void model containing a single box + """ min, max = -10, 10 box = openmc.model.RectangularParallelepiped(min, max, min, max, min, max, boundary_type='vacuum') @@ -17,8 +20,10 @@ def test_source_mesh(): settings.batches = 10 settings.run_mode = 'fixed source' + model = openmc.Model(geometry=geometry, settings=settings) + # define a 2 x 2 x 2 mesh - mesh = openmc.RegularMesh.from_domain(-box, (2, 2, 2)) + mesh = openmc.RegularMesh.from_domain(model.geometry, (2, 2, 2)) strengths = np.ones((2, 2, 2)) @@ -29,34 +34,48 @@ def test_source_mesh(): # out of the problem from there. This demonstrates that # 1) particles are only being sourced within the intented mesh voxel based on source strength # 2) particles are respecting the angle distributions assigned to each voxel - vec = np.ones((3,)) - vec /= -np.linalg.norm(vec) - angle = openmc.stats.Monodirectional(vec) - sources = [openmc.Source(energy=energy, angle=angle, strength=0.0) for _ in range(strengths.size)] + x, y, z = mesh.centroids - sources[0].strength = 1.0 + sources = [] + for idx in mesh.indices: - mesh_source = openmc.MeshSource(mesh, sources=sources) + idx = tuple(i-1 for i in idx) + centroid = np.array((x[idx], y[idx], z[idx])) + print(centroid) + vec = np.sign(centroid, dtype=float) + print(vec) + vec /= np.linalg.norm(vec) + angle = openmc.stats.Monodirectional(vec) + sources.append(openmc.Source(energy=energy, angle=angle, strength=0.0)) - settings.source = mesh_source + for idx in range(mesh.num_mesh_cells): + + for s in sources: + s.strength = 0.0 + + sources[idx].strength = 1.0 + + mesh_source = openmc.MeshSource(mesh, sources=sources) + + model.settings.source = mesh_source - model = openmc.Model(geometry=geometry, settings=settings) - # tally the flux on the mesh - mesh_filter = openmc.MeshFilter(mesh) - tally = openmc.Tally() - tally.filters = [mesh_filter] - tally.scores = ['flux'] + # tally the flux on the mesh + mesh_filter = openmc.MeshFilter(mesh) + tally = openmc.Tally() + tally.filters = [mesh_filter] + tally.scores = ['flux'] - model.tallies = openmc.Tallies([tally]) + model.tallies = openmc.Tallies([tally]) - sp_file = model.run() + sp_file = model.run() - with openmc.StatePoint(sp_file) as sp: - tally_out = sp.get_tally(id=tally.id) + with openmc.StatePoint(sp_file) as sp: + tally_out = sp.get_tally(id=tally.id) + mean = tally_out.mean - assert tally_out.mean[0] != 0 - assert np.all(tally_out.mean[1:] == 0) + assert mean[idx] != 0 + mean[idx] = 0 + assert np.all(mean == 0) - \ No newline at end of file From e5bb4137be6293bae862d873b5c66d0ef366b244 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 24 Sep 2023 22:39:51 -0500 Subject: [PATCH 28/73] Write type and strength in the base source class --- openmc/source.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index 995bb1656bd..a023d139b17 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -61,7 +61,8 @@ def populate_xml_element(self, element): XML element containing source data """ - pass + element.set("type", self.type) + element.set("strength", str(self.strength)) def to_xml_element(self) -> ET.Element: """Return XML representation of the source @@ -73,8 +74,6 @@ def to_xml_element(self) -> ET.Element: """ element = ET.Element("source") - element.set("type", self.type) - element.set("strength", str(self.strength)) self.populate_xml_element(element) return element From 0ede245328b7e076fde155327514e6de35193ec8 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 24 Sep 2023 22:40:22 -0500 Subject: [PATCH 29/73] Refactor of source class constructors s.t. they all accept an XML node --- include/openmc/source.h | 10 ++++--- src/source.cpp | 66 +++++++++++++++++++++++++++-------------- 2 files changed, 49 insertions(+), 27 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 41794c3b31b..c0ee9b55540 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -103,12 +103,12 @@ class IndependentSource : public Source { class FileSource : public Source { public: // Constructors - explicit FileSource(std::string path); - explicit FileSource(const vector& 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 sites_; //!< Source sites from a file }; @@ -120,7 +120,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 @@ -131,6 +131,8 @@ 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 compiled_source_; diff --git a/src/source.cpp b/src/source.cpp index 50c170a1dec..c9e4cb14272 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -52,27 +52,28 @@ vector> external_sources; unique_ptr Source::create(pugi::xml_node node) { - if (check_for_node(node, "file")) { - auto path = get_node_value(node, "file", false, true); - if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { - auto sites = mcpl_source_sites(path); - return make_unique(sites); - } else { - return make_unique(path); - } - } else if (check_for_node(node, "library")) { - // Get shared library path and parameters - auto path = get_node_value(node, "library", false, true); - std::string parameters; - if (check_for_node(node, "parameters")) { - parameters = get_node_value(node, "parameters", false, true); + // if the source type is present, use it to determine the type + // of object to create + if (check_for_node(node, "type")) { + std::string source_type = get_node_value(node, "type"); + if (source_type == "independent") { + return make_unique(node); + } else if (source_type == "file") { + return make_unique(node); + } else if (source_type == "compiled") { + return make_unique(node); + } else if (source_type == "mesh") { + return make_unique(node); } - // Create compiled source - return make_unique(path, parameters); - } else if (get_node_value(node, "type") == "mesh") { - return make_unique(node); } else { - return make_unique(node); + // support legacy source format + if (check_for_node(node, "file")) { + return make_unique(node); + } else if (check_for_node(node, "library")) { + return make_unique(node); + } else { + return make_unique(node); + } } } @@ -297,8 +298,20 @@ SourceSite IndependentSource::sample(uint64_t* seed) const //============================================================================== // FileSource implementation //============================================================================== +FileSource::FileSource(pugi::xml_node node) { + auto path = get_node_value(node, "file", false, true); + if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { + sites_ = mcpl_source_sites(path); + } else { + this->load_sites_from_file(path); + } +} + +FileSource::FileSource(const std::string& path) { + load_sites_from_file(path); +} -FileSource::FileSource(std::string path) +void FileSource::load_sites_from_file(const std::string& path) { // Check if source file exists if (!file_exists(path)) { @@ -335,10 +348,17 @@ SourceSite FileSource::sample(uint64_t* seed) const //============================================================================== // CompiledSourceWrapper implementation //============================================================================== +CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) { + // Get shared library path and parameters + auto path = get_node_value(node, "library", false, true); + std::string parameters; + if (check_for_node(node, "parameters")) { + parameters = get_node_value(node, "parameters", false, true); + } + setup(path, parameters); +} -CompiledSourceWrapper::CompiledSourceWrapper( - std::string path, std::string parameters) -{ +void CompiledSourceWrapper::setup(const std::string& path, const std::string& parameters) { #ifdef HAS_DYNAMIC_LINKING // Open the library shared_library_ = dlopen(path.c_str(), RTLD_LAZY); From bf82136953c4113776b019e0a7ab876187d1940d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 25 Sep 2023 09:49:17 -0500 Subject: [PATCH 30/73] More informative error messages --- include/openmc/mesh.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 93fee30beb0..9e04f61753d 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -252,7 +252,8 @@ class StructuredMesh : public Mesh { //! \param[in] i Direction index virtual double positive_grid_boundary(const MeshIndex& ijk, int i) const { - fatal_error("Shouldn't be here"); + 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 @@ -261,7 +262,8 @@ class StructuredMesh : public Mesh { //! \param[in] i Direction index virtual double negative_grid_boundary(const MeshIndex& ijk, int i) const { - fatal_error("Shouldn't be here"); + 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 From ea083189d9efa1100e3b9a7d45ca3517c3a5e35a Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 25 Sep 2023 13:04:15 -0500 Subject: [PATCH 31/73] Altering some aspects of the MeshSource class - MeshSource.strength is represented by the sum of all provided element sources. It can be scaled by a floating point factor. It has been removed as a parameter of MeshSource as a result. - Element sources are written using the index ordering of the mesh for safety. - Enforcing that MeshBase.sources is a numpy.ndarray of a shape consistent with that of the mesh. --- openmc/source.py | 83 +++++++++++++++++++++++++++++++++++++----------- 1 file changed, 65 insertions(+), 18 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index a023d139b17..649d3f8de37 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -18,7 +18,7 @@ from openmc.stats.multivariate import UnitSphere, Spatial from openmc.stats.univariate import Univariate from ._xml import get_text -from .mesh import MeshBase +from .mesh import MeshBase, StructuredMesh, UnstructuredMesh class SourceBase(ABC): @@ -394,19 +394,36 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': return source class MeshSource(): - """ + """A mesh-based source in which random positions are uniformly sampled + within mesh elements and each element can have independent angle, energy, + and time distributions. The element sampled is chosen based on the relative + strengths of the sources applied to the elements. The strength of the mesh + source as a whole is the sum of all source strengths applied to the elements. + + Parameters + ---------- + mesh : openmc.MeshBase + The mesh on which source sites will be generated. + sources : iterable of openmc.SourceBase + Sources for each element in the mesh. If spatial distributions are set + on any of the source objects, they will be ignored during source site + sampling. + + Attributes + ---------- + mesh : openmc.MeshBase + The mesh on which source sites will be generated. + sources : numpy.ndarray or iterable of openmc.SourceBase + The set of sources to apply to each element. The shape of this array must match + the shape of the mesh with and exception in the case of unstructured mesh, which + allows for application of 1-D array or iterable. """ def __init__(self, mesh: openmc.MeshBase, - strength: Optional[float] = 1.0, - sources: Optional[Iterable] = None): + sources: Optional[Iterable]): self.mesh = mesh - self.strength = strength - self._sources = cv.CheckedList(openmc.Source, 'sources') - - if sources is not None: - self.sources = sources + self.sources = sources @property def mesh(self): @@ -414,7 +431,9 @@ def mesh(self): @property def strength(self): - return self._strength + if self.sources is None: + return 0.0 + return sum(s.strength for s in self.sources.flat) @property def sources(self): @@ -427,9 +446,20 @@ def mesh(self, m): @sources.setter def sources(self, s): - cv.check_iterable_type('mesh sources', s, openmc.SourceBase) + cv.check_iterable_type('mesh sources', s, openmc.SourceBase, max_depth=3) + + s = np.asarray(s) + + if isinstance(self.mesh, StructuredMesh) and s.shape != self.mesh.dimension: + raise ValueError('The shape of the source array' \ + f'({s.shape}) does not match the' \ + f'dimensions of the structured mesh ({self.mesh.dimension})') + elif isinstance(self.mesh, UnstructuredMesh): + if len(s.shape) > 1: + raise ValueError('Sources must be a 1-D array for unstructured mesh') + self._sources = s - for src in self.sources: + for src in self.sources.flat: if src.space is not None: warnings.warn('Some sources on the mesh have ' 'spatial distributions that will ' @@ -437,10 +467,23 @@ def sources(self, s): break @strength.setter - def strength(self, strength): - cv.check_type('source strength', strength, Real) - cv.check_greater_than('source strength', strength, 0.0, True) - self._strength = strength + def strength(self, val): + cv.check_type('mesh source strength', val, Real) + self._set_total_strength(val) + + def _set_total_strength(self, new_strength): + """Scales the element source strengths based on a desired + total mesh strength. + """ + current_strength = self.strength + + for s in self.sources.flat: + s.strength *= new_strength / current_strength + + def normalize_source_strengths(self): + """Update all element source strengths such that they sum to 1.0. + """ + self.set_total_strength(1.0) def to_xml_element(self): """Return XML representation of the mesh source @@ -456,11 +499,15 @@ def to_xml_element(self): element.set('strength', str(self.strength)) # build a mesh spatial distribution - mesh_dist = openmc.MeshSpatial(self.mesh, [s.strength for s in self.sources]) + strengths = [self.sources[i-1, j-1, k-1].strength for i, j, k in self.mesh.indices] + mesh_dist = openmc.MeshSpatial(self.mesh, strengths) mesh_dist.to_xml_element(element) - for source in self.sources: + # write in the order of mesh indices + for idx in self.mesh.indices: + idx = tuple(i - 1 for i in idx) + source = self.sources[idx] src_element = ET.SubElement(element, 'source') source.populate_xml_element(src_element) From 66fdcf864fe37a45c6447e1595051a01938b7965 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 25 Sep 2023 13:09:43 -0500 Subject: [PATCH 32/73] Adding IO format documentation for the mesh source type --- docs/source/io_formats/settings.rst | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 03abe9a3c4f..4f57ccf98ec 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -491,7 +491,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 +665,9 @@ attributes/sub-elements: *Default*: false + :source: + For mesh sources, a set of sources that will be applied to each element. + .. _univariate: Univariate Probability Distributions From eecd0ec6bc1cc27e7dcfe44dd17e5449f8ddc398 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 25 Sep 2023 13:19:43 -0500 Subject: [PATCH 33/73] Improving mesh source test by using indices ordering for setting sources --- tests/regression_tests/source_mesh/test.py | 68 ++++++++++++---------- 1 file changed, 37 insertions(+), 31 deletions(-) diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index e64e1636627..65561823dd0 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -25,8 +25,6 @@ def test_source_mesh(): # define a 2 x 2 x 2 mesh mesh = openmc.RegularMesh.from_domain(model.geometry, (2, 2, 2)) - strengths = np.ones((2, 2, 2)) - energy = openmc.stats.Discrete([1.e6], [1.0]) # create sources with only one non-zero strength for the source in the @@ -37,45 +35,53 @@ def test_source_mesh(): x, y, z = mesh.centroids - sources = [] - for idx in mesh.indices: + sources = np.ndarray((2,2,2), dtype=openmc.SourceBase) + for i, j, k in mesh.indices: + # mesh.indices is currently one-indexed, adjust for Python arrays + idx = (i-1, j-1, k-1) - idx = tuple(i-1 for i in idx) + # get the centroid of the ijk mesh element, set a particle source + # vector based on the centroid = np.array((x[idx], y[idx], z[idx])) - print(centroid) vec = np.sign(centroid, dtype=float) - print(vec) vec /= np.linalg.norm(vec) angle = openmc.stats.Monodirectional(vec) - sources.append(openmc.Source(energy=energy, angle=angle, strength=0.0)) - - for idx in range(mesh.num_mesh_cells): - - for s in sources: - s.strength = 0.0 - - sources[idx].strength = 1.0 - - mesh_source = openmc.MeshSource(mesh, sources=sources) - - model.settings.source = mesh_source - - - # tally the flux on the mesh - mesh_filter = openmc.MeshFilter(mesh) - tally = openmc.Tally() - tally.filters = [mesh_filter] - tally.scores = ['flux'] - - model.tallies = openmc.Tallies([tally]) + sources[idx] = openmc.Source(energy=energy, angle=angle, strength=0.0) + + # create and apply the mesh source + mesh_source = openmc.MeshSource(mesh, sources) + model.settings.source = mesh_source + + # tally the flux on the mesh + mesh_filter = openmc.MeshFilter(mesh) + tally = openmc.Tally() + tally.filters = [mesh_filter] + tally.scores = ['flux'] + + model.tallies = openmc.Tallies([tally]) + + # for each element, set a single-non zero source with particles + # traveling out of the mesh (and geometry) w/o crossing any other + # mesh elements + for i, j, k in mesh.indices: + ijk = (i-1, j-1, k-1) + # zero-out all source strengths and set the strength + # on the element of interest + mesh_source.strength = 0.0 + mesh_source.sources[ijk].strength = 1.0 sp_file = model.run() with openmc.StatePoint(sp_file) as sp: tally_out = sp.get_tally(id=tally.id) - mean = tally_out.mean + mean = tally_out.get_reshaped_data(expand_dims=True).squeeze() - assert mean[idx] != 0 - mean[idx] = 0 + assert mean[ijk] != 0 + mean[ijk] = 0 assert np.all(mean == 0) + # check strength adjustment methods + assert mesh_source.strength == 1.0 + mesh_source.strength = 100.0 + assert mesh_source.strength == 100.0 + From 609f0a1dd6ef725581680f58a623f42cf5ffa353 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 27 Sep 2023 11:05:09 -0500 Subject: [PATCH 34/73] Handling case where all source strengths are zero --- openmc/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index 649d3f8de37..70144faecd0 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -475,7 +475,7 @@ def _set_total_strength(self, new_strength): """Scales the element source strengths based on a desired total mesh strength. """ - current_strength = self.strength + current_strength = self.strength if self.strength != 0.0 else 1.0 for s in self.sources.flat: s.strength *= new_strength / current_strength From 8f11fb309b9fa10b2b6735427320ded137b1149c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 27 Sep 2023 12:48:47 -0500 Subject: [PATCH 35/73] Allowing construction of DiscreteIndex with gsl::span. Adding new constructor for MeshSpatial. --- include/openmc/distribution.h | 7 ++--- include/openmc/distribution_spatial.h | 4 +++ src/distribution.cpp | 14 +++++----- src/distribution_spatial.cpp | 37 +++++++++++++++++---------- 4 files changed, 37 insertions(+), 25 deletions(-) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 65c68ad6457..e70c803de2f 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -7,6 +7,7 @@ #include // for size_t #include "pugixml.hpp" +#include #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 p); - void assign(const double* p, int n); + void assign(gsl::span 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 diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index df314df3a20..bd10a299605 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -104,6 +104,7 @@ class SphericalIndependent : public SpatialDistribution { class MeshSpatial : public SpatialDistribution { public: explicit MeshSpatial(pugi::xml_node node); + explicit MeshSpatial(int32_t mesh_id, gsl::span strengths); //! Sample a position from the distribution //! \param seed Pseudorandom number seed pointer @@ -115,6 +116,9 @@ class MeshSpatial : public SpatialDistribution { //! \return Sampled element index and position within that element std::pair 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(); } diff --git a/src/distribution.cpp b/src/distribution.cpp index 9c76147ee27..1d85b922f94 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -27,17 +27,15 @@ DiscreteIndex::DiscreteIndex(pugi::xml_node node) auto params = get_node_array(node, "parameters"); std::size_t n = params.size() / 2; - assign(params.data() + n, n); + assign({params.data() + n, n}); } -DiscreteIndex::DiscreteIndex(const double* p, int n) -{ - assign(p, n); +DiscreteIndex::DiscreteIndex(gsl::span p) { + assign(p); } -void DiscreteIndex::assign(const double* p, int n) -{ - prob_.assign(p, p + n); +void DiscreteIndex::assign(gsl::span p) { + prob_.assign(p.begin(), p.end()); this->init_alias(); } @@ -126,7 +124,7 @@ Discrete::Discrete(pugi::xml_node node) : di_(node) x_.assign(params.begin(), params.begin() + n); } -Discrete::Discrete(const double* x, const double* p, int n) : di_(p, n) +Discrete::Discrete(const double* x, const double* p, size_t n) : di_({p, n}) { x_.assign(x, x + n); diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 2fcd17b1008..7c2d513d840 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -233,19 +233,9 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) const auto mesh_ptr = model::meshes.at(mesh_idx_).get(); - const auto umesh_ptr = - dynamic_cast(model::meshes.at(mesh_idx_).get()); - if (umesh_ptr) { - // ensure that the unstructured mesh contains only linear tets - for (int bin = 0; bin < mesh_ptr->n_bins(); bin++) { - if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) { - fatal_error( - "Mesh specified for source must contain only linear tetrahedra."); - } - } - } + check_element_types(); - int32_t n_bins = this->n_sources(); + size_t n_bins = this->n_sources(); std::vector strengths(n_bins, 1.0); // Create cdfs for sampling for an element over a mesh @@ -263,11 +253,30 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) if (get_node_value_bool(node, "volume_normalized")) { for (int i = 0; i < n_bins; i++) { - strengths[i] *= mesh()->volume(i); + strengths[i] *= this->mesh()->volume(i); } } - elem_idx_dist_.assign(strengths.data(), n_bins); + elem_idx_dist_.assign(strengths); +} + +MeshSpatial::MeshSpatial(int32_t mesh_idx, gsl::span strengths) : mesh_idx_(mesh_idx) { + check_element_types(); + elem_idx_dist_.assign(strengths); +} + +void MeshSpatial::check_element_types() const { + const auto umesh_ptr = + dynamic_cast(this->mesh()); + if (umesh_ptr) { + // ensure that the unstructured mesh contains only linear tets + for (int bin = 0; bin < umesh_ptr->n_bins(); bin++) { + if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) { + fatal_error( + "Mesh specified for source must contain only linear tetrahedra."); + } + } + } } std::pair MeshSpatial::sample_mesh(uint64_t* seed) const From c481fdfdc9ed2e86a39cf9724775cb6d34807d80 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 27 Sep 2023 12:49:41 -0500 Subject: [PATCH 36/73] MeshSpatial object is not only on the C++ side. --- include/openmc/source.h | 2 +- openmc/source.py | 30 +++++++++--------------------- src/source.cpp | 28 +++++++++++----------------- 3 files changed, 21 insertions(+), 39 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index c0ee9b55540..ee8781d816d 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -167,7 +167,7 @@ class MeshSource : public Source { private: // Data members ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted - double strength_ {1.0}; //!< Source strength + double strength_ {1.0}; //!< Total source strength unique_ptr space_; //!< Mesh spatial vector sources_; //!< Source distributions }; diff --git a/openmc/source.py b/openmc/source.py index 70144faecd0..5d66d9400e6 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -393,7 +393,7 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': return source -class MeshSource(): +class MeshSource(SourceBase): """A mesh-based source in which random positions are uniformly sampled within mesh elements and each element can have independent angle, energy, and time distributions. The element sampled is chosen based on the relative @@ -425,6 +425,10 @@ def __init__(self, self.mesh = mesh self.sources = sources + @property + def type(self): + return "mesh" + @property def mesh(self): return self._mesh @@ -485,34 +489,18 @@ def normalize_source_strengths(self): """ self.set_total_strength(1.0) - def to_xml_element(self): - """Return XML representation of the mesh source - - Returns - ------- - element : xml.etree.ElementTree.Element - XML element containing source data - - """ - element = ET.Element('source') - element.set('type', 'mesh') - element.set('strength', str(self.strength)) - - # build a mesh spatial distribution - strengths = [self.sources[i-1, j-1, k-1].strength for i, j, k in self.mesh.indices] - mesh_dist = openmc.MeshSpatial(self.mesh, strengths) + def populate_xml_element(self, elem): + super().populate_xml_element(elem) - mesh_dist.to_xml_element(element) + elem.set("mesh", str(self.mesh.id)) # write in the order of mesh indices for idx in self.mesh.indices: idx = tuple(i - 1 for i in idx) source = self.sources[idx] - src_element = ET.SubElement(element, 'source') + src_element = ET.SubElement(elem, 'source') source.populate_xml_element(src_element) - return element - @classmethod def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.MeshSource: """ diff --git a/src/source.cpp b/src/source.cpp index c9e4cb14272..ceb60541f36 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -424,33 +424,27 @@ MeshSource::MeshSource(pugi::xml_node node) } } - // read mesh distribution - if (check_for_node(node, "space")) { - space_ = std::make_unique(node.child("space")); - } else { - fatal_error("No spatial distribution found for mesh source."); - } - - MeshSpatial* mesh_spatial = dynamic_cast(space_.get()); - if (!mesh_spatial) - fatal_error("Failed to recast mesh spatital distribution for mesh source."); + int32_t mesh_id = stoi(get_node_value(node, "mesh")); + int32_t mesh_idx = model::mesh_map.at(mesh_id); + const auto& mesh = model::meshes[mesh_idx]; - strength_ = mesh_spatial->total_strength(); - - int32_t n_mesh_sources = mesh_spatial->n_sources(); - - // read all source distributions + std::vector strengths; + // read all source distributions and populate strengths vector for MeshSpatial object for (auto source_node : node.children("source")) { sources_.emplace_back(IndependentSource(source_node)); + strengths.push_back(sources_.back().strength()); } // the number of source distributions should either be one or equal to the // number of mesh elements - if (sources_.size() > 1 && sources_.size() != n_mesh_sources) { + if (sources_.size() > 1 && sources_.size() != mesh->n_bins()) { fatal_error(fmt::format("Incorrect number of source distributions ({}) for " "mesh source with {} elements.", - sources_.size(), n_mesh_sources)); + sources_.size(), mesh->n_bins())); } + + space_ = std::make_unique(mesh_idx, strengths); + strength_ = std::accumulate(strengths.begin(), strengths.end(), 0.0); } SourceSite MeshSource::sample(uint64_t* seed) const From 8c9305e02c08ce2b596918c5540757b7d070d625 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 27 Sep 2023 13:08:24 -0500 Subject: [PATCH 37/73] Expand test to include cylindrical mesh --- tests/regression_tests/source_mesh/test.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index 65561823dd0..2ed768ed681 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -5,8 +5,8 @@ import openmc import openmc.model - -def test_source_mesh(): +@pytest.mark.parametrize('mesh_type', ('rectangular', 'cylindrical')) +def test_source_mesh(mesh_type): """ A void model containing a single box """ @@ -23,7 +23,10 @@ def test_source_mesh(): model = openmc.Model(geometry=geometry, settings=settings) # define a 2 x 2 x 2 mesh - mesh = openmc.RegularMesh.from_domain(model.geometry, (2, 2, 2)) + if mesh_type == 'rectangular': + mesh = openmc.RegularMesh.from_domain(model.geometry, (2, 2, 2)) + elif mesh_type == 'cylindrical': + mesh = openmc.CylindricalMesh.from_domain(model.geometry, (1, 4, 2)) energy = openmc.stats.Discrete([1.e6], [1.0]) @@ -32,10 +35,11 @@ def test_source_mesh(): # out of the problem from there. This demonstrates that # 1) particles are only being sourced within the intented mesh voxel based on source strength # 2) particles are respecting the angle distributions assigned to each voxel - x, y, z = mesh.centroids + if mesh_type == 'cylindrical': + x, y, z = mesh._convert_to_cartesian(mesh.centroids, mesh.origin) - sources = np.ndarray((2,2,2), dtype=openmc.SourceBase) + sources = np.ndarray(mesh.dimension, dtype=openmc.SourceBase) for i, j, k in mesh.indices: # mesh.indices is currently one-indexed, adjust for Python arrays idx = (i-1, j-1, k-1) @@ -65,6 +69,7 @@ def test_source_mesh(): # mesh elements for i, j, k in mesh.indices: ijk = (i-1, j-1, k-1) + print(ijk) # zero-out all source strengths and set the strength # on the element of interest mesh_source.strength = 0.0 @@ -74,7 +79,10 @@ def test_source_mesh(): with openmc.StatePoint(sp_file) as sp: tally_out = sp.get_tally(id=tally.id) - mean = tally_out.get_reshaped_data(expand_dims=True).squeeze() + mean = tally_out.get_reshaped_data(expand_dims=True) + + # remove nuclides and scores axes + mean = mean[..., 0, 0] assert mean[ijk] != 0 mean[ijk] = 0 From d73c002ae7221731c92cfd7babd4e158a21cb0fa Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 28 Sep 2023 13:13:24 -0500 Subject: [PATCH 38/73] Use total strength from DiscreteIndex data member --- include/openmc/source.h | 3 +-- src/source.cpp | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index ee8781d816d..7a2b70241ff 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -156,7 +156,7 @@ class MeshSource : public Source { // Properties ParticleType particle_type() const { return particle_; } - double strength() const override { return strength_; } + double strength() const { return space_->total_strength(); } //!< Total source strength // Accessors const IndependentSource& source(int32_t i) const @@ -167,7 +167,6 @@ class MeshSource : public Source { private: // Data members ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted - double strength_ {1.0}; //!< Total source strength unique_ptr space_; //!< Mesh spatial vector sources_; //!< Source distributions }; diff --git a/src/source.cpp b/src/source.cpp index ceb60541f36..addf93f8ebd 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -444,7 +444,6 @@ MeshSource::MeshSource(pugi::xml_node node) } space_ = std::make_unique(mesh_idx, strengths); - strength_ = std::accumulate(strengths.begin(), strengths.end(), 0.0); } SourceSite MeshSource::sample(uint64_t* seed) const From 5ab09b74b4bc29c8b009b412f4381490e6f39f22 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Oct 2023 23:35:50 -0500 Subject: [PATCH 39/73] Fix in StructuredMesh::sample_element --- src/mesh.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 094b15f6784..264a7d5852d 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -172,11 +172,11 @@ Position StructuredMesh::sample_element( double x_min = negative_grid_boundary(ijk, 0); double x_max = positive_grid_boundary(ijk, 0); - double y_min = negative_grid_boundary(ijk, 1); - double y_max = positive_grid_boundary(ijk, 1); + double y_min = (n_dimension_ >= 2) ? negative_grid_boundary(ijk, 1) : 0.0; + double y_max = (n_dimension_ >= 2) ? positive_grid_boundary(ijk, 1) : 0.0; - double z_min = negative_grid_boundary(ijk, 2); - double z_max = positive_grid_boundary(ijk, 2); + double z_min = (n_dimension_ == 3) ? negative_grid_boundary(ijk, 2) : 0.0; + double z_max = (n_dimension_ == 3) ? positive_grid_boundary(ijk, 2) : 0.0; return {x_min + (x_max - x_min) * prn(seed), y_min + (y_max - y_min) * prn(seed), z_min + (z_max - z_min) * prn(seed)}; From 52d0431e1888af62f730fffef78349a38cf0e441 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 18 Oct 2023 21:40:29 -0500 Subject: [PATCH 40/73] Patching up some I/O methods --- openmc/settings.py | 6 ++++++ openmc/source.py | 15 ++++++--------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 5a45d454eec..64ca0d6fa18 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1878,6 +1878,12 @@ def from_xml_element(cls, elem, meshes=None): Settings object """ + settings_meshes = _read_meshes(elem) + if meshes is None: + meshes = settings_meshes + else: + meshes.update(settings_meshes) + settings = cls() settings._eigenvalue_from_xml_element(elem) settings._run_mode_from_xml_element(elem) diff --git a/openmc/source.py b/openmc/source.py index 5d66d9400e6..7cce3155cdb 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -114,7 +114,7 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.SourceBase: elif source_type == 'file': return FileSource.from_xml_element(elem) elif source_type == 'mesh': - return MeshSource.from_xml_element(elem) + return MeshSource.from_xml_element(elem, meshes) else: raise ValueError(f'Source type {source_type} is not recognized') @@ -421,7 +421,6 @@ class MeshSource(SourceBase): def __init__(self, mesh: openmc.MeshBase, sources: Optional[Iterable]): - self.mesh = mesh self.sources = sources @@ -456,7 +455,7 @@ def sources(self, s): if isinstance(self.mesh, StructuredMesh) and s.shape != self.mesh.dimension: raise ValueError('The shape of the source array' \ - f'({s.shape}) does not match the' \ + f'({s.shape}) does not match the ' \ f'dimensions of the structured mesh ({self.mesh.dimension})') elif isinstance(self.mesh, UnstructuredMesh): if len(s.shape) > 1: @@ -519,15 +518,13 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.MeshSource: openmc.MeshSource MeshSource generated from the XML element """ - mesh_id = int(elem.get('mesh_id')) + mesh_id = int(get_text(elem, 'mesh')) mesh = meshes[mesh_id] - strength = get_text(elem, 'strength') - - sources = [Source.from_xml_element(e) for e in elem.iter('source')] - - return cls(mesh, strength, sources) + sources = [SourceBase.from_xml_element(e) for e in elem.iter('source') if e != elem] + sources = np.asarray(sources).reshape(mesh.dimension, order='F') + return cls(mesh, sources) def Source(*args, **kwargs): From 807fb3f4df5d9a07b66cf38e79e2245264a4c932 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 18 Oct 2023 21:40:49 -0500 Subject: [PATCH 41/73] Test updates --- tests/regression_tests/source_mesh/test.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index 2ed768ed681..cd77ff50a2a 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -35,11 +35,8 @@ def test_source_mesh(mesh_type): # out of the problem from there. This demonstrates that # 1) particles are only being sourced within the intented mesh voxel based on source strength # 2) particles are respecting the angle distributions assigned to each voxel - x, y, z = mesh.centroids - if mesh_type == 'cylindrical': - x, y, z = mesh._convert_to_cartesian(mesh.centroids, mesh.origin) - sources = np.ndarray(mesh.dimension, dtype=openmc.SourceBase) + x, y, z = mesh.centroids for i, j, k in mesh.indices: # mesh.indices is currently one-indexed, adjust for Python arrays idx = (i-1, j-1, k-1) @@ -83,10 +80,9 @@ def test_source_mesh(mesh_type): # remove nuclides and scores axes mean = mean[..., 0, 0] - assert mean[ijk] != 0 mean[ijk] = 0 - assert np.all(mean == 0) + assert np.all(mean == 0), f'Failed on index {ijk} with centroid {mesh.centroids[(..., *ijk)]}' # check strength adjustment methods assert mesh_source.strength == 1.0 From a11949e007077a37452bdf09c4dcaf4b80a9cc0c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 24 Oct 2023 15:38:04 -0500 Subject: [PATCH 42/73] rm print --- tests/regression_tests/source_mesh/test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index cd77ff50a2a..bed0937dc74 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -66,7 +66,6 @@ def test_source_mesh(mesh_type): # mesh elements for i, j, k in mesh.indices: ijk = (i-1, j-1, k-1) - print(ijk) # zero-out all source strengths and set the strength # on the element of interest mesh_source.strength = 0.0 From 2e58d77be193b32cecc30b22d099b676297619af Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 25 Oct 2023 11:39:18 -0500 Subject: [PATCH 43/73] Removing mesh source test from unit test dir --- tests/unit_tests/test_source_mesh.py | 49 ---------------------------- 1 file changed, 49 deletions(-) diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index 925bba1249d..fd8b04148e7 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -201,52 +201,3 @@ def test_roundtrip(run_in_tmpdir, model, request): assert space_in.mesh.id == space_out.mesh.id assert space_in.volume_normalized == space_out.volume_normalized -def test_mesh_source(model): - - - tally = openmc.Tally() - all_cells = list(model.geometry.get_all_cells().values()) - cell_filter = openmc.CellFilter(all_cells) - # energy_filter = openmc.EnergyFilter([0.0, 1e4]) - tally.filters = [cell_filter] - tally.scores = ['flux'] - - model.tallies = [tally] - - - model.settings.particles = 10000 - - - # define a source for the upper right voxel - source = openmc.Source() - source.space = openmc.stats.Box((0, 0, 0), (2, 2, 2)) - model.settings.source = source - - - sp_filename = model.run() - with openmc.StatePoint(sp_filename) as sp: - tally_mean_box = sp.tallies[tally.id].mean - - # change from a standard source to a mesh source with - mesh = openmc.RegularMesh() - mesh.lower_left = (0, 0, 0) - mesh.upper_right = (2, 2, 2) - mesh.dimension = (1, 1, 1) - - ms = openmc.MeshSource(mesh) - - angle = openmc.stats.Isotropic() - src1 = openmc.Source(angle=angle) - - ms.sources = [src1] - - model.settings.source = ms - - model.settings.export_to_xml() - - sp_filename = model.run() - - with openmc.StatePoint(sp_filename) as sp: - tally_mean_mesh = sp.tallies[tally.id].mean - - # np.testing.assert_allclose(tally_mean_box, tally_mean_mesh) From 1879e1603e949e67d4e5e54a939f3f64d21d5894 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Nov 2023 13:22:11 -0600 Subject: [PATCH 44/73] Cleaning up mesh indexing after reordering PR --- tests/regression_tests/source_mesh/test.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index bed0937dc74..4b02a80dd88 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -36,18 +36,19 @@ def test_source_mesh(mesh_type): # 1) particles are only being sourced within the intented mesh voxel based on source strength # 2) particles are respecting the angle distributions assigned to each voxel sources = np.ndarray(mesh.dimension, dtype=openmc.SourceBase) - x, y, z = mesh.centroids + centroids = mesh.centroids + x, y, z = np.swapaxes(mesh.centroids, -1, 0) for i, j, k in mesh.indices: # mesh.indices is currently one-indexed, adjust for Python arrays - idx = (i-1, j-1, k-1) + ijk = (i-1, j-1, k-1) # get the centroid of the ijk mesh element, set a particle source # vector based on the - centroid = np.array((x[idx], y[idx], z[idx])) + centroid = centroids[ijk] vec = np.sign(centroid, dtype=float) vec /= np.linalg.norm(vec) angle = openmc.stats.Monodirectional(vec) - sources[idx] = openmc.Source(energy=energy, angle=angle, strength=0.0) + sources[ijk] = openmc.Source(energy=energy, angle=angle, strength=0.0) # create and apply the mesh source mesh_source = openmc.MeshSource(mesh, sources) @@ -79,9 +80,11 @@ def test_source_mesh(mesh_type): # remove nuclides and scores axes mean = mean[..., 0, 0] + # the mesh elment with a non-zero source strength should have a value assert mean[ijk] != 0 + # all other values should be zero mean[ijk] = 0 - assert np.all(mean == 0), f'Failed on index {ijk} with centroid {mesh.centroids[(..., *ijk)]}' + assert np.all(mean == 0), f'Failed on index {ijk} with centroid {mesh.centroids[ijk]}' # check strength adjustment methods assert mesh_source.strength == 1.0 From a0e2b4777507b2f4bce764fb53272b3aebe21559 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Nov 2023 13:30:49 -0600 Subject: [PATCH 45/73] Add some roundtrip testing --- tests/regression_tests/source_mesh/test.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index 4b02a80dd88..fddd7e58547 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -86,8 +86,17 @@ def test_source_mesh(mesh_type): mean[ijk] = 0 assert np.all(mean == 0), f'Failed on index {ijk} with centroid {mesh.centroids[ijk]}' + # test roundtrip + xml_model = openmc.Model.from_model_xml() + xml_source = xml_model.settings.source[0] + assert isinstance(xml_source, openmc.MeshSource) + assert xml_source.strength == 1.0 + assert isinstance(xml_source.mesh, type(mesh_source.mesh)) + assert xml_source.mesh.dimension == mesh_source.mesh.dimension + assert xml_source.mesh.id == mesh_source.mesh.id + assert len(xml_source.sources) == len(mesh_source.sources) + # check strength adjustment methods assert mesh_source.strength == 1.0 mesh_source.strength = 100.0 assert mesh_source.strength == 100.0 - From 7f220f869591ba7a05afd9d8285fa37973cc71cd Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Nov 2023 14:04:06 -0600 Subject: [PATCH 46/73] Updating format check triggers --- .github/workflows/format-check.yml | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/.github/workflows/format-check.yml b/.github/workflows/format-check.yml index 4914d5e5ea9..36ccece3cc5 100644 --- a/.github/workflows/format-check.yml +++ b/.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: From 827b75570b68c5a6aa8144972bdd7c0f5cba6bee Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Nov 2023 14:06:34 -0600 Subject: [PATCH 47/73] Applying formatter --- include/openmc/mesh.h | 8 ++++++-- include/openmc/source.h | 8 ++++++-- src/distribution.cpp | 6 ++++-- src/distribution_spatial.cpp | 10 ++++++---- src/source.cpp | 32 +++++++++++++++++++------------- 5 files changed, 41 insertions(+), 23 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 9e04f61753d..75de700ff8a 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -252,7 +252,9 @@ class StructuredMesh : public Mesh { //! \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()); + auto msg = + fmt::format("Attempting to call positive_grid_boundary on a {} mesh.", + get_mesh_type()); fatal_error(msg); }; @@ -262,7 +264,9 @@ class StructuredMesh : public Mesh { //! \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()); + auto msg = + fmt::format("Attempting to call negative_grid_boundary on a {} mesh.", + get_mesh_type()); fatal_error(msg); }; diff --git a/include/openmc/source.h b/include/openmc/source.h index 7a2b70241ff..9e64117241b 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -108,7 +108,8 @@ class FileSource : public Source { // Methods SourceSite sample(uint64_t* seed) const override; - void load_sites_from_file(const std::string& path); //!< Load source sites from file + void load_sites_from_file( + const std::string& path); //!< Load source sites from file private: vector sites_; //!< Source sites from a file }; @@ -156,7 +157,10 @@ class MeshSource : public Source { // Properties ParticleType particle_type() const { return particle_; } - double strength() const { return space_->total_strength(); } //!< Total source strength + double strength() const + { + return space_->total_strength(); + } //!< Total source strength // Accessors const IndependentSource& source(int32_t i) const diff --git a/src/distribution.cpp b/src/distribution.cpp index 1d85b922f94..3026630b335 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -30,11 +30,13 @@ DiscreteIndex::DiscreteIndex(pugi::xml_node node) assign({params.data() + n, n}); } -DiscreteIndex::DiscreteIndex(gsl::span p) { +DiscreteIndex::DiscreteIndex(gsl::span p) +{ assign(p); } -void DiscreteIndex::assign(gsl::span p) { +void DiscreteIndex::assign(gsl::span p) +{ prob_.assign(p.begin(), p.end()); this->init_alias(); diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 7c2d513d840..bbbf5b8f5b6 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -260,14 +260,16 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) elem_idx_dist_.assign(strengths); } -MeshSpatial::MeshSpatial(int32_t mesh_idx, gsl::span strengths) : mesh_idx_(mesh_idx) { +MeshSpatial::MeshSpatial(int32_t mesh_idx, gsl::span strengths) + : mesh_idx_(mesh_idx) +{ check_element_types(); elem_idx_dist_.assign(strengths); } -void MeshSpatial::check_element_types() const { - const auto umesh_ptr = - dynamic_cast(this->mesh()); +void MeshSpatial::check_element_types() const +{ + const auto umesh_ptr = dynamic_cast(this->mesh()); if (umesh_ptr) { // ensure that the unstructured mesh contains only linear tets for (int bin = 0; bin < umesh_ptr->n_bins(); bin++) { diff --git a/src/source.cpp b/src/source.cpp index addf93f8ebd..33f4a08f14d 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -59,9 +59,9 @@ unique_ptr Source::create(pugi::xml_node node) if (source_type == "independent") { return make_unique(node); } else if (source_type == "file") { - return make_unique(node); + return make_unique(node); } else if (source_type == "compiled") { - return make_unique(node); + return make_unique(node); } else if (source_type == "mesh") { return make_unique(node); } @@ -298,16 +298,18 @@ SourceSite IndependentSource::sample(uint64_t* seed) const //============================================================================== // FileSource implementation //============================================================================== -FileSource::FileSource(pugi::xml_node node) { +FileSource::FileSource(pugi::xml_node node) +{ auto path = get_node_value(node, "file", false, true); - if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { - sites_ = mcpl_source_sites(path); - } else { - this->load_sites_from_file(path); - } + if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { + sites_ = mcpl_source_sites(path); + } else { + this->load_sites_from_file(path); + } } -FileSource::FileSource(const std::string& path) { +FileSource::FileSource(const std::string& path) +{ load_sites_from_file(path); } @@ -348,8 +350,9 @@ SourceSite FileSource::sample(uint64_t* seed) const //============================================================================== // CompiledSourceWrapper implementation //============================================================================== -CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) { - // Get shared library path and parameters +CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) +{ + // Get shared library path and parameters auto path = get_node_value(node, "library", false, true); std::string parameters; if (check_for_node(node, "parameters")) { @@ -358,7 +361,9 @@ CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) { setup(path, parameters); } -void CompiledSourceWrapper::setup(const std::string& path, const std::string& parameters) { +void CompiledSourceWrapper::setup( + const std::string& path, const std::string& parameters) +{ #ifdef HAS_DYNAMIC_LINKING // Open the library shared_library_ = dlopen(path.c_str(), RTLD_LAZY); @@ -429,7 +434,8 @@ MeshSource::MeshSource(pugi::xml_node node) const auto& mesh = model::meshes[mesh_idx]; std::vector strengths; - // read all source distributions and populate strengths vector for MeshSpatial object + // read all source distributions and populate strengths vector for MeshSpatial + // object for (auto source_node : node.children("source")) { sources_.emplace_back(IndependentSource(source_node)); strengths.push_back(sources_.back().strength()); From 06f7b3aa9ea51972ae57c7c1c8b37dea44f00fff Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Nov 2023 16:15:26 -0600 Subject: [PATCH 48/73] Update openmc/stats/multivariate.py Co-authored-by: Jonathan Shimwell --- openmc/stats/multivariate.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index ea1ff31d0b2..c8b3f1dbcdd 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -714,7 +714,6 @@ def to_xml_element(self, parent_elem=None): else: element = ET.Element('space') - element.set('type', 'mesh') element.set("mesh_id", str(self.mesh.id)) element.set("volume_normalized", str(self.volume_normalized)) From c1f7034927c5d648c84f7c7faf69bf1b2751ac05 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Nov 2023 15:11:20 -0500 Subject: [PATCH 49/73] Remove redundant to_xml_element, small fixes in MeshSource --- docs/source/pythonapi/base.rst | 1 + openmc/source.py | 53 +++++++++++++++------------------- 2 files changed, 24 insertions(+), 30 deletions(-) diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index afecaff5333..478a57609fb 100644 --- a/docs/source/pythonapi/base.rst +++ b/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 diff --git a/openmc/source.py b/openmc/source.py index 7cce3155cdb..789d77686c9 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -313,19 +313,6 @@ def populate_xml_element(self, element): id_elem = ET.SubElement(element, "domain_ids") id_elem.text = ' '.join(str(uid) for uid in self.domain_ids) - def to_xml_element(self) -> ET.Element: - """Return XML representation of the source - - Returns - ------- - element : xml.etree.ElementTree.Element - XML element containing source data - - """ - element = ET.Element("source") - self.populate_xml_element(element) - return element - @classmethod def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': """Generate source from an XML element @@ -393,17 +380,23 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': return source + class MeshSource(SourceBase): - """A mesh-based source in which random positions are uniformly sampled - within mesh elements and each element can have independent angle, energy, - and time distributions. The element sampled is chosen based on the relative - strengths of the sources applied to the elements. The strength of the mesh - source as a whole is the sum of all source strengths applied to the elements. + """A source with a spatial distribution over mesh elements + + This class represents a mesh-based source in which random positions are + uniformly sampled within mesh elements and each element can have independent + angle, energy, and time distributions. The element sampled is chosen based + on the relative strengths of the sources applied to the elements. The + strength of the mesh source as a whole is the sum of all source strengths + applied to the elements. + + .. versionadded:: 0.14.1 Parameters ---------- mesh : openmc.MeshBase - The mesh on which source sites will be generated. + The mesh over which source sites will be generated. sources : iterable of openmc.SourceBase Sources for each element in the mesh. If spatial distributions are set on any of the source objects, they will be ignored during source site @@ -412,15 +405,16 @@ class MeshSource(SourceBase): Attributes ---------- mesh : openmc.MeshBase - The mesh on which source sites will be generated. + The mesh over which source sites will be generated. sources : numpy.ndarray or iterable of openmc.SourceBase - The set of sources to apply to each element. The shape of this array must match - the shape of the mesh with and exception in the case of unstructured mesh, which - allows for application of 1-D array or iterable. + The set of sources to apply to each element. The shape of this array + must match the shape of the mesh with and exception in the case of + unstructured mesh, which allows for application of 1-D array or + iterable. """ def __init__(self, mesh: openmc.MeshBase, - sources: Optional[Iterable]): + sources: Sequence[SourceBase]): self.mesh = mesh self.sources = sources @@ -462,19 +456,18 @@ def sources(self, s): raise ValueError('Sources must be a 1-D array for unstructured mesh') self._sources = s - for src in self.sources.flat: + for src in self._sources.flat: if src.space is not None: - warnings.warn('Some sources on the mesh have ' - 'spatial distributions that will ' - 'be ignored at runtime.') + warnings.warn('Some sources on the mesh have spatial ' + 'distributions that will be ignored at runtime.') break @strength.setter def strength(self, val): cv.check_type('mesh source strength', val, Real) - self._set_total_strength(val) + self.set_total_strength(val) - def _set_total_strength(self, new_strength): + def set_total_strength(self, new_strength): """Scales the element source strengths based on a desired total mesh strength. """ From 95f7c8d3dd2c2d000dd254294e6054b22a90e8fc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Nov 2023 15:11:41 -0500 Subject: [PATCH 50/73] Don't allow MeshSources.sources to be empty --- openmc/source.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index 789d77686c9..b6c4d951b5b 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -428,8 +428,6 @@ def mesh(self): @property def strength(self): - if self.sources is None: - return 0.0 return sum(s.strength for s in self.sources.flat) @property From b22a710edda0ea00b75c4ceaa5e4fb97933ec3c8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Nov 2023 22:17:49 -0500 Subject: [PATCH 51/73] Add type hints to MeshSource --- openmc/source.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index b6c4d951b5b..81865a7bb12 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -419,19 +419,19 @@ def __init__(self, self.sources = sources @property - def type(self): + def type(self) -> str: return "mesh" @property - def mesh(self): + def mesh(self) -> openmc.MeshBase: return self._mesh @property - def strength(self): + def strength(self) -> float: return sum(s.strength for s in self.sources.flat) @property - def sources(self): + def sources(self) -> np.ndarray: return self._sources @mesh.setter @@ -465,7 +465,7 @@ def strength(self, val): cv.check_type('mesh source strength', val, Real) self.set_total_strength(val) - def set_total_strength(self, new_strength): + def set_total_strength(self, new_strength: float): """Scales the element source strengths based on a desired total mesh strength. """ @@ -475,11 +475,10 @@ def set_total_strength(self, new_strength): s.strength *= new_strength / current_strength def normalize_source_strengths(self): - """Update all element source strengths such that they sum to 1.0. - """ + """Update all element source strengths such that they sum to 1.0.""" self.set_total_strength(1.0) - def populate_xml_element(self, elem): + def populate_xml_element(self, elem: ET.Element): super().populate_xml_element(elem) elem.set("mesh", str(self.mesh.id)) From e1a24956dc3ea70913c06492c35b688292c874df Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Nov 2023 10:35:14 -0600 Subject: [PATCH 52/73] Don't require populate_xml_element to call super --- openmc/source.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index 81865a7bb12..6882bb041ba 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -61,8 +61,6 @@ def populate_xml_element(self, element): XML element containing source data """ - element.set("type", self.type) - element.set("strength", str(self.strength)) def to_xml_element(self) -> ET.Element: """Return XML representation of the source @@ -74,6 +72,8 @@ def to_xml_element(self) -> ET.Element: """ element = ET.Element("source") + element.set("type", self.type) + element.set("strength", str(self.strength)) self.populate_xml_element(element) return element @@ -297,7 +297,6 @@ def populate_xml_element(self, element): XML element containing source data """ - super().populate_xml_element(element) element.set("particle", self.particle) if self.space is not None: element.append(self.space.to_xml_element()) @@ -479,8 +478,6 @@ def normalize_source_strengths(self): self.set_total_strength(1.0) def populate_xml_element(self, elem: ET.Element): - super().populate_xml_element(elem) - elem.set("mesh", str(self.mesh.id)) # write in the order of mesh indices @@ -594,8 +591,6 @@ def populate_xml_element(self, element): XML element containing source data """ - super().populate_xml_element(element) - element.set("library", self.library) if self.parameters is not None: @@ -686,8 +681,6 @@ def populate_xml_element(self, element): XML element containing source data """ - super().populate_xml_element(element) - if self.path is not None: element.set("file", self.path) From 1300b012963351b56fc1540275626f5d58e201db Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 9 Nov 2023 11:35:47 -0600 Subject: [PATCH 53/73] Switch to to_xml_element for mesh source element sources --- openmc/source.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index 6882bb041ba..00344a95984 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -483,9 +483,7 @@ def populate_xml_element(self, elem: ET.Element): # write in the order of mesh indices for idx in self.mesh.indices: idx = tuple(i - 1 for i in idx) - source = self.sources[idx] - src_element = ET.SubElement(elem, 'source') - source.populate_xml_element(src_element) + elem.append(self.sources[idx].to_xml_element()) @classmethod def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.MeshSource: From 7eabd7982013f57bd62591962e6f6a065fa09528 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 9 Nov 2023 11:57:05 -0600 Subject: [PATCH 54/73] Allow any source type to be used in the MeshSource (C++). --- include/openmc/source.h | 8 +++--- src/source.cpp | 56 ++++++----------------------------------- 2 files changed, 10 insertions(+), 54 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 9e64117241b..1de85f6e3bc 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -156,23 +156,21 @@ class MeshSource : public Source { SourceSite sample(uint64_t* seed) const override; // Properties - ParticleType particle_type() const { return particle_; } - double strength() const + double strength() const override { return space_->total_strength(); } //!< Total source strength // Accessors - const IndependentSource& source(int32_t i) const + const std::unique_ptr& source(int32_t i) const { return sources_.size() == 1 ? sources_[0] : sources_[i]; } private: // Data members - ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted unique_ptr space_; //!< Mesh spatial - vector sources_; //!< Source distributions + vector> sources_; //!< Source distributions }; //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index 33f4a08f14d..20a46a98c31 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -64,6 +64,8 @@ unique_ptr Source::create(pugi::xml_node node) return make_unique(node); } else if (source_type == "mesh") { return make_unique(node); + } else { + fatal_error(fmt::format("Invalid source type '{}' found.", source_type)); } } else { // support legacy source format @@ -416,19 +418,6 @@ CompiledSourceWrapper::~CompiledSourceWrapper() MeshSource::MeshSource(pugi::xml_node node) { - // Check for particle type - if (check_for_node(node, "particle")) { - auto temp_str = get_node_value(node, "particle", true, true); - if (temp_str == "neutron") { - particle_ = ParticleType::neutron; - } else if (temp_str == "photon") { - particle_ = ParticleType::photon; - settings::photon_transport = true; - } else { - fatal_error(std::string("Unknown source particle type: ") + temp_str); - } - } - int32_t mesh_id = stoi(get_node_value(node, "mesh")); int32_t mesh_idx = model::mesh_map.at(mesh_id); const auto& mesh = model::meshes[mesh_idx]; @@ -437,8 +426,8 @@ MeshSource::MeshSource(pugi::xml_node node) // read all source distributions and populate strengths vector for MeshSpatial // object for (auto source_node : node.children("source")) { - sources_.emplace_back(IndependentSource(source_node)); - strengths.push_back(sources_.back().strength()); + sources_.emplace_back(Source::create(source_node)); + strengths.push_back(sources_.back()->strength()); } // the number of source distributions should either be one or equal to the @@ -454,48 +443,17 @@ MeshSource::MeshSource(pugi::xml_node node) SourceSite MeshSource::sample(uint64_t* seed) const { - SourceSite site; - site.particle = particle_; - - int n_reject = 0; - static int n_accept = 0; // sample location and element from mesh auto mesh_location = space_->sample_mesh(seed); - site.r = mesh_location.second; + Position r_in_element = mesh_location.second; int32_t element = mesh_location.first; - const auto& src = source(element); - - site.u = src.angle()->sample(seed); - - Distribution* e_dist = src.energy(); + SourceSite site = source(element)->sample(seed); - auto p = static_cast(particle_); - // TODO: this is copied code from IndependentSource::sample, refactor - while (true) { - // Sample energy spectrum - site.E = e_dist->sample(seed); - - // Resample if energy falls above maximum particle energy - if (site.E < data::energy_max[p]) - break; - - n_reject++; - if (n_reject >= EXTSRC_REJECT_THRESHOLD && - static_cast(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) { - fatal_error("More than 95% of external source sites sampled were " - "rejected. Please check your external source energy spectrum " - "definition."); - } - } - - site.time = src.time()->sample(seed); - - // Increment number of accepted samples - ++n_accept; + site.r = r_in_element; return site; } From 21405df213afa7856ca55b7c6955d5bbad5c4cb1 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 9 Nov 2023 12:01:11 -0600 Subject: [PATCH 55/73] Add check for source normalization method on mesh sources --- tests/regression_tests/source_mesh/test.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py index fddd7e58547..02eb7f8d742 100644 --- a/tests/regression_tests/source_mesh/test.py +++ b/tests/regression_tests/source_mesh/test.py @@ -100,3 +100,6 @@ def test_source_mesh(mesh_type): assert mesh_source.strength == 1.0 mesh_source.strength = 100.0 assert mesh_source.strength == 100.0 + + mesh_source.normalize_source_strengths() + assert mesh_source.strength == 1.0 \ No newline at end of file From 4d2f2219be09763e315871ca6f9665dd59da8846 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 9 Nov 2023 12:02:01 -0600 Subject: [PATCH 56/73] Only display warning about spatial distributions for IndependentSource objects --- openmc/source.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index 00344a95984..5167c98ee9b 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -454,7 +454,7 @@ def sources(self, s): self._sources = s for src in self._sources.flat: - if src.space is not None: + if isinstance(src, IndependentSource) and src.space is not None: warnings.warn('Some sources on the mesh have spatial ' 'distributions that will be ignored at runtime.') break From ab2c5eb937d90dfd9c7651b9d31131b1ce299d04 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 9 Nov 2023 12:05:24 -0600 Subject: [PATCH 57/73] Move MeshSource test to uni tests directory and add tmpdir fixture --- .../regression_tests/source_mesh/__init__.py | 0 tests/regression_tests/source_mesh/test.py | 105 ------- tests/unit_tests/test_source_mesh.py | 290 ++++++------------ 3 files changed, 96 insertions(+), 299 deletions(-) delete mode 100644 tests/regression_tests/source_mesh/__init__.py delete mode 100644 tests/regression_tests/source_mesh/test.py diff --git a/tests/regression_tests/source_mesh/__init__.py b/tests/regression_tests/source_mesh/__init__.py deleted file mode 100644 index e69de29bb2d..00000000000 diff --git a/tests/regression_tests/source_mesh/test.py b/tests/regression_tests/source_mesh/test.py deleted file mode 100644 index 02eb7f8d742..00000000000 --- a/tests/regression_tests/source_mesh/test.py +++ /dev/null @@ -1,105 +0,0 @@ - -import numpy as np -import pytest - -import openmc -import openmc.model - -@pytest.mark.parametrize('mesh_type', ('rectangular', 'cylindrical')) -def test_source_mesh(mesh_type): - """ - A void model containing a single box - """ - min, max = -10, 10 - box = openmc.model.RectangularParallelepiped(min, max, min, max, min, max, boundary_type='vacuum') - - geometry = openmc.Geometry([openmc.Cell(region=-box)]) - - settings = openmc.Settings() - settings.particles = 100 - settings.batches = 10 - settings.run_mode = 'fixed source' - - model = openmc.Model(geometry=geometry, settings=settings) - - # define a 2 x 2 x 2 mesh - if mesh_type == 'rectangular': - mesh = openmc.RegularMesh.from_domain(model.geometry, (2, 2, 2)) - elif mesh_type == 'cylindrical': - mesh = openmc.CylindricalMesh.from_domain(model.geometry, (1, 4, 2)) - - energy = openmc.stats.Discrete([1.e6], [1.0]) - - # create sources with only one non-zero strength for the source in the - # mesh voxel occupyting the lowest octant. Direct source particles straight - # out of the problem from there. This demonstrates that - # 1) particles are only being sourced within the intented mesh voxel based on source strength - # 2) particles are respecting the angle distributions assigned to each voxel - sources = np.ndarray(mesh.dimension, dtype=openmc.SourceBase) - centroids = mesh.centroids - x, y, z = np.swapaxes(mesh.centroids, -1, 0) - for i, j, k in mesh.indices: - # mesh.indices is currently one-indexed, adjust for Python arrays - ijk = (i-1, j-1, k-1) - - # get the centroid of the ijk mesh element, set a particle source - # vector based on the - centroid = centroids[ijk] - vec = np.sign(centroid, dtype=float) - vec /= np.linalg.norm(vec) - angle = openmc.stats.Monodirectional(vec) - sources[ijk] = openmc.Source(energy=energy, angle=angle, strength=0.0) - - # create and apply the mesh source - mesh_source = openmc.MeshSource(mesh, sources) - model.settings.source = mesh_source - - # tally the flux on the mesh - mesh_filter = openmc.MeshFilter(mesh) - tally = openmc.Tally() - tally.filters = [mesh_filter] - tally.scores = ['flux'] - - model.tallies = openmc.Tallies([tally]) - - # for each element, set a single-non zero source with particles - # traveling out of the mesh (and geometry) w/o crossing any other - # mesh elements - for i, j, k in mesh.indices: - ijk = (i-1, j-1, k-1) - # zero-out all source strengths and set the strength - # on the element of interest - mesh_source.strength = 0.0 - mesh_source.sources[ijk].strength = 1.0 - - sp_file = model.run() - - with openmc.StatePoint(sp_file) as sp: - tally_out = sp.get_tally(id=tally.id) - mean = tally_out.get_reshaped_data(expand_dims=True) - - # remove nuclides and scores axes - mean = mean[..., 0, 0] - # the mesh elment with a non-zero source strength should have a value - assert mean[ijk] != 0 - # all other values should be zero - mean[ijk] = 0 - assert np.all(mean == 0), f'Failed on index {ijk} with centroid {mesh.centroids[ijk]}' - - # test roundtrip - xml_model = openmc.Model.from_model_xml() - xml_source = xml_model.settings.source[0] - assert isinstance(xml_source, openmc.MeshSource) - assert xml_source.strength == 1.0 - assert isinstance(xml_source.mesh, type(mesh_source.mesh)) - assert xml_source.mesh.dimension == mesh_source.mesh.dimension - assert xml_source.mesh.id == mesh_source.mesh.id - assert len(xml_source.sources) == len(mesh_source.sources) - - # check strength adjustment methods - assert mesh_source.strength == 1.0 - mesh_source.strength = 100.0 - assert mesh_source.strength == 100.0 - - mesh_source.normalize_source_strengths() - assert mesh_source.strength == 1.0 \ No newline at end of file diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index fd8b04148e7..eacbb5149fb 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -1,203 +1,105 @@ -from itertools import product -from pathlib import Path -from subprocess import call -import pytest import numpy as np -import openmc -import openmc.lib - -from tests import cdtemp -from tests.regression_tests import config - - -TETS_PER_VOXEL = 12 - -# This test uses a geometry file with cells that match a regular mesh. Each cell -# in the geometry corresponds to 12 tetrahedra in the unstructured mesh file. -@pytest.fixture -def model(): - openmc.reset_auto_ids() - - ### Materials ### - materials = openmc.Materials() - - water_mat = openmc.Material(name="water") - water_mat.add_nuclide("H1", 2.0) - water_mat.add_nuclide("O16", 1.0) - water_mat.set_density("atom/b-cm", 0.07416) - materials.append(water_mat) - - ### Geometry ### - # This test uses a geometry file that resembles a regular mesh. - # 12 tets are used to match each voxel in the geometry. +import pytest - # create a regular mesh that matches the superimposed mesh - regular_mesh = openmc.RegularMesh(mesh_id=10) - regular_mesh.lower_left = (-10, -10, -10) - regular_mesh.dimension = (10, 10, 10) - regular_mesh.width = (2, 2, 2) +import openmc +import openmc.model - root_cell, _ = regular_mesh.build_cells(bc=['vacuum']*6) +@pytest.mark.parametrize('mesh_type', ('rectangular', 'cylindrical')) +def test_source_mesh(run_in_tmpdir, mesh_type): + """ + A void model containing a single box + """ + min, max = -10, 10 + box = openmc.model.RectangularParallelepiped(min, max, min, max, min, max, boundary_type='vacuum') - geometry = openmc.Geometry(root=[root_cell]) + geometry = openmc.Geometry([openmc.Cell(region=-box)]) - ### Settings ### settings = openmc.Settings() - settings.run_mode = 'fixed source' settings.particles = 100 - settings.batches = 2 - - return openmc.model.Model(geometry=geometry, - materials=materials, - settings=settings) - -### Setup test cases ### -param_values = (['libmesh', 'moab'], # mesh libraries - ['uniform', 'manual']) # Element weighting schemes - -test_cases = [] -for i, (lib, schemes) in enumerate(product(*param_values)): - test_cases.append({'library' : lib, - 'source_strengths' : schemes}) - -def ids(params): - """Test naming function for clarity""" - return f"{params['library']}-{params['source_strengths']}" - -@pytest.mark.parametrize("test_cases", test_cases, ids=ids) -def test_unstructured_mesh_sampling(model, request, test_cases): - # skip the test if the library is not enabled - if test_cases['library'] == 'moab' and not openmc.lib._dagmc_enabled(): - pytest.skip("DAGMC (and MOAB) mesh not enabled in this build.") - - if test_cases['library'] == 'libmesh' and not openmc.lib._libmesh_enabled(): - pytest.skip("LibMesh is not enabled in this build.") - - # setup mesh source ### - mesh_filename = Path(request.fspath).parent / "test_mesh_tets.e" - uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_cases['library']) - - # subtract one to account for root cell produced by RegularMesh.build_cells - n_cells = len(model.geometry.get_all_cells()) - 1 - - # set source weights according to test case - if test_cases['source_strengths'] == 'uniform': - vol_norm = True - strengths = None - elif test_cases['source_strengths'] == 'manual': - vol_norm = False - # assign random weights - strengths = np.random.rand(n_cells*TETS_PER_VOXEL) - - # create the spatial distribution based on the mesh - space = openmc.stats.MeshSpatial(uscd_mesh, strengths, vol_norm) - - energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0]) - source = openmc.IndependentSource(space=space, energy=energy) - model.settings.source = source - - with cdtemp([mesh_filename]): - model.export_to_xml() - - n_measurements = 100 - n_samples = 1000 - - cell_counts = np.zeros((n_cells, n_measurements)) - - # This model contains 1000 geometry cells. Each cell is a hex - # corresponding to 12 of the tets. This test runs 1000 samples. This - # results in the following average for each cell - openmc.lib.init([]) - - # perform many sets of samples and track counts for each cell - for m in range(n_measurements): - sites = openmc.lib.sample_external_source(n_samples) - cells = [openmc.lib.find_cell(s.r) for s in sites] - - for c in cells: - # subtract one from index to account for root cell - cell_counts[c[0]._index - 1, m] += 1 - - # make sure particle transport is successful - openmc.lib.run() - openmc.lib.finalize() - - # normalize cell counts to get sampling frequency per particle - cell_counts /= n_samples - - # get the mean and std. dev. of the cell counts - mean = cell_counts.mean(axis=1) - std_dev = cell_counts.std(axis=1) - - if test_cases['source_strengths'] == 'uniform': - exp_vals = np.ones(n_cells) / n_cells - else: - # sum up the source strengths for each tet, these are the expected true mean - # of the sampling frequency for that cell - exp_vals = strengths.reshape(-1, 12).sum(axis=1) / sum(strengths) - - diff = np.abs(mean - exp_vals) - assert((diff < 2*std_dev).sum() / diff.size >= 0.95) - assert((diff < 6*std_dev).sum() / diff.size >= 0.997) - - -def test_strengths_size_failure(request, model): - # setup mesh source ### - mesh_filename = Path(request.fspath).parent / "test_mesh_tets.e" - uscd_mesh = openmc.UnstructuredMesh(mesh_filename, 'libmesh') - - # intentionally incorrectly sized to trigger an error - n_cells = len(model.geometry.get_all_cells()) - strengths = np.random.rand(n_cells*TETS_PER_VOXEL) - - # create the spatial distribution based on the mesh - space = openmc.stats.MeshSpatial(uscd_mesh, strengths) - - energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0]) - source = openmc.IndependentSource(space=space, energy=energy) - model.settings.source = source - - # skip the test if unstructured mesh is not available - if not openmc.lib._libmesh_enabled(): - if openmc.lib._dagmc_enabled(): - source.space.mesh.library = 'moab' - else: - pytest.skip("Unstructured mesh support unavailable.") - - # make sure that an incorrrectly sized strengths array causes a failure - source.space.strengths = source.space.strengths[:-1] - - mesh_filename = Path(request.fspath).parent / source.space.mesh.filename - - with pytest.raises(RuntimeError, match=r'strengths array'), cdtemp([mesh_filename]): - model.export_to_xml() - openmc.run() - -def test_roundtrip(run_in_tmpdir, model, request): - if not openmc.lib._libmesh_enabled() and not openmc.lib._dagmc_enabled(): - pytest.skip("Unstructured mesh is not enabled in this build.") - - mesh_filename = Path(request.fspath).parent / 'test_mesh_tets.e' - ucd_mesh = openmc.UnstructuredMesh(mesh_filename, library='libmesh') - - if not openmc.lib._libmesh_enabled(): - ucd_mesh.library = 'moab' - - n_cells = len(model.geometry.get_all_cells()) - - space_out = openmc.MeshSpatial(ucd_mesh) - space_out.strengths = np.random.rand(n_cells*TETS_PER_VOXEL) - model.settings.source = openmc.IndependentSource(space=space_out) - - # write out the model - model.export_to_xml() - - model_in = openmc.Model.from_xml() - - space_in = model_in.settings.source[0].space - - np.testing.assert_equal(space_out.strengths, space_in.strengths) + settings.batches = 10 + settings.run_mode = 'fixed source' - assert space_in.mesh.id == space_out.mesh.id - assert space_in.volume_normalized == space_out.volume_normalized + model = openmc.Model(geometry=geometry, settings=settings) + + # define a 2 x 2 x 2 mesh + if mesh_type == 'rectangular': + mesh = openmc.RegularMesh.from_domain(model.geometry, (2, 2, 2)) + elif mesh_type == 'cylindrical': + mesh = openmc.CylindricalMesh.from_domain(model.geometry, (1, 4, 2)) + + energy = openmc.stats.Discrete([1.e6], [1.0]) + + # create sources with only one non-zero strength for the source in the + # mesh voxel occupyting the lowest octant. Direct source particles straight + # out of the problem from there. This demonstrates that + # 1) particles are only being sourced within the intented mesh voxel based on source strength + # 2) particles are respecting the angle distributions assigned to each voxel + sources = np.ndarray(mesh.dimension, dtype=openmc.SourceBase) + centroids = mesh.centroids + x, y, z = np.swapaxes(mesh.centroids, -1, 0) + for i, j, k in mesh.indices: + # mesh.indices is currently one-indexed, adjust for Python arrays + ijk = (i-1, j-1, k-1) + + # get the centroid of the ijk mesh element, set a particle source + # vector based on the + centroid = centroids[ijk] + vec = np.sign(centroid, dtype=float) + vec /= np.linalg.norm(vec) + angle = openmc.stats.Monodirectional(vec) + sources[ijk] = openmc.Source(energy=energy, angle=angle, strength=0.0) + + # create and apply the mesh source + mesh_source = openmc.MeshSource(mesh, sources) + model.settings.source = mesh_source + + # tally the flux on the mesh + mesh_filter = openmc.MeshFilter(mesh) + tally = openmc.Tally() + tally.filters = [mesh_filter] + tally.scores = ['flux'] + + model.tallies = openmc.Tallies([tally]) + + # for each element, set a single-non zero source with particles + # traveling out of the mesh (and geometry) w/o crossing any other + # mesh elements + for i, j, k in mesh.indices: + ijk = (i-1, j-1, k-1) + # zero-out all source strengths and set the strength + # on the element of interest + mesh_source.strength = 0.0 + mesh_source.sources[ijk].strength = 1.0 + + sp_file = model.run() + + with openmc.StatePoint(sp_file) as sp: + tally_out = sp.get_tally(id=tally.id) + mean = tally_out.get_reshaped_data(expand_dims=True) + + # remove nuclides and scores axes + mean = mean[..., 0, 0] + # the mesh elment with a non-zero source strength should have a value + assert mean[ijk] != 0 + # all other values should be zero + mean[ijk] = 0 + assert np.all(mean == 0), f'Failed on index {ijk} with centroid {mesh.centroids[ijk]}' + + # test roundtrip + xml_model = openmc.Model.from_model_xml() + xml_source = xml_model.settings.source[0] + assert isinstance(xml_source, openmc.MeshSource) + assert xml_source.strength == 1.0 + assert isinstance(xml_source.mesh, type(mesh_source.mesh)) + assert xml_source.mesh.dimension == mesh_source.mesh.dimension + assert xml_source.mesh.id == mesh_source.mesh.id + assert len(xml_source.sources) == len(mesh_source.sources) + + # check strength adjustment methods + assert mesh_source.strength == 1.0 + mesh_source.strength = 100.0 + assert mesh_source.strength == 100.0 + + mesh_source.normalize_source_strengths() + assert mesh_source.strength == 1.0 \ No newline at end of file From b4369f73d9130fd4b7b4ab473f9b1a7045afd1d5 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 9 Nov 2023 14:16:57 -0600 Subject: [PATCH 58/73] Adding a test for a file source --- openmc/bounding_box.py | 5 +++ tests/unit_tests/test_source_mesh.py | 61 +++++++++++++++++++++++++++- 2 files changed, 65 insertions(+), 1 deletion(-) diff --git a/openmc/bounding_box.py b/openmc/bounding_box.py index b74931d73ae..0d117ae2265 100644 --- a/openmc/bounding_box.py +++ b/openmc/bounding_box.py @@ -95,6 +95,11 @@ def __or__(self, other: BoundingBox) -> BoundingBox: new |= other return new + def __contains__(self, point): + """Check whether or not a point is in the bounding box + """ + return all(point > self.lower_left) and all(point < self.upper_right) + @property def center(self) -> np.ndarray: return (self[0] + self[1]) / 2 diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index eacbb5149fb..82386929d5b 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -3,6 +3,7 @@ import pytest import openmc +import openmc.lib import openmc.model @pytest.mark.parametrize('mesh_type', ('rectangular', 'cylindrical')) @@ -102,4 +103,62 @@ def test_source_mesh(run_in_tmpdir, mesh_type): assert mesh_source.strength == 100.0 mesh_source.normalize_source_strengths() - assert mesh_source.strength == 1.0 \ No newline at end of file + assert mesh_source.strength == 1.0 + + +def test_file_source(run_in_tmpdir): + source_particle = openmc.SourceParticle(r=(0.0, 0.0, 0.0), + u=(0.0, 0.0, 1.0), + E=1e6, + time=10.0) + + openmc.write_source_file([source_particle], 'source.h5') + + file_source = openmc.FileSource('source.h5') + + model = openmc.Model() + + rect_prism = openmc.model.RectangularParallelepiped(-5.0, 5.0, + -5.0, 5.0, + -5.0, 5.0, + boundary_type='vacuum') + mat = openmc.Material() + mat.add_nuclide('H1', 1.0) + + model.geometry = openmc.Geometry([openmc.Cell(fill=mat, region=-rect_prism)]) + model.settings.particles = 1000 + model.settings.batches = 10 + model.settings.inactive = 1 + model.settings.run_mode = 'fixed source' + + mesh = openmc.RegularMesh() + mesh.lower_left = (-1, -2, -3) + mesh.upper_right = (2, 3, 4) + mesh.dimension = (1, 1, 1) + + mesh_source_arr = np.asarray([file_source]).reshape(mesh.dimension) + source = openmc.MeshSource(mesh, mesh_source_arr) + + model.settings.source = source + + model.export_to_model_xml() + + openmc.lib.init() + openmc.lib.simulation_init() + sites = openmc.lib.sample_external_source(10) + openmc.lib.finalize() + + + # the mesh bounds do not contain the point of the lone + # source site in the file source, so it should not appear + # in the set of source sites produced from the mesh source + + # additionally, the source should be located within the mesh + bbox = mesh.bounding_box + for site in sites: + assert site.r != (0, 0, 0) + assert site.E == source_particle.E + assert site.u == source_particle.u + assert site.time == source_particle.time + assert site.r in bbox + From 1f5b99c93e96db8cf0bd3c5e6300f8b7256fe095 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 9 Nov 2023 14:18:29 -0600 Subject: [PATCH 59/73] applying C++ style guide --- include/openmc/source.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 1de85f6e3bc..4a399b6fd5f 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -169,8 +169,8 @@ class MeshSource : public Source { private: // Data members - unique_ptr space_; //!< Mesh spatial - vector> sources_; //!< Source distributions + unique_ptr space_; //!< Mesh spatial + vector> sources_; //!< Source distributions }; //============================================================================== From 1a755c79ed22c93758975e6809c785b1eceba348 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 15 Nov 2023 14:07:35 -0600 Subject: [PATCH 60/73] Adding tmate debug option --- .github/workflows/ci.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ab2121b9c96..2291f40ea87 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -145,6 +145,10 @@ jobs: CTEST_OUTPUT_ON_FAILURE=1 make test -C $GITHUB_WORKSPACE/build/ $GITHUB_WORKSPACE/tools/ci/gha-script.sh + - name: tmate debug + if: ${{ failure() }} + uses: mxschmitt/action-tmate@v3 + - name: after_success shell: bash run: | From 3ab9ff59a50f1dc902d94e8e6cc4009fcffcf1d5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 27 Nov 2023 14:00:34 -0600 Subject: [PATCH 61/73] Fix openmc.lib use in test_source_mesh.py --- .github/workflows/ci.yml | 4 ---- tests/unit_tests/test_source_mesh.py | 3 ++- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2291f40ea87..ab2121b9c96 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -145,10 +145,6 @@ jobs: CTEST_OUTPUT_ON_FAILURE=1 make test -C $GITHUB_WORKSPACE/build/ $GITHUB_WORKSPACE/tools/ci/gha-script.sh - - name: tmate debug - if: ${{ failure() }} - uses: mxschmitt/action-tmate@v3 - - name: after_success shell: bash run: | diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index 82386929d5b..0374343d8a3 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -49,7 +49,7 @@ def test_source_mesh(run_in_tmpdir, mesh_type): vec = np.sign(centroid, dtype=float) vec /= np.linalg.norm(vec) angle = openmc.stats.Monodirectional(vec) - sources[ijk] = openmc.Source(energy=energy, angle=angle, strength=0.0) + sources[ijk] = openmc.IndependentSource(energy=energy, angle=angle, strength=0.0) # create and apply the mesh source mesh_source = openmc.MeshSource(mesh, sources) @@ -146,6 +146,7 @@ def test_file_source(run_in_tmpdir): openmc.lib.init() openmc.lib.simulation_init() sites = openmc.lib.sample_external_source(10) + openmc.lib.simulation_finalize() openmc.lib.finalize() From 71e66300241048cb6ef6707971a08e8675255c9b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 28 Nov 2023 11:11:38 -0600 Subject: [PATCH 62/73] Comment cleanup --- tests/unit_tests/test_source_mesh.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index 0374343d8a3..34befcffa13 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -152,9 +152,8 @@ def test_file_source(run_in_tmpdir): # the mesh bounds do not contain the point of the lone # source site in the file source, so it should not appear - # in the set of source sites produced from the mesh source - - # additionally, the source should be located within the mesh + # in the set of source sites produced from the mesh source. + # Additionally, the source should be located within the mesh bbox = mesh.bounding_box for site in sites: assert site.r != (0, 0, 0) From f0b5a343cf9fa7387bfe084894bd83e9f6356aa4 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 29 Nov 2023 06:15:49 -0600 Subject: [PATCH 63/73] Making public data member placement consistent --- include/openmc/mesh.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 75de700ff8a..00abf73cbc8 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -450,8 +450,6 @@ class CylindricalMesh : public PeriodicStructuredMesh { double volume(const MeshIndex& ijk) const override; - array, 3> grid_; - // grid accessors double r(int i) const { return grid_[0][i]; } double phi(int i) const { return grid_[1][i]; } @@ -459,6 +457,9 @@ class CylindricalMesh : public PeriodicStructuredMesh { int set_grid(); + // Data members + array, 3> grid_; + private: double find_r_crossing( const Position& r, const Direction& u, double l, int shell) const; @@ -514,10 +515,12 @@ class SphericalMesh : public PeriodicStructuredMesh { 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]; } - array, 3> grid_; int set_grid(); + // Data members + array, 3> grid_; + private: double find_r_crossing( const Position& r, const Direction& u, double l, int shell) const; From 9c16efe94846ea2c0d04f78cddad7599964e02a3 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 29 Nov 2023 06:17:51 -0600 Subject: [PATCH 64/73] Back to SourceBase for settings.source CheckedList --- openmc/settings.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 64ca0d6fa18..3ef2514e04f 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -492,7 +492,7 @@ def source(self) -> typing.List[SourceBase]: def source(self, source: typing.Union[SourceBase, typing.Iterable[SourceBase]]): if not isinstance(source, MutableSequence): source = [source] - self._source = cv.CheckedList((SourceBase, MeshSource), 'source distributions', source) + self._source = cv.CheckedList(SourceBase, 'source distributions', source) @property def confidence_intervals(self) -> bool: @@ -1878,11 +1878,10 @@ def from_xml_element(cls, elem, meshes=None): Settings object """ + # read all meshes under the settings node and update settings_meshes = _read_meshes(elem) - if meshes is None: - meshes = settings_meshes - else: - meshes.update(settings_meshes) + meshes = dict() if meshes is None else meshes + meshes.update(settings_meshes) settings = cls() settings._eigenvalue_from_xml_element(elem) From ddbd73cc71d1e93f21e14e431fdb619b26b9005c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 29 Nov 2023 06:18:43 -0600 Subject: [PATCH 65/73] I don't think we need the parent element argument for MeshSpatial --- openmc/stats/multivariate.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index c8b3f1dbcdd..22c3d7ea86b 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -700,7 +700,7 @@ def num_strength_bins(self): raise ValueError('Strengths are not set') return self.strengths.size - def to_xml_element(self, parent_elem=None): + def to_xml_element(self): """Return XML representation of the spatial distribution Returns @@ -709,10 +709,7 @@ def to_xml_element(self, parent_elem=None): XML element containing spatial distribution data """ - if parent_elem is not None: - element = ET.SubElement(parent_elem, 'space') - else: - element = ET.Element('space') + element = ET.Element('space') element.set('type', 'mesh') element.set("mesh_id", str(self.mesh.id)) From 66b6baad1b7d532516446792cec3e707ccb0ce69 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 29 Nov 2023 06:22:22 -0600 Subject: [PATCH 66/73] Replacing old tests for MeshSpatial that I clobbered accidentally --- tests/unit_tests/test_source_mesh.py | 210 ++++++++++++++++++++++++++- 1 file changed, 206 insertions(+), 4 deletions(-) diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index 34befcffa13..079f2b819db 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -1,11 +1,214 @@ +from itertools import product +from pathlib import Path -import numpy as np import pytest - +import numpy as np import openmc import openmc.lib import openmc.model +from tests import cdtemp +from tests.regression_tests import config + + +############### +# MeshSpatial Tests +############### +TETS_PER_VOXEL = 12 + +# This test uses a geometry file with cells that match a regular mesh. Each cell +# in the geometry corresponds to 12 tetrahedra in the unstructured mesh file. +@pytest.fixture +def model(): + openmc.reset_auto_ids() + + ### Materials ### + materials = openmc.Materials() + + water_mat = openmc.Material(name="water") + water_mat.add_nuclide("H1", 2.0) + water_mat.add_nuclide("O16", 1.0) + water_mat.set_density("atom/b-cm", 0.07416) + materials.append(water_mat) + + ### Geometry ### + # This test uses a geometry file that resembles a regular mesh. + # 12 tets are used to match each voxel in the geometry. + + # create a regular mesh that matches the superimposed mesh + regular_mesh = openmc.RegularMesh(mesh_id=10) + regular_mesh.lower_left = (-10, -10, -10) + regular_mesh.dimension = (10, 10, 10) + regular_mesh.width = (2, 2, 2) + + root_cell, _ = regular_mesh.build_cells(bc=['vacuum']*6) + + geometry = openmc.Geometry(root=[root_cell]) + + ### Settings ### + settings = openmc.Settings() + settings.run_mode = 'fixed source' + settings.particles = 100 + settings.batches = 2 + + return openmc.model.Model(geometry=geometry, + materials=materials, + settings=settings) + +### Setup test cases ### +param_values = (['libmesh', 'moab'], # mesh libraries + ['uniform', 'manual']) # Element weighting schemes + +test_cases = [] +for i, (lib, schemes) in enumerate(product(*param_values)): + test_cases.append({'library' : lib, + 'source_strengths' : schemes}) + +def ids(params): + """Test naming function for clarity""" + return f"{params['library']}-{params['source_strengths']}" + +@pytest.mark.parametrize("test_cases", test_cases, ids=ids) +def test_unstructured_mesh_sampling(model, request, test_cases): + # skip the test if the library is not enabled + if test_cases['library'] == 'moab' and not openmc.lib._dagmc_enabled(): + pytest.skip("DAGMC (and MOAB) mesh not enabled in this build.") + + if test_cases['library'] == 'libmesh' and not openmc.lib._libmesh_enabled(): + pytest.skip("LibMesh is not enabled in this build.") + + # setup mesh source ### + mesh_filename = Path(request.fspath).parent / "test_mesh_tets.e" + uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_cases['library']) + + # subtract one to account for root cell produced by RegularMesh.build_cells + n_cells = len(model.geometry.get_all_cells()) - 1 + + # set source weights according to test case + if test_cases['source_strengths'] == 'uniform': + vol_norm = True + strengths = None + elif test_cases['source_strengths'] == 'manual': + vol_norm = False + # assign random weights + strengths = np.random.rand(n_cells*TETS_PER_VOXEL) + + # create the spatial distribution based on the mesh + space = openmc.stats.MeshSpatial(uscd_mesh, strengths, vol_norm) + + energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0]) + source = openmc.IndependentSource(space=space, energy=energy) + model.settings.source = source + + with cdtemp([mesh_filename]): + model.export_to_xml() + + n_measurements = 100 + n_samples = 1000 + + cell_counts = np.zeros((n_cells, n_measurements)) + + # This model contains 1000 geometry cells. Each cell is a hex + # corresponding to 12 of the tets. This test runs 1000 samples. This + # results in the following average for each cell + openmc.lib.init([]) + + # perform many sets of samples and track counts for each cell + for m in range(n_measurements): + sites = openmc.lib.sample_external_source(n_samples) + cells = [openmc.lib.find_cell(s.r) for s in sites] + + for c in cells: + # subtract one from index to account for root cell + cell_counts[c[0]._index - 1, m] += 1 + + # make sure particle transport is successful + openmc.lib.run() + openmc.lib.finalize() + + # normalize cell counts to get sampling frequency per particle + cell_counts /= n_samples + + # get the mean and std. dev. of the cell counts + mean = cell_counts.mean(axis=1) + std_dev = cell_counts.std(axis=1) + + if test_cases['source_strengths'] == 'uniform': + exp_vals = np.ones(n_cells) / n_cells + else: + # sum up the source strengths for each tet, these are the expected true mean + # of the sampling frequency for that cell + exp_vals = strengths.reshape(-1, 12).sum(axis=1) / sum(strengths) + + diff = np.abs(mean - exp_vals) + assert((diff < 2*std_dev).sum() / diff.size >= 0.95) + assert((diff < 6*std_dev).sum() / diff.size >= 0.997) + + +def test_strengths_size_failure(request, model): + # setup mesh source ### + mesh_filename = Path(request.fspath).parent / "test_mesh_tets.e" + uscd_mesh = openmc.UnstructuredMesh(mesh_filename, 'libmesh') + + # intentionally incorrectly sized to trigger an error + n_cells = len(model.geometry.get_all_cells()) + strengths = np.random.rand(n_cells*TETS_PER_VOXEL) + + # create the spatial distribution based on the mesh + space = openmc.stats.MeshSpatial(uscd_mesh, strengths) + + energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0]) + source = openmc.IndependentSource(space=space, energy=energy) + model.settings.source = source + + # skip the test if unstructured mesh is not available + if not openmc.lib._libmesh_enabled(): + if openmc.lib._dagmc_enabled(): + source.space.mesh.library = 'moab' + else: + pytest.skip("Unstructured mesh support unavailable.") + + # make sure that an incorrrectly sized strengths array causes a failure + source.space.strengths = source.space.strengths[:-1] + + mesh_filename = Path(request.fspath).parent / source.space.mesh.filename + + with pytest.raises(RuntimeError, match=r'strengths array'), cdtemp([mesh_filename]): + model.export_to_xml() + openmc.run() + +def test_roundtrip(run_in_tmpdir, model, request): + if not openmc.lib._libmesh_enabled() and not openmc.lib._dagmc_enabled(): + pytest.skip("Unstructured mesh is not enabled in this build.") + + mesh_filename = Path(request.fspath).parent / 'test_mesh_tets.e' + ucd_mesh = openmc.UnstructuredMesh(mesh_filename, library='libmesh') + + if not openmc.lib._libmesh_enabled(): + ucd_mesh.library = 'moab' + + n_cells = len(model.geometry.get_all_cells()) + + space_out = openmc.MeshSpatial(ucd_mesh) + space_out.strengths = np.random.rand(n_cells*TETS_PER_VOXEL) + model.settings.source = openmc.IndependentSource(space=space_out) + + # write out the model + model.export_to_xml() + + model_in = openmc.Model.from_xml() + + space_in = model_in.settings.source[0].space + + np.testing.assert_equal(space_out.strengths, space_in.strengths) + + assert space_in.mesh.id == space_out.mesh.id + assert space_in.volume_normalized == space_out.volume_normalized + + +################### +# MeshSource tests +################### @pytest.mark.parametrize('mesh_type', ('rectangular', 'cylindrical')) def test_source_mesh(run_in_tmpdir, mesh_type): """ @@ -160,5 +363,4 @@ def test_file_source(run_in_tmpdir): assert site.E == source_particle.E assert site.u == source_particle.u assert site.time == source_particle.time - assert site.r in bbox - + assert site.r in bbox \ No newline at end of file From 5f0f958d3ffadd65ffc934af97b39396f1d296aa Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Dec 2023 20:52:46 -0600 Subject: [PATCH 67/73] Update mesh source io_formats documentation --- docs/source/io_formats/settings.rst | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 4f57ccf98ec..9b855b25029 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -473,6 +473,8 @@ pseudo-random number generator. *Default*: 1 +.. _source_element: + -------------------- ```` Element -------------------- @@ -665,8 +667,13 @@ attributes/sub-elements: *Default*: false + :mesh: + For mesh sources, this indicates the ID of the corresponding mesh. + :source: - For mesh sources, a set of sources that will be applied to each element. + 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 + ```` sub-elements should correspond to the number of mesh elements. .. _univariate: From b2f8bd144efc2ad0b75a8af25311e6cb97522e96 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Dec 2023 20:52:59 -0600 Subject: [PATCH 68/73] Utilize UnitSphereDistribution::create in source.cpp --- src/source.cpp | 19 +------------------ 1 file changed, 1 insertion(+), 18 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 20a46a98c31..87a22b74135 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -124,24 +124,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) // Determine external source angular distribution if (check_for_node(node, "angle")) { - // Get pointer to angular distribution - pugi::xml_node node_angle = node.child("angle"); - - // Check for type of angular distribution - std::string type; - if (check_for_node(node_angle, "type")) - type = get_node_value(node_angle, "type", true, true); - if (type == "isotropic") { - angle_ = UPtrAngle {new Isotropic()}; - } else if (type == "monodirectional") { - angle_ = UPtrAngle {new Monodirectional(node_angle)}; - } else if (type == "mu-phi") { - angle_ = UPtrAngle {new PolarAzimuthal(node_angle)}; - } else { - fatal_error(fmt::format( - "Invalid angular distribution for external source: {}", type)); - } - + angle_ = UnitSphereDistribution::create(node.child("angle")); } else { angle_ = UPtrAngle {new Isotropic()}; } From fa6edfde242610f1c0104e058492b4a2440c3712 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Dec 2023 22:56:25 -0600 Subject: [PATCH 69/73] Style, documentation fixed --- include/openmc/mesh.h | 8 +++---- include/openmc/source.h | 5 +--- openmc/bounding_box.py | 3 +-- openmc/settings.py | 2 +- openmc/source.py | 53 +++++++++++++++++++++++++++-------------- src/source.cpp | 8 +++---- 6 files changed, 45 insertions(+), 34 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 00abf73cbc8..dca34b53ee6 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -82,11 +82,11 @@ class Mesh { //! Return a position in the local coordinates of the mesh virtual Position local_coords(const Position& r) const { return r; }; - //! Sample mesh element + //! 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 + //! \param[inout] seed Seed to use for random sampling + //! \param[in] bin Bin value of the mesh element sampled + //! \return sampled position within mesh element virtual Position sample_element(uint64_t* seed, int32_t bin) const = 0; //! Determine which bins were crossed by a particle diff --git a/include/openmc/source.h b/include/openmc/source.h index 4a399b6fd5f..d5ea9bd1d56 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -156,10 +156,7 @@ class MeshSource : public Source { SourceSite sample(uint64_t* seed) const override; // Properties - double strength() const override - { - return space_->total_strength(); - } //!< Total source strength + double strength() const override { return space_->total_strength(); } // Accessors const std::unique_ptr& source(int32_t i) const diff --git a/openmc/bounding_box.py b/openmc/bounding_box.py index 0d117ae2265..e5655cf8cae 100644 --- a/openmc/bounding_box.py +++ b/openmc/bounding_box.py @@ -96,8 +96,7 @@ def __or__(self, other: BoundingBox) -> BoundingBox: return new def __contains__(self, point): - """Check whether or not a point is in the bounding box - """ + """Check whether or not a point is in the bounding box""" return all(point > self.lower_left) and all(point < self.upper_right) @property diff --git a/openmc/settings.py b/openmc/settings.py index 3ef2514e04f..0d287cb9e3b 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1880,7 +1880,7 @@ def from_xml_element(cls, elem, meshes=None): """ # read all meshes under the settings node and update settings_meshes = _read_meshes(elem) - meshes = dict() if meshes is None else meshes + meshes = {} if meshes is None else meshes meshes.update(settings_meshes) settings = cls() diff --git a/openmc/source.py b/openmc/source.py index 5167c98ee9b..de73f13b7b4 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -31,7 +31,7 @@ class SourceBase(ABC): Attributes ---------- - type : {'independent', 'file', 'compiled'} + type : {'independent', 'file', 'compiled', 'mesh'} Indicator of source type. strength : float Strength of the source @@ -78,7 +78,7 @@ def to_xml_element(self) -> ET.Element: return element @classmethod - def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.SourceBase: + def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: """Generate source from an XML element Parameters @@ -291,6 +291,7 @@ def domain_type(self, domain_type): def populate_xml_element(self, element): """Add necessary source information to an XML element + Returns ------- element : lxml.etree._Element @@ -313,7 +314,7 @@ def populate_xml_element(self, element): id_elem.text = ' '.join(str(uid) for uid in self.domain_ids) @classmethod - def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': + def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: """Generate source from an XML element Parameters @@ -410,10 +411,13 @@ class MeshSource(SourceBase): must match the shape of the mesh with and exception in the case of unstructured mesh, which allows for application of 1-D array or iterable. + strength : float + Strength of the source + type : str + Indicator of source type: 'mesh' + """ - def __init__(self, - mesh: openmc.MeshBase, - sources: Sequence[SourceBase]): + def __init__(self, mesh: MeshBase, sources: Sequence[SourceBase]): self.mesh = mesh self.sources = sources @@ -422,7 +426,7 @@ def type(self) -> str: return "mesh" @property - def mesh(self) -> openmc.MeshBase: + def mesh(self) -> MeshBase: return self._mesh @property @@ -435,18 +439,18 @@ def sources(self) -> np.ndarray: @mesh.setter def mesh(self, m): - cv.check_type('source mesh', m, openmc.MeshBase) + cv.check_type('source mesh', m, MeshBase) self._mesh = m @sources.setter def sources(self, s): - cv.check_iterable_type('mesh sources', s, openmc.SourceBase, max_depth=3) + cv.check_iterable_type('mesh sources', s, SourceBase, max_depth=3) s = np.asarray(s) if isinstance(self.mesh, StructuredMesh) and s.shape != self.mesh.dimension: - raise ValueError('The shape of the source array' \ - f'({s.shape}) does not match the ' \ + raise ValueError('The shape of the source array' + f'({s.shape}) does not match the ' f'dimensions of the structured mesh ({self.mesh.dimension})') elif isinstance(self.mesh, UnstructuredMesh): if len(s.shape) > 1: @@ -464,20 +468,33 @@ def strength(self, val): cv.check_type('mesh source strength', val, Real) self.set_total_strength(val) - def set_total_strength(self, new_strength: float): - """Scales the element source strengths based on a desired - total mesh strength. + def set_total_strength(self, strength: float): + """Scales the element source strengths based on a desired total strength. + + Parameters + ---------- + strength : float + Total source strength + """ current_strength = self.strength if self.strength != 0.0 else 1.0 for s in self.sources.flat: - s.strength *= new_strength / current_strength + s.strength *= strength / current_strength def normalize_source_strengths(self): """Update all element source strengths such that they sum to 1.0.""" self.set_total_strength(1.0) def populate_xml_element(self, elem: ET.Element): + """Add necessary source information to an XML element + + Returns + ------- + element : lxml.etree._Element + XML element containing source data + + """ elem.set("mesh", str(self.mesh.id)) # write in the order of mesh indices @@ -486,7 +503,7 @@ def populate_xml_element(self, elem: ET.Element): elem.append(self.sources[idx].to_xml_element()) @classmethod - def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.MeshSource: + def from_xml_element(cls, elem: ET.Element, meshes) -> openmc.MeshSource: """ Generate MeshSource from an XML element @@ -507,7 +524,7 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.MeshSource: mesh = meshes[mesh_id] - sources = [SourceBase.from_xml_element(e) for e in elem.iter('source') if e != elem] + sources = [SourceBase.from_xml_element(e) for e in elem.iterchildren('source')] sources = np.asarray(sources).reshape(mesh.dimension, order='F') return cls(mesh, sources) @@ -595,7 +612,7 @@ def populate_xml_element(self, element): element.set("parameters", self.parameters) @classmethod - def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.CompiledSource: + def from_xml_element(cls, elem: ET.Element) -> openmc.CompiledSource: """Generate a compiled source from an XML element Parameters diff --git a/src/source.cpp b/src/source.cpp index 87a22b74135..bf52b568a26 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -426,17 +426,15 @@ MeshSource::MeshSource(pugi::xml_node node) SourceSite MeshSource::sample(uint64_t* seed) const { - // sample location and element from mesh auto mesh_location = space_->sample_mesh(seed); - Position r_in_element = mesh_location.second; - + // Sample source for the chosen element int32_t element = mesh_location.first; - SourceSite site = source(element)->sample(seed); - site.r = r_in_element; + // Replace spatial position with the one already sampled + site.r = mesh_location.second; return site; } From 8fbd212a68e8621dd079c2c60c763e1cd821692d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Dec 2023 22:56:40 -0600 Subject: [PATCH 70/73] Use uniform_distribution in sample_element methods --- src/mesh.cpp | 31 ++++++++++++++----------------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 264a7d5852d..0a315c1c89c 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -27,6 +27,7 @@ #include "openmc/hdf5_interface.h" #include "openmc/memory.h" #include "openmc/message_passing.h" +#include "openmc/random_dist.h" #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/tallies/filter.h" @@ -1103,12 +1104,12 @@ Position CylindricalMesh::sample_element( double r_min_sq = r_min * r_min; double r_max_sq = r_max * r_max; - double r = sqrt(r_min_sq + (r_max_sq - r_min_sq) * prn(seed)); - double phi = phi_min + (phi_max - phi_min) * prn(seed); - double z = z_min + (z_max - z_min) * prn(seed); + double r = std::sqrt(uniform_distribution(r_min_sq, r_max_sq, seed)); + double phi = uniform_distribution(phi_min, phi_max, seed); + double z = uniform_distribution(z_min, z_max, seed); - double x = r * cos(phi); - double y = r * sin(phi); + double x = r * std::cos(phi); + double y = r * std::sin(phi); return origin_ + Position(x, y, z); } @@ -1386,20 +1387,16 @@ Position SphericalMesh::sample_element( double phi_min = this->phi(ijk[2] - 1); double phi_max = this->phi(ijk[2]); - double u1 = prn(seed); - double u2 = prn(seed); - double u3 = prn(seed); - - double cos_theta = cos(theta_min) + (cos(theta_max) - cos(theta_min)) * u1; - double sin_theta = sin(acos(cos_theta)); - double phi = phi_min + (phi_max - phi_min) * u2; - double r_min_cub = pow(r_min, 3); - double r_max_cub = pow(r_max, 3); + double cos_theta = uniform_distribution(theta_min, theta_max, seed); + double sin_theta = std::sin(std::acos(cos_theta)); + double phi = uniform_distribution(phi_min, phi_max, seed); + double r_min_cub = std::pow(r_min, 3); + double r_max_cub = std::pow(r_max, 3); // might be faster to do rejection here? - double r = cbrt(r_min_cub + (r_max_cub - r_min_cub) * u3); + double r = std::cbrt(uniform_distribution(r_min_cub, r_max_cub, seed)); - double x = r * cos(phi) * sin_theta; - double y = r * sin(phi) * sin_theta; + double x = r * std::cos(phi) * sin_theta; + double y = r * std::sin(phi) * sin_theta; double z = r * cos_theta; return origin_ + Position(x, y, z); From 05902a21da4e6bc82748063c7c0c337b925e8acb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 2 Dec 2023 09:27:59 -0600 Subject: [PATCH 71/73] Small changes in MeshSource tests --- tests/unit_tests/test_source_mesh.py | 62 +++++++++++++--------------- 1 file changed, 28 insertions(+), 34 deletions(-) diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index 079f2b819db..7bec2d0724b 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -5,15 +5,13 @@ import numpy as np import openmc import openmc.lib -import openmc.model from tests import cdtemp -from tests.regression_tests import config -############### +################### # MeshSpatial Tests -############### +################### TETS_PER_VOXEL = 12 # This test uses a geometry file with cells that match a regular mesh. Each cell @@ -51,9 +49,9 @@ def model(): settings.particles = 100 settings.batches = 2 - return openmc.model.Model(geometry=geometry, - materials=materials, - settings=settings) + return openmc.Model(geometry=geometry, + materials=materials, + settings=settings) ### Setup test cases ### param_values = (['libmesh', 'moab'], # mesh libraries @@ -177,6 +175,7 @@ def test_strengths_size_failure(request, model): model.export_to_xml() openmc.run() + def test_roundtrip(run_in_tmpdir, model, request): if not openmc.lib._libmesh_enabled() and not openmc.lib._dagmc_enabled(): pytest.skip("Unstructured mesh is not enabled in this build.") @@ -210,12 +209,13 @@ def test_roundtrip(run_in_tmpdir, model, request): # MeshSource tests ################### @pytest.mark.parametrize('mesh_type', ('rectangular', 'cylindrical')) -def test_source_mesh(run_in_tmpdir, mesh_type): +def test_mesh_source_independent(run_in_tmpdir, mesh_type): """ A void model containing a single box """ min, max = -10, 10 - box = openmc.model.RectangularParallelepiped(min, max, min, max, min, max, boundary_type='vacuum') + box = openmc.model.RectangularParallelepiped( + min, max, min, max, min, max, boundary_type='vacuum') geometry = openmc.Geometry([openmc.Cell(region=-box)]) @@ -234,20 +234,21 @@ def test_source_mesh(run_in_tmpdir, mesh_type): energy = openmc.stats.Discrete([1.e6], [1.0]) - # create sources with only one non-zero strength for the source in the - # mesh voxel occupyting the lowest octant. Direct source particles straight - # out of the problem from there. This demonstrates that - # 1) particles are only being sourced within the intented mesh voxel based on source strength + # create sources with only one non-zero strength for the source in the mesh + # voxel occupying the lowest octant. Direct source particles straight out of + # the problem from there. This demonstrates that + # 1) particles are only being sourced within the intented mesh voxel based + # on source strength # 2) particles are respecting the angle distributions assigned to each voxel - sources = np.ndarray(mesh.dimension, dtype=openmc.SourceBase) + sources = np.empty(mesh.dimension, dtype=openmc.SourceBase) centroids = mesh.centroids x, y, z = np.swapaxes(mesh.centroids, -1, 0) for i, j, k in mesh.indices: # mesh.indices is currently one-indexed, adjust for Python arrays ijk = (i-1, j-1, k-1) - # get the centroid of the ijk mesh element, set a particle source - # vector based on the + # get the centroid of the ijk mesh element and use it to set the + # direction of the source directly out of the problem centroid = centroids[ijk] vec = np.sign(centroid, dtype=float) vec /= np.linalg.norm(vec) @@ -309,29 +310,23 @@ def test_source_mesh(run_in_tmpdir, mesh_type): assert mesh_source.strength == 1.0 -def test_file_source(run_in_tmpdir): - source_particle = openmc.SourceParticle(r=(0.0, 0.0, 0.0), - u=(0.0, 0.0, 1.0), - E=1e6, - time=10.0) - +def test_mesh_source_file(run_in_tmpdir): + # Creating a source file with a single particle + source_particle = openmc.SourceParticle(time=10.0) openmc.write_source_file([source_particle], 'source.h5') - file_source = openmc.FileSource('source.h5') model = openmc.Model() - rect_prism = openmc.model.RectangularParallelepiped(-5.0, 5.0, - -5.0, 5.0, - -5.0, 5.0, - boundary_type='vacuum') + rect_prism = openmc.model.RectangularParallelepiped( + -5.0, 5.0, -5.0, 5.0, -5.0, 5.0, boundary_type='vacuum') + mat = openmc.Material() mat.add_nuclide('H1', 1.0) model.geometry = openmc.Geometry([openmc.Cell(fill=mat, region=-rect_prism)]) model.settings.particles = 1000 model.settings.batches = 10 - model.settings.inactive = 1 model.settings.run_mode = 'fixed source' mesh = openmc.RegularMesh() @@ -352,15 +347,14 @@ def test_file_source(run_in_tmpdir): openmc.lib.simulation_finalize() openmc.lib.finalize() - - # the mesh bounds do not contain the point of the lone - # source site in the file source, so it should not appear - # in the set of source sites produced from the mesh source. - # Additionally, the source should be located within the mesh + # The mesh bounds do not contain the point of the lone source site in the + # file source, so it should not appear in the set of source sites produced + # from the mesh source. Additionally, the source should be located within + # the mesh bbox = mesh.bounding_box for site in sites: assert site.r != (0, 0, 0) assert site.E == source_particle.E assert site.u == source_particle.u assert site.time == source_particle.time - assert site.r in bbox \ No newline at end of file + assert site.r in bbox From 77ceb628d8184f1968f72321cd4251212b8f1b18 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 2 Dec 2023 09:35:21 -0600 Subject: [PATCH 72/73] Swap order of arguments in sample_element --- include/openmc/mesh.h | 18 +++++++++--------- src/distribution_spatial.cpp | 2 +- src/mesh.cpp | 10 +++++----- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index dca34b53ee6..de556321f22 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -84,10 +84,10 @@ class Mesh { //! Sample a position within a mesh element // - //! \param[inout] seed Seed to use for random sampling //! \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(uint64_t* seed, int32_t bin) const = 0; + virtual Position sample_element(int32_t bin, uint64_t* seed) const = 0; //! Determine which bins were crossed by a particle // @@ -184,12 +184,12 @@ class StructuredMesh : public Mesh { } }; - Position sample_element(uint64_t* seed, int32_t bin) const override + Position sample_element(int32_t bin, uint64_t* seed) const override { - return sample_element(seed, get_indices_from_bin(bin)); + return sample_element(get_indices_from_bin(bin), seed); }; - virtual Position sample_element(uint64_t* seed, const MeshIndex& ijk) const; + virtual Position sample_element(const MeshIndex& ijk, uint64_t* seed) const; int get_bin(Position r) const override; @@ -438,7 +438,7 @@ class CylindricalMesh : public PeriodicStructuredMesh { static const std::string mesh_type; - Position sample_element(uint64_t* seed, const MeshIndex& ijk) const override; + 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; @@ -502,7 +502,7 @@ class SphericalMesh : public PeriodicStructuredMesh { static const std::string mesh_type; - Position sample_element(uint64_t* seed, const MeshIndex& ijk) const override; + 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; @@ -664,7 +664,7 @@ class MOABMesh : public UnstructuredMesh { // Overridden Methods - Position sample_element(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& bins, vector& lengths) const override; @@ -832,7 +832,7 @@ class LibMesh : public UnstructuredMesh { void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; - Position sample_element(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; diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index bbbf5b8f5b6..8f34ea6b936 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -285,7 +285,7 @@ std::pair MeshSpatial::sample_mesh(uint64_t* seed) const { // Sample the CDF defined in initialization above int32_t elem_idx = elem_idx_dist_.sample(seed); - return {elem_idx, mesh()->sample_element(seed, elem_idx)}; + return {elem_idx, mesh()->sample_element(elem_idx, seed)}; } Position MeshSpatial::sample(uint64_t* seed) const diff --git a/src/mesh.cpp b/src/mesh.cpp index 0a315c1c89c..ea66552de40 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -167,7 +167,7 @@ xt::xtensor StructuredMesh::get_x_shape() const } Position StructuredMesh::sample_element( - uint64_t* seed, const MeshIndex& ijk) const + const MeshIndex& ijk, uint64_t* seed) const { // lookup the lower/upper bounds for the mesh element double x_min = negative_grid_boundary(ijk, 0); @@ -1091,7 +1091,7 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( } Position CylindricalMesh::sample_element( - uint64_t* seed, const MeshIndex& ijk) const + const MeshIndex& ijk, uint64_t* seed) const { double r_min = this->r(ijk[0] - 1); double r_max = this->r(ijk[0]); @@ -1376,7 +1376,7 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( } Position SphericalMesh::sample_element( - uint64_t* seed, const MeshIndex& ijk) const + const MeshIndex& ijk, uint64_t* seed) const { double r_min = this->r(ijk[0] - 1); double r_max = this->r(ijk[0]); @@ -2249,7 +2249,7 @@ std::string MOABMesh::library() const } // Sample position within a tet for MOAB type tets -Position MOABMesh::sample_element(uint64_t* seed, int32_t bin) const +Position MOABMesh::sample_element(int32_t bin, uint64_t* seed) const { moab::EntityHandle tet_ent = get_ent_handle_from_bin(bin); @@ -2736,7 +2736,7 @@ void LibMesh::initialize() } // Sample position within a tet for LibMesh type tets -Position LibMesh::sample_element(uint64_t* seed, int32_t bin) const +Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const { const auto& elem = get_element_from_bin(bin); // Get tet vertex coordinates from LibMesh From 709b95299af4d85725dbd469e1e473947a98c49d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 2 Dec 2023 09:40:44 -0600 Subject: [PATCH 73/73] Make sure MeshSource writing updates mesh_memo --- openmc/settings.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 0d287cb9e3b..5d0800ac13d 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1072,7 +1072,7 @@ def _create_max_order_subelement(self, root): element = ET.SubElement(root, "max_order") element.text = str(self._max_order) - def _create_source_subelement(self, root): + def _create_source_subelement(self, root, mesh_memo=None): for source in self.source: root.append(source.to_xml_element()) if isinstance(source, IndependentSource) and isinstance(source.space, MeshSpatial): @@ -1083,6 +1083,8 @@ def _create_source_subelement(self, root): path = f"./mesh[@id='{source.mesh.id}']" if root.find(path) is None: root.append(source.mesh.to_xml_element()) + if mesh_memo is not None: + mesh_memo.add(source.mesh.id) def _create_volume_calcs_subelement(self, root): for calc in self.volume_calculations: @@ -1369,7 +1371,8 @@ def _create_weight_windows_subelement(self, root, mesh_memo=None): path = f"./mesh[@id='{ww.mesh.id}']" if root.find(path) is None: root.append(ww.mesh.to_xml_element()) - if mesh_memo is not None: mesh_memo.add(ww.mesh.id) + if mesh_memo is not None: + mesh_memo.add(ww.mesh.id) if self._weight_windows_on is not None: elem = ET.SubElement(root, "weight_windows_on") @@ -1791,7 +1794,7 @@ def to_xml_element(self, mesh_memo=None): self._create_max_write_lost_particles_subelement(element) self._create_generations_per_batch_subelement(element) self._create_keff_trigger_subelement(element) - self._create_source_subelement(element) + self._create_source_subelement(element, mesh_memo) self._create_output_subelement(element) self._create_statepoint_subelement(element) self._create_sourcepoint_subelement(element)