diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index f61e967ca8a..82d20a76302 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -72,7 +72,10 @@ The current version of the statepoint file format is 18.1. **/tallies/meshes/mesh /** -:Datasets: - **type** (*char[]*) -- Type of mesh. +:Attributes: - **id** (*int*) -- ID of the mesh + - **type** (*char[]*) -- Type of mesh. + +:Datasets: - **name** (*char[]*) -- Name of the mesh. - **dimension** (*int*) -- Number of mesh cells in each dimension. - **Regular Mesh Only:** - **lower_left** (*double[]*) -- Coordinates of lower-left corner of diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 78101ab668d..9f29a949acb 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -109,7 +109,7 @@ The ```` element accepts the following sub-elements: prematurely if there are no hits in any bins at the first evalulation. It is the user's responsibility to specify enough particles per batch to get a nonzero score in at least one bin. - + *Default*: False :scores: @@ -329,6 +329,11 @@ If a mesh is desired as a filter for a tally, it must be specified in a separate element with the tag name ````. This element has the following attributes/sub-elements: + :name: + An optional string name to identify the mesh in output files. + + *Default*: "" + :type: The type of mesh. This can be either "regular", "rectilinear", "cylindrical", "spherical", or "unstructured". diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 5edb03003f9..bd86d54fb72 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -132,13 +132,18 @@ class Mesh { int32_t id() const { return id_; } + const std::string& name() const { return name_; } + //! Set the mesh ID void set_id(int32_t id = -1); + //! Write the mesh data to an HDF5 group + void to_hdf5(hid_t group) const; + //! Write mesh data to an HDF5 group // //! \param[in] group HDF5 group - virtual void to_hdf5(hid_t group) const = 0; + virtual void to_hdf5_inner(hid_t group) const = 0; //! Find the mesh lines that intersect an axis-aligned slice plot // @@ -202,7 +207,8 @@ class Mesh { // Data members xt::xtensor lower_left_; //!< Lower-left coordinates of mesh xt::xtensor upper_right_; //!< Upper-right coordinates of mesh - int id_ {-1}; //!< User-specified ID + int id_ {-1}; //!< Mesh ID + std::string name_; //!< User-specified name int n_dimension_ {-1}; //!< Number of dimensions }; @@ -410,7 +416,7 @@ class RegularMesh : public StructuredMesh { std::pair, vector> plot( Position plot_ll, Position plot_ur) const override; - void to_hdf5(hid_t group) const override; + void to_hdf5_inner(hid_t group) const override; //! Get the coordinate for the mesh grid boundary in the positive direction //! @@ -460,7 +466,7 @@ class RectilinearMesh : public StructuredMesh { std::pair, vector> plot( Position plot_ll, Position plot_ur) const override; - void to_hdf5(hid_t group) const override; + void to_hdf5_inner(hid_t group) const override; //! Get the coordinate for the mesh grid boundary in the positive direction //! @@ -506,7 +512,7 @@ class CylindricalMesh : public PeriodicStructuredMesh { std::pair, vector> plot( Position plot_ll, Position plot_ur) const override; - void to_hdf5(hid_t group) const override; + void to_hdf5_inner(hid_t group) const override; double volume(const MeshIndex& ijk) const override; @@ -570,7 +576,7 @@ class SphericalMesh : public PeriodicStructuredMesh { std::pair, vector> plot( Position plot_ll, Position plot_ur) const override; - void to_hdf5(hid_t group) const override; + void to_hdf5_inner(hid_t group) const override; double r(int i) const { return grid_[0][i]; } double theta(int i) const { return grid_[1][i]; } @@ -632,7 +638,7 @@ class UnstructuredMesh : public Mesh { void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; - void to_hdf5(hid_t group) const override; + void to_hdf5_inner(hid_t group) const override; std::string bin_label(int bin) const override; diff --git a/openmc/mesh.py b/openmc/mesh.py index ef5c6c78410..0b6f6b84e4b 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -99,21 +99,40 @@ def from_hdf5(cls, group: h5py.Group): Instance of a MeshBase subclass """ + mesh_type = 'regular' if 'type' not in group.attrs else group.attrs['type'].decode() + mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) + mesh_name = '' if not 'name' in group else group['name'][()].decode() - mesh_type = group['type'][()].decode() if mesh_type == 'regular': - return RegularMesh.from_hdf5(group) + return RegularMesh.from_hdf5(group, mesh_id, mesh_name) elif mesh_type == 'rectilinear': - return RectilinearMesh.from_hdf5(group) + return RectilinearMesh.from_hdf5(group, mesh_id, mesh_name) elif mesh_type == 'cylindrical': - return CylindricalMesh.from_hdf5(group) + return CylindricalMesh.from_hdf5(group, mesh_id, mesh_name) elif mesh_type == 'spherical': - return SphericalMesh.from_hdf5(group) + return SphericalMesh.from_hdf5(group, mesh_id, mesh_name) elif mesh_type == 'unstructured': - return UnstructuredMesh.from_hdf5(group) + return UnstructuredMesh.from_hdf5(group, mesh_id, mesh_name) else: raise ValueError('Unrecognized mesh type: "' + mesh_type + '"') + def to_xml_element(self): + """Return XML representation of the mesh + + Returns + ------- + element : lxml.etree._Element + XML element containing mesh data + + """ + elem = ET.Element("mesh") + + elem.set("id", str(self._id)) + if self.name: + elem.set("name", self.name) + + return elem + @classmethod def from_xml_element(cls, elem: ET.Element): """Generates a mesh from an XML element @@ -132,18 +151,21 @@ def from_xml_element(cls, elem: ET.Element): mesh_type = get_text(elem, 'type') if mesh_type == 'regular' or mesh_type is None: - return RegularMesh.from_xml_element(elem) + mesh = RegularMesh.from_xml_element(elem) elif mesh_type == 'rectilinear': - return RectilinearMesh.from_xml_element(elem) + mesh = RectilinearMesh.from_xml_element(elem) elif mesh_type == 'cylindrical': - return CylindricalMesh.from_xml_element(elem) + mesh = CylindricalMesh.from_xml_element(elem) elif mesh_type == 'spherical': - return SphericalMesh.from_xml_element(elem) + mesh = SphericalMesh.from_xml_element(elem) elif mesh_type == 'unstructured': - return UnstructuredMesh.from_xml_element(elem) + mesh = UnstructuredMesh.from_xml_element(elem) else: raise ValueError(f'Unrecognized mesh type "{mesh_type}" found.') + mesh.name = get_text(elem, 'name', default='') + return mesh + def get_homogenized_materials( self, model: openmc.Model, @@ -791,11 +813,9 @@ def __repr__(self): return string @classmethod - def from_hdf5(cls, group: h5py.Group): - mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) - + def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): # Read and assign mesh properties - mesh = cls(mesh_id) + mesh = cls(mesh_id=mesh_id, name=name) mesh.dimension = group['dimension'][()] mesh.lower_left = group['lower_left'][()] if 'width' in group: @@ -899,9 +919,7 @@ def to_xml_element(self): XML element containing mesh data """ - - element = ET.Element("mesh") - element.set("id", str(self._id)) + element = super().to_xml_element() if self._dimension is not None: subelement = ET.SubElement(element, "dimension") @@ -937,10 +955,6 @@ def from_xml_element(cls, elem: ET.Element): mesh_id = int(get_text(elem, 'id')) mesh = cls(mesh_id=mesh_id) - mesh_type = get_text(elem, 'type') - if mesh_type is not None: - mesh.type = mesh_type - dimension = get_text(elem, 'dimension') if dimension is not None: mesh.dimension = [int(x) for x in dimension.split()] @@ -1235,11 +1249,9 @@ def __repr__(self): return string @classmethod - def from_hdf5(cls, group: h5py.Group): - mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) - + def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): # Read and assign mesh properties - mesh = cls(mesh_id=mesh_id) + mesh = cls(mesh_id=mesh_id, name=name) mesh.x_grid = group['x_grid'][()] mesh.y_grid = group['y_grid'][()] mesh.z_grid = group['z_grid'][()] @@ -1279,8 +1291,7 @@ def to_xml_element(self): """ - element = ET.Element("mesh") - element.set("id", str(self._id)) + element = super().to_xml_element() element.set("type", "rectilinear") subelement = ET.SubElement(element, "x_grid") @@ -1541,12 +1552,11 @@ def get_indices_at_coords( return (r_index, phi_index, z_index) @classmethod - def from_hdf5(cls, group: h5py.Group): - mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) - + def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): # Read and assign mesh properties mesh = cls( mesh_id=mesh_id, + name=name, r_grid = group['r_grid'][()], phi_grid = group['phi_grid'][()], z_grid = group['z_grid'][()], @@ -1647,8 +1657,7 @@ def to_xml_element(self): """ - element = ET.Element("mesh") - element.set("id", str(self._id)) + element = super().to_xml_element() element.set("type", "cylindrical") subelement = ET.SubElement(element, "r_grid") @@ -1926,15 +1935,14 @@ def __repr__(self): return string @classmethod - def from_hdf5(cls, group: h5py.Group): - mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) - + def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): # Read and assign mesh properties mesh = cls( r_grid = group['r_grid'][()], theta_grid = group['theta_grid'][()], phi_grid = group['phi_grid'][()], mesh_id=mesh_id, + name=name ) if 'origin' in group: mesh.origin = group['origin'][()] @@ -1951,8 +1959,7 @@ def to_xml_element(self): """ - element = ET.Element("mesh") - element.set("id", str(self._id)) + element = super().to_xml_element() element.set("type", "spherical") subelement = ET.SubElement(element, "r_grid") @@ -2444,8 +2451,7 @@ def write_data_to_vtk( writer.Write() @classmethod - def from_hdf5(cls, group: h5py.Group): - mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) + def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): filename = group['filename'][()].decode() library = group['library'][()].decode() if 'options' in group.attrs: @@ -2453,7 +2459,7 @@ def from_hdf5(cls, group: h5py.Group): else: options = None - mesh = cls(filename=filename, library=library, mesh_id=mesh_id, options=options) + mesh = cls(filename=filename, library=library, mesh_id=mesh_id, name=name, options=options) mesh._has_statepoint_data = True vol_data = group['volumes'][()] mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],)) @@ -2480,9 +2486,9 @@ def to_xml_element(self): """ - element = ET.Element("mesh") - element.set("id", str(self._id)) + element = super().to_xml_element() element.set("type", "unstructured") + element.set("library", self._library) if self.options is not None: element.set('options', self.options) diff --git a/src/mesh.cpp b/src/mesh.cpp index c572d255a36..2ee616f28dc 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -109,6 +109,8 @@ Mesh::Mesh(pugi::xml_node node) { // Read mesh id id_ = std::stoi(get_node_value(node, "id")); + if (check_for_node(node, "name")) + name_ = get_node_value(node, "name"); } void Mesh::set_id(int32_t id) @@ -236,6 +238,28 @@ vector Mesh::material_volumes( return result; } +void Mesh::to_hdf5(hid_t group) const +{ + // Create group for mesh + std::string group_name = fmt::format("mesh {}", id_); + hid_t mesh_group = create_group(group, group_name.c_str()); + + // Write mesh type + write_attribute(mesh_group, "type", this->get_mesh_type()); + + // Write mesh ID + write_attribute(mesh_group, "id", id_); + + // Write mesh name + write_dataset(mesh_group, "name", name_); + + // Write mesh data + this->to_hdf5_inner(mesh_group); + + // Close group + close_group(mesh_group); +} + //============================================================================== // Structured Mesh implementation //============================================================================== @@ -389,11 +413,8 @@ std::string UnstructuredMesh::bin_label(int bin) const return fmt::format("Mesh Index ({})", bin); }; -void UnstructuredMesh::to_hdf5(hid_t group) const +void UnstructuredMesh::to_hdf5_inner(hid_t mesh_group) const { - hid_t mesh_group = create_group(group, fmt::format("mesh {}", id_)); - - write_dataset(mesh_group, "type", mesh_type); write_dataset(mesh_group, "filename", filename_); write_dataset(mesh_group, "library", this->library()); if (!options_.empty()) { @@ -453,8 +474,6 @@ void UnstructuredMesh::to_hdf5(hid_t group) const write_dataset(mesh_group, "volumes", volumes); write_dataset(mesh_group, "connectivity", connectivity); write_dataset(mesh_group, "element_types", elem_types); - - close_group(mesh_group); } void UnstructuredMesh::set_length_multiplier(double length_multiplier) @@ -948,17 +967,13 @@ std::pair, vector> RegularMesh::plot( return {axis_lines[0], axis_lines[1]}; } -void RegularMesh::to_hdf5(hid_t group) const +void RegularMesh::to_hdf5_inner(hid_t mesh_group) const { - hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); - write_dataset(mesh_group, "type", "regular"); write_dataset(mesh_group, "dimension", get_x_shape()); write_dataset(mesh_group, "lower_left", lower_left_); write_dataset(mesh_group, "upper_right", upper_right_); write_dataset(mesh_group, "width", width_); - - close_group(mesh_group); } xt::xtensor RegularMesh::count_sites( @@ -1138,16 +1153,12 @@ std::pair, vector> RectilinearMesh::plot( return {axis_lines[0], axis_lines[1]}; } -void RectilinearMesh::to_hdf5(hid_t group) const +void RectilinearMesh::to_hdf5_inner(hid_t mesh_group) const { - hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); - write_dataset(mesh_group, "type", "rectilinear"); write_dataset(mesh_group, "x_grid", grid_[0]); write_dataset(mesh_group, "y_grid", grid_[1]); write_dataset(mesh_group, "z_grid", grid_[2]); - - close_group(mesh_group); } double RectilinearMesh::volume(const MeshIndex& ijk) const @@ -1417,17 +1428,13 @@ std::pair, vector> CylindricalMesh::plot( return {axis_lines[0], axis_lines[1]}; } -void CylindricalMesh::to_hdf5(hid_t group) const +void CylindricalMesh::to_hdf5_inner(hid_t mesh_group) const { - hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); - write_dataset(mesh_group, "type", "cylindrical"); write_dataset(mesh_group, "r_grid", grid_[0]); write_dataset(mesh_group, "phi_grid", grid_[1]); write_dataset(mesh_group, "z_grid", grid_[2]); write_dataset(mesh_group, "origin", origin_); - - close_group(mesh_group); } double CylindricalMesh::volume(const MeshIndex& ijk) const @@ -1733,17 +1740,13 @@ std::pair, vector> SphericalMesh::plot( return {axis_lines[0], axis_lines[1]}; } -void SphericalMesh::to_hdf5(hid_t group) const +void SphericalMesh::to_hdf5_inner(hid_t mesh_group) const { - hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); - write_dataset(mesh_group, "type", SphericalMesh::mesh_type); write_dataset(mesh_group, "r_grid", grid_[0]); write_dataset(mesh_group, "theta_grid", grid_[1]); write_dataset(mesh_group, "phi_grid", grid_[2]); write_dataset(mesh_group, "origin", origin_); - - close_group(mesh_group); } double SphericalMesh::volume(const MeshIndex& ijk) const diff --git a/tests/regression_tests/tally_slice_merge/inputs_true.dat b/tests/regression_tests/tally_slice_merge/inputs_true.dat index 356962f8e4f..7159d833a5a 100644 --- a/tests/regression_tests/tally_slice_merge/inputs_true.dat +++ b/tests/regression_tests/tally_slice_merge/inputs_true.dat @@ -308,7 +308,7 @@ - + 2 2 -50.0 -50.0 50.0 50.0 diff --git a/tests/regression_tests/tally_slice_merge/test.py b/tests/regression_tests/tally_slice_merge/test.py index 090e3144838..aec73979b91 100644 --- a/tests/regression_tests/tally_slice_merge/test.py +++ b/tests/regression_tests/tally_slice_merge/test.py @@ -149,6 +149,9 @@ def _get_results(self, hash_output=False): sum2 = mesh_tally.summation(filter_type=openmc.MeshFilter, filter_bins=[(2, 1), (2, 2)]) + mesh = mesh_tally.find_filter(openmc.MeshFilter).mesh + assert mesh.name == 'mesh' + # Merge the mesh tally slices merge_tally = sum1.merge(sum2) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index ed08be816f5..27322165f64 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -357,6 +357,27 @@ def test_CylindricalMesh_get_indices_at_coords(): assert mesh.get_indices_at_coords([102, 199.1, 299]) == (0, 3, 0) # forth angle quadrant +def test_mesh_name_roundtrip(run_in_tmpdir): + + mesh = openmc.RegularMesh() + mesh.name = 'regular-mesh' + mesh.lower_left = (-1, -1, -1) + mesh.width = (1, 1, 1) + mesh.dimension = (1, 1, 1) + + mesh_filter = openmc.MeshFilter(mesh) + tally = openmc.Tally() + tally.filters = [mesh_filter] + tally.scores = ['flux'] + + openmc.Tallies([tally]).export_to_xml() + + xml_tallies = openmc.Tallies.from_xml() + + mesh = xml_tallies[0].find_filter(openmc.MeshFilter).mesh + assert mesh.name == 'regular-mesh' + + def test_umesh_roundtrip(run_in_tmpdir, request): umesh = openmc.UnstructuredMesh(request.path.parent / 'test_mesh_tets.e', 'moab') umesh.output = True