Consider replacing manual /viz topology assembly with PETSc HDF5_VIZ viewer

[lmoresi]

Tracking the cleaner alternative to the XDMF viz-topology fix landed in PR #205.

Background

PR #205 (tiannh7) added a Python+MPI helper _write_mesh_viz_groups that gathers vertex coordinates and cell-to-vertex connectivity across ranks, deduplicates global vertex ids, and writes /viz/geometry/vertices + /viz/topology/cells into the mesh HDF5. This fixes ParaView crashes reading the XDMF generated by write_timestep. ~115 lines of new code, including 3D simplex/quad cell reordering.

The alternative @gthyagi developed

Commit 22dd9ab2 (gthyagi, on a now-deleted branch feature/annulus-spherical-benchmark-debug) addressed the same bug with a much smaller change:

• Add petsc_format: bool = False kwarg to Mesh.write()
• write_timestep calls self.write(mesh_file, petsc_format=False) → uses PETSc.Viewer.Format.HDF5_VIZ → PETSc itself writes /viz/topology/cells natively, no Python-side MPI gather needed
• write_checkpoint keeps HDF5_PETSC (DMPlex restart format)
• Also hardens checkpoint_xdmf: uses integer topology precision from the HDF5 dtype rather than hardcoded Precision=\"8\", and removes the HyperSlab wrapper from variable DataItems (simpler/correct XDMF)

Roughly 10 lines of behavior change vs ~115 in PR #205.

Why we landed #205 anyway

• It was already submitted by a first-time contributor with working tests
• bknight1 confirmed it resolves his ParaView crashes
• The 22dd9ab branch had already been deleted from upstream; reviving it required separate work

What's needed to close this issue

1. Revive 22dd9ab2 as a fresh PR off current development (commit is still reachable via SHA on underworldcode/underworld3)
2. Verify it generates the same XDMF behaviour PR Fix XDMF output to explicitly write cell-to-vertex topology and ensure strict ParaView compatibility #205's tests exercise (test_0005_xdmf_compat.py::test_xdmf_viz_topology_written_correctly and _3d)
3. Run the parallel cases (PR Fix XDMF output to explicitly write cell-to-vertex topology and ensure strict ParaView compatibility #205 only had serial tests)
4. Once landed, _write_mesh_viz_groups and the _local_viz_cell_connectivity / _petsc_numbering_to_global_ids helpers from PR Fix XDMF output to explicitly write cell-to-vertex topology and ensure strict ParaView compatibility #205 can be removed — leaving only PETSc-managed viz topology

If the HDF5_VIZ viewer turns out not to handle some cell topology that the manual Python helper does, this issue can be closed as "alternative not viable" instead.

cc @gthyagi @bknight1

Underworld development team with AI support from Claude Code

[tiannh7]

Thanks for opening this issue and for keeping the PETSc-native route visible.

I agree that PETSc.Viewer.Format.HDF5_VIZ is probably the cleaner long-term direction if it can robustly provide the same /viz/geometry and /viz/topology semantics required by the XDMF output. In principle, letting PETSc/DMPlex manage the visualization topology is preferable to maintaining a Python-side MPI gather / vertex-deduplication / cell-reordering path in UW3.

One thing I think we should test explicitly is high-order / high-degree meshes and curved geometries. XDMF/ParaView compatibility can be fragile when the analysis mesh and visualization mesh are not identical, and some workflows may require a linearized visualization topology even if the underlying FE mesh is higher order.

It would be good to check the HDF5_VIZ path against serial/MPI output, 2D/3D meshes, simplex and quad/hex topologies, higher-degree or curved meshes where relevant, and consistency of node- and cell-centered field sizes in the generated XDMF.

If those cases work, I agree that replacing the manual _write_mesh_viz_groups path with the PETSc-managed output would be preferable. If not, perhaps the current implementation can remain as a fallback or serve as a regression-test baseline.

Happy to help test the revived branch against the ParaView-crash cases from #205.

[gthyagi]

@tiannh7 I'll look into the issue over the weekend and get back to you. I also have a number of uncommitted changes in my current working branch, which I'll clean up and push as well.

@lmoresi and I have already discussed the general direction for how we want to handle output writing, and I'll start working on that implementation too.

At the moment, higher-order mesh variables are working well with the current output-writing framework. Supporting higher-order mesh geometry (i.e., curved geometries) is a different challenge and would require substantial changes to the existing implementation. Given our current priorities, curved geometry support is not high on the roadmap.

That said, if you're interested in implementing it, we'd definitely be happy to help and provide guidance along the way.

[gthyagi]

I revived the PETSc-native alternative in draft PR #218: #218

The branch starts from current development, cherry-picks the two relevant commits (22dd9ab2 and 9f41035e as 90774778 and 47719277), then removes the manual _write_mesh_viz_groups path added in PR #205. The resulting output uses PETSc-native /viz/topology/cells and root /geometry/vertices, with XDMF pointing to both real datasets.

Validation run locally:

• .pixi/envs/amr-dev/bin/python -m pytest -q tests/test_0005_xdmf_compat.py
• .pixi/envs/amr-dev/bin/mpirun -n 2 .pixi/envs/amr-dev/bin/python -m pytest -q --with-mpi tests/parallel/test_0005_xdmf_viz_topology_mpi.py

[gthyagi]

Update: PR #218 has been pushed and its description has been expanded to cover both parts now on the branch:

• PETSc-native XDMF visualization topology: removes the manual _write_mesh_viz_groups() path and uses PETSc-written /viz/topology/cells with XDMF pointing to that dataset.
• Unified timestep/reload output: adds write_timestep(..., petsc_reload=True) so the same timestep variable HDF5 files can contain both XDMF/remap payloads and PETSc DMPlex section/vector reload payloads.

Validation completed locally:

• ./uw build
• tests/test_0003_save_load.py tests/test_0005_xdmf_compat.py: 16 passed
• MPI unified timestep reload test with 2 ranks: passed
• MPI XDMF viz-topology test with 2 ranks: passed in earlier validation
• Spherical Thieulot benchmark at cellsize=1/8, P2/P1, m=-1, 8 MPI ranks: solve completed, checkpoint-only metrics completed from unified files using read_checkpoint(), and serial comparison of read_checkpoint() vs read_timestep() from the same files gave zero velocity and pressure differences.

If we treat curved/isoparametric geometry as out of scope for this issue, then the original #217 request is addressed by PR #218. Curved geometry support should be tracked separately if/when it becomes a priority.