From c19d5dd4b959de9ab2affc08f9ac1e2a4c7bd191 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sat, 12 Sep 2026 19:33:30 +1000 Subject: [PATCH 1/9] Fix DG1 XDMF export with element-local simplex topology Write independent physical vertices and cell-local linear field traces for DG1 triangles and tetrahedra. Use native cell maps and owned-cell PETSc output to preserve discontinuities without coordinate matching or inter-element averaging. Expose continuous and DG1 fields as separate grids in one XDMF file, support scalar/vector/tensor layouts, and retain native checkpoint data unchanged. Reject unsupported DG visualization explicitly while preserving native-only output. Add serial and MPI regression coverage for affine fields, jumps, orientation, ownership, tensor packing, and checkpoint reload. Document supported layouts and ParaView usage. Validated 24 focused serial tests, four DG tests on eight ranks, and actual ParaView 2D/3D reads. --- .../subsystems/checkpointing-system.md | 20 ++- .../discretisation/discretisation_mesh.py | 100 +++++++++++- src/underworld3/function/field_projection.py | 40 +++++ tests/test_0005_xdmf_dg1.py | 142 ++++++++++++++++++ 4 files changed, 298 insertions(+), 4 deletions(-) create mode 100644 tests/test_0005_xdmf_dg1.py diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index c0502b894..c1445ba41 100644 --- a/docs/developer/subsystems/checkpointing-system.md +++ b/docs/developer/subsystems/checkpointing-system.md @@ -30,13 +30,31 @@ Optional payloads are controlled by explicit flags: ### Visualisation and Coordinate Remap +Continuous fields use the standard mesh vertices, and DG0 fields use cell data. +DG1 fields on full-dimensional triangular and tetrahedral meshes use a second +grid named `DG1` in the same XDMF file. Each simplex has its own three or four +physical vertices: the saved linear polynomial is evaluated within that cell, +without averaging traces across shared edges or faces. Interior DG interpolation +nodes are not mistaken for the physical mesh vertices. + +The DG1 visualization arrays (`vertices`, `cells`, `values`) live under `/dg1` +in each variable HDF5 file. Tensor visualization uses a nine-component 3-by-3 +layout (zero-padded in 2D); `/fields` and PETSc reload data retain their native +layout and precision. Open the one `.xdmf` file in ParaView and select the +`domain` or `DG1` block for the corresponding fields. Do not merge coincident +points or apply point-averaging filters if discontinuities must be preserved. + +Higher-degree discontinuous fields, tensor-product DG cells, embedded manifolds, +and integration-point fields are not supported by this DG1 exporter. They raise +an explicit error when visualization is requested; `create_xdmf=False` still +allows native checkpoint output. Parallel export uses owned cells only. + ```python mesh.write_timestep( "output", index=100, outputPath="output", meshVars=[velocity, pressure, temperature], - time=100.0, create_xdmf=True, ) ``` diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 7f045997e..679f025b4 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -4877,11 +4877,31 @@ def write_timestep( exists. If ``True``, write an indexed mesh file for this timestep. create_xdmf Write ParaView/XDMF-compatible datasets and companion XDMF file. + DG1 on full-dimensional triangles/tetrahedra uses a separate grid + with independent vertices per cell, preserving jumps without + smoothing. Visualization-only arrays live under ``/dg1`` in the + variable files; native checkpoint and reload data are unchanged. + Higher-order discontinuous and non-simplex DG visualization are + not supported (use ``create_xdmf=False`` for native-only output). petsc_reload Write PETSc DMPlex section/vector metadata for reload with ``MeshVariable.read_checkpoint()``. """ + if create_xdmf: + for var in meshVars or []: + integration_point = getattr(var, "is_integration_point", False) + if integration_point or (not var.continuous and var.degree > 0): + if ( + var.degree != 1 or not self.isSimplex + or self.dim not in (2, 3) or self.cdim != self.dim + or integration_point + ): + raise NotImplementedError( + "DG XDMF supports degree-one fields on full-dimensional " + "triangle/tetrahedron meshes only; use create_xdmf=False " + "for native-only checkpoints." + ) options = PETSc.Options() options.setValue("viewer_hdf5_sp_output", True) options.setValue("viewer_hdf5_collective", False) @@ -9632,6 +9652,8 @@ def _write_compat_groups(mesh, var, var_h5_path): Uses ``uw.function.write_vertices_to_viewer`` (PETSc interpolation + ViewerHDF5) for continuous variables, and ``uw.function.write_cell_field_to_viewer`` for cell/DG-0 variables. + DG1 uses ``/dg1`` with disconnected simplex vertices and nodal traces, + never the one-value-per-cell compatibility path. PETSc handles all parallel I/O natively. Vertex coordinates are also written to ``/vertex_fields/coordinates`` @@ -9649,8 +9671,9 @@ def _write_compat_groups(mesh, var, var_h5_path): """ import underworld3 as uw - is_cell = (not var.continuous) or (var.degree == 0) - group = "cell_fields" if is_cell else "vertex_fields" + is_dg1 = not var.continuous and var.degree == 1 + is_cell = var.degree == 0 + group = "dg1" if is_dg1 else ("cell_fields" if is_cell else "vertex_fields") # Some PETSc versions (3.21+) write /vertex_fields/ or /cell_fields/ # automatically during var.write(). Remove any pre-existing group so @@ -9668,7 +9691,10 @@ def _write_compat_groups(mesh, var, var_h5_path): var_h5_path, "a", comm=PETSc.COMM_WORLD, ) - if is_cell: + if is_dg1: + from underworld3.function.field_projection import _write_dg1_to_viewer + _write_dg1_to_viewer(var, viewer) + elif is_cell: uw.function.write_cell_field_to_viewer(var, viewer) else: uw.function.write_vertices_to_viewer(var, viewer) @@ -9676,6 +9702,18 @@ def _write_compat_groups(mesh, var, var_h5_path): viewer.destroy() + if is_dg1: + # Only topology is generated on rank zero; field values and vertices + # were written collectively by PETSc in the same owned-cell order. + if uw.mpi.rank == 0: + with h5py.File(var_h5_path, "a") as f: + nvertices = f["dg1/vertices"].shape[0] + ncorners = mesh.dim + 1 + f["dg1"].create_dataset( + "cells", data=numpy.arange(nvertices, dtype=numpy.int64).reshape(-1, ncorners) + ) + uw.mpi.barrier() + def checkpoint_xdmf( filename: str, @@ -9802,6 +9840,13 @@ def checkpoint_xdmf( header += """ ]>""" + dg_vars = [var for var in meshVars if not var.continuous and var.degree == 1] + collection_start = ( + '' + f'" if dg_vars else "" xdmf_start = f""" @@ -9818,6 +9863,7 @@ def checkpoint_xdmf( &MeshData;:/{geomPath}/vertices + {collection_start} + + + &{first.clean_name}_Data;:/dg1/cells + + + + + &{first.clean_name}_Data;:/dg1/vertices + + +""" + for var in dg_vars: + var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" + with h5py.File(var_filename, "r") as f: + shape = f["dg1/values"].shape + if shape[0] != dg_points[0]: + raise ValueError(f"DG1 visualization size mismatch for {var.clean_name}") + components = shape[1] if len(shape) == 2 else 1 + if var.vtype in (uw.VarType.TENSOR, uw.VarType.SYM_TENSOR): + kind = "Tensor" + elif var.vtype == uw.VarType.MATRIX: + kind = "Matrix" + else: + kind = "Scalar" if components == 1 else "Vector" + dimensions = " ".join(str(value) for value in shape) + dg_grid += f""" + + + &{var.clean_name}_Data;:/dg1/values + + +""" + dg_grid += " " xdmf_end = f""" + {dg_grid} + {collection_end} """ diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py index 28fc4809d..91a01b2cc 100644 --- a/src/underworld3/function/field_projection.py +++ b/src/underworld3/function/field_projection.py @@ -350,3 +350,43 @@ def write_cell_field_to_viewer( mesh_var._sync_lvec_to_gvec() data = mesh_var._gvec.array.reshape(-1, nc).copy() _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) + + +def _write_dg1_to_viewer(mesh_var, viewer): + """Write owned simplex cells with independent physical vertices and DG1 traces. + + Coordinate-section cell maps preserve element ownership and node ordering; + no point location, coordinate matching, or inter-element averaging is used. + Interior DG interpolation nodes define an affine polynomial, evaluated at + that same cell's physical vertices. Native checkpoint vectors are untouched. + """ + mesh = mesh_var.mesh + if ( + mesh_var.continuous or mesh_var.degree != 1 or not mesh.isSimplex + or mesh.dim not in (2, 3) or mesh.cdim != mesh.dim + ): + raise NotImplementedError("DG1 XDMF requires a full-dimensional triangle/tetrahedron mesh") + cstart, cend = mesh.dm.getHeightStratum(0) + owned = np.ones(cend - cstart, dtype=bool) + # Serial DMPlex may have an unset point SF; no cells are ghosts there. + if mesh.dm.comm.getSize() > 1: + _, leaves, remote = mesh.dm.getPointSF().getGraph() + if leaves is None: + leaves = np.arange(len(remote)) + leaves = np.asarray(leaves) + owned[leaves[(leaves >= cstart) & (leaves < cend)] - cstart] = False + rows = mesh._cell_node_indices(1, False).reshape(-1, mesh.dim + 1)[owned] + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, mesh.dim + 1)[owned] + corners = mesh._get_coords_for_basis(1, True)[vertex_rows] + # Closure order is arbitrary; give VTK positively oriented simplices. + negative = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) < 0 + corners[negative] = corners[negative][:, [0, 2, 1] if mesh.dim == 2 else [0, 2, 1, 3]] + nodes = mesh_var.coords[rows] + coefficients = mesh_var._lvec.array.reshape(-1, mesh_var.num_components)[rows] + matrix = (nodes[:, 1:] - nodes[:, :1]).transpose(0, 2, 1) + local = np.linalg.solve(matrix, (corners - nodes[:, :1]).transpose(0, 2, 1)) + weights = np.concatenate((1 - local.sum(axis=1, keepdims=True), local), axis=1) + values = np.einsum("cij,cik->cjk", weights, coefficients).reshape(-1, mesh_var.num_components) + values = _repack_tensor_to_paraview(values, mesh_var.vtype, mesh.dim) + _write_vec_to_group(viewer, corners.reshape(-1, mesh.cdim), "vertices", "/dg1", PETSc.COMM_WORLD) + _write_vec_to_group(viewer, values, "values", "/dg1", PETSc.COMM_WORLD) diff --git a/tests/test_0005_xdmf_dg1.py b/tests/test_0005_xdmf_dg1.py new file mode 100644 index 000000000..163aea3ca --- /dev/null +++ b/tests/test_0005_xdmf_dg1.py @@ -0,0 +1,142 @@ +"""DG1 visualization preserves affine fields and jumps, including MPI ownership.""" + +from pathlib import Path +import xml.etree.ElementTree as ET + +import h5py +import numpy as np +import pytest +import underworld3 as uw + + +@pytest.mark.level_1 +@pytest.mark.tier_b +@pytest.mark.parametrize("dim", [2, 3]) +def test_dg1_simplex_output(tmp_path, dim): + directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0,) * dim, + maxCoords=(1.0,) * dim, + cellSize=0.5, + regular=True, + qdegree=3, + ) + scalar = uw.discretisation.MeshVariable("dg_scalar", mesh, 1, degree=1, continuous=False) + tensor = uw.discretisation.MeshVariable( + "dg_tensor", + mesh, + degree=1, + continuous=False, + vtype=uw.VarType.TENSOR, + ) + pressure = uw.discretisation.MeshVariable("pressure", mesh, 1, degree=1) + vector = uw.discretisation.MeshVariable("dg_vector", mesh, dim, degree=1, continuous=False) + symmetric = uw.discretisation.MeshVariable( + "dg_symmetric", mesh, degree=1, continuous=False, vtype=uw.VarType.SYM_TENSOR + ) + rows = mesh._cell_node_indices(1, False).reshape(-1, dim + 1) + coords = scalar.coords + offset = np.floor(coords[rows].mean(axis=1)[:, 0] * 7 + 1e-8) + scalar.array[rows, 0, 0] = 1 + coords[rows, 0] + 2 * coords[rows, 1] + offset[:, None] + tensor.array[:] = 0 + tensor.array[:, 0, 0] = scalar.array[:, 0, 0] + tensor.array[:, 0, 1] = 3 + coords[:, 0] + tensor.array[:, 1, 0] = -2 + coords[:, 1] + tensor.array[:, 1, 1] = 5 + pressure.array[:, 0, 0] = pressure.coords[:, 0] + vector.array[:, 0, :] = coords + symmetric.array[:] = 0 + symmetric.array[:, 0, 0] = 2 + symmetric.array[:, 1, 1] = 3 + symmetric.array[:, 0, 1] = coords[:, 0] + original = np.array(scalar.array) + mesh.write_timestep( + "fields", + index=0, + outputPath=str(directory), + meshVars=[pressure, scalar, tensor, vector, symmetric], + petsc_reload=True, + ) + restored = uw.discretisation.MeshVariable("restored", mesh, 1, degree=1, continuous=False) + restored.read_checkpoint( + str(directory / "fields.mesh.dg_scalar.00000.h5"), data_name="dg_scalar", same_layout=True + ) + np.testing.assert_allclose(restored.array, original, rtol=1e-12, atol=1e-12) + if uw.mpi.rank != 0: + return + with h5py.File(directory / "fields.mesh.dg_scalar.00000.h5", "r") as handle: + points = handle["dg1/vertices"][:] + cells = handle["dg1/cells"][:] + values = handle["dg1/values"][:].reshape(-1) + native = handle["fields/dg_scalar"][:] + assert len(points) == len(cells) * (dim + 1) + assert len(np.unique(cells)) == len(points) + assert native.size == len(points) + with h5py.File(directory / "fields.mesh.00000.h5", "r") as handle: + assert len(cells) == len(handle["viz/topology/cells"]) + corners = points[cells] + assert np.all(np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) > 0) + centers = corners.mean(axis=1) + assert len(np.unique(np.round(centers, 10), axis=0)) == len(cells) + offset = np.repeat(np.floor(centers[:, 0] * 7 + 1e-8), dim + 1) + expected = 1 + points[:, 0] + 2 * points[:, 1] + offset + np.testing.assert_allclose(values, expected, rtol=1e-12, atol=1e-12) + # Repeated positions can have different values: these traces must not merge. + _, inverse = np.unique(np.round(points, 10), axis=0, return_inverse=True) + low = np.full(inverse.max() + 1, np.inf) + high = np.full(inverse.max() + 1, -np.inf) + np.minimum.at(low, inverse, values) + np.maximum.at(high, inverse, values) + assert np.max(high - low) >= 1 + with h5py.File(directory / "fields.mesh.dg_tensor.00000.h5", "r") as handle: + tensor_values = handle["dg1/values"][:] + np.testing.assert_allclose(handle["dg1/vertices"][:], points) + assert tensor_values.shape == (len(points), 9) + np.testing.assert_allclose(tensor_values[:, 0], expected) + np.testing.assert_allclose(tensor_values[:, 1], 3 + points[:, 0]) + np.testing.assert_allclose(tensor_values[:, 3], -2 + points[:, 1], atol=1e-12) + np.testing.assert_allclose(tensor_values[:, 4], 5) + tree = ET.parse(directory / "fields.mesh.00000.xdmf") + with h5py.File(directory / "fields.mesh.dg_vector.00000.h5", "r") as handle: + np.testing.assert_allclose(handle["dg1/values"][:], points, atol=1e-12) + with h5py.File(directory / "fields.mesh.dg_symmetric.00000.h5", "r") as handle: + sym_values = handle["dg1/values"][:] + assert sym_values.shape == (len(points), 9) + np.testing.assert_allclose(sym_values[:, 0], 2) + np.testing.assert_allclose(sym_values[:, 4], 3) + np.testing.assert_allclose(sym_values[:, 1], points[:, 0], atol=1e-12) + np.testing.assert_allclose(sym_values[:, 3], points[:, 0], atol=1e-12) + grids = tree.findall(".//Grid[@GridType='Uniform']") + assert len(grids) == 2 + dg = next(grid for grid in grids if grid.get("Name") == "DG1") + assert {a.get("Name") for a in dg.findall("Attribute")} == { + "dg_scalar", + "dg_tensor", + "dg_vector", + "dg_symmetric", + } + assert all(a.get("Center") == "Node" for a in dg.findall("Attribute")) + for item in tree.findall(".//DataItem[@Format='HDF']"): + filename, dataset = item.text.strip().split(":", 1) + with h5py.File(directory / filename, "r") as handle: + assert tuple(map(int, item.get("Dimensions").split())) == handle[dataset].shape + + +@pytest.mark.level_1 +@pytest.mark.tier_b +@pytest.mark.parametrize("degree", [1, 2]) +def test_unsupported_dg_layout_fails_before_writing(tmp_path, degree): + directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) + mesh = ( + uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + if degree == 1 + else uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) + ) + dg = uw.discretisation.MeshVariable("dg", mesh, 1, degree=degree, continuous=False) + with pytest.raises(NotImplementedError, match="DG.*XDMF"): + mesh.write_timestep("unsupported", index=0, outputPath=str(directory), meshVars=[dg]) + assert not (directory / "unsupported.mesh.00000.h5").exists() + # Unsupported visualization must not prevent native-only checkpoints. + mesh.write_timestep( + "native", index=0, outputPath=str(directory), meshVars=[dg], create_xdmf=False + ) From 8cd6e3827c1a6d0dceb6ccce5fca2df3c7b6bb91 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 18 Sep 2026 15:14:08 +1000 Subject: [PATCH 2/9] Write physical XDMF data while preserving native checkpoints Keep mesh geometry and field checkpoint datasets in the nondimensional solver frame, and create dimensional visualization copies using declared mesh and variable units. Record those units in HDF5 and XDMF metadata. Use nondimensional coordinates for DG1 interpolation before converting its disconnected visualization grid. This removes the need for benchmark-side coordinate wrappers, field-unit mappings, and HDF5 postprocessing. Cover continuous P2, DG0, and DG1 output in serial and MPI tests. Underworld development team with AI support from Claude Code --- .../subsystems/checkpointing-system.md | 9 ++ .../discretisation/discretisation_mesh.py | 134 ++++++++++++++++-- src/underworld3/function/field_projection.py | 53 ++++++- tests/test_0005_xdmf_physical_units.py | 126 ++++++++++++++++ 4 files changed, 308 insertions(+), 14 deletions(-) create mode 100644 tests/test_0005_xdmf_physical_units.py diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index c1445ba41..ba157810b 100644 --- a/docs/developer/subsystems/checkpointing-system.md +++ b/docs/developer/subsystems/checkpointing-system.md @@ -28,6 +28,15 @@ Optional payloads are controlled by explicit flags: | `create_xdmf=True` | XDMF-compatible visualisation datasets and a companion `.xdmf` file | ParaView and other XDMF tools | | `petsc_reload=True` | PETSc DMPlex section/vector metadata | `MeshVariable.read_checkpoint()` | +When nondimensional scaling is active, one output family contains both numeric +frames. Native datasets under `/geometry`, `/fields`, and `/uw_checkpoint` +remain nondimensional for remapping, exact reload, and numerical comparisons. +The XDMF datasets under `/viz/geometry`, `/vertex_fields`, `/cell_fields`, and +`/dg1` are converted during the write to the mesh and variable units declared +in the model. HDF5 dataset attributes and XDMF `Information` elements record +those units. A model script does not need a field-to-unit mapping or an HDF5 +postprocessing pass. + ### Visualisation and Coordinate Remap Continuous fields use the standard mesh vertices, and DG0 fields use cell data. diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 679f025b4..eb43e376c 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -4839,7 +4839,10 @@ def write_timestep( - ``create_xdmf=True`` writes ParaView/XDMF output. Variable files also receive ``/vertex_fields`` or ``/cell_fields`` compatibility groups, - and rank 0 writes the companion ``.xdmf`` file. + and rank 0 writes the companion ``.xdmf`` file. With active + nondimensional scaling, these visualisation datasets use the mesh and + variable units while native ``/geometry`` and ``/fields`` datasets + remain nondimensional for reload. - ``petsc_reload=True`` writes PETSc DMPlex section/local-vector metadata and an in-place global-vector payload into the same per-variable HDF5 files. These files can then be loaded with @@ -4930,6 +4933,9 @@ def write_timestep( mesh_file = output_base_name + f".mesh.{index:05}.h5" self.write(mesh_file) + if create_xdmf: + _write_visualisation_geometry(self, mesh_file) + variables = [] if meshVars is not None: for var in meshVars: @@ -9688,7 +9694,9 @@ def _write_compat_groups(mesh, var, var_h5_path): uw.mpi.barrier() viewer = PETSc.ViewerHDF5().create( - var_h5_path, "a", comm=PETSc.COMM_WORLD, + var_h5_path, + "a", + comm=PETSc.COMM_WORLD, ) if is_dg1: @@ -9714,6 +9722,94 @@ def _write_compat_groups(mesh, var, var_h5_path): ) uw.mpi.barrier() + _write_visualisation_metadata(mesh, var, var_h5_path, group) + + +def _write_visualisation_geometry(mesh, mesh_h5_path): + """Write physical mesh coordinates without changing restart geometry.""" + import h5py + import underworld3 as uw + from underworld3.function.field_projection import ( + _physical_visualisation_enabled, + _physical_visualisation_values, + ) + + if not _physical_visualisation_enabled(mesh.units): + return + + with uw.selective_ranks(0) as should_execute: + if should_execute: + with h5py.File(mesh_h5_path, "a") as handle: + if "viz/geometry" in handle: + del handle["viz/geometry"] + uw.mpi.barrier() + + viewer = PETSc.ViewerHDF5().create(mesh_h5_path, "a", comm=PETSc.COMM_WORLD) + uw.function.write_coordinates_to_viewer( + mesh, + viewer, + group="/viz/geometry", + name="vertices", + ) + viewer.destroy() + + _, unit_label = _physical_visualisation_values(numpy.ones(1), mesh.units) + with uw.selective_ranks(0) as should_execute: + if should_execute: + with h5py.File(mesh_h5_path, "a") as handle: + handle.attrs["checkpoint_units"] = "nondimensional" + handle.attrs["visualisation_units"] = unit_label + handle["geometry/vertices"].attrs["units"] = "nondimensional" + physical = handle["viz/geometry/vertices"] + physical.attrs["units"] = unit_label + uw.mpi.barrier() + + +def _write_visualisation_metadata(mesh, var, var_h5_path, group): + """Describe native and physical datasets after collective output closes.""" + import h5py + import underworld3 as uw + from underworld3.function.field_projection import ( + _physical_visualisation_enabled, + _physical_visualisation_values, + ) + + field_enabled = _physical_visualisation_enabled(var.units) + coordinates_enabled = _physical_visualisation_enabled(mesh.units) + if not (field_enabled or coordinates_enabled): + return + + _, field_units = _physical_visualisation_values(numpy.ones(1), var.units) + _, coordinate_units = _physical_visualisation_values( + numpy.ones(1), mesh.units + ) + if group == "dg1": + field_path = "dg1/values" + coordinate_path = "dg1/vertices" + else: + field_path = f"{group}/{var.clean_name}_{var.clean_name}" + coordinate_path = f"{group}/coordinates" + + with uw.selective_ranks(0) as should_execute: + if should_execute: + with h5py.File(var_h5_path, "a") as handle: + handle.attrs["checkpoint_units"] = "nondimensional" + if field_enabled: + handle.attrs["visualisation_units"] = field_units + native = handle[f"fields/{var.clean_name}"] + native.attrs["units"] = "nondimensional" + native.attrs["physical_units"] = field_units + physical = handle[field_path] + physical.attrs["units"] = field_units + + if coordinates_enabled: + handle["fields/coordinates"].attrs["units"] = "nondimensional" + + if coordinate_path in handle and coordinate_units is not None: + coordinates = handle[coordinate_path] + coordinates.attrs["units"] = coordinate_units + uw.mpi.barrier() + def checkpoint_xdmf( filename: str, @@ -9725,6 +9821,7 @@ def checkpoint_xdmf( import h5py import os import warnings + from xml.sax.saxutils import escape """Create xdmf file for checkpoints""" @@ -9760,6 +9857,7 @@ def checkpoint_xdmf( ) vertices = geom["vertices"] + geometry_units = vertices.attrs.get("units") numVertices = vertices.shape[0] spaceDim = vertices.shape[1] cells = topo["cells"] @@ -9801,6 +9899,17 @@ def checkpoint_xdmf( h5.close() + def units_information(units, indent): + """Return an XDMF Information element for an optional unit label.""" + if units is None: + return "" + if isinstance(units, bytes): + units = units.decode() + value = escape(str(units), {'"': """}) + return f'\n{indent}' + + geometry_information = units_information(geometry_units, " ") + # We only use a subset of the possible cell types if spaceDim == 2: if numCorners == 3: @@ -9875,7 +9984,7 @@ def checkpoint_xdmf( /Xdmf/Domain/DataItem[@Name="vertices"] - + {geometry_information} """ @@ -9883,7 +9992,7 @@ def checkpoint_xdmf( def get_field_info(h5_filename, mesh_var, center): """ - Return (num_items, num_components, dataset_path) for a mesh variable. + Return shape, path, and units for a mesh variable. Prefers vertex/cell compatibility groups, falls back to /fields layout. """ compat_name = f"{mesh_var.clean_name}_{mesh_var.clean_name}" @@ -9898,9 +10007,10 @@ def get_field_info(h5_filename, mesh_var, center): for path in candidates: if path in f: shp = f[path].shape + units = f[path].attrs.get("units") if len(shp) == 1: - return shp[0], 1, path - return shp[0], shp[1], path + return shp[0], 1, path, units + return shp[0], shp[1], path, units raise RuntimeError( f"Could not locate data for variable '{mesh_var.clean_name}' in {h5_filename}" @@ -9917,7 +10027,9 @@ def get_field_info(h5_filename, mesh_var, center): center = "Cell" else: center = "Node" - numItems, numComponents, dataset_path = get_field_info(var_filename, var, center) + numItems, numComponents, dataset_path, field_units = get_field_info( + var_filename, var, center + ) if center == "Node" and numItems != numVertices: warnings.warn( @@ -9955,7 +10067,7 @@ def get_field_info(h5_filename, mesh_var, center): Dimensions="{data_dimensions}" Format="HDF"> &{var.clean_name+"_Data"};:/{dataset_path} - + {units_information(field_units, " ")} """ attributes += var_attribute @@ -10001,6 +10113,7 @@ def get_field_info(h5_filename, mesh_var, center): with h5py.File(first_filename, "r") as f: dg_cells = f["dg1/cells"].shape dg_points = f["dg1/vertices"].shape + dg_geometry_units = f["dg1/vertices"].attrs.get("units") if dg_cells != (numCells, numCorners) or dg_points != (numCells * numCorners, spaceDim): raise ValueError("DG1 visualization topology does not match the checkpoint mesh") dg_grid = f""" @@ -10013,13 +10126,14 @@ def get_field_info(h5_filename, mesh_var, center): &{first.clean_name}_Data;:/dg1/vertices - + {units_information(dg_geometry_units, " ")} """ for var in dg_vars: var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" with h5py.File(var_filename, "r") as f: shape = f["dg1/values"].shape + field_units = f["dg1/values"].attrs.get("units") if shape[0] != dg_points[0]: raise ValueError(f"DG1 visualization size mismatch for {var.clean_name}") components = shape[1] if len(shape) == 2 else 1 @@ -10034,7 +10148,7 @@ def get_field_info(h5_filename, mesh_var, center): &{var.clean_name}_Data;:/dg1/values - + {units_information(field_units, " ")} """ dg_grid += " " diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py index 91a01b2cc..44b6fc57b 100644 --- a/src/underworld3/function/field_projection.py +++ b/src/underworld3/function/field_projection.py @@ -226,6 +226,45 @@ def _write_vec_to_group(viewer, data_array, name, group, comm): vec.destroy() +def _physical_visualisation_enabled(units): + """Return whether declared units require a physical XDMF copy.""" + import underworld3 as uw + + return ( + units is not None + and uw.get_default_model().has_units_active() + and uw.is_nondimensional_scaling_active() + ) + + +def _physical_visualisation_values(data, units): + """Return output values and their declared unit label. + + Solver and checkpoint vectors use model magnitudes. XDMF arrays are a + user-facing boundary, so an active nondimensional model is converted to + the units declared by the mesh or variable before those arrays are written. + """ + import underworld3 as uw + + if units is None: + return data, None + + target_units = uw.units(units).units if isinstance(units, str) else units + unit_label = str(target_units) + if not _physical_visualisation_enabled(units): + return data, unit_label + + dimensionality = dict(target_units.dimensionality) + if not dimensionality: + return data, unit_label + + physical = uw.dimensionalise( + np.asarray(data), + target_dimensionality=dimensionality, + ).to(target_units) + return np.asarray(physical), unit_label + + def write_vertices_to_viewer( mesh_var: "MeshVariable", viewer: "PETSc.ViewerHDF5", @@ -293,6 +332,8 @@ def write_vertices_to_viewer( if is_tensor: data = _repack_tensor_to_paraview(data, mesh_var.vtype, mesh.dim) + data, _ = _physical_visualisation_values(data, mesh_var.units) + _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) @@ -316,7 +357,8 @@ def write_coordinates_to_viewer( Dataset name (default ``coordinates``). """ coord_gvec = mesh.dm.getCoordinates() - coords = coord_gvec.array.reshape(-1, mesh.dim).copy() + coords = coord_gvec.array.reshape(-1, mesh.cdim).copy() + coords, _ = _physical_visualisation_values(coords, mesh.units) _write_vec_to_group(viewer, coords, name, group, PETSc.COMM_WORLD) @@ -349,16 +391,17 @@ def write_cell_field_to_viewer( nc = mesh_var.num_components mesh_var._sync_lvec_to_gvec() data = mesh_var._gvec.array.reshape(-1, nc).copy() + data, _ = _physical_visualisation_values(data, mesh_var.units) _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) def _write_dg1_to_viewer(mesh_var, viewer): - """Write owned simplex cells with independent physical vertices and DG1 traces. + """Write owned simplex cells with independent vertices and DG1 traces. Coordinate-section cell maps preserve element ownership and node ordering; no point location, coordinate matching, or inter-element averaging is used. Interior DG interpolation nodes define an affine polynomial, evaluated at - that same cell's physical vertices. Native checkpoint vectors are untouched. + that same cell's vertices. Native checkpoint vectors are untouched. """ mesh = mesh_var.mesh if ( @@ -381,12 +424,14 @@ def _write_dg1_to_viewer(mesh_var, viewer): # Closure order is arbitrary; give VTK positively oriented simplices. negative = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) < 0 corners[negative] = corners[negative][:, [0, 2, 1] if mesh.dim == 2 else [0, 2, 1, 3]] - nodes = mesh_var.coords[rows] + nodes = mesh_var.coords_nd[rows] coefficients = mesh_var._lvec.array.reshape(-1, mesh_var.num_components)[rows] matrix = (nodes[:, 1:] - nodes[:, :1]).transpose(0, 2, 1) local = np.linalg.solve(matrix, (corners - nodes[:, :1]).transpose(0, 2, 1)) weights = np.concatenate((1 - local.sum(axis=1, keepdims=True), local), axis=1) values = np.einsum("cij,cik->cjk", weights, coefficients).reshape(-1, mesh_var.num_components) values = _repack_tensor_to_paraview(values, mesh_var.vtype, mesh.dim) + corners, _ = _physical_visualisation_values(corners, mesh.units) + values, _ = _physical_visualisation_values(values, mesh_var.units) _write_vec_to_group(viewer, corners.reshape(-1, mesh.cdim), "vertices", "/dg1", PETSc.COMM_WORLD) _write_vec_to_group(viewer, values, "values", "/dg1", PETSc.COMM_WORLD) diff --git a/tests/test_0005_xdmf_physical_units.py b/tests/test_0005_xdmf_physical_units.py new file mode 100644 index 000000000..7d17e8584 --- /dev/null +++ b/tests/test_0005_xdmf_physical_units.py @@ -0,0 +1,126 @@ +"""Physical XDMF output keeps native checkpoint arrays nondimensional.""" + +from pathlib import Path + +import h5py +import numpy as np +import pytest + +import underworld3 as uw + + +def _set_reference_scales(): + """Use exact scales with easy-to-check physical conversions.""" + orchestration_model = uw.get_default_model() + orchestration_model.set_scaling_mode("exact") + orchestration_model.set_reference_quantities( + length=uw.quantity(10, "km"), + velocity=uw.quantity(5, "mm/year"), + pressure=uw.quantity(2, "MPa"), + ) + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_xdmf_uses_declared_physical_units(tmp_path): + """Visualisation copies are physical while restart data stays native.""" + _set_reference_scales() + + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + velocity = uw.discretisation.MeshVariable("velocity", mesh, mesh.dim, degree=2, units="mm/year") + pressure = uw.discretisation.MeshVariable( + "pressure", mesh, 1, degree=0, continuous=False, units="MPa" + ) + velocity.data[:, 0] = 2.0 + velocity.data[:, 1] = 3.0 + pressure.data[:, 0] = 4.0 + + directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) + mesh.write_timestep( + "physical", + index=0, + outputPath=str(directory), + meshVars=[velocity, pressure], + petsc_reload=True, + ) + + with uw.selective_ranks(0) as should_execute: + if not should_execute: + return + + mesh_file = directory / "physical.mesh.00000.h5" + velocity_file = directory / "physical.mesh.velocity.00000.h5" + pressure_file = directory / "physical.mesh.pressure.00000.h5" + + with h5py.File(mesh_file, "r") as handle: + native_coordinates = handle["geometry/vertices"][:] + physical_coordinates = handle["viz/geometry/vertices"][:] + np.testing.assert_allclose(physical_coordinates, native_coordinates * 10.0) + assert handle["geometry/vertices"].attrs["units"] == "nondimensional" + assert handle["viz/geometry/vertices"].attrs["units"] == "kilometer" + + with h5py.File(velocity_file, "r") as handle: + native = handle["fields/velocity"][:].reshape(-1, mesh.dim) + physical = handle["vertex_fields/velocity_velocity"][:].reshape(-1, mesh.dim) + np.testing.assert_allclose(native, [[2.0, 3.0]] * len(native)) + np.testing.assert_allclose(physical, [[10.0, 15.0]] * len(physical)) + np.testing.assert_allclose( + handle["vertex_fields/coordinates"][:].reshape(-1, mesh.dim), + native_coordinates * 10.0, + ) + assert handle["fields/velocity"].attrs["units"] == "nondimensional" + assert handle["vertex_fields/velocity_velocity"].attrs["units"] == "millimeter / year" + assert handle["vertex_fields/coordinates"].attrs["units"] == "kilometer" + np.testing.assert_allclose( + handle["uw_checkpoint/velocity"][:].reshape(-1, mesh.dim), native + ) + + with h5py.File(pressure_file, "r") as handle: + native = handle["fields/pressure"][:].reshape(-1) + physical = handle["cell_fields/pressure_pressure"][:].reshape(-1) + np.testing.assert_allclose(native, 4.0) + np.testing.assert_allclose(physical, 8.0) + assert handle["cell_fields/pressure_pressure"].attrs["units"] == "megapascal" + + xdmf = (directory / "physical.mesh.00000.xdmf").read_text() + assert "&MeshData;:/viz/geometry/vertices" in xdmf + assert 'Name="velocity"' in xdmf + assert 'Information Name="Units" Value="millimeter / year"' in xdmf + assert 'Information Name="Units" Value="kilometer"' in xdmf + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_dg1_xdmf_uses_native_interpolation_and_physical_output(tmp_path): + """DG1 interpolation stays native before its disconnected grid is scaled.""" + _set_reference_scales() + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) + pressure = uw.discretisation.MeshVariable( + "dg_pressure", mesh, 1, degree=1, continuous=False, units="MPa" + ) + pressure.data[:, 0] = 4.0 + + directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) + mesh.write_timestep( + "dg_physical", + index=0, + outputPath=str(directory), + meshVars=[pressure], + petsc_reload=True, + ) + + with uw.selective_ranks(0) as should_execute: + if not should_execute: + return + + field_file = directory / "dg_physical.mesh.dg_pressure.00000.h5" + with h5py.File(field_file, "r") as handle: + np.testing.assert_allclose(handle["fields/dg_pressure"][:], 4.0) + np.testing.assert_allclose(handle["dg1/values"][:], 8.0) + assert np.isclose(handle["dg1/vertices"][:].max(), 10.0) + assert handle["dg1/vertices"].attrs["units"] == "kilometer" + assert handle["dg1/values"].attrs["units"] == "megapascal" + + xdmf = (directory / "dg_physical.mesh.00000.xdmf").read_text() + assert 'Information Name="Units" Value="megapascal"' in xdmf + assert 'Information Name="Units" Value="kilometer"' in xdmf From 22a4a0f30cb1d913e26d2ab2c0996459e9d1172c Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 18 Sep 2026 15:24:30 +1000 Subject: [PATCH 3/9] Propagate mesh units to Surface distance fields Declare signed and unsigned Surface distance variables with the mesh coordinate unit so physical XDMF export scales geometric distance fields without benchmark-side mappings. Extend the physical-output regression to cover Surface.abs_distance while preserving its nondimensional checkpoint values. Underworld development team with AI support from Claude Code --- src/underworld3/meshing/surfaces.py | 2 ++ tests/test_0005_xdmf_physical_units.py | 33 +++++++++++++++++++++++--- 2 files changed, 32 insertions(+), 3 deletions(-) diff --git a/src/underworld3/meshing/surfaces.py b/src/underworld3/meshing/surfaces.py index d62027e28..030e65ead 100644 --- a/src/underworld3/meshing/surfaces.py +++ b/src/underworld3/meshing/surfaces.py @@ -1318,6 +1318,7 @@ def _compute_distance_field(self) -> None: 1, degree=self.mesh.degree, varsymbol=f"d_{{{self._symbol}}}", + units=self.mesh.units, ) # Get mesh coordinates in model (internal) space. @@ -1351,6 +1352,7 @@ def _compute_distance_field(self) -> None: 1, degree=self.mesh.degree, varsymbol=f"|d_{{{self._symbol}}}|", + units=self.mesh.units, ) with uw.synchronised_array_update(): diff --git a/tests/test_0005_xdmf_physical_units.py b/tests/test_0005_xdmf_physical_units.py index 7d17e8584..00c743851 100644 --- a/tests/test_0005_xdmf_physical_units.py +++ b/tests/test_0005_xdmf_physical_units.py @@ -27,10 +27,19 @@ def test_xdmf_uses_declared_physical_units(tmp_path): _set_reference_scales() mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) - velocity = uw.discretisation.MeshVariable("velocity", mesh, mesh.dim, degree=2, units="mm/year") + velocity = uw.discretisation.MeshVariable( + "velocity", mesh, mesh.dim, degree=2, units="mm/year" + ) pressure = uw.discretisation.MeshVariable( "pressure", mesh, 1, degree=0, continuous=False, units="MPa" ) + surface = uw.meshing.Surface( + "unit_test", + mesh, + control_points=uw.quantity([[0.0, 0.0], [10.0, 0.0]], "km"), + ) + surface.discretize() + distance = surface.abs_distance velocity.data[:, 0] = 2.0 velocity.data[:, 1] = 3.0 pressure.data[:, 0] = 4.0 @@ -40,7 +49,7 @@ def test_xdmf_uses_declared_physical_units(tmp_path): "physical", index=0, outputPath=str(directory), - meshVars=[velocity, pressure], + meshVars=[velocity, pressure, distance], petsc_reload=True, ) @@ -51,6 +60,7 @@ def test_xdmf_uses_declared_physical_units(tmp_path): mesh_file = directory / "physical.mesh.00000.h5" velocity_file = directory / "physical.mesh.velocity.00000.h5" pressure_file = directory / "physical.mesh.pressure.00000.h5" + distance_file = directory / "physical.mesh.surf_unit_test_absdistance.00000.h5" with h5py.File(mesh_file, "r") as handle: native_coordinates = handle["geometry/vertices"][:] @@ -69,7 +79,10 @@ def test_xdmf_uses_declared_physical_units(tmp_path): native_coordinates * 10.0, ) assert handle["fields/velocity"].attrs["units"] == "nondimensional" - assert handle["vertex_fields/velocity_velocity"].attrs["units"] == "millimeter / year" + assert ( + handle["vertex_fields/velocity_velocity"].attrs["units"] + == "millimeter / year" + ) assert handle["vertex_fields/coordinates"].attrs["units"] == "kilometer" np.testing.assert_allclose( handle["uw_checkpoint/velocity"][:].reshape(-1, mesh.dim), native @@ -82,6 +95,20 @@ def test_xdmf_uses_declared_physical_units(tmp_path): np.testing.assert_allclose(physical, 8.0) assert handle["cell_fields/pressure_pressure"].attrs["units"] == "megapascal" + with h5py.File(distance_file, "r") as handle: + native = handle["fields/surf_unit_test_absdistance"][:].reshape(-1) + physical = handle[ + "vertex_fields/surf_unit_test_absdistance_surf_unit_test_absdistance" + ][:].reshape(-1) + np.testing.assert_allclose(physical, native * 10.0) + assert distance.units == uw.units("km").units + assert ( + handle[ + "vertex_fields/surf_unit_test_absdistance_surf_unit_test_absdistance" + ].attrs["units"] + == "kilometer" + ) + xdmf = (directory / "physical.mesh.00000.xdmf").read_text() assert "&MeshData;:/viz/geometry/vertices" in xdmf assert 'Name="velocity"' in xdmf From b25d2bad05fa583d433ca51634bf6003c826bd9a Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 18 Sep 2026 17:40:30 +1000 Subject: [PATCH 4/9] Minimize XDMF field duplication Make dimensional /fields datasets authoritative for analysis and direct XDMF rendering. Export continuous P2 triangles with Triangle_6 connectivity and DG1 simplices with exact disconnected-corner basis conversion, while reducing only unsupported continuous P3+ and DG2+ layouts to compact P1 and DG0 visualization arrays.\n\nStore optional native PETSc restart data under /uw_checkpoint, reuse its existing variable vector for same-layout reload, and retain migration-based reload for reconstructed meshes. Teach read_timestep to convert dimensional fields back to the active model frame and reject DG1 corner data that cannot be inverted by nearest-neighbour remapping.\n\nReplace compatibility-group tests with serial and MPI coverage for direct layouts, high-order reductions through P4/DG4, dimensional units, ParaView-facing topology, and checkpoint round trips. Update checkpoint documentation for the new storage contract. --- .../checkpoint-output-and-reload-methods.md | 57 +- .../subsystems/checkpointing-system.md | 64 +- .../discretisation/discretisation_mesh.py | 434 ++++++------ .../discretisation_mesh_variables.py | 101 ++- src/underworld3/function/field_projection.py | 164 ++++- tests/test_0003_save_load.py | 10 +- tests/test_0005_xdmf_compat.py | 637 ++++++------------ tests/test_0005_xdmf_dg1.py | 145 ++-- tests/test_0005_xdmf_physical_units.py | 165 ++--- tests/test_0010_snapshot_disk_format.py | 5 +- 10 files changed, 859 insertions(+), 923 deletions(-) diff --git a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md index 929301078..1acae8383 100644 --- a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md +++ b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md @@ -12,16 +12,34 @@ but new code should use `write_timestep(..., petsc_reload=True)`. ## Standard API `write_timestep()` always writes the mesh file and one HDF5 file per mesh -variable. Mesh-variable files always contain raw coordinate/value datasets under -`/fields`, which are the source data used by `MeshVariable.read_timestep()` for -coordinate/KDTree remapping. +variable. With XDMF enabled, mesh-variable files contain dimensional +coordinates and values under `/fields`. These are the authoritative arrays for +analysis and the source data used by `MeshVariable.read_timestep()` for +coordinate remapping. The two optional payloads are selected with explicit flags: | Flag | Output payload | Reader/use case | | --- | --- | --- | -| `create_xdmf=True` | `/vertex_fields` or `/cell_fields` compatibility datasets plus a companion `.xdmf` file | ParaView/XDMF visualisation | -| `petsc_reload=True` | PETSc DMPlex section/vector metadata under `/topologies/uw_mesh/dms/...` | `MeshVariable.read_checkpoint()` PETSc-native reload | +| `create_xdmf=True` | Dimensional `/fields`, compact high-order reductions when needed, and a companion `.xdmf` file | Analysis and ParaView/XDMF visualisation | +| `petsc_reload=True` | Native PETSc DMPlex section/vector data under `/uw_checkpoint` | `MeshVariable.read_checkpoint()` exact reload | + +The XDMF storage choice follows the finite-element layout: + +| Field layout | XDMF representation | +| --- | --- | +| Continuous P1 | Direct `/fields` node values | +| Continuous P2 triangles | Direct `/fields` values with `Triangle_6` connectivity | +| DG0 | Direct `/fields` cell values | +| DG1 triangles/tetrahedra | Exact element-local values at disconnected corners | +| Continuous P3+ or unsupported P2 | Compact P1 dataset under `/visualization` | +| DG2+ or unsupported DG1 | Compact DG0 dataset under `/visualization` | + +Direct DG1 output is an exact basis conversion for analysis and visualization. +Its disconnected corner coordinates differ from UW3's interior DG1 +interpolation nodes, so use `/uw_checkpoint` with `read_checkpoint()` for exact +DG1 restart. `read_timestep()` performs nearest-neighbour remapping and does not +invert that basis conversion. ### Visualisation And Remap @@ -48,11 +66,12 @@ output.mesh.Pressure.00000.h5 output.mesh.00000.xdmf ``` -The field files contain coordinate/value datasets such as `/fields/` and -`/fields/coordinates`, plus vertex-field datasets for visualisation. Reloading -uses coordinate-based remapping. In practice this means the target variable is -filled by comparing target coordinates to source coordinates, using a KDTree or -similar nearest-neighbour/remap process. +The field files contain `/fields/` and `/fields/coordinates`. P1, P2 +triangles, DG0, and DG1 simplices are visualized directly from those datasets. +Continuous P3+ fields use one compact P1 visualization reduction, and DG2+ +fields use DG0. Reloading with `read_timestep()` compares target coordinates to +the dimensional source coordinates and converts the saved values back to the +active model's nondimensional solver frame. ### Unified Visualisation And PETSc Reload @@ -75,9 +94,9 @@ velocity.read_checkpoint( ) ``` -With both `create_xdmf=True` and `petsc_reload=True`, the same variable file can -be used by `read_timestep()` for coordinate/KDTree remapping and by -`read_checkpoint()` for exact PETSc-native reload. +With both flags enabled, the same variable file contains dimensional `/fields` +for analysis, XDMF, and coordinate remapping plus native `/uw_checkpoint` data +for exact PETSc reload. ### PETSc Reload Without XDMF @@ -94,8 +113,8 @@ mesh.write_timestep( ) ``` -This still writes raw `/fields` datasets, but it does not write -`/vertex_fields`, `/cell_fields`, or a companion `.xdmf` file. +This uses the established native field writer and does not create a companion +`.xdmf` file. The optional checkpoint data are added under `/uw_checkpoint`. Typical PETSc-reload-only files still use the timestep naming convention: @@ -105,12 +124,14 @@ restart.mesh.Velocity.00000.h5 restart.mesh.Pressure.00000.h5 ``` -The variable files contain raw `/fields` datasets and PETSc reload metadata -under `/topologies/uw_mesh/dms//`. +The variable files contain native `/fields` datasets and PETSc reload metadata +under `/uw_checkpoint/topologies/uw_mesh/dms//`. ### Advantages -- Produces XDMF/HDF5 files suitable for visualisation workflows. +- Produces dimensional HDF5 fields suitable for analysis and visualisation. +- Avoids duplicate P1 field arrays when the finite-element layout can be + represented directly. - Can remap data onto a different mesh or a different node layout. - Useful for postprocessing where exact finite-element section identity is not required. diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index ba157810b..3f295cacc 100644 --- a/docs/developer/subsystems/checkpointing-system.md +++ b/docs/developer/subsystems/checkpointing-system.md @@ -18,7 +18,9 @@ including registered meshes, variables, swarms, and Python-side state bearers. `Mesh.write_timestep()` is the standard mesh and mesh-variable output method. It writes one mesh HDF5 file and one HDF5 file per requested mesh variable. -Each variable file always contains `/fields` coordinate/value datasets used by +With `create_xdmf=True`, each variable file contains dimensional +`/fields/coordinates` and `/fields/` datasets. These arrays are the +authoritative analysis output and are also used by `MeshVariable.read_timestep()`. Optional payloads are controlled by explicit flags: @@ -28,35 +30,35 @@ Optional payloads are controlled by explicit flags: | `create_xdmf=True` | XDMF-compatible visualisation datasets and a companion `.xdmf` file | ParaView and other XDMF tools | | `petsc_reload=True` | PETSc DMPlex section/vector metadata | `MeshVariable.read_checkpoint()` | -When nondimensional scaling is active, one output family contains both numeric -frames. Native datasets under `/geometry`, `/fields`, and `/uw_checkpoint` -remain nondimensional for remapping, exact reload, and numerical comparisons. -The XDMF datasets under `/viz/geometry`, `/vertex_fields`, `/cell_fields`, and -`/dg1` are converted during the write to the mesh and variable units declared -in the model. HDF5 dataset attributes and XDMF `Information` elements record -those units. A model script does not need a field-to-unit mapping or an HDF5 -postprocessing pass. +When nondimensional scaling is active, `/fields` is converted during the write +to the mesh and variable units declared in the model. HDF5 attributes and XDMF +`Information` elements record those units. Analysis scripts can therefore read +physical values directly, without maintaining their own conversion table. + +Set `petsc_reload=True` only when an exact restart is needed. It adds the native +nondimensional PETSc payload under `/uw_checkpoint`; the visualization and +analysis arrays remain dimensional. ### Visualisation and Coordinate Remap -Continuous fields use the standard mesh vertices, and DG0 fields use cell data. -DG1 fields on full-dimensional triangular and tetrahedral meshes use a second -grid named `DG1` in the same XDMF file. Each simplex has its own three or four -physical vertices: the saved linear polynomial is evaluated within that cell, -without averaging traces across shared edges or faces. Interior DG interpolation -nodes are not mistaken for the physical mesh vertices. - -The DG1 visualization arrays (`vertices`, `cells`, `values`) live under `/dg1` -in each variable HDF5 file. Tensor visualization uses a nine-component 3-by-3 -layout (zero-padded in 2D); `/fields` and PETSc reload data retain their native -layout and precision. Open the one `.xdmf` file in ParaView and select the -`domain` or `DG1` block for the corresponding fields. Do not merge coincident -points or apply point-averaging filters if discontinuities must be preserved. - -Higher-degree discontinuous fields, tensor-product DG cells, embedded manifolds, -and integration-point fields are not supported by this DG1 exporter. They raise -an explicit error when visualization is requested; `create_xdmf=False` still -allows native checkpoint output. Parallel export uses owned cells only. +XDMF reads P1 and DG0 values directly from `/fields`. Continuous P2 fields on +triangles use XDMF `Triangle_6` connectivity, including the three edge nodes, +so no P1 projection is stored. DG1 fields on full-dimensional triangular and +tetrahedral meshes use disconnected element corners. The element polynomial is +evaluated at each cell's corners, preserving jumps without averaging traces +across shared edges or faces. + +XDMF cannot represent every UW3 finite-element layout directly. Continuous P3+ +fields and unsupported P2 layouts receive one compact P1 dataset under +`/visualization`. DG2+ fields and unsupported DG1 layouts receive one compact +DG0 dataset. Their exact dimensional values remain under `/fields`. Integration +point fields are not supported by this writer. + +Direct DG1 output stores an exact basis conversion at disconnected element +corners. This preserves the field for analysis and visualization, but those +corner coordinates differ from UW3's interior DG1 interpolation nodes. Use the +optional `/uw_checkpoint` payload and `read_checkpoint()` for exact DG1 solver +restart; nearest-neighbour `read_timestep()` is not an inverse basis conversion. ```python mesh.write_timestep( @@ -123,10 +125,10 @@ output/restart.mesh.velocity.00100.h5 output/restart.mesh.pressure.00100.h5 ``` -The variable files contain `/fields` datasets and PETSc reload metadata under -`/topologies/uw_mesh/dms//`. `read_checkpoint()` uses PETSc DMPlex -topology, section, vector, and `PetscSF` metadata. It does not use KDTree -remapping. +The variable files contain PETSc reload metadata and native values under +`/uw_checkpoint/topologies/uw_mesh/dms//`. `read_checkpoint()` uses +PETSc DMPlex topology, section, vector, and `PetscSF` metadata. It does not use +the dimensional `/fields` values or KDTree remapping. ### Unified Visualisation and PETSc Reload diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index eb43e376c..96f6b37b1 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -4832,21 +4832,20 @@ def write_timestep( - one mesh HDF5 file, shared across timesteps unless ``meshUpdates=True`` - one HDF5 file per mesh variable - - raw coordinate/value datasets under ``/fields`` for coordinate-based - reload with ``MeshVariable.read_timestep()`` + - dimensional coordinate/value datasets under ``/fields`` for + analysis and coordinate-based reload with + ``MeshVariable.read_timestep()`` The optional payloads are controlled explicitly: - - ``create_xdmf=True`` writes ParaView/XDMF output. Variable files also - receive ``/vertex_fields`` or ``/cell_fields`` compatibility groups, - and rank 0 writes the companion ``.xdmf`` file. With active - nondimensional scaling, these visualisation datasets use the mesh and - variable units while native ``/geometry`` and ``/fields`` datasets - remain nondimensional for reload. - - ``petsc_reload=True`` writes PETSc DMPlex section/local-vector - metadata and an in-place global-vector payload into the same - per-variable HDF5 files. These files can then be loaded with - ``MeshVariable.read_checkpoint()`` for exact restart. + - ``create_xdmf=True`` writes a companion XDMF file that reads the + dimensional ``/fields`` datasets directly for P1, P2 triangles, + DG0, and DG1 simplices. Continuous P3+ fields receive one compact P1 + visualization dataset; discontinuous DG2+ fields receive DG0. + - ``petsc_reload=True`` additionally writes native nondimensional + PETSc DMPlex section/local-vector data under ``/uw_checkpoint``. + Load that optional payload with ``MeshVariable.read_checkpoint()`` + for an exact solver restart. Common choices are: @@ -4857,9 +4856,9 @@ def write_timestep( - unified visualisation/remap and PETSc reload: ``create_xdmf=True, petsc_reload=True`` - With both flags enabled, the same variable HDF5 file can be used by - ``MeshVariable.read_timestep()`` for coordinate/KDTree remapping and by - ``MeshVariable.read_checkpoint()`` for exact PETSc-native reload. + With both flags enabled, the same variable HDF5 file supports + dimensional analysis, XDMF visualization, coordinate remapping, and + exact PETSc-native restart without duplicate P1 field copies. Parameters ---------- @@ -4879,13 +4878,12 @@ def write_timestep( If ``False``, reuse ``.mesh.00000.h5`` when it already exists. If ``True``, write an indexed mesh file for this timestep. create_xdmf - Write ParaView/XDMF-compatible datasets and companion XDMF file. - DG1 on full-dimensional triangles/tetrahedra uses a separate grid - with independent vertices per cell, preserving jumps without - smoothing. Visualization-only arrays live under ``/dg1`` in the - variable files; native checkpoint and reload data are unchanged. - Higher-order discontinuous and non-simplex DG visualization are - not supported (use ``create_xdmf=False`` for native-only output). + Write ParaView/XDMF-compatible dimensional datasets and the + companion XDMF file. DG1 on full-dimensional triangles/tetrahedra + uses independent vertices per cell, preserving jumps without + smoothing. Higher-order fields use compact P1 or DG0 visualization + reductions while their exact dimensional values remain in + ``/fields``. petsc_reload Write PETSc DMPlex section/vector metadata for reload with ``MeshVariable.read_checkpoint()``. @@ -4893,18 +4891,11 @@ def write_timestep( """ if create_xdmf: for var in meshVars or []: - integration_point = getattr(var, "is_integration_point", False) - if integration_point or (not var.continuous and var.degree > 0): - if ( - var.degree != 1 or not self.isSimplex - or self.dim not in (2, 3) or self.cdim != self.dim - or integration_point - ): - raise NotImplementedError( - "DG XDMF supports degree-one fields on full-dimensional " - "triangle/tetrahedron meshes only; use create_xdmf=False " - "for native-only checkpoints." - ) + if getattr(var, "is_integration_point", False): + raise NotImplementedError( + "Integration-point XDMF output is not supported; use " + "create_xdmf=False for native-only output." + ) options = PETSc.Options() options.setValue("viewer_hdf5_sp_output", True) options.setValue("viewer_hdf5_collective", False) @@ -4936,24 +4927,25 @@ def write_timestep( if create_xdmf: _write_visualisation_geometry(self, mesh_file) - variables = [] + checkpoint_variables = [] if meshVars is not None: for var in meshVars: save_location = output_base_name + f".mesh.{var.clean_name}.{index:05}.h5" var.write(save_location) + if petsc_reload: + self._write_petsc_reload_file(save_location, [var], mode="a") if create_xdmf: - _write_compat_groups(self, var, save_location) - variables.append((var, save_location)) + _write_xdmf_field(self, var, save_location) if swarmVars is not None: for svar in swarmVars: save_location = output_base_name + f".proxy.{svar.clean_name}.{index:05}.h5" svar.write_proxy(save_location) if petsc_reload: - variables.append((svar._meshVar, save_location)) + checkpoint_variables.append((svar._meshVar, save_location)) if petsc_reload: - for var, save_location in variables: + for var, save_location in checkpoint_variables: self._write_petsc_reload_file(save_location, [var], mode="a") if create_xdmf and uw.mpi.rank == 0: @@ -4978,9 +4970,8 @@ def petsc_save_checkpoint( This is a convenience wrapper around ``write_timestep()`` that provides the simpler interface used by earlier Underworld3 code. - Output uses the same per-variable file layout and XDMF generation - (including vertex/cell compatibility groups, field projection, and - tensor repacking) as ``write_timestep()``. + Output uses the same compact dimensional ``/fields`` layout and XDMF + generation as ``write_timestep()``. Parameters ---------- @@ -5046,7 +5037,7 @@ def _write_petsc_reload_variable(self, viewer, var): subdm.destroy() def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): - """Write DMPlex reload metadata and in-place vector payloads.""" + """Write compact DMPlex reload metadata and native local vectors.""" old_dm_name = self.dm.getName() self.dm.setName("uw_mesh") @@ -5055,7 +5046,12 @@ def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): checkpoint_file, mode, comm=PETSc.COMM_WORLD ) viewer.pushFormat(PETSc.Viewer.Format.HDF5_PETSC) + viewer.pushGroup("/uw_checkpoint") try: + # PETSc needs the complete source section to construct the + # migration SF when the checkpoint is read with a different MPI + # ownership ordering. This is metadata only; field values are + # still stored once in each variable's local-vector payload. self.dm.sectionView(viewer, self.dm) for var in variables: @@ -5063,28 +5059,12 @@ def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): uw.mpi.barrier() finally: + viewer.popGroup() viewer.popFormat() viewer.destroy() if old_dm_name is not None: self.dm.setName(old_dm_name) - viewer = PETSc.ViewerHDF5().create( - checkpoint_file, "a", comm=PETSc.COMM_WORLD - ) - try: - viewer.pushGroup("/uw_checkpoint") - for var in variables: - var._sync_lvec_to_gvec() - checkpoint_vec = PETSc.Vec().createWithArray( - var._gvec.array_r, comm=PETSc.COMM_WORLD - ) - checkpoint_vec.setName(var.clean_name) - viewer(checkpoint_vec) - checkpoint_vec.destroy() - viewer.popGroup() - finally: - viewer.destroy() - @timing.routine_timer_decorator def write_checkpoint( self, @@ -5130,8 +5110,8 @@ def write_checkpoint( into one file: ``.checkpoint..h5``. create_xdmf If ``True``, route through ``write_timestep()`` and write XDMF, - vertex/cell compatibility groups, coordinate/KDTree remap data, - and PETSc reload metadata. The output uses the timestep filename + dimensional field/remap data, and PETSc reload metadata. The output + uses the timestep filename convention ``.mesh...h5``. This mode does not support ``unique_id=True`` or ``separate_variable_files=False``. """ @@ -9652,18 +9632,13 @@ def remesh(self, metric_field, verbose=False): return -def _write_compat_groups(mesh, var, var_h5_path): - """Write ``/vertex_fields/`` or ``/cell_fields/`` compatibility groups. - - Uses ``uw.function.write_vertices_to_viewer`` (PETSc interpolation + - ViewerHDF5) for continuous variables, and - ``uw.function.write_cell_field_to_viewer`` for cell/DG-0 variables. - DG1 uses ``/dg1`` with disconnected simplex vertices and nodal traces, - never the one-value-per-cell compatibility path. - PETSc handles all parallel I/O natively. +def _write_xdmf_field(mesh, var, var_h5_path): + """Replace native remap arrays with dimensional field output for XDMF. - Vertex coordinates are also written to ``/vertex_fields/coordinates`` - for XDMF compatibility. + P1, P2 triangles, DG0, and DG1 simplices are represented directly under + ``/fields``. Unsupported higher-order layouts retain their exact physical + values under ``/fields`` and receive one compact visualization reduction: + continuous fields use P1 and discontinuous fields use DG0. Parameters ---------- @@ -9673,24 +9648,43 @@ def _write_compat_groups(mesh, var, var_h5_path): The variable whose data has already been written to *var_h5_path* by ``var.write()`` (so ``var._gvec`` is up-to-date). var_h5_path : str - Path to the HDF5 file (already contains ``/fields/``). + Path to the HDF5 file. Native restart data, when requested, has already + been written under ``/uw_checkpoint``. """ + import h5py import underworld3 as uw + from underworld3.function.field_projection import ( + _physical_visualisation_values, + _write_dg1_to_viewer, + write_field_coordinates_to_viewer, + write_field_to_viewer, + write_p2_triangle_topology_to_viewer, + write_projected_field_to_viewer, + ) - is_dg1 = not var.continuous and var.degree == 1 - is_cell = var.degree == 0 - group = "dg1" if is_dg1 else ("cell_fields" if is_cell else "vertex_fields") - - # Some PETSc versions (3.21+) write /vertex_fields/ or /cell_fields/ - # automatically during var.write(). Remove any pre-existing group so - # that our compat writer can create it afresh (otherwise PETSc error 76 - # on duplicate dataset). - import h5py + direct_p2 = ( + var.continuous + and var.degree == 2 + and mesh.isSimplex + and mesh.dim == 2 + and mesh.cdim == 2 + ) + direct_dg1 = ( + not var.continuous + and var.degree == 1 + and mesh.isSimplex + and mesh.dim in (2, 3) + and mesh.cdim == mesh.dim + ) + needs_projection = (var.continuous and var.degree > 2) or ( + var.continuous and var.degree == 2 and not direct_p2 + ) or (not var.continuous and var.degree > 0 and not direct_dg1) if uw.mpi.rank == 0: with h5py.File(var_h5_path, "a") as f: - if group in f: - del f[group] + for group in ("fields", "vertex_fields", "cell_fields", "dg1", "visualization"): + if group in f: + del f[group] uw.mpi.barrier() viewer = PETSc.ViewerHDF5().create( @@ -9699,30 +9693,56 @@ def _write_compat_groups(mesh, var, var_h5_path): comm=PETSc.COMM_WORLD, ) - if is_dg1: - from underworld3.function.field_projection import _write_dg1_to_viewer - _write_dg1_to_viewer(var, viewer) - elif is_cell: - uw.function.write_cell_field_to_viewer(var, viewer) + if direct_dg1: + _write_dg1_to_viewer( + var, + viewer, + group="/fields", + coordinate_name="coordinates", + value_name=var.clean_name, + repack_tensors=False, + ) else: - uw.function.write_vertices_to_viewer(var, viewer) - uw.function.write_coordinates_to_viewer(mesh, viewer) + write_field_to_viewer(var, viewer, "/fields", var.clean_name) + write_field_coordinates_to_viewer(var, viewer, "/fields") + if direct_p2: + write_p2_triangle_topology_to_viewer(var, viewer, group="/fields") + elif needs_projection: + target_degree = 1 if var.continuous else 0 + write_projected_field_to_viewer( + var, + viewer, + target_degree=target_degree, + continuous=var.continuous, + group="/visualization", + name=var.clean_name, + ) viewer.destroy() - - if is_dg1: - # Only topology is generated on rank zero; field values and vertices - # were written collectively by PETSc in the same owned-cell order. - if uw.mpi.rank == 0: - with h5py.File(var_h5_path, "a") as f: - nvertices = f["dg1/vertices"].shape[0] - ncorners = mesh.dim + 1 - f["dg1"].create_dataset( - "cells", data=numpy.arange(nvertices, dtype=numpy.int64).reshape(-1, ncorners) + _, field_units = _physical_visualisation_values(numpy.ones(1), var.units) + _, coordinate_units = _physical_visualisation_values(numpy.ones(1), mesh.units) + with uw.selective_ranks(0) as should_execute: + if should_execute: + with h5py.File(var_h5_path, "a") as handle: + handle.attrs["storage_frame"] = "physical" + field = handle[f"fields/{var.clean_name}"] + field.attrs["units"] = field_units or "dimensionless" + field.attrs["storage_frame"] = "physical" + field.attrs["degree"] = var.degree + field.attrs["continuous"] = var.continuous + field.attrs["representation"] = ( + "basis_conversion" if direct_dg1 else "exact" ) - uw.mpi.barrier() - - _write_visualisation_metadata(mesh, var, var_h5_path, group) + coordinates = handle["fields/coordinates"] + coordinates.attrs["units"] = coordinate_units or "dimensionless" + coordinates.attrs["storage_frame"] = "physical" + if needs_projection: + projected = handle[f"visualization/{var.clean_name}"] + projected.attrs["units"] = field_units or "dimensionless" + projected.attrs["source_degree"] = var.degree + projected.attrs["visualization_degree"] = 1 if var.continuous else 0 + projected.attrs["representation"] = "projection" + uw.mpi.barrier() def _write_visualisation_geometry(mesh, mesh_h5_path): @@ -9765,52 +9785,6 @@ def _write_visualisation_geometry(mesh, mesh_h5_path): uw.mpi.barrier() -def _write_visualisation_metadata(mesh, var, var_h5_path, group): - """Describe native and physical datasets after collective output closes.""" - import h5py - import underworld3 as uw - from underworld3.function.field_projection import ( - _physical_visualisation_enabled, - _physical_visualisation_values, - ) - - field_enabled = _physical_visualisation_enabled(var.units) - coordinates_enabled = _physical_visualisation_enabled(mesh.units) - if not (field_enabled or coordinates_enabled): - return - - _, field_units = _physical_visualisation_values(numpy.ones(1), var.units) - _, coordinate_units = _physical_visualisation_values( - numpy.ones(1), mesh.units - ) - if group == "dg1": - field_path = "dg1/values" - coordinate_path = "dg1/vertices" - else: - field_path = f"{group}/{var.clean_name}_{var.clean_name}" - coordinate_path = f"{group}/coordinates" - - with uw.selective_ranks(0) as should_execute: - if should_execute: - with h5py.File(var_h5_path, "a") as handle: - handle.attrs["checkpoint_units"] = "nondimensional" - if field_enabled: - handle.attrs["visualisation_units"] = field_units - native = handle[f"fields/{var.clean_name}"] - native.attrs["units"] = "nondimensional" - native.attrs["physical_units"] = field_units - physical = handle[field_path] - physical.attrs["units"] = field_units - - if coordinates_enabled: - handle["fields/coordinates"].attrs["units"] = "nondimensional" - - if coordinate_path in handle and coordinate_units is not None: - coordinates = handle[coordinate_path] - coordinates.attrs["units"] = coordinate_units - uw.mpi.barrier() - - def checkpoint_xdmf( filename: str, meshUpdates: bool = True, @@ -9949,13 +9923,31 @@ def units_information(units, indent): header += """ ]>""" - dg_vars = [var for var in meshVars if not var.continuous and var.degree == 1] + def direct_p2(var): + return ( + var.continuous + and var.degree == 2 + and var.mesh.isSimplex + and var.mesh.dim == 2 + and var.mesh.cdim == 2 + ) + + def direct_dg1(var): + return ( + not var.continuous + and var.degree == 1 + and var.mesh.isSimplex + and var.mesh.dim in (2, 3) + and var.mesh.cdim == var.mesh.dim + ) + + special_vars = [var for var in meshVars if direct_p2(var) or direct_dg1(var)] collection_start = ( '' f'" if dg_vars else "" + collection_end = "" if special_vars else "" xdmf_start = f""" @@ -9990,46 +9982,40 @@ def units_information(units, indent): ## The mesh Var attributes - def get_field_info(h5_filename, mesh_var, center): - """ - Return shape, path, and units for a mesh variable. - Prefers vertex/cell compatibility groups, falls back to /fields layout. - """ - compat_name = f"{mesh_var.clean_name}_{mesh_var.clean_name}" - candidates = [] - - if center == "Cell": - candidates = [f"cell_fields/{compat_name}", f"fields/{mesh_var.clean_name}"] - else: - candidates = [f"vertex_fields/{compat_name}", f"fields/{mesh_var.clean_name}"] - + def get_field_info(h5_filename, dataset_path): + """Return item/component counts and units for one stored field.""" with h5py.File(h5_filename, "r") as f: - for path in candidates: - if path in f: - shp = f[path].shape - units = f[path].attrs.get("units") - if len(shp) == 1: - return shp[0], 1, path, units - return shp[0], shp[1], path, units - - raise RuntimeError( - f"Could not locate data for variable '{mesh_var.clean_name}' in {h5_filename}" - ) + shp = f[dataset_path].shape + units = f[dataset_path].attrs.get("units") + if len(shp) == 1: + return shp[0], 1, units + return shp[0], shp[1], units + + def attribute_kind(var, components): + if var.vtype in (uw.VarType.TENSOR, uw.VarType.SYM_TENSOR) and components == 9: + return "Tensor" + if var.vtype in (uw.VarType.TENSOR, uw.VarType.SYM_TENSOR, uw.VarType.MATRIX): + return "Matrix" + return "Scalar" if components == 1 else "Vector" attributes = "" for var in meshVars: - if not var.continuous and var.degree == 1: + if direct_p2(var) or direct_dg1(var): continue var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" - - # Determine if data is stored on nodes (vertex_fields) or cells (cell_fields) - if not getattr(var, "continuous") or getattr(var, "degree") == 0: + projected = (var.continuous and var.degree > 1) or ( + not var.continuous and var.degree > 0 + ) + dataset_path = ( + f"visualization/{var.clean_name}" + if projected + else f"fields/{var.clean_name}" + ) + if not var.continuous or var.degree == 0: center = "Cell" else: center = "Node" - numItems, numComponents, dataset_path, field_units = get_field_info( - var_filename, var, center - ) + numItems, numComponents, field_units = get_field_info(var_filename, dataset_path) if center == "Node" and numItems != numVertices: warnings.warn( @@ -10044,17 +10030,7 @@ def get_field_info(h5_filename, mesh_var, center): stacklevel=2, ) - # Use variable type when available, but reflect actual stored component count. - if hasattr(var, "vtype") and var.vtype in ( - uw.VarType.TENSOR, - uw.VarType.SYM_TENSOR, - uw.VarType.MATRIX, - ): - variable_type = "Tensor" - elif numComponents == 1: - variable_type = "Scalar" - else: - variable_type = "Vector" + variable_type = attribute_kind(var, numComponents) data_dimensions = f"{numItems}" if numComponents == 1 else f"{numItems} {numComponents}" var_attribute = f""" @@ -10106,55 +10082,43 @@ def get_field_info(h5_filename, mesh_var, center): """ attributes += var_attribute - dg_grid = "" - if dg_vars: - first = dg_vars[0] - first_filename = filename + f".mesh.{first.clean_name}.{index:05}.h5" - with h5py.File(first_filename, "r") as f: - dg_cells = f["dg1/cells"].shape - dg_points = f["dg1/vertices"].shape - dg_geometry_units = f["dg1/vertices"].attrs.get("units") - if dg_cells != (numCells, numCorners) or dg_points != (numCells * numCorners, spaceDim): - raise ValueError("DG1 visualization topology does not match the checkpoint mesh") - dg_grid = f""" - - - - &{first.clean_name}_Data;:/dg1/cells + special_grids = "" + for var in special_vars: + var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" + with h5py.File(var_filename, "r") as f: + cells_shape = f["fields/cells"].shape + points_shape = f["fields/coordinates"].shape + values_shape = f[f"fields/{var.clean_name}"].shape + special_geometry_units = f["fields/coordinates"].attrs.get("units") + field_units = f[f"fields/{var.clean_name}"].attrs.get("units") + special_topology = "Triangle_6" if direct_p2(var) else topology_type + components = values_shape[1] if len(values_shape) == 2 else 1 + kind = attribute_kind(var, components) + dimensions = " ".join(str(value) for value in values_shape) + special_grids += f""" + + + + &{var.clean_name}_Data;:/fields/cells - - &{first.clean_name}_Data;:/dg1/vertices - {units_information(dg_geometry_units, " ")} + + &{var.clean_name}_Data;:/fields/coordinates + {units_information(special_geometry_units, " ")} -""" - for var in dg_vars: - var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" - with h5py.File(var_filename, "r") as f: - shape = f["dg1/values"].shape - field_units = f["dg1/values"].attrs.get("units") - if shape[0] != dg_points[0]: - raise ValueError(f"DG1 visualization size mismatch for {var.clean_name}") - components = shape[1] if len(shape) == 2 else 1 - if var.vtype in (uw.VarType.TENSOR, uw.VarType.SYM_TENSOR): - kind = "Tensor" - elif var.vtype == uw.VarType.MATRIX: - kind = "Matrix" - else: - kind = "Scalar" if components == 1 else "Vector" - dimensions = " ".join(str(value) for value in shape) - dg_grid += f""" - &{var.clean_name}_Data;:/dg1/values + &{var.clean_name}_Data;:/fields/{var.clean_name} {units_information(field_units, " ")} -""" - dg_grid += " " + """ xdmf_end = f""" - {dg_grid} + {special_grids} {collection_end} diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index a2622df28..ebedf35f5 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -1244,6 +1244,7 @@ def read_timestep( varsymbol=r"\cal{S}", ) + field_representation = None if uw.mpi.rank == 0: if verbose: print( @@ -1257,14 +1258,47 @@ def read_timestep( D_src = D_src.reshape(-1, n_components) else: with h5py.File(data_file, "r") as h5f: - X_src = h5f["fields"]["coordinates"][()].reshape(-1, dim) - D_src = h5f["fields"][data_name][()].reshape( - -1, n_components + coordinate_dataset = h5f["fields"]["coordinates"] + field_dataset = h5f["fields"][data_name] + field_representation = field_dataset.attrs.get("representation") + if isinstance(field_representation, bytes): + field_representation = field_representation.decode() + X_src = coordinate_dataset[()].reshape(-1, dim) + D_src = field_dataset[()].reshape(-1, n_components) + storage_frame = field_dataset.attrs.get( + "storage_frame", h5f.attrs.get("storage_frame") ) + if isinstance(storage_frame, bytes): + storage_frame = storage_frame.decode() + if storage_frame == "physical": + coordinate_units = coordinate_dataset.attrs.get("units") + field_units = field_dataset.attrs.get("units") + if isinstance(coordinate_units, bytes): + coordinate_units = coordinate_units.decode() + if isinstance(field_units, bytes): + field_units = field_units.decode() + if coordinate_units and coordinate_units != "dimensionless": + X_src = np.asarray( + uw.non_dimensionalise( + uw.quantity(X_src, coordinate_units) + ) + ) + if field_units and field_units != "dimensionless": + D_src = np.asarray( + uw.non_dimensionalise(uw.quantity(D_src, field_units)) + ) else: X_src = np.empty((0, dim), dtype=np.float64) D_src = np.empty((0, n_components), dtype=np.float64) + field_representation = uw.mpi.comm.bcast(field_representation, root=0) + if field_representation == "basis_conversion": + raise RuntimeError( + "read_timestep cannot invert the element-local DG1 corner " + "basis conversion. Write with petsc_reload=True and use " + "read_checkpoint() for an exact DG1 solver restart." + ) + src_size_before = max(source_swarm.dm.getLocalSize(), 0) source_swarm.add_particles_with_global_coordinates(X_src, migrate=False) source_swarm._invalidate_canonical_data() @@ -1453,11 +1487,12 @@ def read_checkpoint( ): """Load this mesh variable from PETSc reload output. - The default path restores DMPlex section/local-vector data through the + By default, DMPlex section/local-vector data are restored through the topology migration SF, so a mesh reconstructed from its checkpoint may - have a different parallel DOF ordering. Set ``same_layout=True`` only - for an in-place restore onto the exact mesh object that wrote the file; - that path reloads the saved global vector directly. + have a different parallel DOF ordering. With ``same_layout=True``, the + existing PETSc variable vector is loaded directly into the original + layout. Both paths read the same stored values; no duplicate checkpoint + vector is required. This method does not use the coordinate/KDTree remapping provided by ``read_timestep()``. New output should be written with @@ -1472,31 +1507,30 @@ def read_checkpoint( if self._lvec is None: self._set_vec(available=True) - if same_layout: - import h5py + import h5py - if uw.mpi.rank == 0: - with h5py.File(filename, "r") as checkpoint_h5: - has_direct_vector = ( - "uw_checkpoint" in checkpoint_h5 - and data_name in checkpoint_h5["uw_checkpoint"] - ) - else: - has_direct_vector = None - has_direct_vector = uw.mpi.comm.bcast( - has_direct_vector, - root=0, + if uw.mpi.rank == 0: + with h5py.File(filename, "r") as checkpoint_h5: + grouped_checkpoint = "uw_checkpoint/topologies" in checkpoint_h5 + legacy_direct_vector = f"uw_checkpoint/{data_name}" in checkpoint_h5 + else: + grouped_checkpoint = None + legacy_direct_vector = None + grouped_checkpoint = uw.mpi.comm.bcast(grouped_checkpoint, root=0) + legacy_direct_vector = uw.mpi.comm.bcast(legacy_direct_vector, root=0) + + if same_layout and not (grouped_checkpoint or legacy_direct_vector): + raise RuntimeError( + f"{filename} has no direct checkpoint vector for {data_name!r}. " + "Reload it with same_layout=False." ) - if not has_direct_vector: - raise RuntimeError( - f"{filename} has no in-place checkpoint vector for " - f"{data_name!r}. Reload it with same_layout=False." - ) indexset, subdm = self.mesh.dm.createSubDM(self.field_id) sectiondm = self.mesh.dm.clone() viewer = PETSc.ViewerHDF5().create(filename, "r", comm=PETSc.COMM_WORLD) viewer.pushFormat(PETSc.Viewer.Format.HDF5_PETSC) + if grouped_checkpoint and not same_layout: + viewer.pushGroup("/uw_checkpoint") old_mesh_name = self.mesh.dm.getName() old_lvec_name = self._lvec.getName() @@ -1510,16 +1544,15 @@ def read_checkpoint( self._gvec.setName(data_name) if same_layout: - checkpoint_vec = PETSc.Vec().createMPI( - (self._gvec.getLocalSize(), self._gvec.getSize()), - comm=PETSc.COMM_WORLD, + vector_group = ( + f"/uw_checkpoint/topologies/uw_mesh/dms/{data_name}/" + f"vecs/{data_name}" + if grouped_checkpoint + else "/uw_checkpoint" ) - checkpoint_vec.setName(data_name) - viewer.pushGroup("/uw_checkpoint") - checkpoint_vec.load(viewer) + viewer.pushGroup(vector_group) + self._gvec.load(viewer) viewer.popGroup() - self._gvec.array[...] = checkpoint_vec.array_r - checkpoint_vec.destroy() subdm.globalToLocal(self._gvec, self._lvec, addv=False) else: from underworld3.cython.petsc_discretisation import ( @@ -1561,6 +1594,8 @@ def read_checkpoint( self._gvec.setName(old_vec_name) if old_mesh_name is not None: self.mesh.dm.setName(old_mesh_name) + if grouped_checkpoint and not same_layout: + viewer.popGroup() viewer.popFormat() viewer.destroy() sectiondm.destroy() diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py index 44b6fc57b..7a1546a49 100644 --- a/src/underworld3/function/field_projection.py +++ b/src/underworld3/function/field_projection.py @@ -226,6 +226,20 @@ def _write_vec_to_group(viewer, data_array, name, group, comm): vec.destroy() +def _write_index_array_to_group(viewer, data_array, name, group, comm): + """Write a distributed integer connectivity array to an HDF5 group.""" + indices = PETSc.IS().createGeneral( + np.asarray(data_array, dtype=PETSc.IntType).reshape(-1), comm=comm + ) + if data_array.ndim == 2: + indices.setBlockSize(data_array.shape[1]) + indices.setName(name) + viewer.pushGroup(group) + viewer(indices) + viewer.popGroup() + indices.destroy() + + def _physical_visualisation_enabled(units): """Return whether declared units require a physical XDMF copy.""" import underworld3 as uw @@ -337,6 +351,129 @@ def write_vertices_to_viewer( _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) +def write_field_to_viewer( + mesh_var: "MeshVariable", + viewer: "PETSc.ViewerHDF5", + group: str, + name: str, +) -> None: + """Write a variable's owned values in physical units without projection.""" + mesh_var._sync_lvec_to_gvec() + data = mesh_var._gvec.array.reshape(-1, mesh_var.num_components).copy() + data, _ = _physical_visualisation_values(data, mesh_var.units) + _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) + + +def write_field_coordinates_to_viewer( + mesh_var: "MeshVariable", + viewer: "PETSc.ViewerHDF5", + group: str, + name: str = "coordinates", +) -> None: + """Write owned coordinates for a variable's exact finite-element layout.""" + mesh = mesh_var.mesh + coordinate_dm = mesh._basis_coordinate_dm(mesh_var.degree, mesh_var.continuous) + local = coordinate_dm.getLocalVec() + global_vector = coordinate_dm.getGlobalVec() + local.array[...] = np.asarray(mesh_var.coords_nd).reshape(-1) + coordinate_dm.localToGlobal(local, global_vector, addv=False) + coordinates = global_vector.array.reshape(-1, mesh.cdim).copy() + coordinates, _ = _physical_visualisation_values(coordinates, mesh.units) + _write_vec_to_group(viewer, coordinates, name, group, PETSc.COMM_WORLD) + coordinate_dm.restoreGlobalVec(global_vector) + coordinate_dm.restoreLocalVec(local) + coordinate_dm.destroy() + + +def write_projected_field_to_viewer( + mesh_var: "MeshVariable", + viewer: "PETSc.ViewerHDF5", + target_degree: int, + continuous: bool, + group: str, + name: str, +) -> None: + """Write a compact visualization projection for an unsupported layout.""" + data = project_to_degree( + mesh_var, + target_degree=target_degree, + continuous=continuous, + include_ghosts=False, + ) + data, _ = _physical_visualisation_values(data, mesh_var.units) + _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) + + +def write_p2_triangle_topology_to_viewer(mesh_var, viewer, group="/fields"): + """Write VTK-ordered Triangle_6 connectivity for a continuous P2 field.""" + mesh = mesh_var.mesh + if not ( + mesh_var.continuous + and mesh_var.degree == 2 + and mesh.isSimplex + and mesh.dim == 2 + and mesh.cdim == 2 + ): + raise NotImplementedError("direct P2 XDMF currently requires a 2D triangle mesh") + + coordinate_dm = mesh._basis_coordinate_dm(2, True) + local_section = coordinate_dm.getLocalSection() + global_section = coordinate_dm.getGlobalSection() + local_to_global = np.full(len(mesh_var.coords_nd), -1, dtype=PETSc.IntType) + point_start, point_end = local_section.getChart() + for point in range(point_start, point_end): + node_count = local_section.getDof(point) // mesh.cdim + if node_count == 0: + continue + local_offset = local_section.getOffset(point) // mesh.cdim + global_offset = global_section.getOffset(point) + if global_offset < 0: + global_offset = -(global_offset + 1) + global_offset //= mesh.cdim + local_to_global[local_offset : local_offset + node_count] = np.arange( + global_offset, global_offset + node_count, dtype=PETSc.IntType + ) + coordinate_dm.destroy() + + cell_start, cell_end = mesh.dm.getHeightStratum(0) + owned = np.ones(cell_end - cell_start, dtype=bool) + if mesh.dm.comm.getSize() > 1: + _, leaves, remote = mesh.dm.getPointSF().getGraph() + if leaves is None: + leaves = np.arange(len(remote)) + leaves = np.asarray(leaves) + owned[leaves[(leaves >= cell_start) & (leaves < cell_end)] - cell_start] = False + + p2_rows = mesh._cell_node_indices(2, True).reshape(-1, 6)[owned] + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, 3)[owned] + p2_coordinates = mesh_var.coords_nd[p2_rows] + corners = mesh._get_coords_for_basis(1, True)[vertex_rows] + negative = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) < 0 + corners[negative] = corners[negative][:, [0, 2, 1]] + + connectivity = np.empty_like(p2_rows, dtype=PETSc.IntType) + for cell_index, (row, nodes, vertices) in enumerate( + zip(p2_rows, p2_coordinates, corners, strict=True) + ): + targets = np.vstack( + ( + vertices, + 0.5 * (vertices[0] + vertices[1]), + 0.5 * (vertices[1] + vertices[2]), + 0.5 * (vertices[2] + vertices[0]), + ) + ) + order = [np.argmin(np.linalg.norm(nodes - target, axis=1)) for target in targets] + if len(set(order)) != 6: + raise RuntimeError("could not map UW3 P2 nodes to Triangle_6 ordering") + connectivity[cell_index] = local_to_global[row[order]] + if np.any(connectivity < 0): + raise RuntimeError("P2 XDMF connectivity contains an unmapped global node") + _write_index_array_to_group( + viewer, connectivity, "cells", group, PETSc.COMM_WORLD + ) + + def write_coordinates_to_viewer( mesh, viewer: "PETSc.ViewerHDF5", @@ -395,7 +532,14 @@ def write_cell_field_to_viewer( _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) -def _write_dg1_to_viewer(mesh_var, viewer): +def _write_dg1_to_viewer( + mesh_var, + viewer, + group="/dg1", + coordinate_name="vertices", + value_name="values", + repack_tensors=True, +): """Write owned simplex cells with independent vertices and DG1 traces. Coordinate-section cell maps preserve element ownership and node ordering; @@ -430,8 +574,20 @@ def _write_dg1_to_viewer(mesh_var, viewer): local = np.linalg.solve(matrix, (corners - nodes[:, :1]).transpose(0, 2, 1)) weights = np.concatenate((1 - local.sum(axis=1, keepdims=True), local), axis=1) values = np.einsum("cij,cik->cjk", weights, coefficients).reshape(-1, mesh_var.num_components) - values = _repack_tensor_to_paraview(values, mesh_var.vtype, mesh.dim) + if repack_tensors: + values = _repack_tensor_to_paraview(values, mesh_var.vtype, mesh.dim) corners, _ = _physical_visualisation_values(corners, mesh.units) values, _ = _physical_visualisation_values(values, mesh_var.units) - _write_vec_to_group(viewer, corners.reshape(-1, mesh.cdim), "vertices", "/dg1", PETSc.COMM_WORLD) - _write_vec_to_group(viewer, values, "values", "/dg1", PETSc.COMM_WORLD) + corner_rows = corners.reshape(-1, mesh.cdim) + _write_vec_to_group( + viewer, corner_rows, coordinate_name, group, PETSc.COMM_WORLD + ) + _write_vec_to_group(viewer, values, value_name, group, PETSc.COMM_WORLD) + local_count = len(corner_rows) + offset = mesh.dm.comm.tompi4py().exscan(local_count) + if offset is None: + offset = 0 + cells = np.arange(offset, offset + local_count, dtype=PETSc.IntType).reshape( + -1, mesh.dim + 1 + ) + _write_index_array_to_group(viewer, cells, "cells", group, PETSc.COMM_WORLD) diff --git a/tests/test_0003_save_load.py b/tests/test_0003_save_load.py index 0b2d1839b..5b23ef668 100644 --- a/tests/test_0003_save_load.py +++ b/tests/test_0003_save_load.py @@ -172,9 +172,9 @@ def test_timestep_with_petsc_reload_roundtrip(tmp_path): with h5py.File(var_file, "r") as h5f: assert "fields/x" in h5f assert "fields/coordinates" in h5f - assert "vertex_fields/x_x" in h5f - assert "topologies/uw_mesh/dms/x/section" in h5f - assert "topologies/uw_mesh/dms/x/vecs/x" in h5f + assert "vertex_fields" not in h5f + assert "uw_checkpoint/topologies/uw_mesh/dms/x/section" in h5f + assert "uw_checkpoint/topologies/uw_mesh/dms/x/vecs/x" in h5f uw.mpi.barrier() mesh_reloaded = uw.discretisation.Mesh(f"{tmp_path}/unified.mesh.00000.h5") @@ -205,8 +205,8 @@ def test_timestep_with_petsc_reload_roundtrip(tmp_path): assert checkpoint_xdmf_file.is_file() with h5py.File(checkpoint_var_file, "r") as h5f: assert "fields/x" in h5f - assert "vertex_fields/x_x" in h5f - assert "topologies/uw_mesh/dms/x/section" in h5f + assert "vertex_fields" not in h5f + assert "uw_checkpoint/topologies/uw_mesh/dms/x/section" in h5f uw.mpi.barrier() mesh_reloaded = uw.discretisation.Mesh( diff --git a/tests/test_0005_xdmf_compat.py b/tests/test_0005_xdmf_compat.py index 3c7b34d5c..ce5b4d233 100644 --- a/tests/test_0005_xdmf_compat.py +++ b/tests/test_0005_xdmf_compat.py @@ -1,11 +1,8 @@ -"""Test XDMF/HDF5 compatibility groups written by write_timestep. +"""Compact XDMF/HDF5 layouts for continuous and cell fields.""" -Validates that /vertex_fields/ and /cell_fields/ groups are created -correctly, with proper component counts and data values. -""" - -import os +from pathlib import Path import re +import xml.etree.ElementTree as ET import h5py import numpy as np @@ -14,425 +11,229 @@ import underworld3 as uw -# --------------------------------------------------------------------------- -# Helpers -# --------------------------------------------------------------------------- - - -def _check_h5_group_exists(h5_path, group_path): - """Return True if group_path exists in the HDF5 file.""" - with h5py.File(h5_path, "r") as f: - return group_path in f - - -def _read_h5_dataset(h5_path, dataset_path): - """Read and return an HDF5 dataset as a numpy array.""" - with h5py.File(h5_path, "r") as f: - return f[dataset_path][:] - - -def _check_xdmf_refs(xdmf_path, tmp_dir): - """Verify all XDMF entity references point to real HDF5 datasets.""" - with open(xdmf_path, "r") as f: - content = f.read() - - doctype_match = re.search(r"", content, re.DOTALL) - assert doctype_match, "No DOCTYPE entity block found in XDMF file" - - entity_block = doctype_match.group(1) - entities = dict(re.findall(r'', entity_block)) - - # Check both vertex_fields and cell_fields references - refs = re.findall(r"&(\w+);:(/(vertex_fields|cell_fields)/[A-Za-z0-9_]+)", content) - errors = [] - for entity_name, dataset_path, _ in refs: - h5_file = entities.get(entity_name) - if not h5_file: - errors.append(f"Entity {entity_name} not found") - continue - h5_full = os.path.join(tmp_dir, h5_file) - if not os.path.exists(h5_full): - errors.append(f"File {h5_file} not found") - continue - with h5py.File(h5_full, "r") as f: - if dataset_path.lstrip("/") not in f: - errors.append(f"{h5_file}: {dataset_path} missing") - - assert not errors, "XDMF reference errors:\n" + "\n".join(errors) - - -# --------------------------------------------------------------------------- -# Test: P1 scalar + P2 vector (2D) -# --------------------------------------------------------------------------- - - -def test_xdmf_compat_2d(tmp_path): - """write_timestep creates correct compat groups for 2D mesh variables.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - - # P1 scalar, P2 vector - p_var = uw.discretisation.MeshVariable("p", mesh, 1, degree=1) - v_var = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2) - - # Initialise with known values - x, y = mesh.X - p_var.data[:, 0] = mesh._coords[:, 0] # p = x coordinate - v_var.data[:, 0] = 1.0 - v_var.data[:, 1] = 2.0 - - mesh.write_timestep( - "test", index=0, outputPath=str(tmp_path), meshVars=[p_var, v_var] - ) - - # Check XDMF file was created - xdmf_file = os.path.join(str(tmp_path), "test.mesh.00000.xdmf") - assert os.path.exists(xdmf_file), "XDMF file not created" - - # Check P1 scalar: /vertex_fields/p_p should exist with correct shape - p_h5 = os.path.join(str(tmp_path), "test.mesh.p.00000.h5") - assert _check_h5_group_exists(p_h5, "vertex_fields/p_p") - assert _check_h5_group_exists(p_h5, "vertex_fields/coordinates") - - p_compat = _read_h5_dataset(p_h5, "vertex_fields/p_p") - # PETSc writes 1-component scalars as 1D (N,); accept both (N,) and (N,1) - effective_ncomp = p_compat.shape[1] if p_compat.ndim == 2 else 1 - assert effective_ncomp == 1, f"P1 scalar should have 1 component, got {effective_ncomp}" - - # P1 scalar: compat values should match var.data exactly - p_original = _read_h5_dataset(p_h5, "fields/p") - np.testing.assert_allclose( - p_compat.ravel(), p_original.ravel(), atol=1e-10, - err_msg="P1 compat data should match original field data exactly" - ) - - # Check P2 vector: /vertex_fields/v_v should exist with dim components - v_h5 = os.path.join(str(tmp_path), "test.mesh.v.00000.h5") - assert _check_h5_group_exists(v_h5, "vertex_fields/v_v") - - v_compat = _read_h5_dataset(v_h5, "vertex_fields/v_v") - # Standalone Vec writes as 1D — infer components from total size / vertex count - n_verts = mesh._coords.shape[0] - v_total = v_compat.size - v_ncomp = v_total // n_verts if n_verts > 0 else 0 - assert v_ncomp == mesh.dim, ( - f"P2 vector should have {mesh.dim} components, got {v_ncomp}" - ) - - # P2 vector: vertex count should match mesh vertices, not P2 DOFs - coords_compat = _read_h5_dataset(v_h5, "vertex_fields/coordinates") - coords_n_verts = coords_compat.size // mesh.dim - assert coords_n_verts == mesh._coords.shape[0], ( - "Vertex count mismatch between compat coords and mesh" - ) - - # Verify XDMF references are valid - _check_xdmf_refs(xdmf_file, str(tmp_path)) - - del mesh - - -def test_write_timestep_mesh_keeps_restart_and_viz_topology(tmp_path): - """write_timestep mesh output keeps restart data plus XDMF topology.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - s_var = uw.discretisation.MeshVariable("s", mesh, 1, degree=1) - s_var.data[:, 0] = mesh._coords[:, 0] - - mesh.write_timestep( - "viztopo", index=0, outputPath=str(tmp_path), meshVars=[s_var] - ) - - mesh_h5 = os.path.join(str(tmp_path), "viztopo.mesh.00000.h5") - xdmf_file = os.path.join(str(tmp_path), "viztopo.mesh.00000.xdmf") - - with h5py.File(mesh_h5, "r") as h5f: - assert "labels" in h5f - assert "topology/cells" in h5f - assert "topology/cones" in h5f - assert "viz/topology/cells" in h5f - cells = h5f["viz/topology/cells"] - assert len(cells.shape) == 2 - assert cells.shape[1] > 1 - - with open(xdmf_file, "r") as f: - xdmf_text = f.read() - assert "&MeshData;:/viz/topology/cells" in xdmf_text - assert "&MeshData;:/topology/cells" not in xdmf_text - assert 'ItemType="HyperSlab"' not in xdmf_text - - del mesh - - -def test_write_checkpoint_mesh_uses_petsc_topology(tmp_path): - """write_checkpoint keeps PETSc DMPlex topology for restart output.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - s_var = uw.discretisation.MeshVariable("s", mesh, 1, degree=1) - s_var.data[:, 0] = mesh._coords[:, 0] - - with pytest.warns(FutureWarning, match="write_checkpoint\\(\\) is deprecated"): - mesh.write_checkpoint( - "restart", - outputPath=str(tmp_path), - meshUpdates=True, - meshVars=[s_var], - index=0, - ) - - mesh_h5 = os.path.join(str(tmp_path), "restart.mesh.00000.h5") - - with h5py.File(mesh_h5, "r") as h5f: - assert "topologies/uw_mesh/topology" in h5f - assert "viz/topology/cells" not in h5f - - del mesh - - -# --------------------------------------------------------------------------- -# Test: 3D mesh -# --------------------------------------------------------------------------- - - -def test_xdmf_compat_3d(tmp_path): - """write_timestep creates correct compat groups for 3D mesh.""" - - mesh = uw.meshing.StructuredQuadBox( - elementRes=(3, 3, 3), - minCoords=(0.0, 0.0, 0.0), - maxCoords=(1.0, 1.0, 1.0), - ) - - s_var = uw.discretisation.MeshVariable("s", mesh, 1, degree=1) - s_var.data[:, 0] = mesh._coords[:, 2] # s = z - - mesh.write_timestep( - "test3d", index=0, outputPath=str(tmp_path), meshVars=[s_var] - ) - - s_h5 = os.path.join(str(tmp_path), "test3d.mesh.s.00000.h5") - assert _check_h5_group_exists(s_h5, "vertex_fields/s_s") - - s_compat = _read_h5_dataset(s_h5, "vertex_fields/s_s") - s_original = _read_h5_dataset(s_h5, "fields/s") - np.testing.assert_allclose(s_compat.ravel(), s_original.ravel(), atol=1e-10) - - # Verify XDMF - xdmf_file = os.path.join(str(tmp_path), "test3d.mesh.00000.xdmf") - _check_xdmf_refs(xdmf_file, str(tmp_path)) - - mesh_h5 = os.path.join(str(tmp_path), "test3d.mesh.00000.h5") - with h5py.File(mesh_h5, "r") as h5f: - assert "viz/topology/cells" in h5f - cells = h5f["viz/topology/cells"] - assert len(cells.shape) == 2 - assert cells.shape[1] > 1 - - with open(xdmf_file, "r") as f: - xdmf_text = f.read() - assert "&MeshData;:/viz/topology/cells" in xdmf_text - assert "&MeshData;:/topology/cells" not in xdmf_text - assert f'Dimensions="{s_compat.shape[0]}"' in xdmf_text +_OLD_GROUPS = ("vertex_fields", "cell_fields", "dg1") - del mesh +def _shared_path(tmp_path): + return Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) -# --------------------------------------------------------------------------- -# Test: Cell (discontinuous) variable -# --------------------------------------------------------------------------- +def _assert_xdmf_references_exist(xdmf_path): + """Check every HDF reference and advertised shape in an XDMF file.""" + text = xdmf_path.read_text() + entities = dict(re.findall(r'', text)) + tree = ET.parse(xdmf_path) + for item in tree.findall(".//DataItem[@Format='HDF']"): + reference, dataset = item.text.strip().split(":", 1) + entity = reference.removeprefix("&").removesuffix(";") + h5_name = entities.get(entity, reference) + with h5py.File(xdmf_path.parent / h5_name, "r") as handle: + assert dataset in handle + dimensions = item.get("Dimensions") + if dimensions: + assert tuple(map(int, dimensions.split())) == handle[dataset].shape -def test_xdmf_compat_cell_variable(tmp_path): - """Cell variables go to /cell_fields/ group.""" - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) +def _assert_no_compatibility_copies(handle): + for group in _OLD_GROUPS: + assert group not in handle - # Discontinuous (cell-centred) variable - c_var = uw.discretisation.MeshVariable( - "c", mesh, 1, degree=0, continuous=False - ) - c_var.data[:, 0] = 42.0 - - mesh.write_timestep( - "testcell", index=0, outputPath=str(tmp_path), meshVars=[c_var] - ) - - c_h5 = os.path.join(str(tmp_path), "testcell.mesh.c.00000.h5") - assert _check_h5_group_exists(c_h5, "cell_fields/c_c"), ( - "/cell_fields/c_c group should exist for discontinuous variable" - ) - - c_compat = _read_h5_dataset(c_h5, "cell_fields/c_c") - c_effective_ncomp = c_compat.shape[1] if c_compat.ndim == 2 else 1 - assert c_effective_ncomp == 1 - - del mesh - - -# --------------------------------------------------------------------------- -# Test: read_timestep round-trip still works (checkpoint integrity) -# --------------------------------------------------------------------------- - - -def test_xdmf_checkpoint_roundtrip(tmp_path): - """Compat groups don't corrupt the /fields/ data used by read_timestep.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_direct_p1_p2_and_dg0_use_fields_only(tmp_path): + """P1, triangular P2, and DG0 are visualized directly from /fields.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.4, regular=True, qdegree=3) p1 = uw.discretisation.MeshVariable("p1", mesh, 1, degree=1) - p1.data[:, 0] = mesh._coords[:, 0] + mesh._coords[:, 1] - + p2 = uw.discretisation.MeshVariable("p2", mesh, 2, degree=2) + dg0 = uw.discretisation.MeshVariable("dg0", mesh, 1, degree=0, continuous=False) + + p1.array[:, 0, 0] = p1.coords[:, 0] + 2.0 * p1.coords[:, 1] + p2.array[:, 0, 0] = p2.coords[:, 0] + p2.array[:, 0, 1] = p2.coords[:, 1] + dg0.array[:, 0, 0] = 7.0 + mesh.write_timestep("compact", 0, outputPath=str(directory), meshVars=[p1, p2, dg0]) + + if uw.mpi.rank == 0: + p1_file = directory / "compact.mesh.p1.00000.h5" + p2_file = directory / "compact.mesh.p2.00000.h5" + dg0_file = directory / "compact.mesh.dg0.00000.h5" + with h5py.File(p1_file, "r") as handle: + assert set(handle) == {"fields"} + assert set(handle["fields"]) == {"coordinates", "p1"} + assert "visualization" not in handle + _assert_no_compatibility_copies(handle) + with h5py.File(dg0_file, "r") as handle: + assert set(handle) == {"fields"} + assert handle["fields/dg0"].shape[0] == handle["fields/coordinates"].shape[0] + assert "visualization" not in handle + _assert_no_compatibility_copies(handle) + with h5py.File(p2_file, "r") as handle: + assert set(handle) == {"fields"} + coordinates = handle["fields/coordinates"][:] + cells = handle["fields/cells"][:] + values = handle["fields/p2"][:] + assert cells.shape[1] == 6 + assert values.shape == coordinates.shape + np.testing.assert_allclose(values, coordinates, atol=1.0e-12) + points = coordinates[cells] + np.testing.assert_allclose(points[:, 3], 0.5 * (points[:, 0] + points[:, 1])) + np.testing.assert_allclose(points[:, 4], 0.5 * (points[:, 1] + points[:, 2])) + np.testing.assert_allclose(points[:, 5], 0.5 * (points[:, 2] + points[:, 0])) + _assert_no_compatibility_copies(handle) + + xdmf = directory / "compact.mesh.00000.xdmf" + text = xdmf.read_text() + assert 'TopologyType="Triangle_6"' in text + assert "&p1_Data;:/fields/p1" in text + assert "&p2_Data;:/fields/p2" in text + assert "&dg0_Data;:/fields/dg0" in text + _assert_xdmf_references_exist(xdmf) + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_higher_order_fields_keep_exact_data_and_one_compact_reduction(tmp_path): + """P3+ uses P1 and DG2+ uses DG0 only for the XDMF view.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.4, regular=True, qdegree=5) + continuous = [ + uw.discretisation.MeshVariable(f"p{degree}", mesh, 1, degree=degree) for degree in (3, 4) + ] + discontinuous = [ + uw.discretisation.MeshVariable(f"dg{degree}", mesh, 1, degree=degree, continuous=False) + for degree in (2, 3, 4) + ] + for field in continuous: + field.array[:, 0, 0] = 2.0 + for field in discontinuous: + field.array[:, 0, 0] = 3.0 + mesh.write_timestep("high", 0, outputPath=str(directory), meshVars=continuous + discontinuous) + + if uw.mpi.rank == 0: + mesh_file = directory / "high.mesh.00000.h5" + with h5py.File(mesh_file, "r") as handle: + vertex_count = handle["geometry/vertices"].shape[0] + cell_count = handle["viz/topology/cells"].shape[0] + for field in continuous: + with h5py.File(directory / f"high.mesh.{field.clean_name}.00000.h5", "r") as handle: + assert handle[f"fields/{field.clean_name}"].shape[0] > vertex_count + assert handle[f"visualization/{field.clean_name}"].shape[0] == vertex_count + np.testing.assert_allclose(handle[f"visualization/{field.clean_name}"][:], 2.0) + assert ( + handle[f"visualization/{field.clean_name}"].attrs["visualization_degree"] == 1 + ) + _assert_no_compatibility_copies(handle) + for field in discontinuous: + with h5py.File(directory / f"high.mesh.{field.clean_name}.00000.h5", "r") as handle: + assert handle[f"fields/{field.clean_name}"].shape[0] > cell_count + assert handle[f"visualization/{field.clean_name}"].shape[0] == cell_count + np.testing.assert_allclose(handle[f"visualization/{field.clean_name}"][:], 3.0) + assert ( + handle[f"visualization/{field.clean_name}"].attrs["visualization_degree"] == 0 + ) + _assert_no_compatibility_copies(handle) + + xdmf = directory / "high.mesh.00000.xdmf" + text = xdmf.read_text() + for field in continuous + discontinuous: + assert f"&{field.clean_name}_Data;:/visualization/{field.clean_name}" in text + _assert_xdmf_references_exist(xdmf) + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_compact_tensor_keeps_native_component_count(tmp_path): + """The authoritative field avoids a second nine-component tensor copy.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + tensor = uw.discretisation.MeshVariable("tensor", mesh, degree=1, vtype=uw.VarType.TENSOR) + tensor.array[:, 0, 0] = 1.0 + tensor.array[:, 0, 1] = 2.0 + tensor.array[:, 1, 0] = 3.0 + tensor.array[:, 1, 1] = 4.0 + mesh.write_timestep("tensor", 0, outputPath=str(directory), meshVars=[tensor]) + + if uw.mpi.rank == 0: + with h5py.File(directory / "tensor.mesh.tensor.00000.h5", "r") as handle: + assert handle["fields/tensor"].shape[1] == 4 + assert "visualization" not in handle + _assert_no_compatibility_copies(handle) + text = (directory / "tensor.mesh.00000.xdmf").read_text() + assert 'Type="Matrix"' in text + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_timestep_and_optional_checkpoint_roundtrip(tmp_path): + """Dimensional fields support remap; /uw_checkpoint supports exact reload.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.4, regular=True) + source = uw.discretisation.MeshVariable("source", mesh, 1, degree=2) + source.array[:, 0, 0] = source.coords[:, 0] + 2.0 * source.coords[:, 1] + same_layout_expected = np.array(source.array) mesh.write_timestep( - "roundtrip", index=0, outputPath=str(tmp_path), meshVars=[p1] - ) - - # Read back - p1_check = uw.discretisation.MeshVariable("p1check", mesh, 1, degree=1) - p1_check.read_timestep("roundtrip", "p1", 0, outputPath=str(tmp_path)) - - np.testing.assert_allclose(p1.data, p1_check.data, atol=1e-10) - - del mesh - - -# --------------------------------------------------------------------------- -# Test: create_xdmf=False skips compat groups -# --------------------------------------------------------------------------- - - -def test_no_xdmf_when_disabled(tmp_path): - """create_xdmf=False should not create compat groups or XDMF file.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - s = uw.discretisation.MeshVariable("s", mesh, 1, degree=1) - s.data[:, 0] = 1.0 # initialise so gvec exists - - mesh.write_timestep( - "noxdmf", index=0, outputPath=str(tmp_path), - meshVars=[s], create_xdmf=False, - ) - - # XDMF file should not be generated. - # Note: we do NOT check HDF5 groups because some PETSc versions (3.21+) - # write /vertex_fields/ automatically during var.write() — that is PETSc - # behaviour, not ours. - xdmf_file = os.path.join(str(tmp_path), "noxdmf.mesh.00000.xdmf") - assert not os.path.exists(xdmf_file), "XDMF file should not exist" - - del mesh - - -# --------------------------------------------------------------------------- -# Test: Tensor variable repacking -# --------------------------------------------------------------------------- - - -def test_tensor_variable_repacking(tmp_path): - """Tensor variables are repacked to 9 components (3x3) in vertex_fields.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - - # 2D full tensor (4 components internally) - T = uw.discretisation.MeshVariable( - "T", mesh, mesh.dim, degree=1, vtype=uw.VarType.TENSOR - ) - T.data[:, 0] = 1.0 # xx - T.data[:, 1] = 0.1 # xy - T.data[:, 2] = 0.2 # yx - T.data[:, 3] = 2.0 # yy - - mesh.write_timestep( - "tensor", index=0, outputPath=str(tmp_path), meshVars=[T] - ) - - t_h5 = os.path.join(str(tmp_path), "tensor.mesh.T.00000.h5") - assert _check_h5_group_exists(t_h5, "vertex_fields/T_T") - - t_compat = _read_h5_dataset(t_h5, "vertex_fields/T_T") - n_verts = mesh._coords.shape[0] - t_ncomp = t_compat.size // n_verts - assert t_ncomp == 9, f"Tensor should be repacked to 9 components, got {t_ncomp}" - - # Verify XDMF - xdmf_file = os.path.join(str(tmp_path), "tensor.mesh.00000.xdmf") - _check_xdmf_refs(xdmf_file, str(tmp_path)) - - del mesh - - -# --------------------------------------------------------------------------- -# Test: Valid /viz/topology connectivity generation -# --------------------------------------------------------------------------- - - -def test_xdmf_viz_topology_written_correctly(tmp_path): - """write_timestep uses PETSc-native viz topology and valid geometry.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(3, 3)) - - # Write just the mesh (no vars needed to test mesh topology) - mesh.write_timestep("test_topo", index=0, outputPath=str(tmp_path)) - - mesh_h5 = os.path.join(str(tmp_path), "test_topo.mesh.00000.h5") - assert _check_h5_group_exists(mesh_h5, "viz/topology/cells"), ( - "Mesh HDF5 must contain the /viz/topology/cells dataset" - ) - assert _check_h5_group_exists(mesh_h5, "geometry/vertices"), ( - "Mesh HDF5 must contain the /geometry/vertices dataset" - ) - - # Validate cell-to-vertex connectivity bounds - cells = _read_h5_dataset(mesh_h5, "viz/topology/cells") - vertices = _read_h5_dataset(mesh_h5, "geometry/vertices") - - num_vertices = vertices.shape[0] - assert cells.max() < num_vertices, ( - f"Invalid topology: cells max ({cells.max()}) must be < numVertices ({num_vertices})" - ) - assert cells.min() >= 0, "Invalid topology: cells indices cannot be negative" - - # Validate XDMF actually points to the viz group, not the raw DMPlex group - xdmf_file = os.path.join(str(tmp_path), "test_topo.mesh.00000.xdmf") - assert os.path.exists(xdmf_file), "XDMF file should exist" - - with open(xdmf_file, "r") as f: - xdmf_content = f.read() - - assert "/viz/topology/cells" in xdmf_content, ( - "XDMF file should explicitly point to /viz/topology/cells" - ) - assert "/geometry/vertices" in xdmf_content, ( - "XDMF file should explicitly point to /geometry/vertices" - ) - - del mesh - - -def test_xdmf_viz_topology_3d_written_correctly(tmp_path): - """write_timestep uses PETSc-native viz topology for 3D meshes.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2, 2)) - - # Write just the mesh - mesh.write_timestep("test_topo_3d", index=0, outputPath=str(tmp_path)) - - mesh_h5 = os.path.join(str(tmp_path), "test_topo_3d.mesh.00000.h5") - assert _check_h5_group_exists(mesh_h5, "viz/topology/cells"), ( - "3D Mesh HDF5 must contain the /viz/topology/cells dataset" - ) - - # Validate connectivity - cells = _read_h5_dataset(mesh_h5, "viz/topology/cells") - vertices = _read_h5_dataset(mesh_h5, "geometry/vertices") - - num_vertices = vertices.shape[0] - assert cells.max() < num_vertices, "Invalid 3D topology bounds" - - xdmf_file = os.path.join(str(tmp_path), "test_topo_3d.mesh.00000.xdmf") - assert os.path.exists(xdmf_file) - with open(xdmf_file, "r") as f: - xdmf_content = f.read() - assert "/viz/topology/cells" in xdmf_content - assert "/geometry/vertices" in xdmf_content - - del mesh + "roundtrip", 0, outputPath=str(directory), meshVars=[source], petsc_reload=True + ) + + field_file = directory / "roundtrip.mesh.source.00000.h5" + source.array[:, 0, 0] = 0.0 + source.read_checkpoint(str(field_file), data_name="source", same_layout=True) + np.testing.assert_allclose(source.array, same_layout_expected, atol=1.0e-12) + + reloaded_mesh = uw.discretisation.Mesh(str(directory / "roundtrip.mesh.00000.h5")) + restored = uw.discretisation.MeshVariable("source", reloaded_mesh, 1, degree=2) + restored.read_checkpoint(str(field_file), data_name="source") + expected = restored.coords[:, 0] + 2.0 * restored.coords[:, 1] + np.testing.assert_allclose(restored.array[:, 0, 0], expected, atol=1.0e-12) + + remapped = uw.discretisation.MeshVariable("remapped", mesh, 1, degree=2) + remapped.read_timestep("roundtrip", "source", 0, outputPath=str(directory)) + expected = remapped.coords[:, 0] + 2.0 * remapped.coords[:, 1] + np.testing.assert_allclose(remapped.array[:, 0, 0], expected, atol=1.0e-12) + + if uw.mpi.rank == 0: + with h5py.File(field_file, "r") as handle: + assert set(handle) == {"fields", "uw_checkpoint"} + assert "uw_checkpoint/topologies/uw_mesh/dms/source/vecs/source/source" in handle + assert set(handle["uw_checkpoint/topologies/uw_mesh/dms"]) == { + "source", + "uw_mesh", + } + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_create_xdmf_false_preserves_native_writer(tmp_path): + """Disabling XDMF leaves the established native field file untouched.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + field = uw.discretisation.MeshVariable("field", mesh, 1, degree=1) + field.array[:, 0, 0] = 1.0 + mesh.write_timestep("native", 0, outputPath=str(directory), meshVars=[field], create_xdmf=False) + if uw.mpi.rank == 0: + assert not (directory / "native.mesh.00000.xdmf").exists() + with h5py.File(directory / "native.mesh.field.00000.h5", "r") as handle: + assert "fields/field" in handle + assert "visualization" not in handle + + +@pytest.mark.level_1 +@pytest.mark.tier_b +@pytest.mark.parametrize("dim", [2, 3]) +def test_mesh_xdmf_topology_is_direct_and_bounded(tmp_path, dim): + """The base mesh XDMF uses valid direct cell-to-vertex connectivity.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2,) * dim) + mesh.write_timestep("topology", 0, outputPath=str(directory)) + if uw.mpi.rank == 0: + with h5py.File(directory / "topology.mesh.00000.h5", "r") as handle: + cells = handle["viz/topology/cells"][:] + vertices = handle["geometry/vertices"][:] + assert cells.ndim == 2 + assert cells.min() >= 0 + assert cells.max() < len(vertices) + text = (directory / "topology.mesh.00000.xdmf").read_text() + assert "&MeshData;:/viz/topology/cells" in text + assert "&MeshData;:/topology/cells" not in text diff --git a/tests/test_0005_xdmf_dg1.py b/tests/test_0005_xdmf_dg1.py index 163aea3ca..b99f03a91 100644 --- a/tests/test_0005_xdmf_dg1.py +++ b/tests/test_0005_xdmf_dg1.py @@ -1,4 +1,4 @@ -"""DG1 visualization preserves affine fields and jumps, including MPI ownership.""" +"""Direct DG1 XDMF output preserves element-local affine traces.""" from pathlib import Path import xml.etree.ElementTree as ET @@ -6,6 +6,7 @@ import h5py import numpy as np import pytest + import underworld3 as uw @@ -23,120 +24,112 @@ def test_dg1_simplex_output(tmp_path, dim): ) scalar = uw.discretisation.MeshVariable("dg_scalar", mesh, 1, degree=1, continuous=False) tensor = uw.discretisation.MeshVariable( - "dg_tensor", - mesh, - degree=1, - continuous=False, - vtype=uw.VarType.TENSOR, + "dg_tensor", mesh, degree=1, continuous=False, vtype=uw.VarType.TENSOR ) pressure = uw.discretisation.MeshVariable("pressure", mesh, 1, degree=1) - vector = uw.discretisation.MeshVariable("dg_vector", mesh, dim, degree=1, continuous=False) - symmetric = uw.discretisation.MeshVariable( - "dg_symmetric", mesh, degree=1, continuous=False, vtype=uw.VarType.SYM_TENSOR - ) rows = mesh._cell_node_indices(1, False).reshape(-1, dim + 1) coords = scalar.coords - offset = np.floor(coords[rows].mean(axis=1)[:, 0] * 7 + 1e-8) - scalar.array[rows, 0, 0] = 1 + coords[rows, 0] + 2 * coords[rows, 1] + offset[:, None] + offsets = np.floor(coords[rows].mean(axis=1)[:, 0] * 7 + 1.0e-8) + scalar.array[rows, 0, 0] = 1 + coords[rows, 0] + 2 * coords[rows, 1] + offsets[:, None] tensor.array[:] = 0 tensor.array[:, 0, 0] = scalar.array[:, 0, 0] tensor.array[:, 0, 1] = 3 + coords[:, 0] tensor.array[:, 1, 0] = -2 + coords[:, 1] tensor.array[:, 1, 1] = 5 pressure.array[:, 0, 0] = pressure.coords[:, 0] - vector.array[:, 0, :] = coords - symmetric.array[:] = 0 - symmetric.array[:, 0, 0] = 2 - symmetric.array[:, 1, 1] = 3 - symmetric.array[:, 0, 1] = coords[:, 0] - original = np.array(scalar.array) mesh.write_timestep( "fields", - index=0, + 0, outputPath=str(directory), - meshVars=[pressure, scalar, tensor, vector, symmetric], + meshVars=[pressure, scalar, tensor], petsc_reload=True, ) - restored = uw.discretisation.MeshVariable("restored", mesh, 1, degree=1, continuous=False) + reloaded_mesh = uw.discretisation.Mesh(str(directory / "fields.mesh.00000.h5")) + restored = uw.discretisation.MeshVariable( + "dg_scalar", reloaded_mesh, 1, degree=1, continuous=False + ) restored.read_checkpoint( - str(directory / "fields.mesh.dg_scalar.00000.h5"), data_name="dg_scalar", same_layout=True + str(directory / "fields.mesh.dg_scalar.00000.h5"), data_name="dg_scalar" + ) + restored_rows = reloaded_mesh._cell_node_indices(1, False).reshape(-1, dim + 1) + restored_coords = restored.coords + restored_offsets = np.floor(restored_coords[restored_rows].mean(axis=1)[:, 0] * 7 + 1.0e-8) + expected_restored = ( + 1 + + restored_coords[restored_rows, 0] + + 2 * restored_coords[restored_rows, 1] + + restored_offsets[:, None] ) - np.testing.assert_allclose(restored.array, original, rtol=1e-12, atol=1e-12) + np.testing.assert_allclose( + restored.array[restored_rows, 0, 0], + expected_restored, + rtol=1.0e-12, + atol=1.0e-12, + ) + + remapped = uw.discretisation.MeshVariable("remapped", mesh, 1, degree=1, continuous=False) + with pytest.raises(RuntimeError, match="DG1 corner basis conversion"): + remapped.read_timestep("fields", "dg_scalar", 0, outputPath=str(directory)) + if uw.mpi.rank != 0: return with h5py.File(directory / "fields.mesh.dg_scalar.00000.h5", "r") as handle: - points = handle["dg1/vertices"][:] - cells = handle["dg1/cells"][:] - values = handle["dg1/values"][:].reshape(-1) - native = handle["fields/dg_scalar"][:] + points = handle["fields/coordinates"][:] + cells = handle["fields/cells"][:] + values = handle["fields/dg_scalar"][:].reshape(-1) + assert handle["fields/dg_scalar"].attrs["representation"] == "basis_conversion" assert len(points) == len(cells) * (dim + 1) assert len(np.unique(cells)) == len(points) - assert native.size == len(points) + assert "dg1" not in handle + assert "visualization" not in handle with h5py.File(directory / "fields.mesh.00000.h5", "r") as handle: assert len(cells) == len(handle["viz/topology/cells"]) + corners = points[cells] - assert np.all(np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) > 0) + determinants = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) + assert np.all(determinants > 0) centers = corners.mean(axis=1) - assert len(np.unique(np.round(centers, 10), axis=0)) == len(cells) - offset = np.repeat(np.floor(centers[:, 0] * 7 + 1e-8), dim + 1) - expected = 1 + points[:, 0] + 2 * points[:, 1] + offset - np.testing.assert_allclose(values, expected, rtol=1e-12, atol=1e-12) - # Repeated positions can have different values: these traces must not merge. + offsets = np.repeat(np.floor(centers[:, 0] * 7 + 1.0e-8), dim + 1) + expected = 1 + points[:, 0] + 2 * points[:, 1] + offsets + np.testing.assert_allclose(values, expected, rtol=1.0e-12, atol=1.0e-12) + _, inverse = np.unique(np.round(points, 10), axis=0, return_inverse=True) low = np.full(inverse.max() + 1, np.inf) high = np.full(inverse.max() + 1, -np.inf) np.minimum.at(low, inverse, values) np.maximum.at(high, inverse, values) assert np.max(high - low) >= 1 + with h5py.File(directory / "fields.mesh.dg_tensor.00000.h5", "r") as handle: - tensor_values = handle["dg1/values"][:] - np.testing.assert_allclose(handle["dg1/vertices"][:], points) - assert tensor_values.shape == (len(points), 9) + tensor_values = handle["fields/dg_tensor"][:] + assert tensor_values.shape == (len(points), dim * dim) np.testing.assert_allclose(tensor_values[:, 0], expected) np.testing.assert_allclose(tensor_values[:, 1], 3 + points[:, 0]) - np.testing.assert_allclose(tensor_values[:, 3], -2 + points[:, 1], atol=1e-12) - np.testing.assert_allclose(tensor_values[:, 4], 5) + np.testing.assert_allclose(tensor_values[:, dim], -2 + points[:, 1], atol=1.0e-12) + np.testing.assert_allclose(tensor_values[:, dim + 1], 5) + tree = ET.parse(directory / "fields.mesh.00000.xdmf") - with h5py.File(directory / "fields.mesh.dg_vector.00000.h5", "r") as handle: - np.testing.assert_allclose(handle["dg1/values"][:], points, atol=1e-12) - with h5py.File(directory / "fields.mesh.dg_symmetric.00000.h5", "r") as handle: - sym_values = handle["dg1/values"][:] - assert sym_values.shape == (len(points), 9) - np.testing.assert_allclose(sym_values[:, 0], 2) - np.testing.assert_allclose(sym_values[:, 4], 3) - np.testing.assert_allclose(sym_values[:, 1], points[:, 0], atol=1e-12) - np.testing.assert_allclose(sym_values[:, 3], points[:, 0], atol=1e-12) grids = tree.findall(".//Grid[@GridType='Uniform']") - assert len(grids) == 2 - dg = next(grid for grid in grids if grid.get("Name") == "DG1") - assert {a.get("Name") for a in dg.findall("Attribute")} == { - "dg_scalar", - "dg_tensor", - "dg_vector", - "dg_symmetric", - } - assert all(a.get("Center") == "Node" for a in dg.findall("Attribute")) - for item in tree.findall(".//DataItem[@Format='HDF']"): - filename, dataset = item.text.strip().split(":", 1) - with h5py.File(directory / filename, "r") as handle: - assert tuple(map(int, item.get("Dimensions").split())) == handle[dataset].shape + assert {grid.get("Name") for grid in grids} == {"domain", "dg_scalar", "dg_tensor"} + for name in ("dg_scalar", "dg_tensor"): + grid = next(grid for grid in grids if grid.get("Name") == name) + attribute = grid.find("Attribute") + assert attribute.get("Name") == name + assert attribute.get("Center") == "Node" @pytest.mark.level_1 @pytest.mark.tier_b -@pytest.mark.parametrize("degree", [1, 2]) -def test_unsupported_dg_layout_fails_before_writing(tmp_path, degree): +def test_unsupported_dg_layout_uses_dg0_visualization_reduction(tmp_path): + """Non-simplex DG1 remains exact in /fields and gets one DG0 XDMF array.""" directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) - mesh = ( - uw.meshing.StructuredQuadBox(elementRes=(2, 2)) - if degree == 1 - else uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) - ) - dg = uw.discretisation.MeshVariable("dg", mesh, 1, degree=degree, continuous=False) - with pytest.raises(NotImplementedError, match="DG.*XDMF"): - mesh.write_timestep("unsupported", index=0, outputPath=str(directory), meshVars=[dg]) - assert not (directory / "unsupported.mesh.00000.h5").exists() - # Unsupported visualization must not prevent native-only checkpoints. - mesh.write_timestep( - "native", index=0, outputPath=str(directory), meshVars=[dg], create_xdmf=False - ) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + dg = uw.discretisation.MeshVariable("dg", mesh, 1, degree=1, continuous=False) + dg.array[:, 0, 0] = 4.0 + mesh.write_timestep("structured", 0, outputPath=str(directory), meshVars=[dg]) + if uw.mpi.rank == 0: + with h5py.File(directory / "structured.mesh.dg.00000.h5", "r") as handle: + assert "fields/dg" in handle + assert "visualization/dg" in handle + np.testing.assert_allclose(handle["visualization/dg"][:], 4.0) + assert handle["visualization/dg"].attrs["visualization_degree"] == 0 diff --git a/tests/test_0005_xdmf_physical_units.py b/tests/test_0005_xdmf_physical_units.py index 00c743851..7e84dae1f 100644 --- a/tests/test_0005_xdmf_physical_units.py +++ b/tests/test_0005_xdmf_physical_units.py @@ -1,4 +1,4 @@ -"""Physical XDMF output keeps native checkpoint arrays nondimensional.""" +"""Dimensional /fields output and optional native PETSc checkpoints.""" from pathlib import Path @@ -10,10 +10,9 @@ def _set_reference_scales(): - """Use exact scales with easy-to-check physical conversions.""" - orchestration_model = uw.get_default_model() - orchestration_model.set_scaling_mode("exact") - orchestration_model.set_reference_quantities( + model = uw.get_default_model() + model.set_scaling_mode("exact") + model.set_reference_quantities( length=uw.quantity(10, "km"), velocity=uw.quantity(5, "mm/year"), pressure=uw.quantity(2, "MPa"), @@ -22,132 +21,96 @@ def _set_reference_scales(): @pytest.mark.level_1 @pytest.mark.tier_b -def test_xdmf_uses_declared_physical_units(tmp_path): - """Visualisation copies are physical while restart data stays native.""" +def test_fields_are_dimensional_and_checkpoint_is_native(tmp_path): + """Analysis sees physical values while exact reload restores solver values.""" _set_reference_scales() - - mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) - velocity = uw.discretisation.MeshVariable( - "velocity", mesh, mesh.dim, degree=2, units="mm/year" - ) + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) + velocity = uw.discretisation.MeshVariable("velocity", mesh, mesh.dim, degree=2, units="mm/year") pressure = uw.discretisation.MeshVariable( "pressure", mesh, 1, degree=0, continuous=False, units="MPa" ) - surface = uw.meshing.Surface( - "unit_test", - mesh, - control_points=uw.quantity([[0.0, 0.0], [10.0, 0.0]], "km"), - ) - surface.discretize() - distance = surface.abs_distance - velocity.data[:, 0] = 2.0 - velocity.data[:, 1] = 3.0 - pressure.data[:, 0] = 4.0 + velocity.array[:, 0, 0] = uw.quantity(10.0, "mm/year") + velocity.array[:, 0, 1] = uw.quantity(15.0, "mm/year") + pressure.array[:, 0, 0] = uw.quantity(8.0, "MPa") + expected_velocity = np.array(velocity.array) + expected_pressure = np.array(pressure.array) directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) mesh.write_timestep( "physical", - index=0, + 0, outputPath=str(directory), - meshVars=[velocity, pressure, distance], + meshVars=[velocity, pressure], petsc_reload=True, ) - with uw.selective_ranks(0) as should_execute: - if not should_execute: - return - - mesh_file = directory / "physical.mesh.00000.h5" velocity_file = directory / "physical.mesh.velocity.00000.h5" pressure_file = directory / "physical.mesh.pressure.00000.h5" - distance_file = directory / "physical.mesh.surf_unit_test_absdistance.00000.h5" + velocity_restart = uw.discretisation.MeshVariable( + "velocity_restart", mesh, mesh.dim, degree=2, units="mm/year" + ) + pressure_restart = uw.discretisation.MeshVariable( + "pressure_restart", mesh, 1, degree=0, continuous=False, units="MPa" + ) + velocity_restart.read_checkpoint(str(velocity_file), data_name="velocity") + pressure_restart.read_checkpoint(str(pressure_file), data_name="pressure") + np.testing.assert_allclose(np.array(velocity_restart.array), expected_velocity) + np.testing.assert_allclose(np.array(pressure_restart.array), expected_pressure) + + velocity_remap = uw.discretisation.MeshVariable( + "velocity_remap", mesh, mesh.dim, degree=2, units="mm/year" + ) + velocity_remap.read_timestep("physical", "velocity", 0, outputPath=str(directory)) + np.testing.assert_allclose(np.array(velocity_remap.array), expected_velocity) - with h5py.File(mesh_file, "r") as handle: + if uw.mpi.rank != 0: + return + with h5py.File(directory / "physical.mesh.00000.h5", "r") as handle: native_coordinates = handle["geometry/vertices"][:] physical_coordinates = handle["viz/geometry/vertices"][:] np.testing.assert_allclose(physical_coordinates, native_coordinates * 10.0) - assert handle["geometry/vertices"].attrs["units"] == "nondimensional" assert handle["viz/geometry/vertices"].attrs["units"] == "kilometer" with h5py.File(velocity_file, "r") as handle: - native = handle["fields/velocity"][:].reshape(-1, mesh.dim) - physical = handle["vertex_fields/velocity_velocity"][:].reshape(-1, mesh.dim) - np.testing.assert_allclose(native, [[2.0, 3.0]] * len(native)) - np.testing.assert_allclose(physical, [[10.0, 15.0]] * len(physical)) - np.testing.assert_allclose( - handle["vertex_fields/coordinates"][:].reshape(-1, mesh.dim), - native_coordinates * 10.0, - ) - assert handle["fields/velocity"].attrs["units"] == "nondimensional" - assert ( - handle["vertex_fields/velocity_velocity"].attrs["units"] - == "millimeter / year" - ) - assert handle["vertex_fields/coordinates"].attrs["units"] == "kilometer" - np.testing.assert_allclose( - handle["uw_checkpoint/velocity"][:].reshape(-1, mesh.dim), native - ) + physical_velocity = handle["fields/velocity"][:] + np.testing.assert_allclose(physical_velocity[:, 0], 10.0) + np.testing.assert_allclose(physical_velocity[:, 1], 15.0) + assert np.isclose(handle["fields/coordinates"][:].max(), 10.0) + assert handle["fields/velocity"].attrs["units"] == "millimeter / year" + assert handle["fields/coordinates"].attrs["units"] == "kilometer" + assert "visualization" not in handle + assert "vertex_fields" not in handle + assert "uw_checkpoint/topologies/uw_mesh/dms/velocity/vecs/velocity/velocity" in handle with h5py.File(pressure_file, "r") as handle: - native = handle["fields/pressure"][:].reshape(-1) - physical = handle["cell_fields/pressure_pressure"][:].reshape(-1) - np.testing.assert_allclose(native, 4.0) - np.testing.assert_allclose(physical, 8.0) - assert handle["cell_fields/pressure_pressure"].attrs["units"] == "megapascal" - - with h5py.File(distance_file, "r") as handle: - native = handle["fields/surf_unit_test_absdistance"][:].reshape(-1) - physical = handle[ - "vertex_fields/surf_unit_test_absdistance_surf_unit_test_absdistance" - ][:].reshape(-1) - np.testing.assert_allclose(physical, native * 10.0) - assert distance.units == uw.units("km").units - assert ( - handle[ - "vertex_fields/surf_unit_test_absdistance_surf_unit_test_absdistance" - ].attrs["units"] - == "kilometer" - ) - - xdmf = (directory / "physical.mesh.00000.xdmf").read_text() - assert "&MeshData;:/viz/geometry/vertices" in xdmf - assert 'Name="velocity"' in xdmf - assert 'Information Name="Units" Value="millimeter / year"' in xdmf - assert 'Information Name="Units" Value="kilometer"' in xdmf + np.testing.assert_allclose(handle["fields/pressure"][:], 8.0) + assert handle["fields/pressure"].attrs["units"] == "megapascal" + assert "cell_fields" not in handle + + text = (directory / "physical.mesh.00000.xdmf").read_text() + assert "&velocity_Data;:/fields/velocity" in text + assert "&pressure_Data;:/fields/pressure" in text + assert 'Information Name="Units" Value="millimeter / year"' in text + assert 'Information Name="Units" Value="megapascal"' in text + assert 'Information Name="Units" Value="kilometer"' in text @pytest.mark.level_1 @pytest.mark.tier_b -def test_dg1_xdmf_uses_native_interpolation_and_physical_output(tmp_path): - """DG1 interpolation stays native before its disconnected grid is scaled.""" +def test_dg1_fields_are_physical_at_element_corners(tmp_path): + """DG1 basis conversion and units are stored once under /fields.""" _set_reference_scales() mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) pressure = uw.discretisation.MeshVariable( "dg_pressure", mesh, 1, degree=1, continuous=False, units="MPa" ) - pressure.data[:, 0] = 4.0 - + pressure.array[:, 0, 0] = uw.quantity(8.0, "MPa") directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) - mesh.write_timestep( - "dg_physical", - index=0, - outputPath=str(directory), - meshVars=[pressure], - petsc_reload=True, - ) - - with uw.selective_ranks(0) as should_execute: - if not should_execute: - return - - field_file = directory / "dg_physical.mesh.dg_pressure.00000.h5" - with h5py.File(field_file, "r") as handle: - np.testing.assert_allclose(handle["fields/dg_pressure"][:], 4.0) - np.testing.assert_allclose(handle["dg1/values"][:], 8.0) - assert np.isclose(handle["dg1/vertices"][:].max(), 10.0) - assert handle["dg1/vertices"].attrs["units"] == "kilometer" - assert handle["dg1/values"].attrs["units"] == "megapascal" - - xdmf = (directory / "dg_physical.mesh.00000.xdmf").read_text() - assert 'Information Name="Units" Value="megapascal"' in xdmf - assert 'Information Name="Units" Value="kilometer"' in xdmf + mesh.write_timestep("dg", 0, outputPath=str(directory), meshVars=[pressure]) + if uw.mpi.rank == 0: + with h5py.File(directory / "dg.mesh.dg_pressure.00000.h5", "r") as handle: + np.testing.assert_allclose(handle["fields/dg_pressure"][:], 8.0) + assert np.isclose(handle["fields/coordinates"][:].max(), 10.0) + assert handle["fields/dg_pressure"].attrs["units"] == "megapascal" + assert handle["fields/coordinates"].attrs["units"] == "kilometer" + assert "dg1" not in handle diff --git a/tests/test_0010_snapshot_disk_format.py b/tests/test_0010_snapshot_disk_format.py index 0b016e246..773ef74f8 100644 --- a/tests/test_0010_snapshot_disk_format.py +++ b/tests/test_0010_snapshot_disk_format.py @@ -221,8 +221,9 @@ def test_write_snapshot_produces_wrapper_and_bulk_dir(tmp_path): if filename.endswith(f".{variable_name}.00000.h5") ) with h5py.File(os.path.join(bulk, variable_file), "r") as h5: - assert variable_name in h5["topologies"]["uw_mesh"]["dms"] - assert variable_name in h5["uw_checkpoint"] + dms = h5["uw_checkpoint/topologies/uw_mesh/dms"] + assert variable_name in dms + assert f"{variable_name}/vecs/{variable_name}/{variable_name}" in dms def test_snapshot_bulk_filenames_do_not_expand_loaded_mesh_name(tmp_path): From eb689b4273c2f7ae5a4968a7a8ff5905c9188de5 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 18 Sep 2026 17:58:17 +1000 Subject: [PATCH 5/9] Preserve native DG1 fields for timestep reload Keep dimensional DG1 interpolation coordinates and values under /fields so read_timestep can recover solver data exactly. Store the disconnected-corner basis conversion separately under /visualization for XDMF, where it preserves element geometry and discontinuous traces.\n\nRemove the DG1 read_timestep rejection and add serial and two-rank MPI coverage for coordinate reload, PETSc checkpoint reload, dimensional conversion, and ParaView rendering. Update the checkpoint documentation to explain why DG1 requires separate reload and visualization representations. --- .../checkpoint-output-and-reload-methods.md | 22 +++---- .../subsystems/checkpointing-system.md | 18 +++--- .../discretisation/discretisation_mesh.py | 64 ++++++++++++------- .../discretisation_mesh_variables.py | 12 ---- tests/test_0005_xdmf_dg1.py | 20 +++--- tests/test_0005_xdmf_physical_units.py | 19 +++++- 6 files changed, 89 insertions(+), 66 deletions(-) diff --git a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md index 1acae8383..3c45fcbd3 100644 --- a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md +++ b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md @@ -31,15 +31,14 @@ The XDMF storage choice follows the finite-element layout: | Continuous P1 | Direct `/fields` node values | | Continuous P2 triangles | Direct `/fields` values with `Triangle_6` connectivity | | DG0 | Direct `/fields` cell values | -| DG1 triangles/tetrahedra | Exact element-local values at disconnected corners | +| DG1 triangles/tetrahedra | Native `/fields` for reload plus exact disconnected-corner `/visualization` data | | Continuous P3+ or unsupported P2 | Compact P1 dataset under `/visualization` | | DG2+ or unsupported DG1 | Compact DG0 dataset under `/visualization` | -Direct DG1 output is an exact basis conversion for analysis and visualization. -Its disconnected corner coordinates differ from UW3's interior DG1 -interpolation nodes, so use `/uw_checkpoint` with `read_checkpoint()` for exact -DG1 restart. `read_timestep()` performs nearest-neighbour remapping and does not -invert that basis conversion. +DG1 retains native interpolation coordinates and values under `/fields`, so +`read_timestep()` can reload the solver field. XDMF reads an additional exact +basis conversion at disconnected element corners from `/visualization` because +the native interior DG1 points do not describe the full element geometry. ### Visualisation And Remap @@ -67,11 +66,12 @@ output.mesh.00000.xdmf ``` The field files contain `/fields/` and `/fields/coordinates`. P1, P2 -triangles, DG0, and DG1 simplices are visualized directly from those datasets. -Continuous P3+ fields use one compact P1 visualization reduction, and DG2+ -fields use DG0. Reloading with `read_timestep()` compares target coordinates to -the dimensional source coordinates and converts the saved values back to the -active model's nondimensional solver frame. +triangles, and DG0 are visualized directly from those datasets. DG1 uses its +additional exact disconnected-corner representation. Continuous P3+ fields use +one compact P1 visualization reduction, and DG2+ fields use DG0. Reloading with +`read_timestep()` compares target coordinates to the dimensional source +coordinates and converts the saved values back to the active model's +nondimensional solver frame. ### Unified Visualisation And PETSc Reload diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index 3f295cacc..8b345e8e9 100644 --- a/docs/developer/subsystems/checkpointing-system.md +++ b/docs/developer/subsystems/checkpointing-system.md @@ -43,10 +43,11 @@ analysis arrays remain dimensional. XDMF reads P1 and DG0 values directly from `/fields`. Continuous P2 fields on triangles use XDMF `Triangle_6` connectivity, including the three edge nodes, -so no P1 projection is stored. DG1 fields on full-dimensional triangular and -tetrahedral meshes use disconnected element corners. The element polynomial is -evaluated at each cell's corners, preserving jumps without averaging traces -across shared edges or faces. +so no P1 projection is stored. DG1 keeps native interpolation coordinates and +values under `/fields`, allowing `read_timestep()` to recover the solver field. +For XDMF, the same element polynomial is evaluated at disconnected cell corners +under `/visualization`, preserving jumps without averaging traces across shared +edges or faces. XDMF cannot represent every UW3 finite-element layout directly. Continuous P3+ fields and unsupported P2 layouts receive one compact P1 dataset under @@ -54,11 +55,10 @@ fields and unsupported P2 layouts receive one compact P1 dataset under DG0 dataset. Their exact dimensional values remain under `/fields`. Integration point fields are not supported by this writer. -Direct DG1 output stores an exact basis conversion at disconnected element -corners. This preserves the field for analysis and visualization, but those -corner coordinates differ from UW3's interior DG1 interpolation nodes. Use the -optional `/uw_checkpoint` payload and `read_checkpoint()` for exact DG1 solver -restart; nearest-neighbour `read_timestep()` is not an inverse basis conversion. +DG1 therefore has two representations because their coordinate sets serve +different purposes: native `/fields` for analysis and coordinate reload, and an +exact disconnected-corner `/visualization` representation for XDMF. Add +`/uw_checkpoint` when PETSc-native restart is also required. ```python mesh.write_timestep( diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 96f6b37b1..e9cb31763 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -4839,9 +4839,11 @@ def write_timestep( The optional payloads are controlled explicitly: - ``create_xdmf=True`` writes a companion XDMF file that reads the - dimensional ``/fields`` datasets directly for P1, P2 triangles, - DG0, and DG1 simplices. Continuous P3+ fields receive one compact P1 - visualization dataset; discontinuous DG2+ fields receive DG0. + dimensional ``/fields`` datasets directly for P1, P2 triangles, and + DG0. DG1 retains native ``/fields`` values for reload and adds an + exact disconnected-corner visualization. Continuous P3+ fields + receive one compact P1 visualization dataset; discontinuous DG2+ + fields receive DG0. - ``petsc_reload=True`` additionally writes native nondimensional PETSc DMPlex section/local-vector data under ``/uw_checkpoint``. Load that optional payload with ``MeshVariable.read_checkpoint()`` @@ -9635,10 +9637,13 @@ def remesh(self, metric_field, verbose=False): def _write_xdmf_field(mesh, var, var_h5_path): """Replace native remap arrays with dimensional field output for XDMF. - P1, P2 triangles, DG0, and DG1 simplices are represented directly under - ``/fields``. Unsupported higher-order layouts retain their exact physical - values under ``/fields`` and receive one compact visualization reduction: - continuous fields use P1 and discontinuous fields use DG0. + P1, P2 triangles, and DG0 are represented directly under ``/fields``. + DG1 keeps its native interpolation values under ``/fields`` for exact + coordinate reload and receives an exact disconnected-corner representation + under ``/visualization``. Unsupported higher-order layouts retain their + exact physical values under ``/fields`` and receive one compact + visualization reduction: continuous fields use P1 and discontinuous fields + use DG0. Parameters ---------- @@ -9693,18 +9698,18 @@ def _write_xdmf_field(mesh, var, var_h5_path): comm=PETSc.COMM_WORLD, ) + write_field_to_viewer(var, viewer, "/fields", var.clean_name) + write_field_coordinates_to_viewer(var, viewer, "/fields") if direct_dg1: _write_dg1_to_viewer( var, viewer, - group="/fields", + group="/visualization", coordinate_name="coordinates", value_name=var.clean_name, repack_tensors=False, ) else: - write_field_to_viewer(var, viewer, "/fields", var.clean_name) - write_field_coordinates_to_viewer(var, viewer, "/fields") if direct_p2: write_p2_triangle_topology_to_viewer(var, viewer, group="/fields") elif needs_projection: @@ -9730,18 +9735,26 @@ def _write_xdmf_field(mesh, var, var_h5_path): field.attrs["storage_frame"] = "physical" field.attrs["degree"] = var.degree field.attrs["continuous"] = var.continuous - field.attrs["representation"] = ( - "basis_conversion" if direct_dg1 else "exact" - ) + field.attrs["representation"] = "exact" coordinates = handle["fields/coordinates"] coordinates.attrs["units"] = coordinate_units or "dimensionless" coordinates.attrs["storage_frame"] = "physical" - if needs_projection: + if direct_dg1 or needs_projection: projected = handle[f"visualization/{var.clean_name}"] projected.attrs["units"] = field_units or "dimensionless" projected.attrs["source_degree"] = var.degree - projected.attrs["visualization_degree"] = 1 if var.continuous else 0 - projected.attrs["representation"] = "projection" + projected.attrs["visualization_degree"] = ( + 1 if var.continuous or direct_dg1 else 0 + ) + projected.attrs["representation"] = ( + "basis_conversion" if direct_dg1 else "projection" + ) + if direct_dg1: + visual_coordinates = handle["visualization/coordinates"] + visual_coordinates.attrs["units"] = ( + coordinate_units or "dimensionless" + ) + visual_coordinates.attrs["storage_frame"] = "physical" uw.mpi.barrier() @@ -10085,12 +10098,15 @@ def attribute_kind(var, components): special_grids = "" for var in special_vars: var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" + storage_group = "fields" if direct_p2(var) else "visualization" with h5py.File(var_filename, "r") as f: - cells_shape = f["fields/cells"].shape - points_shape = f["fields/coordinates"].shape - values_shape = f[f"fields/{var.clean_name}"].shape - special_geometry_units = f["fields/coordinates"].attrs.get("units") - field_units = f[f"fields/{var.clean_name}"].attrs.get("units") + cells_shape = f[f"{storage_group}/cells"].shape + points_shape = f[f"{storage_group}/coordinates"].shape + values_shape = f[f"{storage_group}/{var.clean_name}"].shape + special_geometry_units = f[f"{storage_group}/coordinates"].attrs.get( + "units" + ) + field_units = f[f"{storage_group}/{var.clean_name}"].attrs.get("units") special_topology = "Triangle_6" if direct_p2(var) else topology_type components = values_shape[1] if len(values_shape) == 2 else 1 kind = attribute_kind(var, components) @@ -10101,18 +10117,18 @@ def attribute_kind(var, components): - &{var.clean_name}_Data;:/fields/cells + &{var.clean_name}_Data;:/{storage_group}/cells - &{var.clean_name}_Data;:/fields/coordinates + &{var.clean_name}_Data;:/{storage_group}/coordinates {units_information(special_geometry_units, " ")} - &{var.clean_name}_Data;:/fields/{var.clean_name} + &{var.clean_name}_Data;:/{storage_group}/{var.clean_name} {units_information(field_units, " ")} """ diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index ebedf35f5..40b69f51f 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -1244,7 +1244,6 @@ def read_timestep( varsymbol=r"\cal{S}", ) - field_representation = None if uw.mpi.rank == 0: if verbose: print( @@ -1260,9 +1259,6 @@ def read_timestep( with h5py.File(data_file, "r") as h5f: coordinate_dataset = h5f["fields"]["coordinates"] field_dataset = h5f["fields"][data_name] - field_representation = field_dataset.attrs.get("representation") - if isinstance(field_representation, bytes): - field_representation = field_representation.decode() X_src = coordinate_dataset[()].reshape(-1, dim) D_src = field_dataset[()].reshape(-1, n_components) storage_frame = field_dataset.attrs.get( @@ -1291,14 +1287,6 @@ def read_timestep( X_src = np.empty((0, dim), dtype=np.float64) D_src = np.empty((0, n_components), dtype=np.float64) - field_representation = uw.mpi.comm.bcast(field_representation, root=0) - if field_representation == "basis_conversion": - raise RuntimeError( - "read_timestep cannot invert the element-local DG1 corner " - "basis conversion. Write with petsc_reload=True and use " - "read_checkpoint() for an exact DG1 solver restart." - ) - src_size_before = max(source_swarm.dm.getLocalSize(), 0) source_swarm.add_particles_with_global_coordinates(X_src, migrate=False) source_swarm._invalidate_canonical_data() diff --git a/tests/test_0005_xdmf_dg1.py b/tests/test_0005_xdmf_dg1.py index b99f03a91..69e7ed219 100644 --- a/tests/test_0005_xdmf_dg1.py +++ b/tests/test_0005_xdmf_dg1.py @@ -68,20 +68,23 @@ def test_dg1_simplex_output(tmp_path, dim): ) remapped = uw.discretisation.MeshVariable("remapped", mesh, 1, degree=1, continuous=False) - with pytest.raises(RuntimeError, match="DG1 corner basis conversion"): - remapped.read_timestep("fields", "dg_scalar", 0, outputPath=str(directory)) + remapped.read_timestep("fields", "dg_scalar", 0, outputPath=str(directory)) + np.testing.assert_allclose(remapped.array, scalar.array, rtol=1.0e-12, atol=1.0e-12) if uw.mpi.rank != 0: return with h5py.File(directory / "fields.mesh.dg_scalar.00000.h5", "r") as handle: - points = handle["fields/coordinates"][:] - cells = handle["fields/cells"][:] - values = handle["fields/dg_scalar"][:].reshape(-1) - assert handle["fields/dg_scalar"].attrs["representation"] == "basis_conversion" + native_points = handle["fields/coordinates"][:] + native_values = handle["fields/dg_scalar"][:].reshape(-1) + points = handle["visualization/coordinates"][:] + cells = handle["visualization/cells"][:] + values = handle["visualization/dg_scalar"][:].reshape(-1) + assert handle["fields/dg_scalar"].attrs["representation"] == "exact" + assert handle["visualization/dg_scalar"].attrs["representation"] == "basis_conversion" + assert len(native_points) == len(native_values) assert len(points) == len(cells) * (dim + 1) assert len(np.unique(cells)) == len(points) assert "dg1" not in handle - assert "visualization" not in handle with h5py.File(directory / "fields.mesh.00000.h5", "r") as handle: assert len(cells) == len(handle["viz/topology/cells"]) @@ -101,7 +104,7 @@ def test_dg1_simplex_output(tmp_path, dim): assert np.max(high - low) >= 1 with h5py.File(directory / "fields.mesh.dg_tensor.00000.h5", "r") as handle: - tensor_values = handle["fields/dg_tensor"][:] + tensor_values = handle["visualization/dg_tensor"][:] assert tensor_values.shape == (len(points), dim * dim) np.testing.assert_allclose(tensor_values[:, 0], expected) np.testing.assert_allclose(tensor_values[:, 1], 3 + points[:, 0]) @@ -116,6 +119,7 @@ def test_dg1_simplex_output(tmp_path, dim): attribute = grid.find("Attribute") assert attribute.get("Name") == name assert attribute.get("Center") == "Node" + assert f"/visualization/{name}" in attribute.find("DataItem").text @pytest.mark.level_1 diff --git a/tests/test_0005_xdmf_physical_units.py b/tests/test_0005_xdmf_physical_units.py index 7e84dae1f..48e8ac4c3 100644 --- a/tests/test_0005_xdmf_physical_units.py +++ b/tests/test_0005_xdmf_physical_units.py @@ -98,7 +98,7 @@ def test_fields_are_dimensional_and_checkpoint_is_native(tmp_path): @pytest.mark.level_1 @pytest.mark.tier_b def test_dg1_fields_are_physical_at_element_corners(tmp_path): - """DG1 basis conversion and units are stored once under /fields.""" + """DG1 native reload data and corner visualization use physical units.""" _set_reference_scales() mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) pressure = uw.discretisation.MeshVariable( @@ -107,10 +107,25 @@ def test_dg1_fields_are_physical_at_element_corners(tmp_path): pressure.array[:, 0, 0] = uw.quantity(8.0, "MPa") directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) mesh.write_timestep("dg", 0, outputPath=str(directory), meshVars=[pressure]) + + remapped = uw.discretisation.MeshVariable( + "dg_pressure_remapped", + mesh, + 1, + degree=1, + continuous=False, + units="MPa", + ) + remapped.read_timestep("dg", "dg_pressure", 0, outputPath=str(directory)) + np.testing.assert_allclose(np.array(remapped.array), np.array(pressure.array)) + if uw.mpi.rank == 0: with h5py.File(directory / "dg.mesh.dg_pressure.00000.h5", "r") as handle: np.testing.assert_allclose(handle["fields/dg_pressure"][:], 8.0) - assert np.isclose(handle["fields/coordinates"][:].max(), 10.0) + np.testing.assert_allclose(handle["visualization/dg_pressure"][:], 8.0) + assert np.isclose(handle["visualization/coordinates"][:].max(), 10.0) assert handle["fields/dg_pressure"].attrs["units"] == "megapascal" assert handle["fields/coordinates"].attrs["units"] == "kilometer" + assert handle["visualization/dg_pressure"].attrs["units"] == "megapascal" + assert handle["visualization/coordinates"].attrs["units"] == "kilometer" assert "dg1" not in handle From 917b6c03d96fa01311353f26e6a359c3a35c3767 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 18 Sep 2026 18:36:54 +1000 Subject: [PATCH 6/9] Add exact P2 tetra XDMF and compact restart output Write continuous P2 tetrahedra as VTK-ordered Tetrahedron_10 grids so ParaView can render the native quadratic field without a P1 reduction. Keep field coordinates and values in their native order so coordinate-based reload remains independent of visualization connectivity. For create_xdmf=False with petsc_reload=True, write only the DMPlex checkpoint payload instead of first duplicating values under /fields. Store one owned global vector and reload it directly for the same layout or through PETSc's migration SF for reconstructed layouts and different MPI decompositions. Add serial and MPI coverage for quadratic tetrahedra, checkpoint-only file structure, same-layout reload, and migrated reload. --- .../checkpoint-output-and-reload-methods.md | 28 +++---- .../subsystems/checkpointing-system.md | 23 +++--- .../discretisation/discretisation_mesh.py | 79 +++++++++++-------- .../discretisation_mesh_variables.py | 32 ++++++-- src/underworld3/function/field_projection.py | 38 +++++---- tests/test_0005_xdmf_compat.py | 76 ++++++++++++++++++ 6 files changed, 196 insertions(+), 80 deletions(-) diff --git a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md index 3c45fcbd3..9948de09c 100644 --- a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md +++ b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md @@ -11,8 +11,8 @@ but new code should use `write_timestep(..., petsc_reload=True)`. ## Standard API -`write_timestep()` always writes the mesh file and one HDF5 file per mesh -variable. With XDMF enabled, mesh-variable files contain dimensional +`write_timestep()` writes the mesh file and one HDF5 file per mesh variable. +With XDMF enabled, mesh-variable files contain dimensional coordinates and values under `/fields`. These are the authoritative arrays for analysis and the source data used by `MeshVariable.read_timestep()` for coordinate remapping. @@ -22,17 +22,17 @@ The two optional payloads are selected with explicit flags: | Flag | Output payload | Reader/use case | | --- | --- | --- | | `create_xdmf=True` | Dimensional `/fields`, compact high-order reductions when needed, and a companion `.xdmf` file | Analysis and ParaView/XDMF visualisation | -| `petsc_reload=True` | Native PETSc DMPlex section/vector data under `/uw_checkpoint` | `MeshVariable.read_checkpoint()` exact reload | +| `petsc_reload=True` | Native PETSc DMPlex section/global-vector data under `/uw_checkpoint` | `MeshVariable.read_checkpoint()` exact reload | The XDMF storage choice follows the finite-element layout: | Field layout | XDMF representation | | --- | --- | | Continuous P1 | Direct `/fields` node values | -| Continuous P2 triangles | Direct `/fields` values with `Triangle_6` connectivity | +| Continuous P2 triangles/tetrahedra | Direct `/fields` values with `Triangle_6`/`Tetrahedron_10` connectivity | | DG0 | Direct `/fields` cell values | | DG1 triangles/tetrahedra | Native `/fields` for reload plus exact disconnected-corner `/visualization` data | -| Continuous P3+ or unsupported P2 | Compact P1 dataset under `/visualization` | +| Continuous P3+ or unsupported P2 layout | Compact P1 dataset under `/visualization` | | DG2+ or unsupported DG1 | Compact DG0 dataset under `/visualization` | DG1 retains native interpolation coordinates and values under `/fields`, so @@ -66,11 +66,11 @@ output.mesh.00000.xdmf ``` The field files contain `/fields/` and `/fields/coordinates`. P1, P2 -triangles, and DG0 are visualized directly from those datasets. DG1 uses its -additional exact disconnected-corner representation. Continuous P3+ fields use -one compact P1 visualization reduction, and DG2+ fields use DG0. Reloading with -`read_timestep()` compares target coordinates to the dimensional source -coordinates and converts the saved values back to the active model's +triangles/tetrahedra, and DG0 are visualized directly from those datasets. DG1 +uses its additional exact disconnected-corner representation. Continuous P3+ +fields use one compact P1 visualization reduction, and DG2+ fields use DG0. +Reloading with `read_timestep()` compares target coordinates to the dimensional +source coordinates and converts the saved values back to the active model's nondimensional solver frame. ### Unified Visualisation And PETSc Reload @@ -113,8 +113,8 @@ mesh.write_timestep( ) ``` -This uses the established native field writer and does not create a companion -`.xdmf` file. The optional checkpoint data are added under `/uw_checkpoint`. +This writes the native checkpoint directly and does not create a companion +`.xdmf` file or duplicate `/fields` coordinate/value arrays. Typical PETSc-reload-only files still use the timestep naming convention: @@ -124,8 +124,8 @@ restart.mesh.Velocity.00000.h5 restart.mesh.Pressure.00000.h5 ``` -The variable files contain native `/fields` datasets and PETSc reload metadata -under `/uw_checkpoint/topologies/uw_mesh/dms//`. +Each variable file contains only PETSc reload metadata and one native global +vector under `/uw_checkpoint/topologies/uw_mesh/dms//`. ### Advantages diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index 8b345e8e9..c3e25db31 100644 --- a/docs/developer/subsystems/checkpointing-system.md +++ b/docs/developer/subsystems/checkpointing-system.md @@ -28,7 +28,7 @@ Optional payloads are controlled by explicit flags: | Flag | Payload | Reader / use | | --- | --- | --- | | `create_xdmf=True` | XDMF-compatible visualisation datasets and a companion `.xdmf` file | ParaView and other XDMF tools | -| `petsc_reload=True` | PETSc DMPlex section/vector metadata | `MeshVariable.read_checkpoint()` | +| `petsc_reload=True` | PETSc DMPlex section/global-vector metadata | `MeshVariable.read_checkpoint()` | When nondimensional scaling is active, `/fields` is converted during the write to the mesh and variable units declared in the model. HDF5 attributes and XDMF @@ -42,12 +42,12 @@ analysis arrays remain dimensional. ### Visualisation and Coordinate Remap XDMF reads P1 and DG0 values directly from `/fields`. Continuous P2 fields on -triangles use XDMF `Triangle_6` connectivity, including the three edge nodes, -so no P1 projection is stored. DG1 keeps native interpolation coordinates and -values under `/fields`, allowing `read_timestep()` to recover the solver field. -For XDMF, the same element polynomial is evaluated at disconnected cell corners -under `/visualization`, preserving jumps without averaging traces across shared -edges or faces. +triangles and tetrahedra use XDMF `Triangle_6` and `Tetrahedron_10` +connectivity, including their edge nodes, so no P1 projection is stored. DG1 +keeps native interpolation coordinates and values under `/fields`, allowing +`read_timestep()` to recover the solver field. For XDMF, the same element +polynomial is evaluated at disconnected cell corners under `/visualization`, +preserving jumps without averaging traces across shared edges or faces. XDMF cannot represent every UW3 finite-element layout directly. Continuous P3+ fields and unsupported P2 layouts receive one compact P1 dataset under @@ -125,10 +125,11 @@ output/restart.mesh.velocity.00100.h5 output/restart.mesh.pressure.00100.h5 ``` -The variable files contain PETSc reload metadata and native values under -`/uw_checkpoint/topologies/uw_mesh/dms//`. `read_checkpoint()` uses -PETSc DMPlex topology, section, vector, and `PetscSF` metadata. It does not use -the dimensional `/fields` values or KDTree remapping. +The variable files contain PETSc reload metadata and one native global vector +under `/uw_checkpoint/topologies/uw_mesh/dms//`. `read_checkpoint()` +uses PETSc DMPlex topology, section, vector, and `PetscSF` metadata. It does not +use dimensional `/fields` values or KDTree remapping. Restart-only output does +not write `/fields`, so the native values are stored only once. ### Unified Visualisation and PETSc Reload diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index e9cb31763..8cb211caf 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -4832,28 +4832,27 @@ def write_timestep( - one mesh HDF5 file, shared across timesteps unless ``meshUpdates=True`` - one HDF5 file per mesh variable - - dimensional coordinate/value datasets under ``/fields`` for - analysis and coordinate-based reload with - ``MeshVariable.read_timestep()`` - - The optional payloads are controlled explicitly: + The variable-file payloads are controlled explicitly: - ``create_xdmf=True`` writes a companion XDMF file that reads the - dimensional ``/fields`` datasets directly for P1, P2 triangles, and - DG0. DG1 retains native ``/fields`` values for reload and adds an - exact disconnected-corner visualization. Continuous P3+ fields - receive one compact P1 visualization dataset; discontinuous DG2+ - fields receive DG0. - - ``petsc_reload=True`` additionally writes native nondimensional - PETSc DMPlex section/local-vector data under ``/uw_checkpoint``. + dimensional ``/fields`` datasets directly for P1, P2 triangles and + tetrahedra, and DG0. DG1 retains native ``/fields`` values for reload + and adds an exact disconnected-corner visualization. Continuous P3+ + fields receive one compact P1 visualization dataset; discontinuous + DG2+ fields receive DG0. + - ``petsc_reload=True`` writes native nondimensional PETSc DMPlex + section/global-vector data under ``/uw_checkpoint``. Load that optional payload with ``MeshVariable.read_checkpoint()`` - for an exact solver restart. + for an exact solver restart. When XDMF is disabled, this is the only + variable payload; no duplicate ``/fields`` values are written. + - With both flags disabled, the established low-level native + ``/fields`` output is retained for compatibility. Common choices are: - visualisation/remap only: ``create_xdmf=True, petsc_reload=False`` - - PETSc-native reload only: + - PETSc-native reload only (no duplicate ``/fields`` values): ``create_xdmf=False, petsc_reload=True`` - unified visualisation/remap and PETSc reload: ``create_xdmf=True, petsc_reload=True`` @@ -4883,9 +4882,10 @@ def write_timestep( Write ParaView/XDMF-compatible dimensional datasets and the companion XDMF file. DG1 on full-dimensional triangles/tetrahedra uses independent vertices per cell, preserving jumps without - smoothing. Higher-order fields use compact P1 or DG0 visualization - reductions while their exact dimensional values remain in - ``/fields``. + smoothing. Continuous P2 triangles and tetrahedra use their exact + quadratic connectivity. Higher-order fields use compact P1 or DG0 + visualization reductions while their exact dimensional values + remain in ``/fields``. petsc_reload Write PETSc DMPlex section/vector metadata for reload with ``MeshVariable.read_checkpoint()``. @@ -4933,11 +4933,19 @@ def write_timestep( if meshVars is not None: for var in meshVars: save_location = output_base_name + f".mesh.{var.clean_name}.{index:05}.h5" - var.write(save_location) - if petsc_reload: - self._write_petsc_reload_file(save_location, [var], mode="a") if create_xdmf: + var.write(save_location) + if petsc_reload: + self._write_petsc_reload_file(save_location, [var], mode="a") _write_xdmf_field(self, var, save_location) + elif petsc_reload: + # Exact-restart-only output needs the PETSc section/vector + # payload, not a second native /fields copy of the values. + self._write_petsc_reload_file(save_location, [var], mode="w") + else: + # Preserve the established low-level native output when no + # optional visualization or restart payload was requested. + var.write(save_location) if swarmVars is not None: for svar in swarmVars: @@ -5024,22 +5032,26 @@ def _write_petsc_reload_variable(self, viewer, var): if var._lvec is None: var._set_vec(available=True) + var._sync_lvec_to_gvec() iset, subdm = self.dm.createSubDM(var.field_id) subdm.setName(var.clean_name) old_lvec_name = var._lvec.getName() + old_gvec_name = var._gvec.getName() try: var._lvec.setName(var.clean_name) + var._gvec.setName(var.clean_name) self.dm.sectionView(viewer, subdm) - self.dm.localVectorView(viewer, subdm, var._lvec) + self.dm.globalVectorView(viewer, subdm, var._gvec) finally: var._lvec.setName(old_lvec_name) + var._gvec.setName(old_gvec_name) iset.destroy() subdm.destroy() def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): - """Write compact DMPlex reload metadata and native local vectors.""" + """Write compact DMPlex reload metadata and native global vectors.""" old_dm_name = self.dm.getName() self.dm.setName("uw_mesh") @@ -5052,8 +5064,8 @@ def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): try: # PETSc needs the complete source section to construct the # migration SF when the checkpoint is read with a different MPI - # ownership ordering. This is metadata only; field values are - # still stored once in each variable's local-vector payload. + # ownership ordering. This is metadata only; owned field values are + # stored once in each variable's global-vector payload. self.dm.sectionView(viewer, self.dm) for var in variables: @@ -9637,7 +9649,7 @@ def remesh(self, metric_field, verbose=False): def _write_xdmf_field(mesh, var, var_h5_path): """Replace native remap arrays with dimensional field output for XDMF. - P1, P2 triangles, and DG0 are represented directly under ``/fields``. + P1, P2 triangles/tetrahedra, and DG0 are represented directly under ``/fields``. DG1 keeps its native interpolation values under ``/fields`` for exact coordinate reload and receives an exact disconnected-corner representation under ``/visualization``. Unsupported higher-order layouts retain their @@ -9663,7 +9675,7 @@ def _write_xdmf_field(mesh, var, var_h5_path): _write_dg1_to_viewer, write_field_coordinates_to_viewer, write_field_to_viewer, - write_p2_triangle_topology_to_viewer, + write_p2_simplex_topology_to_viewer, write_projected_field_to_viewer, ) @@ -9671,8 +9683,8 @@ def _write_xdmf_field(mesh, var, var_h5_path): var.continuous and var.degree == 2 and mesh.isSimplex - and mesh.dim == 2 - and mesh.cdim == 2 + and mesh.dim in (2, 3) + and mesh.cdim == mesh.dim ) direct_dg1 = ( not var.continuous @@ -9711,7 +9723,7 @@ def _write_xdmf_field(mesh, var, var_h5_path): ) else: if direct_p2: - write_p2_triangle_topology_to_viewer(var, viewer, group="/fields") + write_p2_simplex_topology_to_viewer(var, viewer, group="/fields") elif needs_projection: target_degree = 1 if var.continuous else 0 write_projected_field_to_viewer( @@ -9941,8 +9953,8 @@ def direct_p2(var): var.continuous and var.degree == 2 and var.mesh.isSimplex - and var.mesh.dim == 2 - and var.mesh.cdim == 2 + and var.mesh.dim in (2, 3) + and var.mesh.cdim == var.mesh.dim ) def direct_dg1(var): @@ -10107,7 +10119,10 @@ def attribute_kind(var, components): "units" ) field_units = f[f"{storage_group}/{var.clean_name}"].attrs.get("units") - special_topology = "Triangle_6" if direct_p2(var) else topology_type + if direct_p2(var): + special_topology = "Triangle_6" if var.mesh.dim == 2 else "Tetrahedron_10" + else: + special_topology = topology_type components = values_shape[1] if len(values_shape) == 2 else 1 kind = attribute_kind(var, components) dimensions = " ".join(str(value) for value in values_shape) diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 40b69f51f..f0b08eac5 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -1475,7 +1475,7 @@ def read_checkpoint( ): """Load this mesh variable from PETSc reload output. - By default, DMPlex section/local-vector data are restored through the + By default, DMPlex section/global-vector data are restored through the topology migration SF, so a mesh reconstructed from its checkpoint may have a different parallel DOF ordering. With ``same_layout=True``, the existing PETSc variable vector is loaded directly into the original @@ -1538,18 +1538,31 @@ def read_checkpoint( if grouped_checkpoint else "/uw_checkpoint" ) + # A DM-associated Vec ignores the viewer group and redirects + # HDF5 I/O through /fields. A plain Vec reads the requested + # /uw_checkpoint dataset when /fields is intentionally absent. + checkpoint_vec = PETSc.Vec().createMPI( + (self._gvec.getLocalSize(), PETSc.DECIDE), + comm=PETSc.COMM_WORLD, + ) + checkpoint_vec.setName(data_name) viewer.pushGroup(vector_group) - self._gvec.load(viewer) + checkpoint_vec.load(viewer) viewer.popGroup() + self._gvec.array[...] = checkpoint_vec.array_r subdm.globalToLocal(self._gvec, self._lvec, addv=False) + checkpoint_vec.destroy() else: - from underworld3.cython.petsc_discretisation import ( - petsc_dmplex_load_local_vector, + global_sf, local_sf = self.mesh.dm.sectionLoad( + viewer, sectiondm, self.mesh.sf ) - - loaded_lvec = petsc_dmplex_load_local_vector( - self.mesh.dm, viewer, sectiondm, self.mesh.sf, data_name + loaded_gvec = sectiondm.createGlobalVec() + loaded_gvec.setName(data_name) + self.mesh.dm.globalVectorLoad( + viewer, sectiondm, global_sf, loaded_gvec ) + loaded_lvec = sectiondm.createLocalVec() + sectiondm.globalToLocal(loaded_gvec, loaded_lvec, addv=False) source_section = sectiondm.getSection() target_section = subdm.getSection() @@ -1576,6 +1589,11 @@ def read_checkpoint( ) loaded_lvec.destroy() + loaded_gvec.destroy() + if global_sf is not None: + global_sf.destroy() + if local_sf is not None: + local_sf.destroy() self._sync_lvec_to_gvec() finally: self._lvec.setName(old_lvec_name) diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py index 7a1546a49..bfd5b8b83 100644 --- a/src/underworld3/function/field_projection.py +++ b/src/underworld3/function/field_projection.py @@ -404,17 +404,19 @@ def write_projected_field_to_viewer( _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) -def write_p2_triangle_topology_to_viewer(mesh_var, viewer, group="/fields"): - """Write VTK-ordered Triangle_6 connectivity for a continuous P2 field.""" +def write_p2_simplex_topology_to_viewer(mesh_var, viewer, group="/fields"): + """Write VTK-ordered quadratic triangle or tetrahedron connectivity.""" mesh = mesh_var.mesh if not ( mesh_var.continuous and mesh_var.degree == 2 and mesh.isSimplex - and mesh.dim == 2 - and mesh.cdim == 2 + and mesh.dim in (2, 3) + and mesh.cdim == mesh.dim ): - raise NotImplementedError("direct P2 XDMF currently requires a 2D triangle mesh") + raise NotImplementedError( + "direct P2 XDMF requires a full-dimensional triangle or tetrahedron mesh" + ) coordinate_dm = mesh._basis_coordinate_dm(2, True) local_section = coordinate_dm.getLocalSection() @@ -444,28 +446,32 @@ def write_p2_triangle_topology_to_viewer(mesh_var, viewer, group="/fields"): leaves = np.asarray(leaves) owned[leaves[(leaves >= cell_start) & (leaves < cell_end)] - cell_start] = False - p2_rows = mesh._cell_node_indices(2, True).reshape(-1, 6)[owned] - vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, 3)[owned] + corner_count = mesh.dim + 1 + node_count = 6 if mesh.dim == 2 else 10 + p2_rows = mesh._cell_node_indices(2, True).reshape(-1, node_count)[owned] + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, corner_count)[owned] p2_coordinates = mesh_var.coords_nd[p2_rows] corners = mesh._get_coords_for_basis(1, True)[vertex_rows] negative = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) < 0 - corners[negative] = corners[negative][:, [0, 2, 1]] + if mesh.dim == 2: + corners[negative] = corners[negative][:, [0, 2, 1]] + edge_pairs = ((0, 1), (1, 2), (2, 0)) + topology_name = "Triangle_6" + else: + corners[negative] = corners[negative][:, [0, 2, 1, 3]] + edge_pairs = ((0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)) + topology_name = "Tetrahedron_10" connectivity = np.empty_like(p2_rows, dtype=PETSc.IntType) for cell_index, (row, nodes, vertices) in enumerate( zip(p2_rows, p2_coordinates, corners, strict=True) ): targets = np.vstack( - ( - vertices, - 0.5 * (vertices[0] + vertices[1]), - 0.5 * (vertices[1] + vertices[2]), - 0.5 * (vertices[2] + vertices[0]), - ) + (vertices, *(0.5 * (vertices[a] + vertices[b]) for a, b in edge_pairs)) ) order = [np.argmin(np.linalg.norm(nodes - target, axis=1)) for target in targets] - if len(set(order)) != 6: - raise RuntimeError("could not map UW3 P2 nodes to Triangle_6 ordering") + if len(set(order)) != node_count: + raise RuntimeError(f"could not map UW3 P2 nodes to {topology_name} ordering") connectivity[cell_index] = local_to_global[row[order]] if np.any(connectivity < 0): raise RuntimeError("P2 XDMF connectivity contains an unmapped global node") diff --git a/tests/test_0005_xdmf_compat.py b/tests/test_0005_xdmf_compat.py index ce5b4d233..af8a5bbbb 100644 --- a/tests/test_0005_xdmf_compat.py +++ b/tests/test_0005_xdmf_compat.py @@ -92,6 +92,50 @@ def test_direct_p1_p2_and_dg0_use_fields_only(tmp_path): _assert_xdmf_references_exist(xdmf) +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_direct_p2_tetrahedron_uses_exact_tetrahedron_10(tmp_path): + """A 3D continuous P2 field uses exact VTK-ordered tetrahedral data.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), + maxCoords=(1.0, 1.0, 1.0), + cellSize=0.7, + regular=True, + qdegree=3, + ) + p2 = uw.discretisation.MeshVariable("p2_tet", mesh, 3, degree=2) + p2.array[:, 0, :] = p2.coords + mesh.write_timestep("tet", 0, outputPath=str(directory), meshVars=[p2]) + + restored = uw.discretisation.MeshVariable("restored", mesh, 3, degree=2) + restored.read_timestep("tet", "p2_tet", 0, outputPath=str(directory)) + np.testing.assert_allclose(restored.array, p2.array, atol=1.0e-12) + + if uw.mpi.rank == 0: + field_file = directory / "tet.mesh.p2_tet.00000.h5" + with h5py.File(field_file, "r") as handle: + coordinates = handle["fields/coordinates"][:] + values = handle["fields/p2_tet"][:] + cells = handle["fields/cells"][:] + assert cells.shape[1] == 10 + assert values.shape == coordinates.shape + np.testing.assert_allclose(values, coordinates, atol=1.0e-12) + points = coordinates[cells] + edge_pairs = ((0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)) + for local_node, (a, b) in enumerate(edge_pairs, start=4): + np.testing.assert_allclose( + points[:, local_node], 0.5 * (points[:, a] + points[:, b]) + ) + assert "visualization" not in handle + _assert_no_compatibility_copies(handle) + + xdmf = directory / "tet.mesh.00000.xdmf" + assert 'TopologyType="Tetrahedron_10"' in xdmf.read_text() + assert "&p2_tet_Data;:/fields/p2_tet" in xdmf.read_text() + _assert_xdmf_references_exist(xdmf) + + @pytest.mark.level_1 @pytest.mark.tier_b def test_higher_order_fields_keep_exact_data_and_one_compact_reduction(tmp_path): @@ -219,6 +263,38 @@ def test_create_xdmf_false_preserves_native_writer(tmp_path): assert "visualization" not in handle +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_petsc_reload_only_omits_duplicate_fields(tmp_path): + """Restart-only output stores the native PETSc payload exactly once.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + field = uw.discretisation.MeshVariable("field", mesh, 1, degree=2) + field.array[:, 0, 0] = field.coords[:, 0] + 2.0 * field.coords[:, 1] + expected = np.asarray(field.array).copy() + mesh.write_timestep( + "restart", + 0, + outputPath=str(directory), + meshVars=[field], + create_xdmf=False, + petsc_reload=True, + ) + + field.array[...] = 0.0 + field.read_checkpoint( + str(directory / "restart.mesh.field.00000.h5"), + data_name="field", + same_layout=True, + ) + np.testing.assert_allclose(field.array, expected) + + if uw.mpi.rank == 0: + assert not (directory / "restart.mesh.00000.xdmf").exists() + with h5py.File(directory / "restart.mesh.field.00000.h5", "r") as handle: + assert set(handle) == {"uw_checkpoint"} + + @pytest.mark.level_1 @pytest.mark.tier_b @pytest.mark.parametrize("dim", [2, 3]) From 36c01399e74724c78fef3c725927503c27de3342 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 18 Sep 2026 18:44:35 +1000 Subject: [PATCH 7/9] Rename native restart payload to restart/petsc Write PETSc-native DMPlex restart data under /restart/petsc so the HDF5 hierarchy states both the payload purpose and storage backend. Keep read compatibility with grouped and direct-vector files that use the former /uw_checkpoint path, and cover the compatibility path in serial and MPI tests. --- .../checkpoint-output-and-reload-methods.md | 9 +++-- .../subsystems/checkpointing-system.md | 9 +++-- .../discretisation/discretisation_mesh.py | 6 ++-- .../discretisation_mesh_variables.py | 33 +++++++++++++------ tests/test_0003_save_load.py | 6 ++-- tests/test_0005_xdmf_compat.py | 25 +++++++++++--- tests/test_0005_xdmf_physical_units.py | 2 +- tests/test_0010_snapshot_disk_format.py | 2 +- 8 files changed, 63 insertions(+), 29 deletions(-) diff --git a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md index 9948de09c..2a07b02fc 100644 --- a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md +++ b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md @@ -22,7 +22,7 @@ The two optional payloads are selected with explicit flags: | Flag | Output payload | Reader/use case | | --- | --- | --- | | `create_xdmf=True` | Dimensional `/fields`, compact high-order reductions when needed, and a companion `.xdmf` file | Analysis and ParaView/XDMF visualisation | -| `petsc_reload=True` | Native PETSc DMPlex section/global-vector data under `/uw_checkpoint` | `MeshVariable.read_checkpoint()` exact reload | +| `petsc_reload=True` | Native PETSc DMPlex section/global-vector data under `/restart/petsc` | `MeshVariable.read_checkpoint()` exact reload | The XDMF storage choice follows the finite-element layout: @@ -95,7 +95,7 @@ velocity.read_checkpoint( ``` With both flags enabled, the same variable file contains dimensional `/fields` -for analysis, XDMF, and coordinate remapping plus native `/uw_checkpoint` data +for analysis, XDMF, and coordinate remapping plus native `/restart/petsc` data for exact PETSc reload. ### PETSc Reload Without XDMF @@ -125,7 +125,10 @@ restart.mesh.Pressure.00000.h5 ``` Each variable file contains only PETSc reload metadata and one native global -vector under `/uw_checkpoint/topologies/uw_mesh/dms//`. +vector under `/restart/petsc/topologies/uw_mesh/dms//`. + +`read_checkpoint()` also recognizes the former `/uw_checkpoint` group so +existing restart files remain readable. ### Advantages diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index c3e25db31..49b0a6c40 100644 --- a/docs/developer/subsystems/checkpointing-system.md +++ b/docs/developer/subsystems/checkpointing-system.md @@ -36,7 +36,7 @@ to the mesh and variable units declared in the model. HDF5 attributes and XDMF physical values directly, without maintaining their own conversion table. Set `petsc_reload=True` only when an exact restart is needed. It adds the native -nondimensional PETSc payload under `/uw_checkpoint`; the visualization and +nondimensional PETSc payload under `/restart/petsc`; the visualization and analysis arrays remain dimensional. ### Visualisation and Coordinate Remap @@ -58,7 +58,7 @@ point fields are not supported by this writer. DG1 therefore has two representations because their coordinate sets serve different purposes: native `/fields` for analysis and coordinate reload, and an exact disconnected-corner `/visualization` representation for XDMF. Add -`/uw_checkpoint` when PETSc-native restart is also required. +`/restart/petsc` when PETSc-native restart is also required. ```python mesh.write_timestep( @@ -126,11 +126,14 @@ output/restart.mesh.pressure.00100.h5 ``` The variable files contain PETSc reload metadata and one native global vector -under `/uw_checkpoint/topologies/uw_mesh/dms//`. `read_checkpoint()` +under `/restart/petsc/topologies/uw_mesh/dms//`. `read_checkpoint()` uses PETSc DMPlex topology, section, vector, and `PetscSF` metadata. It does not use dimensional `/fields` values or KDTree remapping. Restart-only output does not write `/fields`, so the native values are stored only once. +The reader also accepts the former `/uw_checkpoint` group for compatibility +with existing restart files. + ### Unified Visualisation and PETSc Reload Set both flags when one output family should support ParaView, coordinate diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 8cb211caf..d44e12c17 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -4841,7 +4841,7 @@ def write_timestep( fields receive one compact P1 visualization dataset; discontinuous DG2+ fields receive DG0. - ``petsc_reload=True`` writes native nondimensional PETSc DMPlex - section/global-vector data under ``/uw_checkpoint``. + section/global-vector data under ``/restart/petsc``. Load that optional payload with ``MeshVariable.read_checkpoint()`` for an exact solver restart. When XDMF is disabled, this is the only variable payload; no duplicate ``/fields`` values are written. @@ -5060,7 +5060,7 @@ def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): checkpoint_file, mode, comm=PETSc.COMM_WORLD ) viewer.pushFormat(PETSc.Viewer.Format.HDF5_PETSC) - viewer.pushGroup("/uw_checkpoint") + viewer.pushGroup("/restart/petsc") try: # PETSc needs the complete source section to construct the # migration SF when the checkpoint is read with a different MPI @@ -9666,7 +9666,7 @@ def _write_xdmf_field(mesh, var, var_h5_path): by ``var.write()`` (so ``var._gvec`` is up-to-date). var_h5_path : str Path to the HDF5 file. Native restart data, when requested, has already - been written under ``/uw_checkpoint``. + been written under ``/restart/petsc``. """ import h5py import underworld3 as uw diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index f0b08eac5..42381421a 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -1499,13 +1499,26 @@ def read_checkpoint( if uw.mpi.rank == 0: with h5py.File(filename, "r") as checkpoint_h5: - grouped_checkpoint = "uw_checkpoint/topologies" in checkpoint_h5 - legacy_direct_vector = f"uw_checkpoint/{data_name}" in checkpoint_h5 + if "restart/petsc/topologies" in checkpoint_h5: + checkpoint_group = "/restart/petsc" + elif "uw_checkpoint/topologies" in checkpoint_h5: + checkpoint_group = "/uw_checkpoint" + else: + checkpoint_group = None + + if f"restart/petsc/{data_name}" in checkpoint_h5: + direct_vector_group = "/restart/petsc" + elif f"uw_checkpoint/{data_name}" in checkpoint_h5: + direct_vector_group = "/uw_checkpoint" + else: + direct_vector_group = None else: - grouped_checkpoint = None - legacy_direct_vector = None - grouped_checkpoint = uw.mpi.comm.bcast(grouped_checkpoint, root=0) - legacy_direct_vector = uw.mpi.comm.bcast(legacy_direct_vector, root=0) + checkpoint_group = None + direct_vector_group = None + checkpoint_group = uw.mpi.comm.bcast(checkpoint_group, root=0) + direct_vector_group = uw.mpi.comm.bcast(direct_vector_group, root=0) + grouped_checkpoint = checkpoint_group is not None + legacy_direct_vector = direct_vector_group is not None if same_layout and not (grouped_checkpoint or legacy_direct_vector): raise RuntimeError( @@ -1518,7 +1531,7 @@ def read_checkpoint( viewer = PETSc.ViewerHDF5().create(filename, "r", comm=PETSc.COMM_WORLD) viewer.pushFormat(PETSc.Viewer.Format.HDF5_PETSC) if grouped_checkpoint and not same_layout: - viewer.pushGroup("/uw_checkpoint") + viewer.pushGroup(checkpoint_group) old_mesh_name = self.mesh.dm.getName() old_lvec_name = self._lvec.getName() @@ -1533,14 +1546,14 @@ def read_checkpoint( if same_layout: vector_group = ( - f"/uw_checkpoint/topologies/uw_mesh/dms/{data_name}/" + f"{checkpoint_group}/topologies/uw_mesh/dms/{data_name}/" f"vecs/{data_name}" if grouped_checkpoint - else "/uw_checkpoint" + else direct_vector_group ) # A DM-associated Vec ignores the viewer group and redirects # HDF5 I/O through /fields. A plain Vec reads the requested - # /uw_checkpoint dataset when /fields is intentionally absent. + # restart dataset when /fields is intentionally absent. checkpoint_vec = PETSc.Vec().createMPI( (self._gvec.getLocalSize(), PETSc.DECIDE), comm=PETSc.COMM_WORLD, diff --git a/tests/test_0003_save_load.py b/tests/test_0003_save_load.py index 5b23ef668..f2f2c98fd 100644 --- a/tests/test_0003_save_load.py +++ b/tests/test_0003_save_load.py @@ -173,8 +173,8 @@ def test_timestep_with_petsc_reload_roundtrip(tmp_path): assert "fields/x" in h5f assert "fields/coordinates" in h5f assert "vertex_fields" not in h5f - assert "uw_checkpoint/topologies/uw_mesh/dms/x/section" in h5f - assert "uw_checkpoint/topologies/uw_mesh/dms/x/vecs/x" in h5f + assert "restart/petsc/topologies/uw_mesh/dms/x/section" in h5f + assert "restart/petsc/topologies/uw_mesh/dms/x/vecs/x" in h5f uw.mpi.barrier() mesh_reloaded = uw.discretisation.Mesh(f"{tmp_path}/unified.mesh.00000.h5") @@ -206,7 +206,7 @@ def test_timestep_with_petsc_reload_roundtrip(tmp_path): with h5py.File(checkpoint_var_file, "r") as h5f: assert "fields/x" in h5f assert "vertex_fields" not in h5f - assert "uw_checkpoint/topologies/uw_mesh/dms/x/section" in h5f + assert "restart/petsc/topologies/uw_mesh/dms/x/section" in h5f uw.mpi.barrier() mesh_reloaded = uw.discretisation.Mesh( diff --git a/tests/test_0005_xdmf_compat.py b/tests/test_0005_xdmf_compat.py index af8a5bbbb..2bdf32319 100644 --- a/tests/test_0005_xdmf_compat.py +++ b/tests/test_0005_xdmf_compat.py @@ -211,7 +211,7 @@ def test_compact_tensor_keeps_native_component_count(tmp_path): @pytest.mark.level_1 @pytest.mark.tier_b def test_timestep_and_optional_checkpoint_roundtrip(tmp_path): - """Dimensional fields support remap; /uw_checkpoint supports exact reload.""" + """Dimensional fields support remap; /restart/petsc supports exact reload.""" directory = _shared_path(tmp_path) mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.4, regular=True) source = uw.discretisation.MeshVariable("source", mesh, 1, degree=2) @@ -239,9 +239,9 @@ def test_timestep_and_optional_checkpoint_roundtrip(tmp_path): if uw.mpi.rank == 0: with h5py.File(field_file, "r") as handle: - assert set(handle) == {"fields", "uw_checkpoint"} - assert "uw_checkpoint/topologies/uw_mesh/dms/source/vecs/source/source" in handle - assert set(handle["uw_checkpoint/topologies/uw_mesh/dms"]) == { + assert set(handle) == {"fields", "restart"} + assert "restart/petsc/topologies/uw_mesh/dms/source/vecs/source/source" in handle + assert set(handle["restart/petsc/topologies/uw_mesh/dms"]) == { "source", "uw_mesh", } @@ -292,7 +292,22 @@ def test_petsc_reload_only_omits_duplicate_fields(tmp_path): if uw.mpi.rank == 0: assert not (directory / "restart.mesh.00000.xdmf").exists() with h5py.File(directory / "restart.mesh.field.00000.h5", "r") as handle: - assert set(handle) == {"uw_checkpoint"} + assert set(handle) == {"restart"} + assert set(handle["restart"]) == {"petsc"} + + # Existing files that used the former group name remain readable. + with h5py.File(directory / "restart.mesh.field.00000.h5", "a") as handle: + handle.move("restart/petsc", "uw_checkpoint") + del handle["restart"] + uw.mpi.barrier() + + field.array[...] = 0.0 + field.read_checkpoint( + str(directory / "restart.mesh.field.00000.h5"), + data_name="field", + same_layout=True, + ) + np.testing.assert_allclose(field.array, expected) @pytest.mark.level_1 diff --git a/tests/test_0005_xdmf_physical_units.py b/tests/test_0005_xdmf_physical_units.py index 48e8ac4c3..52b38e7a5 100644 --- a/tests/test_0005_xdmf_physical_units.py +++ b/tests/test_0005_xdmf_physical_units.py @@ -80,7 +80,7 @@ def test_fields_are_dimensional_and_checkpoint_is_native(tmp_path): assert handle["fields/coordinates"].attrs["units"] == "kilometer" assert "visualization" not in handle assert "vertex_fields" not in handle - assert "uw_checkpoint/topologies/uw_mesh/dms/velocity/vecs/velocity/velocity" in handle + assert "restart/petsc/topologies/uw_mesh/dms/velocity/vecs/velocity/velocity" in handle with h5py.File(pressure_file, "r") as handle: np.testing.assert_allclose(handle["fields/pressure"][:], 8.0) diff --git a/tests/test_0010_snapshot_disk_format.py b/tests/test_0010_snapshot_disk_format.py index 773ef74f8..518c13431 100644 --- a/tests/test_0010_snapshot_disk_format.py +++ b/tests/test_0010_snapshot_disk_format.py @@ -221,7 +221,7 @@ def test_write_snapshot_produces_wrapper_and_bulk_dir(tmp_path): if filename.endswith(f".{variable_name}.00000.h5") ) with h5py.File(os.path.join(bulk, variable_file), "r") as h5: - dms = h5["uw_checkpoint/topologies/uw_mesh/dms"] + dms = h5["restart/petsc/topologies/uw_mesh/dms"] assert variable_name in dms assert f"{variable_name}/vecs/{variable_name}/{variable_name}" in dms From 29067306bb335acc08685a89a3cb54bb4421a8bc Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 18 Sep 2026 18:56:12 +1000 Subject: [PATCH 8/9] Fix XDMF3 mesh references for ParaView Inline the base-grid topology and geometry HDF5 DataItems in generated XDMF instead of relying on XML Reference nodes that ParaView 6's XDMF3 reader rejects. Keep the shared mesh DataItems for compatibility while making both the XDMF3 and legacy readers load exact P2 and DG1 output. Add a regression assertion that generated compact XDMF contains no XML-reference DataItems. --- src/underworld3/discretisation/discretisation_mesh.py | 11 +++++++---- tests/test_0005_xdmf_compat.py | 4 ++++ 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index d44e12c17..2751838c6 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -9994,13 +9994,16 @@ def direct_dg1(var): - - /Xdmf/Domain/DataItem[@Name="cells"] + + &MeshData;:/{topoPath}/cells - - /Xdmf/Domain/DataItem[@Name="vertices"] + + &MeshData;:/{geomPath}/vertices {geometry_information} """ diff --git a/tests/test_0005_xdmf_compat.py b/tests/test_0005_xdmf_compat.py index 2bdf32319..fae982db6 100644 --- a/tests/test_0005_xdmf_compat.py +++ b/tests/test_0005_xdmf_compat.py @@ -85,6 +85,10 @@ def test_direct_p1_p2_and_dg0_use_fields_only(tmp_path): xdmf = directory / "compact.mesh.00000.xdmf" text = xdmf.read_text() + # ParaView's XDMF3 reader does not accept XML-reference DataItems + # without a Format attribute. Inline the two mesh HDF references so + # both its XDMF3 reader and the legacy XDMF reader can open the file. + assert 'Reference="XML"' not in text assert 'TopologyType="Triangle_6"' in text assert "&p1_Data;:/fields/p1" in text assert "&p2_Data;:/fields/p2" in text From e4d79d71e1241383be6a48c2fa429da14a33d8e5 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 18 Sep 2026 20:01:47 +1000 Subject: [PATCH 9/9] Share compact DG1 geometry across XDMF fields Store simplex DG1 corner geometry and connectivity once in the mesh HDF5 file, while retaining only each variable's exact dimensional corner values under /fields. Emit one XDMF DG1 grid with all compatible variables to remove per-variable visualization geometry duplication.\n\nTeach read_timestep to identify the compact corner-nodal representation and reconstruct native DG1 values cell by cell. Match complete cell geometry rather than coincident vertices so discontinuous traces remain distinct and reload remains exact across MPI repartitioning.\n\nUpdate DG1 and physical-unit tests to cover shared geometry, compact field files, the absence of optional restart data, XDMF references, and exact remapping. --- .../discretisation/discretisation_mesh.py | 171 +++++++--- .../discretisation_mesh_variables.py | 305 +++++++++++++++++- src/underworld3/function/field_projection.py | 37 ++- tests/test_0005_xdmf_dg1.py | 35 +- tests/test_0005_xdmf_physical_units.py | 17 +- 5 files changed, 478 insertions(+), 87 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 2751838c6..b1089c084 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -4836,10 +4836,11 @@ def write_timestep( - ``create_xdmf=True`` writes a companion XDMF file that reads the dimensional ``/fields`` datasets directly for P1, P2 triangles and - tetrahedra, and DG0. DG1 retains native ``/fields`` values for reload - and adds an exact disconnected-corner visualization. Continuous P3+ - fields receive one compact P1 visualization dataset; discontinuous - DG2+ fields receive DG0. + tetrahedra, DG0, and simplex DG1. All simplex DG1 variables share one + disconnected-corner geometry in the mesh file; each variable stores + one exact corner-nodal value array. Continuous P3+ fields receive one + compact P1 visualization dataset; discontinuous DG2+ fields receive + DG0. - ``petsc_reload=True`` writes native nondimensional PETSc DMPlex section/global-vector data under ``/restart/petsc``. Load that optional payload with ``MeshVariable.read_checkpoint()`` @@ -4937,7 +4938,7 @@ def write_timestep( var.write(save_location) if petsc_reload: self._write_petsc_reload_file(save_location, [var], mode="a") - _write_xdmf_field(self, var, save_location) + _write_xdmf_field(self, var, save_location, mesh_file) elif petsc_reload: # Exact-restart-only output needs the PETSc section/vector # payload, not a second native /fields copy of the values. @@ -9646,16 +9647,15 @@ def remesh(self, metric_field, verbose=False): return -def _write_xdmf_field(mesh, var, var_h5_path): +def _write_xdmf_field(mesh, var, var_h5_path, mesh_h5_path): """Replace native remap arrays with dimensional field output for XDMF. - P1, P2 triangles/tetrahedra, and DG0 are represented directly under ``/fields``. - DG1 keeps its native interpolation values under ``/fields`` for exact - coordinate reload and receives an exact disconnected-corner representation - under ``/visualization``. Unsupported higher-order layouts retain their - exact physical values under ``/fields`` and receive one compact - visualization reduction: continuous fields use P1 and discontinuous fields - use DG0. + P1, P2 triangles/tetrahedra, DG0, and simplex DG1 are represented directly + under ``/fields``. DG1 uses exact element-corner values and shares one + disconnected geometry under ``/viz/dg1`` in the mesh file. Unsupported + higher-order layouts retain exact physical values under ``/fields`` and + receive one compact visualization reduction: continuous fields use P1 and + discontinuous fields use DG0. Parameters ---------- @@ -9667,12 +9667,16 @@ def _write_xdmf_field(mesh, var, var_h5_path): var_h5_path : str Path to the HDF5 file. Native restart data, when requested, has already been written under ``/restart/petsc``. + mesh_h5_path : str + Path to the shared mesh HDF5 file that owns DG1 visualization geometry. """ import h5py import underworld3 as uw from underworld3.function.field_projection import ( + _dg1_corner_data, _physical_visualisation_values, - _write_dg1_to_viewer, + _write_index_array_to_group, + _write_vec_to_group, write_field_coordinates_to_viewer, write_field_to_viewer, write_p2_simplex_topology_to_viewer, @@ -9710,18 +9714,16 @@ def _write_xdmf_field(mesh, var, var_h5_path): comm=PETSc.COMM_WORLD, ) - write_field_to_viewer(var, viewer, "/fields", var.clean_name) - write_field_coordinates_to_viewer(var, viewer, "/fields") if direct_dg1: - _write_dg1_to_viewer( - var, - viewer, - group="/visualization", - coordinate_name="coordinates", - value_name=var.clean_name, - repack_tensors=False, + corner_rows, corner_values, corner_cells = _dg1_corner_data( + var, repack_tensors=False + ) + _write_vec_to_group( + viewer, corner_values, var.clean_name, "/fields", PETSc.COMM_WORLD ) else: + write_field_to_viewer(var, viewer, "/fields", var.clean_name) + write_field_coordinates_to_viewer(var, viewer, "/fields") if direct_p2: write_p2_simplex_topology_to_viewer(var, viewer, group="/fields") elif needs_projection: @@ -9736,6 +9738,36 @@ def _write_xdmf_field(mesh, var, var_h5_path): ) viewer.destroy() + + if direct_dg1: + if uw.mpi.rank == 0: + with h5py.File(mesh_h5_path, "r") as mesh_handle: + has_dg1_geometry = "viz/dg1/coordinates" in mesh_handle + else: + has_dg1_geometry = None + has_dg1_geometry = uw.mpi.comm.bcast(has_dg1_geometry, root=0) + if not has_dg1_geometry: + mesh_viewer = PETSc.ViewerHDF5().create( + mesh_h5_path, + "a", + comm=PETSc.COMM_WORLD, + ) + _write_vec_to_group( + mesh_viewer, + corner_rows, + "coordinates", + "/viz/dg1", + PETSc.COMM_WORLD, + ) + _write_index_array_to_group( + mesh_viewer, + corner_cells, + "cells", + "/viz/dg1", + PETSc.COMM_WORLD, + ) + mesh_viewer.destroy() + _, field_units = _physical_visualisation_values(numpy.ones(1), var.units) _, coordinate_units = _physical_visualisation_values(numpy.ones(1), mesh.units) with uw.selective_ranks(0) as should_execute: @@ -9747,26 +9779,30 @@ def _write_xdmf_field(mesh, var, var_h5_path): field.attrs["storage_frame"] = "physical" field.attrs["degree"] = var.degree field.attrs["continuous"] = var.continuous - field.attrs["representation"] = "exact" - coordinates = handle["fields/coordinates"] - coordinates.attrs["units"] = coordinate_units or "dimensionless" - coordinates.attrs["storage_frame"] = "physical" - if direct_dg1 or needs_projection: + field.attrs["representation"] = ( + "dg1_corner_nodal" if direct_dg1 else "exact" + ) + if direct_dg1: + handle.attrs["dg1_mesh_file"] = os.path.basename(mesh_h5_path) + handle.attrs["dg1_geometry_path"] = "/viz/dg1" + else: + coordinates = handle["fields/coordinates"] + coordinates.attrs["units"] = coordinate_units or "dimensionless" + coordinates.attrs["storage_frame"] = "physical" + if needs_projection: projected = handle[f"visualization/{var.clean_name}"] projected.attrs["units"] = field_units or "dimensionless" projected.attrs["source_degree"] = var.degree - projected.attrs["visualization_degree"] = ( - 1 if var.continuous or direct_dg1 else 0 - ) - projected.attrs["representation"] = ( - "basis_conversion" if direct_dg1 else "projection" - ) - if direct_dg1: - visual_coordinates = handle["visualization/coordinates"] - visual_coordinates.attrs["units"] = ( - coordinate_units or "dimensionless" + projected.attrs["visualization_degree"] = 1 if var.continuous else 0 + projected.attrs["representation"] = "projection" + if direct_dg1: + with h5py.File(mesh_h5_path, "a") as mesh_handle: + coordinates = mesh_handle["viz/dg1/coordinates"] + coordinates.attrs["units"] = coordinate_units or "dimensionless" + coordinates.attrs["storage_frame"] = "physical" + mesh_handle["viz/dg1/cells"].attrs["representation"] = ( + "disconnected_simplex" ) - visual_coordinates.attrs["storage_frame"] = "physical" uw.mpi.barrier() @@ -9966,7 +10002,9 @@ def direct_dg1(var): and var.mesh.cdim == var.mesh.dim ) - special_vars = [var for var in meshVars if direct_p2(var) or direct_dg1(var)] + p2_vars = [var for var in meshVars if direct_p2(var)] + dg1_vars = [var for var in meshVars if direct_dg1(var)] + special_vars = p2_vars + dg1_vars collection_start = ( '' f'""" + + if dg1_vars: + with h5py.File(mesh_filename, "r") as mesh_handle: + dg1_cells = mesh_handle["viz/dg1/cells"] + dg1_coordinates = mesh_handle["viz/dg1/coordinates"] + dg1_cells_shape = dg1_cells.shape + dg1_points_shape = dg1_coordinates.shape + dg1_topology_precision = dg1_cells.dtype.itemsize + dg1_geometry_units = dg1_coordinates.attrs.get("units") + + dg1_attributes = "" + for var in dg1_vars: + var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" + with h5py.File(var_filename, "r") as field_handle: + values = field_handle[f"fields/{var.clean_name}"] + values_shape = values.shape + field_units = values.attrs.get("units") + components = values_shape[1] if len(values_shape) == 2 else 1 + kind = attribute_kind(var, components) + dimensions = " ".join(str(value) for value in values_shape) + dg1_attributes += f""" + + + &{var.clean_name}_Data;:/fields/{var.clean_name} + {units_information(field_units, " ")} + """ + + special_grids += f""" + + + + &MeshData;:/viz/dg1/cells + + + + + &MeshData;:/viz/dg1/coordinates + {units_information(dg1_geometry_units, " ")} + +{dg1_attributes} + """ xdmf_end = f""" {special_grids} diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 42381421a..e849cc16c 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -1164,10 +1164,12 @@ def read_timestep( verbose=False, ): """ - Read a mesh variable from ``Mesh.write_timestep()`` output using the - coordinate-remap path. The saved mesh and the live mesh may have - different sizes or decompositions; values are matched to the live mesh - nodes by nearest-neighbour KDTree interpolation. + Read a mesh variable from ``Mesh.write_timestep()`` output. Standard + nodal fields use coordinate-remap interpolation. Compact simplex DG1 + fields use their exact element-corner representation: complete cells + are matched by geometry and the local affine polynomial is evaluated + at the live DG1 nodes. The latter supports MPI repartitioning without + averaging the discontinuous traces at shared vertices. This is the flexible remap reader. It is distinct from ``read_checkpoint()``, which loads PETSc DMPlex section/vector metadata @@ -1232,6 +1234,24 @@ def read_timestep( n_components = self.num_components dim = self.mesh.dim + if not is_v1_1: + if uw.mpi.rank == 0: + with h5py.File(data_file, "r") as h5f: + representation = h5f[f"fields/{data_name}"].attrs.get( + "representation" + ) + if isinstance(representation, bytes): + representation = representation.decode() + else: + representation = None + representation = uw.mpi.comm.bcast(representation, root=0) + if representation == "dg1_corner_nodal": + return self._read_dg1_corner_timestep( + data_file, + data_name, + verbose=verbose, + ) + # ---- Phase 1: source swarm carries saved (coord, value) pairs ---- source_swarm = uw.swarm.Swarm(self.mesh) saved = uw.swarm.SwarmVariable( @@ -1435,6 +1455,283 @@ def read_timestep( return + def _read_dg1_corner_timestep(self, data_file, data_name, verbose=False): + """Reload one exact simplex DG1 field from disconnected corner values. + + A coordinate-only nearest-neighbour lookup cannot distinguish the + separate traces of cells meeting at one vertex. This reader routes + complete source and destination cells by centroid, matches their full + corner geometry, and evaluates each source cell's affine polynomial at + the live variable's native DG1 nodes. The cell match remains exact when + the same mesh is repartitioned across a different MPI rank count. + """ + import h5py + import numpy as np + from mpi4py import MPI + + mesh = self.mesh + dim = mesh.dim + corner_count = dim + 1 + n_components = self.num_components + if ( + self.continuous + or self.degree != 1 + or not mesh.isSimplex + or dim not in (2, 3) + or mesh.cdim != dim + ): + raise RuntimeError( + "dg1_corner_nodal data require a discontinuous degree-one " + "variable on a full-dimensional triangle or tetrahedron mesh" + ) + + if uw.mpi.rank == 0: + with h5py.File(data_file, "r") as field_handle: + field = field_handle[f"fields/{data_name}"] + values = field[()].reshape(-1, n_components) + field_units = field.attrs.get("units") + mesh_file = field_handle.attrs.get("dg1_mesh_file") + geometry_path = field_handle.attrs.get( + "dg1_geometry_path", "/viz/dg1" + ) + for name, value in ( + ("field units", field_units), + ("mesh file", mesh_file), + ("geometry path", geometry_path), + ): + if isinstance(value, bytes): + if name == "field units": + field_units = value.decode() + elif name == "mesh file": + mesh_file = value.decode() + else: + geometry_path = value.decode() + if not mesh_file: + raise RuntimeError( + f"{data_file} does not identify its shared DG1 mesh geometry" + ) + mesh_path = os.path.join(os.path.dirname(data_file), mesh_file) + with h5py.File(mesh_path, "r") as mesh_handle: + coordinates_dataset = mesh_handle[ + f"{geometry_path.strip('/')}/coordinates" + ] + coordinates = coordinates_dataset[()].reshape(-1, dim) + coordinate_units = coordinates_dataset.attrs.get("units") + cells = mesh_handle[f"{geometry_path.strip('/')}/cells"][()].reshape( + -1, corner_count + ) + if isinstance(coordinate_units, bytes): + coordinate_units = coordinate_units.decode() + if coordinate_units and coordinate_units != "dimensionless": + coordinates = np.asarray( + uw.non_dimensionalise(uw.quantity(coordinates, coordinate_units)) + ) + if field_units and field_units != "dimensionless": + values = np.asarray( + uw.non_dimensionalise(uw.quantity(values, field_units)) + ) + source_corners = coordinates[cells] + source_values = values[cells] + source_centroids = np.round(source_corners.mean(axis=1), decimals=12) + else: + source_corners = np.empty((0, corner_count, dim), dtype=np.float64) + source_values = np.empty( + (0, corner_count, n_components), dtype=np.float64 + ) + source_centroids = np.empty((0, dim), dtype=np.float64) + + source_swarm = uw.swarm.Swarm(mesh) + saved_corners = uw.swarm.SwarmVariable( + "_dg1_source_corners", + source_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * dim), + dtype=float, + _proxy=False, + varsymbol=r"\cal{C}_s", + ) + saved_values = uw.swarm.SwarmVariable( + "_dg1_source_values", + source_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * n_components), + dtype=float, + _proxy=False, + varsymbol=r"\cal{V}_s", + ) + source_before = max(source_swarm.dm.getLocalSize(), 0) + source_swarm.add_particles_with_global_coordinates( + source_centroids, migrate=False + ) + source_swarm._invalidate_canonical_data() + saved_corners.array[source_before:, 0, :] = source_corners.reshape( + -1, corner_count * dim + ) + saved_values.array[source_before:, 0, :] = source_values.reshape( + -1, corner_count * n_components + ) + source_swarm._route_by_nearest_centroid() + + landed_corners = np.asarray(saved_corners.array)[:, 0, :].reshape( + -1, corner_count, dim + ) + landed_values = np.asarray(saved_values.array)[:, 0, :].reshape( + -1, corner_count, n_components + ) + + cell_start, cell_end = mesh.dm.getHeightStratum(0) + owned = np.ones(cell_end - cell_start, dtype=bool) + if mesh.dm.comm.getSize() > 1: + _, leaves, remote = mesh.dm.getPointSF().getGraph() + if leaves is None: + leaves = np.arange(len(remote)) + leaves = np.asarray(leaves) + owned[leaves[(leaves >= cell_start) & (leaves < cell_end)] - cell_start] = ( + False + ) + native_rows = mesh._cell_node_indices(1, False).reshape(-1, corner_count)[owned] + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, corner_count)[owned] + target_nodes = np.asarray(self.coords_nd)[native_rows] + target_corners = mesh._get_coords_for_basis(1, True)[vertex_rows] + target_centroids = np.round(target_corners.mean(axis=1), decimals=12) + n_owned = len(native_rows) + + query_swarm = uw.swarm.Swarm(mesh) + origin_rank = uw.swarm.SwarmVariable( + "_dg1_origin_rank", + query_swarm, + vtype=uw.VarType.SCALAR, + dtype=int, + _proxy=False, + varsymbol=r"\cal{R}_d", + ) + origin_index = uw.swarm.SwarmVariable( + "_dg1_origin_index", + query_swarm, + vtype=uw.VarType.SCALAR, + dtype=int, + _proxy=False, + varsymbol=r"\cal{I}_d", + ) + query_nodes = uw.swarm.SwarmVariable( + "_dg1_query_nodes", + query_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * dim), + dtype=float, + _proxy=False, + varsymbol=r"\cal{X}_d", + ) + query_corners = uw.swarm.SwarmVariable( + "_dg1_query_corners", + query_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * dim), + dtype=float, + _proxy=False, + varsymbol=r"\cal{C}_d", + ) + result = uw.swarm.SwarmVariable( + "_dg1_result", + query_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * n_components), + dtype=float, + _proxy=False, + varsymbol=r"\cal{D}_d", + ) + query_before = max(query_swarm.dm.getLocalSize(), 0) + query_swarm.add_particles_with_global_coordinates( + target_centroids, migrate=False + ) + query_swarm._invalidate_canonical_data() + origin_rank.array[query_before:, 0, 0] = uw.mpi.rank + origin_index.array[query_before:, 0, 0] = np.arange(n_owned) + query_nodes.array[query_before:, 0, :] = target_nodes.reshape( + -1, corner_count * dim + ) + query_corners.array[query_before:, 0, :] = target_corners.reshape( + -1, corner_count * dim + ) + query_swarm._route_by_nearest_centroid() + + def cell_key(corners): + rounded = np.round(np.asarray(corners), decimals=12) + order = np.lexsort( + tuple(rounded[:, axis] for axis in range(dim - 1, -1, -1)) + ) + return tuple(rounded[order].reshape(-1)) + + source_by_cell = {} + duplicate_source = 0 + for source_index, corners in enumerate(landed_corners): + key = cell_key(corners) + if key in source_by_cell: + duplicate_source += 1 + source_by_cell[key] = source_index + + local_query_nodes = np.asarray(query_nodes.array)[:, 0, :].reshape( + -1, corner_count, dim + ) + local_query_corners = np.asarray(query_corners.array)[:, 0, :].reshape( + -1, corner_count, dim + ) + missing = 0 + for query_index, (query_cell_corners, nodes) in enumerate( + zip(local_query_corners, local_query_nodes, strict=True) + ): + key = cell_key(query_cell_corners) + matched_index = source_by_cell.get(key) + if matched_index is None: + missing += 1 + continue + corners = landed_corners[matched_index] + coefficients = landed_values[matched_index] + matrix = (corners[1:] - corners[:1]).T + local = np.linalg.solve(matrix, (nodes - corners[:1]).T).T + weights = np.column_stack((1.0 - local.sum(axis=1), local)) + result.array[query_index, 0, :] = (weights @ coefficients).reshape(-1) + + failures = uw.mpi.comm.allreduce(missing + duplicate_source, op=MPI.SUM) + if failures: + raise RuntimeError( + "DG1 corner reload could not match every destination cell to " + "one source cell. The cell-aware path supports identical mesh " + "geometry across arbitrary MPI repartitioning." + ) + + query_swarm._rank_var.array[...] = origin_rank.array[...] + query_swarm.dm.migrate(remove_sent_points=True) + uw.mpi.barrier() + query_swarm._invalidate_canonical_data() + + returned_indices = origin_index.array[:, 0, 0].astype(int) + returned = np.asarray(result.array)[:, 0, :].reshape( + -1, corner_count, n_components + ) + cell_values = np.empty((n_owned, corner_count, n_components)) + cell_values[returned_indices] = returned + self.data[...] = 0.0 + self.data[native_rows.reshape(-1), :] = cell_values.reshape( + -1, n_components + ) + + indexset, subdm = mesh.dm.createSubDM(self.field_id) + subdm.localToGlobal(self._lvec, self._gvec, addv=False) + subdm.globalToLocal(self._gvec, self._lvec, addv=False) + indexset.destroy() + subdm.destroy() + if mesh._lvec is not None: + mesh._lvec.destroy() + mesh._lvec = None + mesh._stale_lvec = True + if verbose and uw.mpi.rank == 0: + print( + f"Reloaded DG1 field {data_name!r} from exact corner-nodal data", + flush=True, + ) + return + @timing.routine_timer_decorator def load_from_h5_plex_vector( self, diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py index bfd5b8b83..4ed0c713f 100644 --- a/src/underworld3/function/field_projection.py +++ b/src/underworld3/function/field_projection.py @@ -538,20 +538,14 @@ def write_cell_field_to_viewer( _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) -def _write_dg1_to_viewer( - mesh_var, - viewer, - group="/dg1", - coordinate_name="vertices", - value_name="values", - repack_tensors=True, -): - """Write owned simplex cells with independent vertices and DG1 traces. +def _dg1_corner_data(mesh_var, repack_tensors=False): + """Return owned disconnected-cell geometry and exact DG1 corner values. Coordinate-section cell maps preserve element ownership and node ordering; no point location, coordinate matching, or inter-element averaging is used. Interior DG interpolation nodes define an affine polynomial, evaluated at - that same cell's vertices. Native checkpoint vectors are untouched. + that same cell's vertices. The corner values are an exact alternative basis + for the element-local linear polynomial, not a lower-order projection. """ mesh = mesh_var.mesh if ( @@ -585,10 +579,6 @@ def _write_dg1_to_viewer( corners, _ = _physical_visualisation_values(corners, mesh.units) values, _ = _physical_visualisation_values(values, mesh_var.units) corner_rows = corners.reshape(-1, mesh.cdim) - _write_vec_to_group( - viewer, corner_rows, coordinate_name, group, PETSc.COMM_WORLD - ) - _write_vec_to_group(viewer, values, value_name, group, PETSc.COMM_WORLD) local_count = len(corner_rows) offset = mesh.dm.comm.tompi4py().exscan(local_count) if offset is None: @@ -596,4 +586,23 @@ def _write_dg1_to_viewer( cells = np.arange(offset, offset + local_count, dtype=PETSc.IntType).reshape( -1, mesh.dim + 1 ) + return corner_rows, values, cells + + +def _write_dg1_to_viewer( + mesh_var, + viewer, + group="/dg1", + coordinate_name="vertices", + value_name="values", + repack_tensors=True, +): + """Write owned simplex cells with independent vertices and DG1 traces.""" + corner_rows, values, cells = _dg1_corner_data( + mesh_var, repack_tensors=repack_tensors + ) + _write_vec_to_group( + viewer, corner_rows, coordinate_name, group, PETSc.COMM_WORLD + ) + _write_vec_to_group(viewer, values, value_name, group, PETSc.COMM_WORLD) _write_index_array_to_group(viewer, cells, "cells", group, PETSc.COMM_WORLD) diff --git a/tests/test_0005_xdmf_dg1.py b/tests/test_0005_xdmf_dg1.py index 69e7ed219..51708b38b 100644 --- a/tests/test_0005_xdmf_dg1.py +++ b/tests/test_0005_xdmf_dg1.py @@ -74,19 +74,18 @@ def test_dg1_simplex_output(tmp_path, dim): if uw.mpi.rank != 0: return with h5py.File(directory / "fields.mesh.dg_scalar.00000.h5", "r") as handle: - native_points = handle["fields/coordinates"][:] - native_values = handle["fields/dg_scalar"][:].reshape(-1) - points = handle["visualization/coordinates"][:] - cells = handle["visualization/cells"][:] - values = handle["visualization/dg_scalar"][:].reshape(-1) - assert handle["fields/dg_scalar"].attrs["representation"] == "exact" - assert handle["visualization/dg_scalar"].attrs["representation"] == "basis_conversion" - assert len(native_points) == len(native_values) - assert len(points) == len(cells) * (dim + 1) - assert len(np.unique(cells)) == len(points) + values = handle["fields/dg_scalar"][:].reshape(-1) + assert handle["fields/dg_scalar"].attrs["representation"] == "dg1_corner_nodal" + assert set(handle["fields"]) == {"dg_scalar"} + assert "visualization" not in handle assert "dg1" not in handle with h5py.File(directory / "fields.mesh.00000.h5", "r") as handle: + points = handle["viz/dg1/coordinates"][:] + cells = handle["viz/dg1/cells"][:] assert len(cells) == len(handle["viz/topology/cells"]) + assert len(points) == len(cells) * (dim + 1) + assert len(values) == len(points) + assert len(np.unique(cells)) == len(points) corners = points[cells] determinants = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) @@ -104,7 +103,7 @@ def test_dg1_simplex_output(tmp_path, dim): assert np.max(high - low) >= 1 with h5py.File(directory / "fields.mesh.dg_tensor.00000.h5", "r") as handle: - tensor_values = handle["visualization/dg_tensor"][:] + tensor_values = handle["fields/dg_tensor"][:] assert tensor_values.shape == (len(points), dim * dim) np.testing.assert_allclose(tensor_values[:, 0], expected) np.testing.assert_allclose(tensor_values[:, 1], 3 + points[:, 0]) @@ -113,13 +112,15 @@ def test_dg1_simplex_output(tmp_path, dim): tree = ET.parse(directory / "fields.mesh.00000.xdmf") grids = tree.findall(".//Grid[@GridType='Uniform']") - assert {grid.get("Name") for grid in grids} == {"domain", "dg_scalar", "dg_tensor"} - for name in ("dg_scalar", "dg_tensor"): - grid = next(grid for grid in grids if grid.get("Name") == name) - attribute = grid.find("Attribute") - assert attribute.get("Name") == name + assert {grid.get("Name") for grid in grids} == {"domain", "DG1"} + grid = next(grid for grid in grids if grid.get("Name") == "DG1") + assert grid.find("Topology/DataItem").text.strip().endswith("/viz/dg1/cells") + assert grid.find("Geometry/DataItem").text.strip().endswith("/viz/dg1/coordinates") + attributes = {attribute.get("Name"): attribute for attribute in grid.findall("Attribute")} + assert set(attributes) == {"dg_scalar", "dg_tensor"} + for name, attribute in attributes.items(): assert attribute.get("Center") == "Node" - assert f"/visualization/{name}" in attribute.find("DataItem").text + assert f"/fields/{name}" in attribute.find("DataItem").text @pytest.mark.level_1 diff --git a/tests/test_0005_xdmf_physical_units.py b/tests/test_0005_xdmf_physical_units.py index 52b38e7a5..3a3230c6d 100644 --- a/tests/test_0005_xdmf_physical_units.py +++ b/tests/test_0005_xdmf_physical_units.py @@ -98,7 +98,7 @@ def test_fields_are_dimensional_and_checkpoint_is_native(tmp_path): @pytest.mark.level_1 @pytest.mark.tier_b def test_dg1_fields_are_physical_at_element_corners(tmp_path): - """DG1 native reload data and corner visualization use physical units.""" + """DG1 uses one physical corner field and shared mesh geometry.""" _set_reference_scales() mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) pressure = uw.discretisation.MeshVariable( @@ -122,10 +122,15 @@ def test_dg1_fields_are_physical_at_element_corners(tmp_path): if uw.mpi.rank == 0: with h5py.File(directory / "dg.mesh.dg_pressure.00000.h5", "r") as handle: np.testing.assert_allclose(handle["fields/dg_pressure"][:], 8.0) - np.testing.assert_allclose(handle["visualization/dg_pressure"][:], 8.0) - assert np.isclose(handle["visualization/coordinates"][:].max(), 10.0) assert handle["fields/dg_pressure"].attrs["units"] == "megapascal" - assert handle["fields/coordinates"].attrs["units"] == "kilometer" - assert handle["visualization/dg_pressure"].attrs["units"] == "megapascal" - assert handle["visualization/coordinates"].attrs["units"] == "kilometer" + assert ( + handle["fields/dg_pressure"].attrs["representation"] + == "dg1_corner_nodal" + ) + assert set(handle["fields"]) == {"dg_pressure"} + assert "visualization" not in handle assert "dg1" not in handle + assert "restart" not in handle + with h5py.File(directory / "dg.mesh.00000.h5", "r") as handle: + assert np.isclose(handle["viz/dg1/coordinates"][:].max(), 10.0) + assert handle["viz/dg1/coordinates"].attrs["units"] == "kilometer"