Skip to content

Commit

Permalink
Write and read mesh name attribute (openmc-dev#3221)
Browse files Browse the repository at this point in the history
  • Loading branch information
pshriwise authored Dec 16, 2024
1 parent de8132a commit 775c415
Show file tree
Hide file tree
Showing 8 changed files with 126 additions and 79 deletions.
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 @@ -2444,16 +2451,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 @@ -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)
Expand Down
Loading

0 comments on commit 775c415

Please sign in to comment.