diff --git a/src/adios4dolfinx/backends/adios2/backend.py b/src/adios4dolfinx/backends/adios2/backend.py index 947fc8e..1c9ca2f 100644 --- a/src/adios4dolfinx/backends/adios2/backend.py +++ b/src/adios4dolfinx/backends/adios2/backend.py @@ -31,6 +31,8 @@ def get_default_backend_args(arguments: dict[str, Any] | None) -> dict[str, Any] args = arguments or {} if "engine" not in args.keys(): args["engine"] = "BP4" + if "legacy" not in args.keys(): + args["legacy"] = False # Only used for legacy HDF5 meshtags return args @@ -71,7 +73,7 @@ def write_attributes( filename=filename, mode=adios2.Mode.Append, io_name="AttributeWriter", - **backend_args, + engine=backend_args["engine"], ) as adios_file: adios_file.file.BeginStep() @@ -105,7 +107,7 @@ def read_attributes( adios=adios, filename=filename, mode=adios2.Mode.Read, - **backend_args, + engine=backend_args["engine"], io_name="AttributesReader", ) as adios_file: adios_file.file.BeginStep() @@ -143,7 +145,7 @@ def read_timestamps( adios=adios, filename=filename, mode=adios2.Mode.Read, - **backend_args, + engine=backend_args["engine"], io_name="TimestepReader", ) as adios_file: time_name = f"{function_name}_time" @@ -195,7 +197,8 @@ def write_mesh( filename=filename, mode=mode, comm=comm, - **backend_args, + engine=backend_args["engine"], + io_name=backend_args["io_name"], ) as adios_file: adios_file.file.BeginStep() # Write geometry @@ -473,9 +476,10 @@ def read_meshtags_data( """ adios = adios2.ADIOS(comm) - backend_args = {} if backend_args is None else backend_args + backend_args = get_default_backend_args(backend_args) io_name = backend_args.get("io_name", "MeshTagsReader") - engine = backend_args.get("engine", "BP4") + engine = backend_args["engine"] + legacy = backend_args["legacy"] with ADIOSFile( adios=adios, filename=filename, @@ -483,63 +487,115 @@ def read_meshtags_data( engine=engine, io_name=io_name, ) as adios_file: - # Get mesh cell type - dim_attr_name = f"{name}_dim" - step = 0 - for i in range(adios_file.file.Steps()): - adios_file.file.BeginStep() - if dim_attr_name in adios_file.io.AvailableAttributes().keys(): - step = i - break - adios_file.file.EndStep() - if dim_attr_name not in adios_file.io.AvailableAttributes().keys(): - raise KeyError(f"{dim_attr_name} not found in {filename}") + if not legacy: + # Get mesh cell type + dim_attr_name = f"{name}_dim" + step = 0 + for i in range(adios_file.file.Steps()): + adios_file.file.BeginStep() + if dim_attr_name in adios_file.io.AvailableAttributes().keys(): + step = i + break + adios_file.file.EndStep() + if dim_attr_name not in adios_file.io.AvailableAttributes().keys(): + raise KeyError(f"{dim_attr_name} not found in {filename}") - m_dim = adios_file.io.InquireAttribute(dim_attr_name) - dim = int(m_dim.Data()[0]) + m_dim = adios_file.io.InquireAttribute(dim_attr_name) + dim = int(m_dim.Data()[0]) - # Get mesh tags entites - topology_name = f"{name}_topology" - for i in range(step, adios_file.file.Steps()): - if i > step: - adios_file.file.BeginStep() - if topology_name in adios_file.io.AvailableVariables().keys(): - break - adios_file.file.EndStep() - if topology_name not in adios_file.io.AvailableVariables().keys(): - raise KeyError(f"{topology_name} not found in {filename}") + # Get mesh tags entites + topology_name = f"{name}_topology" + for i in range(step, adios_file.file.Steps()): + if i > step: + adios_file.file.BeginStep() + if topology_name in adios_file.io.AvailableVariables().keys(): + break + adios_file.file.EndStep() + if topology_name not in adios_file.io.AvailableVariables().keys(): + raise KeyError(f"{topology_name} not found in {filename}") + + topology = adios_file.io.InquireVariable(topology_name) + top_shape = topology.Shape() + topology_range = compute_local_range(comm, top_shape[0]) + + topology.SetSelection( + [ + [topology_range[0], 0], + [topology_range[1] - topology_range[0], top_shape[1]], + ] + ) + mesh_entities = np.empty( + (topology_range[1] - topology_range[0], top_shape[1]), dtype=np.int64 + ) + adios_file.file.Get(topology, mesh_entities, adios2.Mode.Deferred) + + # Get mesh tags values + values_name = f"{name}_values" + if values_name not in adios_file.io.AvailableVariables().keys(): + raise KeyError(f"{values_name} not found") + + values = adios_file.io.InquireVariable(values_name) + val_shape = values.Shape() + assert val_shape[0] == top_shape[0] + values.SetSelection([[topology_range[0]], [topology_range[1] - topology_range[0]]]) + tag_values = np.empty( + (topology_range[1] - topology_range[0]), dtype=values.Type().strip("_t") + ) + adios_file.file.Get(values, tag_values, adios2.Mode.Deferred) - topology = adios_file.io.InquireVariable(topology_name) - top_shape = topology.Shape() - topology_range = compute_local_range(comm, top_shape[0]) + adios_file.file.PerformGets() + adios_file.file.EndStep() + else: + # Get mesh cell type + dim_attr_name = f"{name}_dim" + assert adios_file.file.Steps() == 0 + if (ct_key := f"/{name}/topology/celltype") in adios_file.io.AvailableAttributes(): + cell_type = adios_file.io.InquireAttribute(ct_key) + else: + raise ValueError(f"Celltype not found for meshtags {name} in {filename}.") + dim = dolfinx.mesh.cell_dim(dolfinx.mesh.to_type(cell_type.DataString()[0])) + + # Get mesh tags entites + if (top_key := f"/{name}/topology") in adios_file.io.AvailableVariables(): + topology = adios_file.io.InquireVariable(top_key) + else: + raise ValueError(f"Topology not found for meshtags {name} in {filename}.") + + top_shape = topology.Shape() + topology_range = compute_local_range(comm, top_shape[0]) + + topology.SetSelection( + [ + [topology_range[0], 0], + [topology_range[1] - topology_range[0], top_shape[1]], + ] + ) + mesh_entities = np.empty( + (topology_range[1] - topology_range[0], top_shape[1]), dtype=np.int64 + ) + adios_file.file.Get(topology, mesh_entities, adios2.Mode.Deferred) - topology.SetSelection( - [ - [topology_range[0], 0], - [topology_range[1] - topology_range[0], top_shape[1]], - ] - ) - mesh_entities = np.empty( - (topology_range[1] - topology_range[0], top_shape[1]), dtype=np.int64 - ) - adios_file.file.Get(topology, mesh_entities, adios2.Mode.Deferred) + # Get mesh tags values + if (val_key := f"/{name}/values") in adios_file.io.AvailableVariables(): + values = adios_file.io.InquireVariable(val_key) + else: + raise ValueError(f"Values not found for meshtags {name} in {filename}.") - # Get mesh tags values - values_name = f"{name}_values" - if values_name not in adios_file.io.AvailableVariables().keys(): - raise KeyError(f"{values_name} not found") + val_shape = values.Shape() + assert val_shape[0] == top_shape[0] - values = adios_file.io.InquireVariable(values_name) - val_shape = values.Shape() - assert val_shape[0] == top_shape[0] - values.SetSelection([[topology_range[0]], [topology_range[1] - topology_range[0]]]) - tag_values = np.empty((topology_range[1] - topology_range[0]), dtype=np.int32) - adios_file.file.Get(values, tag_values, adios2.Mode.Deferred) + values.SetSelection([[topology_range[0]], [topology_range[1] - topology_range[0]]]) + tag_values = np.empty( + (topology_range[1] - topology_range[0]), dtype=values.Type().strip("_t") + ) + adios_file.file.Get(values, tag_values, adios2.Mode.Deferred) - adios_file.file.PerformGets() - adios_file.file.EndStep() + adios_file.file.PerformGets() + adios_file.file.EndStep() - return MeshTagsData(name=name, values=tag_values, indices=mesh_entities, dim=dim) + return MeshTagsData( + name=name, values=tag_values.astype(np.int32), indices=mesh_entities, dim=dim + ) def read_dofmap( diff --git a/src/adios4dolfinx/backends/h5py/backend.py b/src/adios4dolfinx/backends/h5py/backend.py index 7fd77ff..0ba6a64 100644 --- a/src/adios4dolfinx/backends/h5py/backend.py +++ b/src/adios4dolfinx/backends/h5py/backend.py @@ -13,6 +13,7 @@ from mpi4py import MPI import dolfinx +import h5py import numpy as np import numpy.typing as npt from dolfinx.graph import adjacencylist @@ -23,6 +24,10 @@ read_mode = ReadMode.parallel +# try: +# except ModuleNotFoundError: +# raise ModuleNotFoundError("This backend requires h5py to be installed.") + @contextlib.contextmanager def h5pyfile(h5name, filemode="r", force_serial: bool = False, comm=None): @@ -35,8 +40,6 @@ def h5pyfile(h5name, filemode="r", force_serial: bool = False, comm=None): comm: The MPI communicator """ - import h5py - if comm is None: comm = MPI.COMM_WORLD @@ -58,13 +61,7 @@ def h5pyfile(h5name, filemode="r", force_serial: bool = False, comm=None): def get_default_backend_args(arguments: dict[str, Any] | None) -> dict[str, Any]: - args = arguments or {} - - if arguments: - # Currently no default arguments for h5py backend - # TODO: Pehaps we would like to make this into a warning instead? - raise RuntimeError("Unexpected backend arguments to h5py backend") - + args = arguments or {"legacy": False} # If meshtags is read from legacy return args @@ -384,20 +381,32 @@ def read_meshtags_data( Internal data structure for the mesh tags read from file """ backend_args = get_default_backend_args(backend_args) + legacy = backend_args["legacy"] with h5pyfile(filename, filemode="r", comm=comm, force_serial=False) as h5file: - if "mesh" not in h5file.keys(): - raise KeyError("No mesh found") - mesh = h5file["mesh"] - if "tags" not in mesh.keys(): - raise KeyError("Could not find 'tags' in file, are you sure this is a checkpoint?") - tags = mesh["tags"] - if name not in tags.keys(): - raise KeyError(f"Could not find {name} in '/mesh/tags/' in {filename}") - tag = tags[name] - - dim = tag.attrs["dim"] - topology = tag["Topology"] - values = tag["Values"] + if legacy: + if name not in h5file.keys(): + raise RuntimeError(f"MeshTag {name} not found in {filename}.") + mesh = h5file[name] + topology = mesh["topology"] + cell_type = topology.attrs["celltype"] + if isinstance(cell_type, np.bytes_): + cell_type = cell_type.decode("utf-8") + dim = dolfinx.mesh.cell_dim(dolfinx.mesh.to_type(cell_type)) + values = mesh["values"] + else: + if "mesh" not in h5file.keys(): + raise KeyError("No mesh found") + mesh = h5file["mesh"] + if "tags" not in mesh.keys(): + raise KeyError("Could not find 'tags' in file, are you sure this is a checkpoint?") + tags = mesh["tags"] + if name not in tags.keys(): + raise KeyError(f"Could not find {name} in '/mesh/tags/' in {filename}") + tag = tags[name] + + dim = tag.attrs["dim"] + topology = tag["Topology"] + values = tag["Values"] num_entities_global = topology.shape[0] topology_range = compute_local_range(comm, num_entities_global) indices = topology[slice(*topology_range), :].astype(np.int64) diff --git a/src/adios4dolfinx/backends/pyvista/backend.py b/src/adios4dolfinx/backends/pyvista/backend.py index ba3a418..25b46fb 100644 --- a/src/adios4dolfinx/backends/pyvista/backend.py +++ b/src/adios4dolfinx/backends/pyvista/backend.py @@ -10,7 +10,7 @@ try: import pyvista except ImportError: - raise ModuleNotFoundError("This module requires pyvista") + raise ModuleNotFoundError("The PyVista-backend requires pyvista") from pathlib import Path from mpi4py import MPI @@ -138,7 +138,11 @@ def read_mesh_data( grid = in_data elif isinstance(in_data, pyvista.core.composite.MultiBlock): # To handle multiblock like pvd - pyvista._VTK_SNAKE_CASE_STATE = "allow" + if hasattr(pyvista, "_VTK_SNAKE_CASE_STATE"): + pyvista._VTK_SNAKE_CASE_STATE = "allow" + else: + # Compatibility with 0.47 + pyvista.core.vtk_snake_case._state = "allow" number_of_blocks = in_data.number_of_blocks assert number_of_blocks == 1 b0 = in_data.get_block(0) @@ -203,7 +207,11 @@ def read_point_data( grid = in_data elif isinstance(in_data, pyvista.core.composite.MultiBlock): # To handle multiblock like pvd - pyvista._VTK_SNAKE_CASE_STATE = "allow" + if hasattr(pyvista, "_VTK_SNAKE_CASE_STATE"): + pyvista._VTK_SNAKE_CASE_STATE = "allow" + else: + # Compatibility with 0.47 + pyvista.core.vtk_snake_case._state = "allow" number_of_blocks = in_data.number_of_blocks assert number_of_blocks == 1 b0 = in_data.get_block(0) @@ -217,10 +225,10 @@ def read_point_data( else: num_components = dataset.shape[1] if np.issubdtype(dataset.dtype, np.integer): - gtype = in_data.points.dtype + gtype = grid.points.dtype dataset = dataset.astype(gtype) else: - gtype = in_data.dtype + gtype = dataset.dtype num_components, gtype = comm.bcast((num_components, gtype), root=0) local_range_start = 0 else: @@ -246,7 +254,11 @@ def read_cell_data( grid = in_data elif isinstance(in_data, pyvista.core.composite.MultiBlock): # To handle multiblock like pvd - pyvista._VTK_SNAKE_CASE_STATE = "allow" + if hasattr(pyvista, "_VTK_SNAKE_CASE_STATE"): + pyvista._VTK_SNAKE_CASE_STATE = "allow" + else: + # Compatibility with 0.47 + pyvista.core.vtk_snake_case._state = "allow" number_of_blocks = in_data.number_of_blocks assert number_of_blocks == 1 b0 = in_data.get_block(0) @@ -351,7 +363,11 @@ def read_function_names( grid = in_data elif isinstance(in_data, pyvista.core.composite.MultiBlock): # To handle multiblock like pvd - pyvista._VTK_SNAKE_CASE_STATE = "allow" + if hasattr(pyvista, "_VTK_SNAKE_CASE_STATE"): + pyvista._VTK_SNAKE_CASE_STATE = "allow" + else: + # Compatibility with 0.47 + pyvista.core.vtk_snake_case._state = "allow" number_of_blocks = in_data.number_of_blocks assert number_of_blocks == 1 b0 = in_data.get_block(0) diff --git a/src/adios4dolfinx/backends/xdmf/backend.py b/src/adios4dolfinx/backends/xdmf/backend.py index 98a84f7..362ada8 100644 --- a/src/adios4dolfinx/backends/xdmf/backend.py +++ b/src/adios4dolfinx/backends/xdmf/backend.py @@ -2,7 +2,6 @@ Module that uses DOLFINx/H%py to import XDMF files. """ -import contextlib from pathlib import Path from typing import Any from xml.etree import ElementTree @@ -18,43 +17,11 @@ from adios4dolfinx.utils import check_file_exists, compute_local_range from .. import FileMode, ReadMode +from ..h5py.backend import h5pyfile read_mode = ReadMode.parallel -@contextlib.contextmanager -def h5pyfile(h5name, filemode="r", force_serial: bool = False, comm=None): - """Context manager for opening an HDF5 file with h5py. - - Args: - h5name: The name of the HDF5 file. - filemode: The file mode. - force_serial: Force serial access to the file. - comm: The MPI communicator - - """ - import h5py - - if comm is None: - comm = MPI.COMM_WORLD - - if h5py.h5.get_config().mpi and not force_serial: - h5file = h5py.File(h5name, filemode, driver="mpio", comm=comm) - else: - if comm.size > 1 and not force_serial: - raise ValueError( - f"h5py is not installed with MPI support, while using {comm.size} processes.", - "If you really want to do this, turn on the `force_serial` flag.", - ) - h5file = h5py.File(h5name, filemode) - - try: - yield h5file - finally: - if h5file.id: - h5file.close() - - def extract_function_names_and_timesteps(filename: Path | str) -> dict[str, list[str]]: tree = ElementTree.parse(filename) root = tree.getroot() diff --git a/src/adios4dolfinx/comm_helpers.py b/src/adios4dolfinx/comm_helpers.py index c24ea7c..6a4e2f3 100644 --- a/src/adios4dolfinx/comm_helpers.py +++ b/src/adios4dolfinx/comm_helpers.py @@ -228,7 +228,11 @@ def send_dofs_and_recv_values( # Send sizes to create data structures for receiving from NeighAlltoAllv recv_size = np.zeros_like(source, dtype=np.int32) + recv_size.resize(max(len(recv_size), 1)) # Minimal resize to work with ompi + dest_size.resize(max(len(dest_size), 1)) # Mininal resize to work with ompi dofmap_to_values.Neighbor_alltoall(dest_size, recv_size) + dest_size.resize(len(dest)) + recv_size.resize(len(source)) # Send input dofs to processes holding input array inc_dofs = np.zeros(sum(recv_size), dtype=np.int64) diff --git a/src/adios4dolfinx/utils.py b/src/adios4dolfinx/utils.py index 82b738c..ac1f2bd 100644 --- a/src/adios4dolfinx/utils.py +++ b/src/adios4dolfinx/utils.py @@ -27,56 +27,9 @@ "unroll_dofmap", "compute_insert_position", "unroll_insert_position", - "get_hdf5_version", ] -def get_hdf5_version(): - """Get the HDF5 library version found on the system. - - Note: This function first tries to get the version via h5py. - If h5py is not installed, it attempts to load the HDF5 shared library - directly using ctypes to query the version. - """ - import ctypes - from ctypes import byref, c_uint - - from packaging.version import parse as _v - - try: - import h5py - - return _v(h5py.version.hdf5_version) - except ImportError: - try: - # Try to load default HDF5 library names - # If adios2 is already imported, the symbols might be globally available - try: - libhdf5 = ctypes.CDLL("libhdf5.so") - except OSError: - try: - libhdf5 = ctypes.CDLL("libhdf5.dylib") # MacOS - except OSError: - print("\n[Method 2] Could not load libhdf5 shared library directly.") - libhdf5 = None - - if libhdf5: - major = c_uint() - minor = c_uint() - rel = c_uint() - - # Call the C function - libhdf5.H5get_libversion(byref(major), byref(minor), byref(rel)) - return _v(f"{major.value}.{minor.value}.{rel.value}") - except Exception: - pass - - except Exception: - pass - - raise RuntimeError("Failed to get HDF5 version") - - def check_file_exists(filename: Path | str): """Check if file exists.""" if not Path(filename).exists(): diff --git a/tests/create_legacy_data.py b/tests/create_legacy_data.py index 74b8b7d..e2aacdf 100644 --- a/tests/create_legacy_data.py +++ b/tests/create_legacy_data.py @@ -40,12 +40,20 @@ def create_reference_data( tt = dolfin.Function(dolfin.FunctionSpace(mesh, "DG", 0)) tt.interpolate(dolfin.Expression(("x[0]+x[1]"), degree=0)) + ff = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 1) + + class Boundary(dolfin.SubDomain): + def inside(self, x, on_boundary): + return on_boundary + + Boundary().mark(ff, 2) + with dolfin.HDF5File(mesh.mpi_comm(), str(h5_file), "w") as hdf: hdf.write(mesh, mesh_name) hdf.write(v0, function_name) hdf.write(w0, function_name_vec) hdf.write(tt, function_name + "DG") - + hdf.write(ff, "facet_tags") with dolfin.XDMFFile(mesh.mpi_comm(), str(xdmf_file)) as xdmf: xdmf.write(mesh) xdmf.write_checkpoint(v0, function_name, 0, dolfin.XDMFFile.Encoding.HDF5, append=True) diff --git a/tests/test_legacy_readers.py b/tests/test_legacy_readers.py index 29762bc..2eb10e9 100644 --- a/tests/test_legacy_readers.py +++ b/tests/test_legacy_readers.py @@ -23,6 +23,7 @@ read_function_names, read_mesh, read_mesh_from_legacy_h5, + read_meshtags, read_point_data, ) @@ -81,6 +82,26 @@ def test_legacy_function(backend): if not path.exists(): pytest.skip(f"{path} does not exist") mesh = read_mesh_from_legacy_h5(path, comm, "/mesh", backend=backend) + ff = read_meshtags( + path, mesh, "facet_tags", backend_args={"legacy": True, "engine": "HDF5"}, backend=backend + ) + fdim = mesh.topology.dim - 1 + assert ff.dim == fdim + boundary_facets = ff.find(2) + interior_facets = ff.find(1) + mesh.topology.create_connectivity(fdim, fdim + 1) + num_facets_local = ( + mesh.topology.index_map(fdim).size_local + mesh.topology.index_map(fdim).num_ghosts + ) + true_exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology) + assert len(boundary_facets) == len(true_exterior_facets) + assert np.isin(boundary_facets, true_exterior_facets).all() + true_interior = np.ones(num_facets_local, dtype=np.int8) + true_interior[true_exterior_facets] = 0 + interior_marker = np.flatnonzero(true_interior) + assert len(interior_marker) == len(interior_facets) + assert np.isin(interior_facets, interior_marker).all() + V = dolfinx.fem.functionspace(mesh, ("DG", 2)) u = ufl.TrialFunction(V) v = ufl.TestFunction(V) diff --git a/tests/test_mesh_writer.py b/tests/test_mesh_writer.py index 52df0df..83bdb2b 100644 --- a/tests/test_mesh_writer.py +++ b/tests/test_mesh_writer.py @@ -7,7 +7,53 @@ import pytest import ufl -from adios4dolfinx import FileMode, read_mesh, utils, write_mesh +from adios4dolfinx import FileMode, read_mesh, write_mesh + + +def get_hdf5_version(): + """Get the HDF5 library version found on the system. + + Note: This function first tries to get the version via h5py. + If h5py is not installed, it attempts to load the HDF5 shared library + directly using ctypes to query the version. + """ + import ctypes + from ctypes import byref, c_uint + + from packaging.version import parse as _v + + try: + import h5py + + return _v(h5py.version.hdf5_version) + except ImportError: + try: + # Try to load default HDF5 library names + # If adios2 is already imported, the symbols might be globally available + try: + libhdf5 = ctypes.CDLL("libhdf5.so") + except OSError: + try: + libhdf5 = ctypes.CDLL("libhdf5.dylib") # MacOS + except OSError: + print("\n[Method 2] Could not load libhdf5 shared library directly.") + libhdf5 = None + + if libhdf5: + major = c_uint() + minor = c_uint() + rel = c_uint() + + # Call the C function + libhdf5.H5get_libversion(byref(major), byref(minor), byref(rel)) + return _v(f"{major.value}.{minor.value}.{rel.value}") + except Exception: + pass + + except Exception: + pass + + raise RuntimeError("Failed to get HDF5 version") @pytest.mark.parametrize( @@ -26,7 +72,7 @@ def test_mesh_read_writer(backend, encoder, suffix, ghost_mode, tmp_path, store_partition): N = 7 - if backend == "adios2" and encoder == "HDF5" and utils.get_hdf5_version().major >= 2: + if backend == "adios2" and encoder == "HDF5" and get_hdf5_version().major >= 2: pytest.skip("HDF5 version >= 2 is not supported due to ADIOS2 limitations.") try: importlib.import_module(backend)