Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Write and read mesh name attribute #3221

Merged
merged 11 commits into from
Dec 16, 2024
5 changes: 4 additions & 1 deletion docs/source/io_formats/statepoint.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
7 changes: 6 additions & 1 deletion docs/source/io_formats/tallies.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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".
Expand Down
20 changes: 13 additions & 7 deletions include/openmc/mesh.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
//
Expand Down Expand Up @@ -202,7 +207,8 @@ class Mesh {
// 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
};

Expand Down Expand Up @@ -410,7 +416,7 @@ class RegularMesh : public StructuredMesh {
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
//!
Expand Down Expand Up @@ -460,7 +466,7 @@ class RectilinearMesh : public StructuredMesh {
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
//!
Expand Down Expand Up @@ -506,7 +512,7 @@ class CylindricalMesh : public PeriodicStructuredMesh {
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;

Expand Down Expand Up @@ -570,7 +576,7 @@ class SphericalMesh : public PeriodicStructuredMesh {
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]; }
Expand Down Expand Up @@ -632,7 +638,7 @@ class UnstructuredMesh : public Mesh {
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;

Expand Down
92 changes: 49 additions & 43 deletions openmc/mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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()]
Expand Down Expand Up @@ -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'][()]
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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'][()],
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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'][()]
Expand All @@ -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")
Expand Down Expand Up @@ -2442,16 +2449,15 @@ 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:
options = group.attrs['options'].decode()
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],))
Expand All @@ -2478,9 +2484,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)
Expand Down
Loading
Loading