Write and read mesh name attribute (#3221)

This commit is contained in:
Patrick Shriwise 2024-12-16 08:40:02 -06:00 committed by GitHub
parent de8132a5a4
commit 775c415122
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 126 additions and 79 deletions

View file

@ -72,7 +72,10 @@ The current version of the statepoint file format is 18.1.
**/tallies/meshes/mesh <uid>/**
: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

View file

@ -109,7 +109,7 @@ The ``<tally>`` 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 ``<mesh>``. 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".

View file

@ -132,13 +132,18 @@ public:
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 @@ public:
// Data members
xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
xt::xtensor<double, 1> 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 @@ public:
std::pair<vector<double>, vector<double>> 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 @@ public:
std::pair<vector<double>, vector<double>> 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 @@ public:
std::pair<vector<double>, vector<double>> 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 @@ public:
std::pair<vector<double>, vector<double>> 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 @@ public:
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& 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;

View file

@ -99,21 +99,40 @@ class MeshBase(IDManagerMixin, ABC):
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 @@ class MeshBase(IDManagerMixin, ABC):
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 @@ class RegularMesh(StructuredMesh):
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 @@ class RegularMesh(StructuredMesh):
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 @@ class RegularMesh(StructuredMesh):
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 @@ class RectilinearMesh(StructuredMesh):
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 @@ class RectilinearMesh(StructuredMesh):
"""
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 @@ class CylindricalMesh(StructuredMesh):
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 @@ class CylindricalMesh(StructuredMesh):
"""
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 @@ class SphericalMesh(StructuredMesh):
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 @@ class SphericalMesh(StructuredMesh):
"""
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 @@ class UnstructuredMesh(MeshBase):
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 @@ class UnstructuredMesh(MeshBase):
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 @@ class UnstructuredMesh(MeshBase):
"""
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)

View file

@ -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::MaterialVolume> 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<double>, vector<double>> 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<double, 1> RegularMesh::count_sites(
@ -1138,16 +1153,12 @@ std::pair<vector<double>, vector<double>> 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<double>, vector<double>> 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<double>, vector<double>> 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

View file

@ -308,7 +308,7 @@
</settings>
<tallies>
<!--mesh-->
<mesh id="1">
<mesh id="1" name="mesh">
<dimension>2 2</dimension>
<lower_left>-50.0 -50.0</lower_left>
<upper_right>50.0 50.0</upper_right>

View file

@ -149,6 +149,9 @@ class TallySliceMergeTestHarness(PyAPITestHarness):
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)

View file

@ -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