read_timestep mis-reads DG0 cell fields as a constant (spatial mean) for gmsh/units/multi-var checkpoints

[lmoresi]

Summary

MeshVariable.read_timestep() mis-reads DG0 (degree=0, continuous=False) cell fields from some checkpoints: the returned field is a near-constant ≈ the spatial mean instead of the spatially-varying per-cell values. The HDF5 file itself contains the correct data — only the read path is wrong.

Found while reloading the strain-rate-invariant (edot, DG0) from the Spiegelman plasticity benchmark for post-processing/rendering.

Evidence (benchmark edot field)

Raw HDF5 contents are correct, but read_timestep collapses to the mean:

# raw h5 (correct):
cell_fields/dotvarepsilon_dotvarepsilon  min=2.68e-4  max=1.97e-1   # varies ~700x
# read_timestep -> constant:
edot.data  min=1.15e-2  max=1.15e-2   (== spatial mean)

Reproduced from both a serial (np1) and a parallel (np4) checkpoint, so it is not a parallel-only issue.

Minimal reproducer that WORKS (for contrast)

A single DG0 variable on a UnstructuredSimplexBox, with and without refinement, round-trips correctly (READ std == WROTE std):

import numpy as np, underworld3 as uw, os
D = "/tmp/dg0_repro"; os.makedirs(D, exist_ok=True)
mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0,0), maxCoords=(1,1),
                                         cellSize=0.2, refinement=1)
v = uw.discretisation.MeshVariable("phi", mesh, 1, degree=0, continuous=False)
v.data[:,0] = np.sin(3*v.coords[:,0]) + v.coords[:,1]
print("WROTE std", v.data[:,0].std())                    # 0.40
mesh.petsc_save_checkpoint(index=0, meshVars=[v], outputPath=D+"/")
m2 = uw.discretisation.Mesh(f"{D}/checkpoint.mesh.00000.h5")
w = uw.discretisation.MeshVariable("phi", m2, 1, degree=0, continuous=False)
w.read_timestep("checkpoint", "phi", 0, outputPath=D+"/")
print("READ  std", w.data[:,0].std())                    # 0.40 (OK)

So the bug is not universal to DG0 — it is triggered by something in the failing context that the minimal case lacks.

What differs in the failing (benchmark) context

The benchmark checkpoint differs from the working minimal case in several ways; I was not able to isolate the single trigger:

• mesh imported via discretisation._from_gmsh(...) (not a uw.meshing.* constructor);
• a uw.Model scaling is active (dimensional/unit-aware coordinates);
• a multi-variable checkpoint with mixed degrees — meshVars=[edot(DG0), visc(DG0), stress(DG0), p_soln(P1), v_soln(P2)].

The "read == spatial mean" symptom strongly suggests a continuous/averaging projection code path is being taken when the DG0 field is read back (or a cell→point→cell remap that averages), rather than a direct cell-data load — but only for these meshes.

Impact

DG0 fields (material index, strain-rate invariant, viscosity, stress) cannot be reliably reloaded from these checkpoints for restart or post-processing. Workaround: read cell_fields/<name>_<name> directly from the HDF5, or render in-process from the live variable.

Repro scripts

~/+Simulations/jac_unwrap_probe/ style probes and the Spiegelman benchmark
(docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py) demonstrate it.

---

Reported via Claude Code while validating the yield-homotopy / FMG work on the Spiegelman benchmark.