Refine XDMF compat groups: field projection, standalone Vec I/O, tensor repacking by lmoresi · Pull Request #72 · underworldcode/underworld3

lmoresi · 2026-03-05T11:15:28Z

Summary

Builds on @gthyagi's work in PR #69, which added the create_xdmf documentation
and established the design for ParaView-compatible XDMF output from write_timestep().
This PR implements the compatibility group writer with refinements:

New field_projection utility (uw.function.field_projection): general-purpose
degree transfer using PETSc createInterpolation, replacing the KDTree-based remap.
project_to_vertices(var) returns a numpy array; write_vertices_to_viewer() and
write_cell_field_to_viewer() write directly to an open PETSc ViewerHDF5.
Standalone Vec writing: DM-associated Vecs ignore pushGroup() (the DMPlex HDF5
writer redirects to /fields/). The new helpers create plain Vecs with
setBlockSize(ncomp) so pushGroup is respected and datasets are written as (N, ncomp).
Tensor repacking: converts UW3's internal _data_layout ordering to ParaView's
9-component row-major 3x3 format for the visualisation copy in /vertex_fields/.
Corrects the 2D non-symmetric tensor mapping from the original implementation
([xx, yy, xy, yx] assumed vs actual [xx, xy, yx, yy] per _data_layout).
Checkpoint data in /fields/ remains in UW3's native ordering.
write_timestep defaults to create_xdmf=True, matching the previous
petsc_save_checkpoint behaviour.
Tests for 2D/3D vertex fields, cell fields, checkpoint round-trip integrity,
and create_xdmf=False opt-out.

Test plan

All 786 existing tests pass (0 failures, 176 skipped, 20 xfailed)
New test_0005_xdmf_compat.py covers: P1 scalar, P2 vector, 3D mesh,
cell variable, checkpoint round-trip, xdmf-disabled mode
Manual verification: open XDMF in ParaView, confirm tensor display

Underworld development team with AI support from Claude Code

…d of KDTree remap Replace _write_vertex_cell_groups (274 lines, chunked KDTree nearest-neighbour remap from HDF5 file data) with _write_compat_groups (~80 lines) that uses uw.function.evaluate for exact FE interpolation at vertex positions. Strategy by variable type: - P1 continuous: direct copy from _gvec (DOFs = vertices) - P2+ continuous: uw.function.evaluate(var.sym, vertex_coords) - Cell/DG-0: direct copy to /cell_fields/ Also adds _repack_tensor_to_paraview with correct _data_layout mappings (fixes 2D TENSOR: [xx,xy,yx,yy] not [xx,yy,xy,yx]), changes create_xdmf default to True, removes vertex_cell_chunk_size parameter, and adds test_0005_xdmf_compat.py with 5 tests covering 2D/3D/cell/roundtrip/disabled. Underworld development team with AI support from Claude Code

… creation project_to_degree() and project_to_vertices() use a transient scratch DM with PETSc createInterpolation to project/prolong MeshVariable data between polynomial degrees. No MeshVariable is created, the mesh DM is never modified, and all scratch objects are destroyed before returning. Use cases: P2→P1 for visualisation/XDMF, P1→P2 prolongation for initialisation, vertex value extraction from any higher-order field. Underworld development team with AI support from Claude Code

… I/O Replace _write_compat_groups implementation (uw.function.evaluate + gather to rank 0 + h5py) with write_vertices_to_viewer / write_cell_field_to_viewer from field_projection module. These use PETSc createInterpolation for P2->P1 projection and standalone Vecs (no DM) for ViewerHDF5 writes, so pushGroup is respected and parallel I/O is handled natively by PETSc. Key insight: DM-associated Vecs ignore pushGroup because the DMPlex HDF5 writer internally pushes /fields/. Standalone Vecs honour the group stack. Also adds stale build cache documentation to CLAUDE.md and fixes tests to handle PETSc's 1D dataset shape for single-component fields. Underworld development team with AI support from Claude Code

… output write_vertices_to_viewer now repacks TENSOR and SYM_TENSOR variables from UW3's _data_layout ordering to ParaView's row-major 3x3 (9 components). Checkpoint data in /fields/ is unchanged — only the /vertex_fields/ copy is repacked. Corrects the 2D TENSOR component mapping from the original PR #69 which assumed [xx, yy, xy, yx]; the actual _data_layout is [xx, xy, yx, yy] (row-major). Also sets Vec block size in _write_vec_to_group so HDF5 datasets are written as 2D (N, ncomp) instead of flat 1D. Updates create_xdmf.md with tensor ordering details and PETSc I/O notes. Underworld development team with AI support from Claude Code

Copilot

Pull request overview

This PR replaces the KDTree-based remap approach from PR #69 with a PETSc createInterpolation-based field projection for writing ParaView-compatible XDMF output from write_timestep(). It introduces a new field_projection module for degree transfer, standalone Vec I/O helpers to work around DMPlex HDF5 writer limitations, and tensor repacking with corrected 2D component ordering.

Changes:

Added src/underworld3/function/field_projection.py with project_to_degree, project_to_vertices, write_vertices_to_viewer, write_cell_field_to_viewer, and tensor repacking utilities, replacing the removed KDTree-based _write_vertex_cell_groups function.
Changed write_timestep to default create_xdmf=True and simplified the compat group writer to a ~40-line _write_compat_groups function that delegates to the new field projection helpers.
Added tests/test_0005_xdmf_compat.py with tests for 2D/3D vertex fields, cell fields, checkpoint round-trip, and XDMF opt-out; updated docs/beginner/create_xdmf.md and CLAUDE.md accordingly.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 5 comments.

Show a summary per file

File	Description
`src/underworld3/function/field_projection.py`	New module with PETSc-based field projection, tensor repacking, and standalone Vec HDF5 writing utilities
`src/underworld3/function/__init__.py`	Exports new field projection functions as public API
`src/underworld3/discretisation/discretisation_mesh.py`	Replaces `_write_vertex_cell_groups` with simplified `_write_compat_groups`; changes `create_xdmf` default to `True`; removes `vertex_cell_chunk_size` parameter
`tests/test_0005_xdmf_compat.py`	New tests for XDMF compat group generation: 2D, 3D, cell variables, round-trip, and disabled mode
`docs/beginner/create_xdmf.md`	Updated documentation to reflect new PETSc interpolation approach and corrected tensor layout
`CLAUDE.md`	Updated build commands and added stale build cache troubleshooting

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-03-05T11:24:58Z

src/underworld3/function/field_projection.py

+    result = l_dst.array.reshape(-1, nc).copy()
+
+    # --- cleanup ---
+    for obj in (interp, g_src, g_dst, l_dst, dm_scratch, fe_target):


iset and subdm_src returned by createSubDM are never destroyed. The codebase convention (e.g., discretisation_mesh.py:1341-1342, discretisation_mesh_variables.py:685-686) is to always call .destroy() on both the index set and the sub-DM. Add iset.destroy() and subdm_src.destroy() to the cleanup block to prevent PETSc resource leaks.

Suggested change

for obj in (interp, g_src, g_dst, l_dst, dm_scratch, fe_target):

for obj in (interp, g_src, g_dst, l_dst, dm_scratch, fe_target, subdm_src, iset):

Copilot · 2026-03-05T11:24:58Z

src/underworld3/function/field_projection.py

+
+        data = g_dst.array.reshape(-1, nc).copy()
+
+        for obj in (interp, g_src, g_dst, dm_scratch, fe_target):


iset and subdm_src returned by createSubDM are never destroyed here either. Same issue as in project_to_degree — both the index set and the sub-DM should be added to the cleanup loop. See discretisation_mesh.py:1341-1342 for the established pattern.

Suggested change

for obj in (interp, g_src, g_dst, dm_scratch, fe_target):

for obj in (interp, g_src, g_dst, dm_scratch, fe_target, subdm_src, iset):

Copilot · 2026-03-05T11:24:59Z

src/underworld3/function/field_projection.py

+        # Build scratch DM at degree 1, interpolate, extract array
+        options = PETSc.Options()
+        prefix = "_fieldproj_"
+        options.setValue(f"{prefix}petscspace_degree", 1)
+        options.setValue(f"{prefix}petscdualspace_lagrange_continuity", True)
+        options.setValue(f"{prefix}petscdualspace_lagrange_node_endpoints", False)
+
+        fe_target = PETSc.FE().createDefault(
+            mesh.dim, nc, mesh.isSimplex, mesh.qdegree, prefix, PETSc.COMM_SELF,
+        )
+
+        dm_scratch = mesh.dm.clone()
+        dm_scratch.addField(fe_target)
+        dm_scratch.createDS()
+
+        iset, subdm_src = mesh.dm.createSubDM(mesh_var.field_id)
+        interp, _scale = subdm_src.createInterpolation(dm_scratch)
+
+        g_src = subdm_src.createGlobalVec()
+        subdm_src.localToGlobal(mesh_var._lvec, g_src)
+
+        g_dst = dm_scratch.createGlobalVec()
+        interp.mult(g_src, g_dst)
+
+        data = g_dst.array.reshape(-1, nc).copy()
+
+        for obj in (interp, g_src, g_dst, dm_scratch, fe_target):
+            obj.destroy()
+        if _scale is not None:
+            _scale.destroy()
+


The interpolation logic (lines 271-300) duplicates the implementation already in project_to_degree() (lines 57-101). The non-P1 branch could be simplified to data = project_to_vertices(mesh_var), avoiding the duplicated scratch-DM setup, interpolation, and cleanup code. This would reduce maintenance burden and keep the interpolation logic in a single place.

Suggested change

# Build scratch DM at degree 1, interpolate, extract array

options = PETSc.Options()

prefix = "_fieldproj_"

options.setValue(f"{prefix}petscspace_degree", 1)

options.setValue(f"{prefix}petscdualspace_lagrange_continuity", True)

options.setValue(f"{prefix}petscdualspace_lagrange_node_endpoints", False)

fe_target = PETSc.FE().createDefault(

mesh.dim, nc, mesh.isSimplex, mesh.qdegree, prefix, PETSc.COMM_SELF,

)

dm_scratch = mesh.dm.clone()

dm_scratch.addField(fe_target)

dm_scratch.createDS()

iset, subdm_src = mesh.dm.createSubDM(mesh_var.field_id)

interp, _scale = subdm_src.createInterpolation(dm_scratch)

g_src = subdm_src.createGlobalVec()

subdm_src.localToGlobal(mesh_var._lvec, g_src)

g_dst = dm_scratch.createGlobalVec()

interp.mult(g_src, g_dst)

data = g_dst.array.reshape(-1, nc).copy()

for obj in (interp, g_src, g_dst, dm_scratch, fe_target):

obj.destroy()

if _scale is not None:

_scale.destroy()

# Use shared interpolation logic to obtain vertex (P1) values.

data = project_to_vertices(mesh_var)

Copilot · 2026-03-05T11:24:59Z

tests/test_0005_xdmf_compat.py

+    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


The PR description highlights tensor repacking as a key feature (including a correction to the 2D tensor mapping from PR #69), but no tests cover tensor or symmetric tensor variables. The _repack_tensor_to_paraview function has 6 distinct branches. Consider adding at least one test with a TENSOR or SYM_TENSOR variable to validate the repacking logic and prevent regressions.

Suggested change

del mesh

del mesh

def test_tensor_variable_repacking_creates_vertex_field(tmp_path):

"""Tensor MeshVariable should produce a repacked tensor dataset in vertex_fields."""

# 2D mesh so tensor repacking must handle 2D->3D mapping for ParaView.

mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4))

# Full tensor in 2D has 4 independent components.

T = uw.discretisation.MeshVariable(

"T", mesh, 4, degree=1, vtype=uw.VarType.TENSOR

)

# Initialise tensor components with a simple pattern based on coordinates

coords = mesh._coords

T.data[:, 0] = coords[:, 0] # T_xx

T.data[:, 1] = coords[:, 1] # T_yy

T.data[:, 2] = coords[:, 0] + coords[:, 1] # T_xy

T.data[:, 3] = 1.0 # T_yx (constant)

mesh.write_timestep(

"tensor", index=0, outputPath=str(tmp_path), meshVars=[T]

)

# Open the compat HDF5 file for the tensor variable.

t_h5 = os.path.join(str(tmp_path), "tensor.mesh.T.00000.h5")

assert os.path.exists(t_h5), "Tensor HDF5 file was not created"

with h5py.File(t_h5, "r") as f:

assert "vertex_fields" in f, "vertex_fields group missing for tensor variable"

vf = f["vertex_fields"]

# Find the dataset corresponding to the tensor variable.

tensor_ds = None

for name, ds in vf.items():

if "T" in name:

tensor_ds = ds[()]

break

assert tensor_ds is not None, "No tensor dataset found in vertex_fields"

# Expect tensors to be repacked into 9 components (3x3) for ParaView.

assert tensor_ds.ndim == 2, "Tensor vertex field should be 2D (npoints, ncomp)"

assert tensor_ds.shape[1] == 9, (

f"Expected repacked tensor with 9 components, got {tensor_ds.shape[1]}"

)

del mesh

Copilot · 2026-03-05T11:24:59Z

src/underworld3/discretisation/discretisation_mesh.py

+        from underworld3.function.field_projection import _write_vec_to_group

-            source_tree = uw.kdtree.KDTree(source_coordinates)
-            nearest_distance, nearest_index = source_tree.query(
-                local_vertices, k=1, sqr_dists=False
-            )
-            nearest_distance = np.asarray(nearest_distance).reshape(-1)
-            nearest_index = np.asarray(nearest_index).reshape(-1)
-            nearest_distance_sq = nearest_distance * nearest_distance
-            improved = nearest_distance_sq < local_best_distance_sq
-
-            if np.any(improved):
-                local_best_distance_sq[improved] = nearest_distance_sq[improved]
-                local_best_values[improved, :] = source_values[nearest_index[improved], :]
-
-        if local_best_values is None:
-            # Degenerate case: empty source
-            local_best_values = np.zeros((local_vertex_count, 1), dtype=field_dataset_dtype)
-
-        if is_tensor:
-            local_component_count = local_best_values.shape[1]
-            tensor_values = np.zeros((local_vertex_count, 9), dtype=field_dataset_dtype)
-
-            if dim == 2 and local_component_count == 4:
-                tensor_values[:, 0] = local_best_values[:, 0]
-                tensor_values[:, 4] = local_best_values[:, 1]
-                tensor_values[:, 1] = local_best_values[:, 2]
-                tensor_values[:, 3] = local_best_values[:, 3]
-            elif dim == 2 and local_component_count == 3:
-                tensor_values[:, 0] = local_best_values[:, 0]
-                tensor_values[:, 4] = local_best_values[:, 1]
-                tensor_values[:, 1] = local_best_values[:, 2]
-                tensor_values[:, 3] = local_best_values[:, 2]
-            elif dim == 3 and local_component_count == 9:
-                tensor_values[:, :] = local_best_values[:, :]
-            elif dim == 3 and local_component_count == 6:
-                tensor_values[:, 0] = local_best_values[:, 0]
-                tensor_values[:, 4] = local_best_values[:, 1]
-                tensor_values[:, 8] = local_best_values[:, 2]
-                tensor_values[:, 1] = local_best_values[:, 3]
-                tensor_values[:, 3] = local_best_values[:, 3]
-                tensor_values[:, 5] = local_best_values[:, 4]
-                tensor_values[:, 7] = local_best_values[:, 4]
-                tensor_values[:, 2] = local_best_values[:, 5]
-                tensor_values[:, 6] = local_best_values[:, 5]
-            else:
-                tensor_values = None
-
-            if tensor_values is not None:
-                local_best_values = tensor_values
-
-        output_components = local_best_values.shape[1]
-        vertex_group = variable_file.require_group("vertex_fields")
-        if rank == 0:
-            if "coordinates" in vertex_group:
-                del vertex_group["coordinates"]
-            if output_dataset_name in vertex_group:
-                del vertex_group[output_dataset_name]
-        if parallel_h5:
-            comm.Barrier()
-
-        coordinates_output = vertex_group.require_dataset(
-            "coordinates", shape=(num_vertices, dim), dtype=local_vertices.dtype
-        )
-        values_output = vertex_group.require_dataset(
-            output_dataset_name,
-            shape=(num_vertices, output_components),
-            dtype=field_dataset_dtype,
-        )
-        coordinates_output[local_vertex_start:local_vertex_stop, :] = local_vertices
-        values_output[local_vertex_start:local_vertex_stop, :] = local_best_values
+        coord_gvec = mesh.dm.getCoordinates()
+        coords = coord_gvec.array.reshape(-1, mesh.dim).copy()
+        _write_vec_to_group(viewer, coords, "coordinates",
+                            "/vertex_fields", PETSc.COMM_WORLD)


Importing the private helper _write_vec_to_group directly from the module (from underworld3.function.field_projection import _write_vec_to_group) bypasses the public API. This creates a coupling to an internal implementation detail. Consider either making _write_vec_to_group part of the public exports in __init__.py, or encapsulating the coordinate-writing logic within a public function in field_projection.py (e.g., a write_coordinates_to_viewer helper).

lmoresi · 2026-03-05T11:28:33Z

Follow-up to @gthyagi's PR for xdmf in write_timestep. This is just simple things to do with memory efficiency / reducing file writing. There is a low-overhead projection solver that is added to provide consistent sampling of higher-order fields onto the mesh nodes.

- Fix test_no_xdmf_when_disabled: check for specific dataset (vertex_fields/s_s) rather than group existence, since some PETSc versions create /vertex_fields/ as a side effect of DM writes - Fix PETSc resource leak: destroy iset and subdm_src from createSubDM() in both project_to_degree and write_vertices_to_viewer - Add write_coordinates_to_viewer() public API to replace private _write_vec_to_group import in discretisation_mesh.py - Add test_tensor_variable_repacking: validates 2D tensor is repacked to 9 components in /vertex_fields/ - Add comment explaining why write_vertices_to_viewer cannot reuse project_to_vertices (global vs local vector data for MPI writing) Underworld development team with AI support from Claude Code

PETSc 3.21 (conda) automatically writes /vertex_fields/ during var.write(), which conflicts with our compat group writer (error 76 on duplicate dataset, and false test failures). Fix: delete any pre-existing compat group before writing our own. This ensures consistent tensor repacking and correct dataset shapes regardless of PETSc version. Simplify test_no_xdmf_when_disabled to check XDMF file absence only — HDF5 group contents are PETSc-version-dependent. Underworld development team with AI support from Claude Code

gthyagi · 2026-03-10T04:34:44Z

@lmoresi I tested this branch and the outputs look correct. Thanks for the projection function for a more parallel-efficient approach and for fixing the PyVista tensor repacking. Ready to merge into development.

lmoresi added 4 commits March 5, 2026 22:13

Copilot AI review requested due to automatic review settings March 5, 2026 11:15

Copilot started reviewing on behalf of lmoresi March 5, 2026 11:16 View session

Copilot AI reviewed Mar 5, 2026

View reviewed changes

lmoresi assigned gthyagi Mar 5, 2026

lmoresi force-pushed the feature/xdmf-compat-refinement branch from 3cf8bdc to 96431a6 Compare March 5, 2026 20:49

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Refine XDMF compat groups: field projection, standalone Vec I/O, tensor repacking#72

Refine XDMF compat groups: field projection, standalone Vec I/O, tensor repacking#72
lmoresi wants to merge 6 commits intodevelopmentfrom
feature/xdmf-compat-refinement

lmoresi commented Mar 5, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Copilot AI Mar 5, 2026

Uh oh!

Copilot AI Mar 5, 2026

Uh oh!

Copilot AI Mar 5, 2026

Uh oh!

Copilot AI Mar 5, 2026

Uh oh!

Copilot AI Mar 5, 2026

Uh oh!

lmoresi commented Mar 5, 2026

Uh oh!

gthyagi commented Mar 10, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

	for obj in (interp, g_src, g_dst, l_dst, dm_scratch, fe_target):
	for obj in (interp, g_src, g_dst, l_dst, dm_scratch, fe_target, subdm_src, iset):


		data = g_dst.array.reshape(-1, nc).copy()

		for obj in (interp, g_src, g_dst, dm_scratch, fe_target):

-    del mesh
+    del mesh
+def test_tensor_variable_repacking_creates_vertex_field(tmp_path):
+    """Tensor MeshVariable should produce a repacked tensor dataset in vertex_fields."""
+    # 2D mesh so tensor repacking must handle 2D->3D mapping for ParaView.
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4))
+    # Full tensor in 2D has 4 independent components.
+    T = uw.discretisation.MeshVariable(
+        "T", mesh, 4, degree=1, vtype=uw.VarType.TENSOR
+    )
+    # Initialise tensor components with a simple pattern based on coordinates
+    coords = mesh._coords
+    T.data[:, 0] = coords[:, 0]              # T_xx
+    T.data[:, 1] = coords[:, 1]              # T_yy
+    T.data[:, 2] = coords[:, 0] + coords[:, 1]  # T_xy
+    T.data[:, 3] = 1.0                       # T_yx (constant)
+    mesh.write_timestep(
+        "tensor", index=0, outputPath=str(tmp_path), meshVars=[T]
+    )
+    # Open the compat HDF5 file for the tensor variable.
+    t_h5 = os.path.join(str(tmp_path), "tensor.mesh.T.00000.h5")
+    assert os.path.exists(t_h5), "Tensor HDF5 file was not created"
+    with h5py.File(t_h5, "r") as f:
+        assert "vertex_fields" in f, "vertex_fields group missing for tensor variable"
+        vf = f["vertex_fields"]
+        # Find the dataset corresponding to the tensor variable.
+        tensor_ds = None
+        for name, ds in vf.items():
+            if "T" in name:
+                tensor_ds = ds[()]
+                break
+        assert tensor_ds is not None, "No tensor dataset found in vertex_fields"
+        # Expect tensors to be repacked into 9 components (3x3) for ParaView.
+        assert tensor_ds.ndim == 2, "Tensor vertex field should be 2D (npoints, ncomp)"
+        assert tensor_ds.shape[1] == 9, (
+            f"Expected repacked tensor with 9 components, got {tensor_ds.shape[1]}"
+        )
+    del mesh

Conversation

lmoresi commented Mar 5, 2026

Summary

Test plan

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Mar 5, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 5, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 5, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 5, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Mar 5, 2026

Choose a reason for hiding this comment

Uh oh!

lmoresi commented Mar 5, 2026

Uh oh!

gthyagi commented Mar 10, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants