feat(vtkhdf compression): implemented compression, bmarks and docs.

benchmarks/report_vtkhdf_compression.py

@@ -0,0 +1,137 @@
+"""Summarize rank subdirectories produced by vtkhdf_compression.py.
+
+Usage: python benchmarks/report_vtkhdf_compression.py output/vtkhdf-compression
+"""
+
+import argparse
+import csv
+import json
+import os
+import statistics
+from pathlib import Path
+
+os.environ.setdefault("MPLCONFIGDIR", "/tmp/io4dolfinx-matplotlib")
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+MODES = ["uncompressed", "gzip", "gzip_shuffle"]
+LABELS = ["No compression", "Gzip level 4", "Gzip level 4 + shuffle"]
+COLORS = ["#697586", "#3178b8", "#c45525"]
+
+
+def summarize(root: Path):
+    reports = [json.loads(path.read_text()) for path in sorted(root.glob("ranks-*/results.json"))]
+    if not reports:
+        raise ValueError(f"No rank reports in {root}")
+    raw = [row for report in reports for row in report["results"]]
+    sizes = sorted({row["n"] for row in raw})
+    ranks = sorted({row["ranks"] for row in raw})
+    rows = []
+    for n in sizes:
+        for rank in ranks:
+            for mode in MODES:
+                samples = [
+                    r for r in raw if r["n"] == n and r["ranks"] == rank and r["mode"] == mode
+                ]
+                if not samples:
+                    raise ValueError(f"Missing n={n}, ranks={rank}, mode={mode}")
+                row = dict(n=n, cells=2 * n * n, ranks=rank, mode=mode, samples=len(samples))
+                for metric in (
+                    "bytes",
+                    "write_seconds",
+                    "read_seconds",
+                    "sync_seconds",
+                    "mesh_write_seconds",
+                    "point_write_seconds",
+                    "cell_write_seconds",
+                    "mesh_read_seconds",
+                    "field_read_seconds",
+                ):
+                    values = sorted(r[metric] for r in samples)
+                    row[metric] = statistics.median(values)
+                    if metric in ("write_seconds", "read_seconds"):
+                        quartiles = statistics.quantiles(values, n=4, method="inclusive")
+                        row[f"{metric}_p25"] = quartiles[0]
+                        row[f"{metric}_p75"] = quartiles[2]
+                        row[f"{metric}_min"] = min(values)
+                        row[f"{metric}_max"] = max(values)
+                rows.append(row)
+    for row in rows:
+        baseline = next(
+            r
+            for r in rows
+            if r["n"] == row["n"] and r["ranks"] == row["ranks"] and r["mode"] == "uncompressed"
+        )
+        row["size_ratio"] = row["bytes"] / baseline["bytes"]
+        row["write_time_ratio"] = row["write_seconds"] / baseline["write_seconds"]
+        row["read_time_ratio"] = row["read_seconds"] / baseline["read_seconds"]
+    with (root / "summary.csv").open("w", newline="") as stream:
+        writer = csv.DictWriter(stream, fieldnames=list(rows[0]))
+        writer.writeheader()
+        writer.writerows(rows)
+    (root / "all_results.json").write_text(json.dumps(reports, indent=2) + "\n")
+
+    plt.rcParams.update({"font.size": 10, "axes.spines.top": False, "axes.spines.right": False})
+    fig, axes = plt.subplots(len(sizes), 3, figsize=(15, 3.5 * len(sizes)), squeeze=False)
+    metrics = [
+        ("write_seconds", "Write time incl. fsync (ms)", 1000),
+        ("read_seconds", "Read time, warm cache (ms)", 1000),
+        ("bytes", "File size (MB, decimal)", 1e-6),
+    ]
+    for i, n in enumerate(sizes):
+        for j, (metric, label, scale) in enumerate(metrics):
+            ax = axes[i, j]
+            for mode, legend, color in zip(MODES, LABELS, COLORS):
+                selected = [
+                    next(
+                        r for r in rows if r["n"] == n and r["ranks"] == rank and r["mode"] == mode
+                    )
+                    for rank in ranks
+                ]
+                y = [r[metric] * scale for r in selected]
+                error = None
+                if metric != "bytes":
+                    error = [
+                        [(r[metric] - r[f"{metric}_p25"]) * scale for r in selected],
+                        [(r[f"{metric}_p75"] - r[metric]) * scale for r in selected],
+                    ]
+                ax.errorbar(ranks, y, yerr=error, marker="o", capsize=3, color=color, label=legend)
+            ax.set_xscale("log", base=2)
+            ax.set_xticks(ranks, [str(rank) for rank in ranks])
+            ax.set_xlabel("MPI ranks")
+            ax.set_ylabel(label)
+            ax.set_title(f"{2 * n * n:,} triangles ({n} × {n} mesh)")
+            ax.grid(True, alpha=0.2)
+            ax.set_ylim(bottom=0)
+    handles, labels = axes[0, 0].get_legend_handles_labels()
+    fig.legend(handles, labels, loc="upper center", ncol=3, bbox_to_anchor=(0.5, 0.965))
+    fig.suptitle("VTKHDF: compression, shuffle, and MPI scaling", fontsize=17, y=0.995)
+    fig.text(
+        0.5,
+        0.005,
+        "Median of 5 measured runs; bars show the middle 50%. One warm-up per mode. "
+        "Smooth scalar/vector fields; automatic chunks; ext4 storage.",
+        ha="center",
+        fontsize=10,
+    )
+    fig.tight_layout(rect=(0, 0.025, 1, 0.94))
+    fig.savefig(root / "comparison.png", dpi=160)
+    fig.savefig(root / "comparison.svg")
+    plt.close(fig)
+    for row in rows:
+        print(
+            f"n={row['n']:3} ranks={row['ranks']} {row['mode']:14} "
+            f"{row['bytes'] / 1e6:7.3f} MB "
+            f"write={row['write_seconds'] * 1000:8.2f} ms "
+            f"read={row['read_seconds'] * 1000:8.2f} ms "
+            f"write/base={row['write_time_ratio']:.2f}"
+        )
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("root", type=Path)
+    summarize(parser.parse_args().root)

benchmarks/vtkhdf_compression.py

@@ -0,0 +1,143 @@
+"""Compare uncompressed VTKHDF, gzip, and gzip+shuffle under MPI.
+
+Run from the repository root, for example:
+    mpiexec -n 2 python benchmarks/vtkhdf_compression.py --output /tmp/vtkhdf-mpi
+
+Uses smooth scalar/vector fields on triangular unit-square meshes. Timings include
+open/close, mesh repartitioning, and field interpolation, and report the slowest
+rank. Reads use warm filesystem caches; writes close the file without an explicit
+fsync unless --sync is supplied. Results describe this workload and filesystem,
+not general speed guarantees.
+Each mode is warmed up before measurement; measured mode order rotates each round.
+Output files are overwritten between repetitions to bound disk usage.
+"""
+
+import argparse
+import json
+import os
+import platform
+from pathlib import Path
+
+from mpi4py import MPI
+
+import dolfinx
+import h5py
+import numpy as np
+from dolfinx.fem import Function, functionspace
+from dolfinx.mesh import create_unit_square
+
+import io4dolfinx
+
+MODES = {
+    "uncompressed": None,
+    "gzip": {"compression": "gzip", "compression_opts": 4, "shuffle": False},
+    "gzip_shuffle": {"compression": "gzip", "compression_opts": 4, "shuffle": True},
+}
+
+
+def measure(filename, mesh, point, cell, options, sync):
+    comm = mesh.comm
+    comm.barrier()
+    start = MPI.Wtime()
+    io4dolfinx.write_mesh(filename, mesh, backend="vtkhdf", backend_args=options)
+    mesh_written = MPI.Wtime()
+    io4dolfinx.write_point_data(
+        filename, point, time=0.0, mode=io4dolfinx.FileMode.append, backend_args=options
+    )
+    point_written = MPI.Wtime()
+    io4dolfinx.write_cell_data(
+        filename, cell, time=0.0, mode=io4dolfinx.FileMode.append, backend_args=options
+    )
+    cell_written = MPI.Wtime()
+    if sync:
+        comm.barrier()
+        if comm.rank == 0:
+            with filename.open("rb") as stream:
+                os.fsync(stream.fileno())
+        comm.barrier()
+    write_finished = MPI.Wtime()
+    write_times = {
+        "write_seconds": write_finished - start,
+        "mesh_write_seconds": mesh_written - start,
+        "point_write_seconds": point_written - mesh_written,
+        "cell_write_seconds": cell_written - point_written,
+        "sync_seconds": write_finished - cell_written,
+    }
+    # Reduce outside the timed region; retain the critical-path time across ranks.
+    timings = {key: comm.allreduce(value, op=MPI.MAX) for key, value in write_times.items()}
+    comm.barrier()
+    start = MPI.Wtime()
+    restored = io4dolfinx.read_mesh(filename, comm, backend="vtkhdf")
+    mesh_read = MPI.Wtime()
+    io4dolfinx.read_point_data(filename, point.name, restored, 0.0, backend="vtkhdf")
+    io4dolfinx.read_cell_data(filename, cell.name, restored, 0.0, backend="vtkhdf")
+    fields_read = MPI.Wtime()
+    read_times = {
+        "read_seconds": fields_read - start,
+        "mesh_read_seconds": mesh_read - start,
+        "field_read_seconds": fields_read - mesh_read,
+    }
+    timings.update({key: comm.allreduce(value, op=MPI.MAX) for key, value in read_times.items()})
+    return timings
+
+
+def benchmark(output: Path, sizes: list[int], repeats: int, warmups: int, sync: bool):
+    comm = MPI.COMM_WORLD
+    if comm.rank == 0:
+        output.mkdir(parents=True, exist_ok=True)
+    comm.barrier()
+    results = []
+    for n in sizes:
+        mesh = create_unit_square(comm, n, n)
+        point = Function(functionspace(mesh, ("Lagrange", 1, (3,))), name="velocity")
+        point.interpolate(lambda x: np.vstack((np.sin(x[0]), np.cos(x[1]), x[0] * x[1])))
+        cell = Function(functionspace(mesh, ("DG", 0)), name="pressure")
+        cell.interpolate(lambda x: x[0] + 2 * x[1])
+        for repeat in range(-warmups, repeats):
+            modes = list(MODES)
+            shift = max(repeat, 0) % len(modes)
+            modes = modes[shift:] + modes[:shift]
+            for mode in modes:
+                filename = output / f"n{n}-{mode}.vtkhdf"
+                timings = measure(filename, mesh, point, cell, MODES[mode], sync)
+                if repeat >= 0 and comm.rank == 0:
+                    row = dict(
+                        n=n,
+                        cells=2 * n * n,
+                        ranks=comm.size,
+                        mode=mode,
+                        repeat=repeat,
+                        bytes=filename.stat().st_size,
+                        **timings,
+                    )
+                    results.append(row)
+                    print(json.dumps(row), flush=True)
+        if comm.rank == 0:
+            # Save incrementally so completed sizes survive an interrupted run.
+            report = dict(
+                hostname=platform.node(),
+                platform=platform.platform(),
+                h5py=h5py.__version__,
+                hdf5=h5py.version.hdf5_version,
+                dolfinx=dolfinx.__version__,
+                mpi=MPI.Get_library_version().strip(),
+                repeats=repeats,
+                warmups=warmups,
+                fsync=sync,
+                modes=MODES,
+                results=results,
+            )
+            (output / "results.json").write_text(json.dumps(report, indent=2) + "\n")
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--output", type=Path, required=True)
+    parser.add_argument("--sizes", type=int, nargs="+", default=[16, 256, 768])
+    parser.add_argument("--repeats", type=int, default=5)
+    parser.add_argument("--warmups", type=int, default=1)
+    parser.add_argument("--sync", action="store_true", help="Include fsync after all writes")
+    args = parser.parse_args()
+    if any(n <= 0 for n in args.sizes) or args.repeats <= 0 or args.warmups < 0:
+        parser.error("sizes and repeats must be positive; warmups must be nonnegative")
+    benchmark(args.output, args.sizes, args.repeats, args.warmups, args.sync)

docs/backends/vtkhdf.rst

@@ -6,6 +6,93 @@ VTKHDF-backend
 Relies on MPI compatible :code:`h5py`, similar to the  :ref:`H5PY backend<h5pybackend>`.
 See that backend for detailed installation instructions.
 
+Lossless compression
+~~~~~~~~~~~~~~~~~~~~
+
+Compression is opt-in through ``backend_args``. The public writer signatures are
+unchanged. Reuse the same dictionary on each call to compress both the mesh and
+its point or cell data:
+
+.. code-block:: python
+
+    from mpi4py import MPI
+    from dolfinx.fem import Function, functionspace
+    from dolfinx.mesh import create_unit_square
+    import io4dolfinx
+
+    mesh = create_unit_square(MPI.COMM_WORLD, 16, 16)
+    u = Function(functionspace(mesh, ("Lagrange", 1)), name="u")
+    u.interpolate(lambda x: x[0] + x[1])
+    q = Function(functionspace(mesh, ("DG", 0)), name="q")
+    q.interpolate(lambda x: x[0] - x[1])
+
+    storage_args = {
+        "compression": "gzip",
+        "compression_opts": 4,
+        "shuffle": True,
+    }
+    io4dolfinx.write_mesh(
+        "results.vtkhdf", mesh, backend="vtkhdf", backend_args=storage_args
+    )
+    io4dolfinx.write_point_data(
+        "results.vtkhdf", u, time=0.0, mode=io4dolfinx.FileMode.append,
+        backend="vtkhdf", backend_args=storage_args,
+    )
+    io4dolfinx.write_cell_data(
+        "results.vtkhdf", q, time=0.0, mode=io4dolfinx.FileMode.append,
+        backend="vtkhdf", backend_args=storage_args,
+    )
+
+Supported options are ``compression`` (``"gzip"`` or ``None``),
+``compression_opts`` (an integer from 0 to 9, requiring ``compression="gzip"``;
+the gzip default is 4), and ``shuffle`` (a boolean, default ``False``).
+Shuffle is lossless and can improve gzip compression. Other compression filters,
+lossy scale-offset, checksums, and explicit chunk shapes are not supported by this
+interface. Omitting these options creates uncompressed datasets.
+
+Coordinates, connectivity, topology offsets, cell types, and point/cell values
+are compressed independently. Small bookkeeping arrays, including timestep
+metadata, remain uncompressed. Each dataset uses h5py automatic chunking;
+``chunks=True`` is accepted, but a fixed chunk shape is not. A chunk contains
+elements of one dataset, not a combination of mesh and field arrays. See the
+`h5py chunking guidance <https://docs.h5py.org/en/stable/high/dataset.html#chunked-storage>`_.
+
+Settings apply when a dataset is first created. On append, omitted options retain
+that dataset's existing settings; explicitly conflicting settings raise
+``ValueError`` before datasets are resized. Newly created fields use the options
+of their own write call and do not inherit settings from ``write_mesh``.
+Changing an existing dataset's filters requires rewriting it to a new file.
+
+MPI writes to filtered datasets are collective automatically. Every rank in the
+communicator must call the writer, including ranks with no local data, and must
+pass identical filter options. Parallel compression requires MPI-enabled h5py
+and HDF5 with parallel filtered-write support (introduced in HDF5 1.10.2).
+Invalid options and unavailable filters raise errors rather than silently
+disabling compression. Readers do not need compression arguments.
+
+Compression reduces storage for compressible data, but CPU work, chunk sharing,
+and MPI coordination can increase write time, particularly for small datasets.
+There is no automatic size threshold or speed guarantee. Benchmark representative
+data and rank counts; the repository's ``benchmarks/vtkhdf_compression.py`` reports
+file sizes and mesh/field write and read times for uncompressed output, gzip,
+and gzip with shuffle. For example, from the repository root:
+
+.. code-block:: bash
+
+    OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 mpiexec -n 4 \
+        python benchmarks/vtkhdf_compression.py \
+        --output output/compression/ranks-4 --sizes 16 256 768 --repeats 5 --sync
+
+Each mode is warmed up before measurement and the measured order rotates between
+repetitions. ``--sync`` includes a final filesystem ``fsync`` in write time;
+reads use warm caches. The workload contains smooth scalar/vector fields, so
+compression ratios do not predict those of noisy data. To combine runs in
+``ranks-*`` directories into CSV and charts, run
+``python benchmarks/report_vtkhdf_compression.py output/compression`` in an
+environment with matplotlib. Compression enables collective writes, so timing
+differences include the change of I/O mode as well as the filter itself. See the
+`HDF5 parallel compression guidance <https://portal.hdfgroup.org/documentation/hdf5/latest/_par_compr.html>`_.
+
 .. automodule:: io4dolfinx.backends.vtkhdf.backend
     :members:
     :exclude-members: read_attributes, read_timestamps, write_attributes, write_mesh, write_meshtags, read_meshtags_data, read_dofmap, read_dofs, read_cell_perms, write_function, read_legacy_mesh, snapshot_checkpoint, read_hdf5_array