mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 14:35:27 -04:00
Adding vtkhdf option to write vtk data (#3252)
Co-authored-by: shimwell <mail@jshimwell.com> Co-authored-by: Jonathan Shimwell <drshimwell@gmail.com> Co-authored-by: rherrero-pf <156206440+rherrero-pf@users.noreply.github.com> Co-authored-by: Patrick Shriwise <pshriwise@gmail.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Jon Shimwell
GitHub
parent
bd76fc0566
commit
e5c7d0ca88
4 changed files with 377 additions and 18 deletions
166
openmc/mesh.py
166
openmc/mesh.py
|
|
@ -5,6 +5,7 @@ from collections.abc import Iterable, Sequence, Mapping
|
|||
from functools import wraps
|
||||
from math import pi, sqrt, atan2
|
||||
from numbers import Integral, Real
|
||||
from pathlib import Path
|
||||
from typing import Protocol
|
||||
|
||||
import h5py
|
||||
|
|
@ -2443,6 +2444,7 @@ class UnstructuredMesh(MeshBase):
|
|||
_UNSUPPORTED_ELEM = -1
|
||||
_LINEAR_TET = 0
|
||||
_LINEAR_HEX = 1
|
||||
_VTK_TETRA = 10
|
||||
|
||||
def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None,
|
||||
name: str = '', length_multiplier: float = 1.0,
|
||||
|
|
@ -2652,7 +2654,8 @@ class UnstructuredMesh(MeshBase):
|
|||
warnings.warn(
|
||||
"The 'UnstructuredMesh.write_vtk_mesh' method has been renamed "
|
||||
"to 'write_data_to_vtk' and will be removed in a future version "
|
||||
" of OpenMC.", FutureWarning
|
||||
" of OpenMC.",
|
||||
FutureWarning,
|
||||
)
|
||||
self.write_data_to_vtk(**kwargs)
|
||||
|
||||
|
|
@ -2670,9 +2673,10 @@ class UnstructuredMesh(MeshBase):
|
|||
Parameters
|
||||
----------
|
||||
filename : str or pathlib.Path
|
||||
Name of the VTK file to write. If the filename ends in '.vtu' then a
|
||||
binary VTU format file will be written, if the filename ends in
|
||||
'.vtk' then a legacy VTK file will be written.
|
||||
Name of the VTK file to write. If the filename ends in '.vtkhdf'
|
||||
then a VTKHDF format file will be written. If the filename ends in
|
||||
'.vtu' then a binary VTU format file will be written. If the
|
||||
filename ends in '.vtk' then a legacy VTK file will be written.
|
||||
datasets : dict
|
||||
Dictionary whose keys are the data labels and values are numpy
|
||||
appropriately sized arrays of the data
|
||||
|
|
@ -2680,6 +2684,35 @@ class UnstructuredMesh(MeshBase):
|
|||
Whether or not to normalize the data by the volume of the mesh
|
||||
elements
|
||||
"""
|
||||
|
||||
if Path(filename).suffix == ".vtkhdf":
|
||||
|
||||
self._write_data_to_vtk_hdf5_format(
|
||||
filename=filename,
|
||||
datasets=datasets,
|
||||
volume_normalization=volume_normalization,
|
||||
)
|
||||
|
||||
elif Path(filename).suffix == ".vtk" or Path(filename).suffix == ".vtu":
|
||||
|
||||
self._write_data_to_vtk_ascii_format(
|
||||
filename=filename,
|
||||
datasets=datasets,
|
||||
volume_normalization=volume_normalization,
|
||||
)
|
||||
|
||||
else:
|
||||
raise ValueError(
|
||||
"Unsupported file extension, The filename must end with "
|
||||
"'.vtkhdf', '.vtu' or '.vtk'"
|
||||
)
|
||||
|
||||
def _write_data_to_vtk_ascii_format(
|
||||
self,
|
||||
filename: PathLike | None = None,
|
||||
datasets: dict | None = None,
|
||||
volume_normalization: bool = True,
|
||||
):
|
||||
from vtkmodules.util import numpy_support
|
||||
from vtkmodules import vtkCommonCore
|
||||
from vtkmodules import vtkCommonDataModel
|
||||
|
|
@ -2687,9 +2720,7 @@ class UnstructuredMesh(MeshBase):
|
|||
from vtkmodules import vtkIOXML
|
||||
|
||||
if self.connectivity is None or self.vertices is None:
|
||||
raise RuntimeError(
|
||||
"This mesh has not been loaded from a statepoint file."
|
||||
)
|
||||
raise RuntimeError("This mesh has not been loaded from a statepoint file.")
|
||||
|
||||
if filename is None:
|
||||
filename = f"mesh_{self.id}.vtk"
|
||||
|
|
@ -2771,29 +2802,128 @@ class UnstructuredMesh(MeshBase):
|
|||
|
||||
writer.Write()
|
||||
|
||||
def _write_data_to_vtk_hdf5_format(
|
||||
self,
|
||||
filename: PathLike | None = None,
|
||||
datasets: dict | None = None,
|
||||
volume_normalization: bool = True,
|
||||
):
|
||||
def append_dataset(dset, array):
|
||||
"""Convenience function to append data to an HDF5 dataset"""
|
||||
origLen = dset.shape[0]
|
||||
dset.resize(origLen + array.shape[0], axis=0)
|
||||
dset[origLen:] = array
|
||||
|
||||
if self.library != "moab":
|
||||
raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")
|
||||
|
||||
# the self.connectivity contains arrays of length 8 to support hex
|
||||
# elements as well, in the case of tetrahedra mesh elements, the
|
||||
# last 4 values are -1 and are removed
|
||||
trimmed_connectivity = []
|
||||
for cell in self.connectivity:
|
||||
# Find the index of the first -1 value, if any
|
||||
first_negative_index = np.where(cell == -1)[0]
|
||||
if first_negative_index.size > 0:
|
||||
# Slice the array up to the first -1 value
|
||||
trimmed_connectivity.append(cell[: first_negative_index[0]])
|
||||
else:
|
||||
# No -1 values, append the whole cell
|
||||
trimmed_connectivity.append(cell)
|
||||
trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten()
|
||||
|
||||
# MOAB meshes supports tet elements only so we know it has 4 points per cell
|
||||
points_per_cell = 4
|
||||
|
||||
# offsets are the indices of the first point of each cell in the array of points
|
||||
offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell)
|
||||
|
||||
for name, data in datasets.items():
|
||||
if data.shape != self.dimension:
|
||||
raise ValueError(
|
||||
f'Cannot apply dataset "{name}" with '
|
||||
f"shape {data.shape} to mesh {self.id} "
|
||||
f"with dimensions {self.dimension}"
|
||||
)
|
||||
|
||||
with h5py.File(filename, "w") as f:
|
||||
|
||||
root = f.create_group("VTKHDF")
|
||||
vtk_file_format_version = (2, 1)
|
||||
root.attrs["Version"] = vtk_file_format_version
|
||||
ascii_type = "UnstructuredGrid".encode("ascii")
|
||||
root.attrs.create(
|
||||
"Type",
|
||||
ascii_type,
|
||||
dtype=h5py.string_dtype("ascii", len(ascii_type)),
|
||||
)
|
||||
|
||||
# create hdf5 file structure
|
||||
root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8")
|
||||
root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8")
|
||||
root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f")
|
||||
root.create_dataset(
|
||||
"NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8"
|
||||
)
|
||||
root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8")
|
||||
root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8")
|
||||
root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8")
|
||||
|
||||
append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)]))
|
||||
append_dataset(root["Points"], self.vertices)
|
||||
append_dataset(
|
||||
root["NumberOfConnectivityIds"],
|
||||
np.array([len(trimmed_connectivity)]),
|
||||
)
|
||||
append_dataset(root["Connectivity"], trimmed_connectivity)
|
||||
append_dataset(root["NumberOfCells"], np.array([self.n_elements]))
|
||||
append_dataset(root["Offsets"], offsets)
|
||||
|
||||
append_dataset(
|
||||
root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8")
|
||||
)
|
||||
|
||||
cell_data_group = root.create_group("CellData")
|
||||
|
||||
for name, data in datasets.items():
|
||||
|
||||
cell_data_group.create_dataset(
|
||||
name, (0,), maxshape=(None,), dtype="float64", chunks=True
|
||||
)
|
||||
|
||||
if volume_normalization:
|
||||
data /= self.volumes
|
||||
append_dataset(cell_data_group[name], data)
|
||||
|
||||
@classmethod
|
||||
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:
|
||||
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, name=name, 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'][()]
|
||||
vol_data = group["volumes"][()]
|
||||
mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
|
||||
mesh.n_elements = mesh.volumes.size
|
||||
|
||||
vertices = group['vertices'][()]
|
||||
vertices = group["vertices"][()]
|
||||
mesh._vertices = vertices.reshape((-1, 3))
|
||||
connectivity = group['connectivity'][()]
|
||||
connectivity = group["connectivity"][()]
|
||||
mesh._connectivity = connectivity.reshape((-1, 8))
|
||||
mesh._element_types = group['element_types'][()]
|
||||
mesh._element_types = group["element_types"][()]
|
||||
|
||||
if 'length_multiplier' in group:
|
||||
mesh.length_multiplier = group['length_multiplier'][()]
|
||||
if "length_multiplier" in group:
|
||||
mesh.length_multiplier = group["length_multiplier"][()]
|
||||
|
||||
return mesh
|
||||
|
||||
|
|
@ -2812,7 +2942,7 @@ class UnstructuredMesh(MeshBase):
|
|||
|
||||
element.set("library", self._library)
|
||||
if self.options is not None:
|
||||
element.set('options', self.options)
|
||||
element.set("options", self.options)
|
||||
subelement = ET.SubElement(element, "filename")
|
||||
subelement.text = str(self.filename)
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,159 @@
|
|||
# vtk DataFile Version 2.0
|
||||
made_with_cad_to_dagmc_package, Created by Gmsh 4.12.1
|
||||
ASCII
|
||||
DATASET UNSTRUCTURED_GRID
|
||||
POINTS 14 double
|
||||
-0.5 -0.5 0.5
|
||||
-0.5 -0.5 -0.5
|
||||
-0.5 0.5 0.5
|
||||
-0.5 0.5 -0.5
|
||||
0.5 -0.5 0.5
|
||||
0.5 -0.5 -0.5
|
||||
0.5 0.5 0.5
|
||||
0.5 0.5 -0.5
|
||||
-0.5 0 0
|
||||
0.5 0 0
|
||||
0 -0.5 0
|
||||
0 0.5 0
|
||||
0 0 -0.5
|
||||
0 0 0.5
|
||||
|
||||
CELLS 68 268
|
||||
1 0
|
||||
1 1
|
||||
1 2
|
||||
1 3
|
||||
1 4
|
||||
1 5
|
||||
1 6
|
||||
1 7
|
||||
2 1 0
|
||||
2 0 2
|
||||
2 3 2
|
||||
2 1 3
|
||||
2 5 4
|
||||
2 4 6
|
||||
2 7 6
|
||||
2 5 7
|
||||
2 1 5
|
||||
2 0 4
|
||||
2 3 7
|
||||
2 2 6
|
||||
3 1 0 8
|
||||
3 0 2 8
|
||||
3 3 1 8
|
||||
3 2 3 8
|
||||
3 5 9 4
|
||||
3 4 9 6
|
||||
3 7 9 5
|
||||
3 6 9 7
|
||||
3 0 1 10
|
||||
3 4 0 10
|
||||
3 1 5 10
|
||||
3 5 4 10
|
||||
3 2 11 3
|
||||
3 6 11 2
|
||||
3 3 11 7
|
||||
3 7 11 6
|
||||
3 1 3 12
|
||||
3 5 1 12
|
||||
3 3 7 12
|
||||
3 7 5 12
|
||||
3 0 13 2
|
||||
3 4 13 0
|
||||
3 2 13 6
|
||||
3 6 13 4
|
||||
4 13 8 12 10
|
||||
4 11 8 12 13
|
||||
4 12 11 13 9
|
||||
4 10 12 13 9
|
||||
4 12 3 11 7
|
||||
4 13 2 8 0
|
||||
4 8 12 1 3
|
||||
4 11 3 8 2
|
||||
4 10 8 1 0
|
||||
4 0 10 13 4
|
||||
4 1 12 10 5
|
||||
4 11 2 13 6
|
||||
4 4 9 13 6
|
||||
4 6 9 11 7
|
||||
4 10 9 4 5
|
||||
4 7 9 12 5
|
||||
4 3 8 12 11
|
||||
4 1 12 8 10
|
||||
4 8 11 2 13
|
||||
4 10 13 8 0
|
||||
4 13 10 9 4
|
||||
4 11 13 9 6
|
||||
4 12 11 9 7
|
||||
4 9 10 12 5
|
||||
|
||||
CELL_TYPES 68
|
||||
1
|
||||
1
|
||||
1
|
||||
1
|
||||
1
|
||||
1
|
||||
1
|
||||
1
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
3
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
5
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
10
|
||||
|
|
@ -2,6 +2,7 @@ from math import pi
|
|||
from tempfile import TemporaryDirectory
|
||||
from pathlib import Path
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
from scipy.stats import chi2
|
||||
import pytest
|
||||
|
|
@ -486,6 +487,74 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type):
|
|||
with pytest.raises(ValueError, match='Cannot apply dataset "mean"'):
|
||||
simple_umesh.write_data_to_vtk(datasets={'mean': ref_data[:-2]}, filename=filename)
|
||||
|
||||
|
||||
@pytest.mark.skipif(not openmc.lib._dagmc_enabled(), reason="DAGMC not enabled.")
|
||||
def test_write_vtkhdf(request, run_in_tmpdir):
|
||||
"""Performs a minimal UnstructuredMesh simulation, reads in the resulting
|
||||
statepoint file and writes the mesh data to vtk and vtkhdf files. It is
|
||||
necessary to read in the unstructured mesh from a statepoint file to ensure
|
||||
it has all the required attributes
|
||||
"""
|
||||
model = openmc.Model()
|
||||
|
||||
surf1 = openmc.Sphere(r=1000.0, boundary_type="vacuum")
|
||||
cell1 = openmc.Cell(region=-surf1)
|
||||
model.geometry = openmc.Geometry([cell1])
|
||||
|
||||
umesh = openmc.UnstructuredMesh(
|
||||
request.path.parent / "test_mesh_dagmc_tets.vtk",
|
||||
"moab",
|
||||
mesh_id = 1
|
||||
)
|
||||
mesh_filter = openmc.MeshFilter(umesh)
|
||||
|
||||
# Create flux mesh tally to score alpha production
|
||||
mesh_tally = openmc.Tally(name="test_tally")
|
||||
mesh_tally.filters = [mesh_filter]
|
||||
mesh_tally.scores = ["flux"]
|
||||
|
||||
model.tallies = [mesh_tally]
|
||||
|
||||
model.settings.run_mode = "fixed source"
|
||||
model.settings.batches = 2
|
||||
model.settings.particles = 10
|
||||
|
||||
statepoint_file = model.run()
|
||||
|
||||
with openmc.StatePoint(statepoint_file) as statepoint:
|
||||
my_tally = statepoint.get_tally(name="test_tally")
|
||||
|
||||
umesh_from_sp = statepoint.meshes[umesh.id]
|
||||
|
||||
datasets={
|
||||
"mean": my_tally.mean.flatten(),
|
||||
"std_dev": my_tally.std_dev.flatten()
|
||||
}
|
||||
|
||||
umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtkhdf")
|
||||
umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtk")
|
||||
|
||||
with pytest.raises(ValueError, match="Unsupported file extension"):
|
||||
# Supported file extensions are vtk or vtkhdf, not hdf5, so this should raise an error
|
||||
umesh_from_sp.write_data_to_vtk(
|
||||
datasets=datasets,
|
||||
filename="test_mesh.hdf5",
|
||||
)
|
||||
with pytest.raises(ValueError, match="Cannot apply dataset"):
|
||||
# The shape of the data should match the shape of the mesh, so this should raise an error
|
||||
umesh_from_sp.write_data_to_vtk(
|
||||
datasets={'incorrectly_shaped_data': np.array(([1,2,3]))},
|
||||
filename="test_mesh_incorrect_shape.vtkhdf",
|
||||
)
|
||||
|
||||
assert Path("test_mesh.vtk").exists()
|
||||
assert Path("test_mesh.vtkhdf").exists()
|
||||
|
||||
# just ensure we can open the file without error
|
||||
with h5py.File("test_mesh.vtkhdf", "r"):
|
||||
...
|
||||
|
||||
|
||||
def test_mesh_get_homogenized_materials():
|
||||
"""Test the get_homogenized_materials method"""
|
||||
# Simple model with 1 cm of Fe56 next to 1 cm of H1
|
||||
|
|
|
|||
1
tests/unit_tests/test_mesh_dagmc_tets.vtk
Symbolic link
1
tests/unit_tests/test_mesh_dagmc_tets.vtk
Symbolic link
|
|
@ -0,0 +1 @@
|
|||
../regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk
|
||||
Loading…
Add table
Add a link
Reference in a new issue