diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index c0502b89..c1445ba4 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 7f045997..679f025b 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 28fc4809..91a01b2c 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 00000000..163aea3c --- /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 + )