From bd60956d20c9e538179597d0c6e8029cbea2963c Mon Sep 17 00:00:00 2001 From: NybergWISC Date: Wed, 10 Aug 2022 14:42:57 -0500 Subject: [PATCH] Second round of comments. volume method is now implemented from Mesh instead of UnstructuredMesh --- include/openmc/distribution_spatial.h | 10 ++++---- include/openmc/mesh.h | 19 +++++++++++---- openmc/settings.py | 13 +++++++---- openmc/stats/multivariate.py | 8 +++---- src/distribution_spatial.cpp | 23 +++++++++---------- src/mesh.cpp | 8 ++++++- src/source.cpp | 2 +- .../unstructured_mesh/source_sampling/test.py | 21 ++++++++--------- 8 files changed, 60 insertions(+), 44 deletions(-) diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 8864e1139..6d771fd6b 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -99,9 +99,9 @@ private: //! Distribution of points within a mesh //============================================================================== -class MeshIndependent : public SpatialDistribution { +class MeshSpatial : public SpatialDistribution { public: - explicit MeshIndependent(pugi::xml_node node); + explicit MeshSpatial(pugi::xml_node node); //! Sample a position from the distribution //! \param seed Pseudorandom number seed pointer @@ -109,13 +109,13 @@ public: Position sample(uint64_t* seed) const; private: - UnstructuredMesh* umesh_ptr_; - int32_t mesh_map_idx_; + Mesh* mesh_ptr_; + int32_t mesh_idx_; std::string sample_scheme_; double total_strength_; std::vector mesh_CDF_; std::vector mesh_strengths_; - int64_t tot_bins_; + int32_t tot_bins_; }; //============================================================================== diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index a4104f75a..8015d2e4c 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -79,7 +79,13 @@ public: // //! \param[in] seed Seed to use for random sampling //! \param[out] r Position within tet - virtual Position sample(uint64_t* seed, int32_t tet_bin) const=0; + virtual Position sample(uint64_t* seed, int32_t bin) const=0; + + //! Get the volume of a mesh bin + // + //! \param[in] bin Bin to return the volume for + //! \return Volume of the bin + virtual double volume(int bin) const = 0; //! Determine which bins were crossed by a particle // @@ -166,7 +172,9 @@ public: } }; - Position sample(uint64_t* seed, int32_t tet_bin) const override; + Position sample(uint64_t* seed, int32_t bin) const override; + + double volume(int bin) const override; int get_bin(Position r) const override; @@ -469,13 +477,14 @@ public: virtual std::string get_mesh_type() const override; // Overridden Methods - // TODO Position sample(uint64_t* seed) const=0; void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; void to_hdf5(hid_t group) const override; + virtual double volume(int bin) const = 0; + std::string bin_label(int bin) const override; // Methods @@ -565,7 +574,7 @@ public: // Overridden Methods - Position sample(uint64_t* seed, int32_t tet_bin) const override; + Position sample(uint64_t* seed, int32_t bin) const override; void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; @@ -729,7 +738,7 @@ public: void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; - Position sample(uint64_t* seed, int32_t tet_bin) const override; + Position sample(uint64_t* seed, int32_t bin) const override; int get_bin(Position r) const override; diff --git a/openmc/settings.py b/openmc/settings.py index 045cccd5c..b8b1ec2b4 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -11,7 +11,7 @@ from typing import Optional from xml.etree import ElementTree as ET import openmc.checkvalue as cv -from openmc.stats.multivariate import MeshIndependent +from openmc.stats.multivariate import MeshSpatial from . import RegularMesh, Source, VolumeCalculation, WeightWindows from ._xml import clean_indentation, get_text, reorder_attributes @@ -977,7 +977,7 @@ class Settings: def _create_source_subelement(self, root): for source in self.source: root.append(source.to_xml_element()) - if isinstance(source.space, MeshIndependent): + if isinstance(source.space, MeshSpatial): path = f"./mesh[@id='{source.space.mesh.id}']" if root.find(path) is None: root.append(source.space.mesh.to_xml_element()) @@ -1328,13 +1328,16 @@ class Settings: def _source_from_xml_element(self, root): for elem in root.findall('source'): - self.source.append(Source.from_xml_element(elem)) - if isinstance(Source.from_xml_element(elem), MeshIndependent): + src = Source.from_xml_element(elem) + if isinstance(src.space, MeshSpatial): mesh_id = int(get_text(elem, 'mesh')) path = f"./mesh[@id='{mesh_id}']" mesh_elem = root.find(path) if mesh_elem is not None: - self.source.space.mesh = MeshBase.from_xml_element(mesh_elem) + src.space.mesh = MeshBase.from_xml_element(mesh_elem) + else: + raise RuntimeError('No mesh was specified for the mesh source') + self.source.append(src) def _volume_calcs_from_xml_element(self, root): volume_elems = root.findall("volume_calc") diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 0df540ddd..9aa03b2d7 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -271,7 +271,7 @@ class Spatial(ABC): elif distribution == 'spherical': return SphericalIndependent.from_xml_element(elem) elif distribution == 'mesh': - return MeshIndependent.from_xml_element(elem) + return MeshSpatial.from_xml_element(elem) elif distribution == 'box' or distribution == 'fission': return Box.from_xml_element(elem) elif distribution == 'point': @@ -619,10 +619,10 @@ class CylindricalIndependent(Spatial): origin = [float(x) for x in elem.get('origin').split()] return cls(r, phi, z, origin=origin) -class MeshIndependent(Spatial): +class MeshSpatial(Spatial): """Spatial distribution for a mesh. - This distribution allows one to specify a mesh to sample over and the scheme to sample over the entire mesh. + This distribution specifies a mesh to sample over, chooses whether it will be volume normalized, and can set the source strengths. .. versionadded:: 0.13 @@ -716,7 +716,7 @@ class MeshIndependent(Spatial): Returns ------- - openmc.stats.MeshIndependent + openmc.stats.MeshSpatial Spatial distribution generated from XML element """ diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 0ac24e39a..6d00739b0 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -183,24 +183,23 @@ Position SphericalIndependent::sample(uint64_t* seed) const } //============================================================================== -// MeshIndependent implementation +// MeshSpatial implementation //============================================================================== -MeshIndependent::MeshIndependent(pugi::xml_node node) +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 int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id")); // Get pointer to spatial distribution - mesh_map_idx_ = model::mesh_map.at(mesh_id); - const auto& mesh_ptr = model::meshes[mesh_map_idx_]; + mesh_idx_ = model::mesh_map.at(mesh_id); - // Check whether mesh pointer points to an unstructured mesh - umesh_ptr_ = dynamic_cast(mesh_ptr.get()); - if (!umesh_ptr_) {fatal_error("Mesh passed to spatial distribution is not an unstructured mesh object"); } + // Check whether mesh pointer points to a mesh + mesh_ptr_ = dynamic_cast(model::meshes[mesh_idx_].get()); + if (!mesh_ptr_) {fatal_error("Mesh passed to spatial distribution is not a mesh object"); } // Initialize arrays for CDF creation - tot_bins_ = umesh_ptr_->n_bins(); + tot_bins_ = mesh_ptr_->n_bins(); std::vector strengths = {}; strengths.resize(tot_bins_); double temp_total_strength = 0.0; @@ -219,7 +218,7 @@ MeshIndependent::MeshIndependent(pugi::xml_node node) fatal_error("The size of the strengths array from the xml file does not equal the number of elements in the mesh."); } if (get_node_value_bool(node, "volume_normalized")){ for (int i = 0; i < tot_bins_; i++){ - strengths[i] = mesh_strengths_[i]*umesh_ptr_->volume(i); + strengths[i] = mesh_strengths_[i]*mesh_ptr_->volume(i); } } else if (!get_node_value_bool(node, "volume_normalized")){ for (int i = 0; i < tot_bins_; i++){ @@ -228,7 +227,7 @@ MeshIndependent::MeshIndependent(pugi::xml_node node) } } else if (get_node_value_bool(node, "volume_normalized")){ for (int i = 0; ivolume(i); + strengths[i] = mesh_ptr_->volume(i); } } else { for (int i = 0; isample(seed, tet_bin); + return mesh_ptr_->sample(seed, tet_bin); } diff --git a/src/mesh.cpp b/src/mesh.cpp index b37c238f1..453ce572b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -358,8 +358,14 @@ 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(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 diff --git a/src/source.cpp b/src/source.cpp index 2ceff9964..2fdcc319f 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -95,7 +95,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) } else if (type == "spherical") { space_ = UPtrSpace {new SphericalIndependent(node_space)}; } else if (type =="mesh"){ - space_ = UPtrSpace {new MeshIndependent(node_space)}; + space_ = UPtrSpace {new MeshSpatial(node_space)}; } else if (type == "box") { space_ = UPtrSpace {new SpatialBox(node_space)}; } else if (type == "fission") { diff --git a/tests/regression_tests/unstructured_mesh/source_sampling/test.py b/tests/regression_tests/unstructured_mesh/source_sampling/test.py index dfbe74d05..6f8d6fded 100644 --- a/tests/regression_tests/unstructured_mesh/source_sampling/test.py +++ b/tests/regression_tests/unstructured_mesh/source_sampling/test.py @@ -34,15 +34,14 @@ class UnstructuredMeshSourceTest(PyAPITestHarness): # shown below average_in_hex = 10.0 - call(['../../../../scripts/openmc-track-combine', '-o', 'tracks.h5'] + - glob.glob('tracks_p*.h5')) - - # loop over the tracks and get data + # Load in tracks + openmc.Tracks.combine(glob.glob('tracks_p*.h5')) tracks = openmc.Tracks(filepath='tracks.h5') tracks_born = np.empty((len(tracks), 1)) instances = np.zeros(1000) + # loop over the tracks and get data for i in range(0, len(tracks)): tracks_born[i] = tracks[i].particle_tracks[0].states['cell_id'][0] instances[int(tracks_born[i])-1] += 1 @@ -113,9 +112,9 @@ def test_unstructured_mesh(test_opts): for i in range(0,dimen+1): coord = -10.0 + i * size_hex - surfaces[i][0] = openmc.XPlane(coord, name="X plane at "+str(coord)) - surfaces[i][1] = openmc.YPlane(coord, name="Y plane at "+str(coord)) - surfaces[i][2] = openmc.ZPlane(coord, name="Z plane at "+str(coord)) + surfaces[i][0] = openmc.XPlane(coord, name=f"X plane at {coord}") + surfaces[i][1] = openmc.YPlane(coord, name=f"Y plane at {coord}") + surfaces[i][2] = openmc.ZPlane(coord, name=f"Z plane at {coord}") surfaces[i][0].boundary_type = 'vacuum' surfaces[i][1].boundary_type = 'vacuum' @@ -126,8 +125,8 @@ def test_unstructured_mesh(test_opts): for i in range(0,dimen): cell[i][j][k] = openmc.Cell(name=("x = {}, y = {}, z = {}".format(i,j,k))) cell[i][j][k].region = +surfaces[i][0] & -surfaces[i+1][0] & \ - +surfaces[j][1] & -surfaces[j+1][1] & \ - +surfaces[k][2] & -surfaces[k+1][2] + +surfaces[j][1] & -surfaces[j+1][1] & \ + +surfaces[k][2] & -surfaces[k+1][2] cell[i][j][k].fill = None universe.add_cell(cell[i][j][k]) @@ -158,13 +157,13 @@ def test_unstructured_mesh(test_opts): # source setup if test_opts['schemes'] == 'volume': - space = openmc.stats.MeshIndependent(volume_normalized=True, mesh=uscd_mesh) + space = openmc.stats.MeshSpatial(volume_normalized=True, mesh=uscd_mesh) elif test_opts['schemes'] == 'file': array = np.zeros(12000) for i in range(0, 12): array[i] = 10 array[i+12] = 2 - space = openmc.stats.MeshIndependent(volume_normalized=False, strengths=array, mesh=uscd_mesh) + space = openmc.stats.MeshSpatial(volume_normalized=False, strengths=array, mesh=uscd_mesh) energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0]) source = openmc.Source(space=space, energy=energy)