ADC-376 export hyqmom15 campaign to ParaView (.vti + .pvd)

hyqmom15/campaigns/README.md

@@ -44,8 +44,9 @@ threads, wall-clock, dt range, AMR flag, commits, host), and feeds a
 runtime, status). Post-process into exploitable artefacts:
 
 ```bash
-python3 hyqmom15/campaigns/export_h5.py    out/campaign   # -> out/campaign/h5/<case>.h5
-python3 hyqmom15/campaigns/make_rapport.py out/campaign   # -> out/campaign/rapport.md
+python3 hyqmom15/campaigns/export_h5.py     out/campaign  # -> out/campaign/h5/<case>.h5
+python3 hyqmom15/campaigns/make_rapport.py  out/campaign  # -> out/campaign/rapport.md
+python3 hyqmom15/campaigns/to_paraview.py   out/campaign  # -> out/campaign/paraview/<case>.pvd
 python3 hyqmom15/plots/diagnostics_plots.py out/campaign  # -> out/campaign/figures/ (matplotlib)
 ```
 
@@ -56,6 +57,16 @@ realizability verdict, conservation, positivity, symmetry, figure list, HDF5 ref
 plus a synthesis table and the optional speedup. Figures come from `hyqmom15/plots`
 (ADC-377/384).
 
+**ParaView** comes two ways. The run already writes the **native adc_cpp VTK output**
+(`adc.System.write(format="vtk")`): one `step_NNNNNN.vti` (ImageData, CellData
+`mom_<moment>` + `phi`) per snapshot in each case dir, opened directly by ParaView /
+VisIt -- open the `step_*.vti` series for the animation. `to_paraview.py` is an
+optional **enriched** export: from the npz it adds physics-ready cell fields (density,
+ux/uy + speed, the realizability margin `lam_min`) and a `<case>.pvd` time collection
+(real `t` values) under `paraview/`. The runs are single-rank, so each snapshot is one
+grid; parallel pieces (`.pvti`) would only apply to an MPI-distributed run. The enriched
+exporter has no VTK/pyvista dependency (hand-written VTK XML).
+
 ## Matlab speedup baseline (Octave, run where the Matlab source lives)
 
 ```bash

hyqmom15/campaigns/check_campaign.py

@@ -26,6 +26,7 @@
 import octave_matlab as om  # noqa: E402
 import romeo_rie_mom2d as campaign  # noqa: E402
 import synthesis  # noqa: E402
+import to_paraview  # noqa: E402
 
 _MOMENTS = ["M00", "M10", "M20", "M30", "M40", "M01", "M11", "M21", "M31",
             "M02", "M12", "M22", "M03", "M13", "M04"]
@@ -137,8 +138,22 @@ def check_dry_run():
     return "dry-run OK (3 cases -> run_meta.json each + synthesis.md, no adc)"
 
 
+def check_paraview():
+    with tempfile.TemporaryDirectory() as d:
+        root = pathlib.Path(d)
+        _synth_campaign(root, cases=("dicotron",))
+        pvd = to_paraview.export_case(root / "dicotron", root / "paraview")
+        assert pvd and pvd.exists() and pvd.suffix == ".pvd", "no .pvd written"
+        vti = sorted((root / "paraview").glob("dicotron_*.vti"))
+        assert vti, "no .vti written"
+        head = vti[0].read_text()
+        assert 'type="ImageData"' in head and 'Name="density"' in head, "malformed .vti"
+        assert pvd.read_text().count("<DataSet") == len(vti), "pvd/vti count mismatch"
+    return "to_paraview OK (.pvd time series + ImageData .vti per snapshot)"
+
+
 CHECKS = [check_tables, check_run_meta, check_octave_source, check_dry_run,
-          check_make_rapport, check_export_h5]
+          check_make_rapport, check_export_h5, check_paraview]
 
 
 def main() -> int:

hyqmom15/campaigns/export_h5.py

@@ -89,7 +89,7 @@ def main(argv=None) -> int:
     out_dir = pathlib.Path(args.out) if args.out else root / "h5"
     out_dir.mkdir(parents=True, exist_ok=True)
     case_dirs = [root / args.case] if args.case else sorted(
-        d for d in root.iterdir() if d.is_dir() and d.name not in ("figures", "h5"))
+        d for d in root.iterdir() if d.is_dir() and d.name not in ("figures", "h5", "paraview"))
     n = 0
     for cd in case_dirs:
         path = export_case(cd, out_dir / ("%s.h5" % cd.name))

hyqmom15/campaigns/make_rapport.py

@@ -103,7 +103,7 @@ def build_rapport(campaign_dir, matlab_times=None):
     root = pathlib.Path(campaign_dir)
     figures_dir, h5_dir = root / "figures", root / "h5"
     case_dirs = sorted(d for d in root.iterdir()
-                       if d.is_dir() and d.name not in ("figures", "h5"))
+                       if d.is_dir() and d.name not in ("figures", "h5", "paraview"))
     lines = ["# hyqmom15 ROMEO campaign report", "",
              "Per-case detailed analysis (config, realizability, conservation, "
              "symmetry, figures, HDF5 export). D/Dmax is a clarified convention, "

hyqmom15/campaigns/romeo_rie_mom2d.py

@@ -85,6 +85,9 @@ def _run_and_record(sim, case, case_dir, n_snapshots, dt_fn, smoke):
     while True:
         if t >= next_snap_t - 1e-15 or t >= tmax:
             sim.write(str(case_dir / "step"), format="npz", step=snap_k)
+            # native ParaView output (adc_cpp System.write VTK): one ImageData .vti per
+            # snapshot, CellData = mom_<moment> + phi, opened directly by ParaView/VisIt.
+            sim.write(str(case_dir / "step"), format="vtk", step=snap_k)
             u = np.array(sim.get_state("mom"))
             realiz.append(summarize(*field_realizability(u)))
             try:
@@ -118,7 +121,7 @@ def _run_and_record(sim, case, case_dir, n_snapshots, dt_fn, smoke):
     }
 
 
-def _run_diocotron(n, case_dir, n_snapshots, smoke):
+def _run_diocotron(n, case_dir, n_snapshots, smoke, projection=False):
     import run_diocotron_periodic as drv
     from matlab_ref import compute_dt
     case = drv.periodic_case(n)
@@ -128,11 +131,13 @@ def _run_diocotron(n, case_dir, n_snapshots, smoke):
     probe.set_state("mom", u0)
     probe.solve_fields()
     vmax = case.cfl * case.dx / probe.step_cfl(case.cfl)
-    sim = drv.build_periodic_sim(n, rho_bg=rho_bg)
+    sim = drv.build_periodic_sim(n, rho_bg=rho_bg, projection=projection)
     sim.set_state("mom", u0)
     sim.solve_fields()
     rec = _run_and_record(sim, case, case_dir, n_snapshots, lambda t: compute_dt(vmax, case, t), smoke)
-    return case, _solver_config(case, "hll"), rec
+    solver = _solver_config(case, "hll")
+    solver["projection"] = projection  # native relaxation15 projector (build_projection)
+    return case, solver, rec
 
 
 def _run_wave(driver, case_builder, ic, build, n, case_dir, n_snapshots, smoke):
@@ -187,10 +192,14 @@ def _run_constant(n, case_dir, n_snapshots, smoke):
     return case, _solver_config(case, "hll"), rec
 
 
-def run_one(case_name, n, case_dir, n_snapshots, smoke):
-    """Run a single case end-to-end; returns (case, solver_config, record)."""
+def run_one(case_name, n, case_dir, n_snapshots, smoke, projection=False):
+    """Run a single case end-to-end; returns (case, solver_config, record).
+
+    projection enables the native relaxation15 realizability projector; wired for
+    dicotron (the case that loses realizability without it), ignored otherwise.
+    """
     if case_name == "dicotron":
-        return _run_diocotron(n, case_dir, n_snapshots, smoke)
+        return _run_diocotron(n, case_dir, n_snapshots, smoke, projection)
     if case_name == "fluid_wave":
         return _run_fluid(n, case_dir, n_snapshots, smoke)
     if case_name == "constant":
@@ -214,7 +223,7 @@ def _provenance(case_name, params, solver, rec, threads, host, ts):
         host=host, timestamp=ts)
 
 
-def run_campaign(out_dir, cases, smoke, dry_run, threads, matlab_times=None):
+def run_campaign(out_dir, cases, smoke, dry_run, threads, matlab_times=None, projection=False):
     """Run the campaign; write per-case snapshots + run_meta.json and a synthesis table."""
     import datetime
     import socket
@@ -238,7 +247,7 @@ def run_campaign(out_dir, cases, smoke, dry_run, threads, matlab_times=None):
             params, solver, rec = _params(case), _solver_config(case, "hll"), {
                 "Np": n, "n_steps": 0, "status": "dry-run", "wall_clock_s": None}
         else:
-            case, solver, rec = run_one(name, n, case_dir, DEFAULT_SNAPSHOTS, smoke)
+            case, solver, rec = run_one(name, n, case_dir, DEFAULT_SNAPSHOTS, smoke, projection)
             params = _params(case)
             adc_times[name] = rec.get("wall_clock_s")
         meta = _provenance(name, params, solver, rec, threads, host, ts)
@@ -265,6 +274,8 @@ def main(argv=None) -> int:
     grp.add_argument("--smoke", action="store_true", help="reduced Np + capped steps")
     grp.add_argument("--full", action="store_true", help="full Matlab resolution")
     p.add_argument("--dry-run", action="store_true", help="no adc run; structure + meta only")
+    p.add_argument("--projection", action="store_true",
+                   help="enable the native relaxation15 realizability projector (wired for dicotron)")
     args = p.parse_args(argv)
 
     cases = args.cases.split(",") if args.cases else list(CASE_NP)
@@ -276,7 +287,8 @@ def main(argv=None) -> int:
         p.error("specify --smoke, --full, or --dry-run (no silent default)")
     smoke = args.smoke  # --full -> full resolution; --dry-run ignores it (no sim)
     run_campaign(args.out, cases, smoke=smoke, dry_run=args.dry_run,
-                 threads=args.threads, matlab_times=args.matlab_times)
+                 threads=args.threads, matlab_times=args.matlab_times,
+                 projection=args.projection)
     return 0
 
 

hyqmom15/campaigns/romeo_rie_mom2d.sbatch

@@ -49,4 +49,6 @@ $PY hyqmom15/campaigns/romeo_rie_mom2d.py --full --out "$OUT" --threads "$THREAD
 $PY hyqmom15/campaigns/export_h5.py "$OUT"
 # 3) per-case analysis report (config + realizability + figures + h5)
 $PY hyqmom15/campaigns/make_rapport.py "$OUT"
-echo "campaign + HDF5 + rapport.md written under: $OUT"
+# 4) ParaView time series per case (<case>.pvd + .vti) -- open the .pvd in ParaView
+$PY hyqmom15/campaigns/to_paraview.py "$OUT"
+echo "campaign + HDF5 + rapport.md + ParaView (.pvd) written under: $OUT"

hyqmom15/campaigns/to_paraview.py

@@ -0,0 +1,153 @@
+#!/usr/bin/env python3
+"""Enriched ParaView export for a hyqmom15 campaign case (ADC-376).
+
+The run already writes the native ``adc.System.write(format="vtk")`` output (one
+``.vti`` per snapshot, raw moments + phi). This is the ENRICHED complement, built
+from the npz snapshots on the uniform Cartesian grid [-0.5, 0.5]^2: one
+``<case>_NNNN.vti`` (ImageData, cell data) per snapshot plus a ``<case>.pvd``
+collection that ParaView opens as an animated time series with real ``t`` values.
+Fields: density (M00), velocities ux/uy + speed, the potential phi (when present),
+the realizability margin lam_min, and the 15 raw moments.
+
+The runs are single-rank (no MPI), so each snapshot is one grid -> a single .vti
+per step is correct; parallel pieces (.pvti) would only apply to MPI-distributed
+runs. Pure Python (hand-written VTK XML, base64 binary); no VTK/pyvista needed.
+
+Usage:
+    python3 hyqmom15/campaigns/to_paraview.py <campaign_dir> [--case NAME] [--out DIR]
+"""
+from __future__ import annotations
+
+import argparse
+import base64
+import pathlib
+import struct
+import sys
+
+import numpy as np
+
+HERE = pathlib.Path(__file__).resolve().parent
+sys.path.insert(0, str(HERE))
+sys.path.insert(0, str(HERE.parent))
+sys.path.insert(0, str(HERE.parent / "plots"))
+
+from snapshots import load_case  # noqa: E402
+
+try:
+    from diagnostics import field_realizability  # noqa: E402
+    _HAVE_REALIZ = True
+except Exception:  # noqa: BLE001
+    _HAVE_REALIZ = False
+
+
+def _b64(arr):
+    """Encode a numpy array as a VTK inline-binary payload (UInt32 length + LE f8)."""
+    data = np.ascontiguousarray(arr, dtype="<f8").tobytes()
+    return base64.b64encode(struct.pack("<I", len(data)) + data).decode("ascii")
+
+
+def _fields(snap):
+    """Return an ordered dict of cell-data field name -> flat (ny*nx,) array."""
+    m = snap.moments  # (15, ny, nx)
+    names = list(snap.names)
+    idx = {n: k for k, n in enumerate(names)}
+    rho = m[idx.get("M00", 0)]
+    ok = np.abs(rho) > 1e-30  # vacuum cells read as zero velocity, not the raw moment
+    safe = np.where(ok, rho, 1.0)
+    ux = np.where(ok, m[idx["M10"]] / safe, 0.0) if "M10" in idx else np.zeros_like(rho)
+    uy = np.where(ok, m[idx["M01"]] / safe, 0.0) if "M01" in idx else np.zeros_like(rho)
+    fields = {
+        "density": rho, "ux": ux, "uy": uy, "speed": np.hypot(ux, uy),
+    }
+    if snap.phi is not None:
+        fields["phi"] = snap.phi
+    if _HAVE_REALIZ:
+        try:
+            fields["lam_min"] = field_realizability(m)[0]["lam_min"]
+        except Exception:  # noqa: BLE001
+            pass
+    for n in names:
+        fields[n] = m[idx[n]]
+    return {k: np.asarray(v, dtype=float).ravel() for k, v in fields.items()}
+
+
+def write_vti(path, ny, nx, fields):
+    """Write one ImageData .vti (cell data) on the [-0.5,0.5]^2 uniform grid."""
+    dx, dy = 1.0 / nx, 1.0 / ny
+    head = (
+        '<?xml version="1.0"?>\n'
+        '<VTKFile type="ImageData" version="1.0" byte_order="LittleEndian" '
+        'header_type="UInt32">\n'
+        '  <ImageData WholeExtent="0 %d 0 %d 0 0" Origin="-0.5 -0.5 0" '
+        'Spacing="%.17g %.17g 1">\n'
+        '    <Piece Extent="0 %d 0 %d 0 0">\n'
+        '      <CellData Scalars="density">\n' % (nx, ny, dx, dy, nx, ny)
+    )
+    body = "".join(
+        '        <DataArray type="Float64" Name="%s" format="binary">%s</DataArray>\n'
+        % (name, _b64(arr)) for name, arr in fields.items()
+    )
+    tail = "      </CellData>\n    </Piece>\n  </ImageData>\n</VTKFile>\n"
+    pathlib.Path(path).write_text(head + body + tail, encoding="ascii")
+
+
+def write_pvd(path, entries):
+    """Write a ParaView .pvd collection: (timestep, relative .vti file) pairs."""
+    rows = "".join(
+        '    <DataSet timestep="%.17g" group="" part="0" file="%s"/>\n' % (t, f)
+        for t, f in entries
+    )
+    pathlib.Path(path).write_text(
+        '<?xml version="1.0"?>\n'
+        '<VTKFile type="Collection" version="1.0" byte_order="LittleEndian">\n'
+        "  <Collection>\n%s  </Collection>\n</VTKFile>\n" % rows,
+        encoding="ascii")
+
+
+def export_case(case_dir, out_dir):
+    """Convert one case directory to ``<case>.pvd`` + ``<case>_NNNN.vti``; return the .pvd path."""
+    snaps = load_case(case_dir)
+    if not snaps:
+        return None
+    case = pathlib.Path(case_dir).name
+    out_dir = pathlib.Path(out_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    entries = []
+    for k, s in enumerate(snaps):
+        ny, nx = s.density.shape
+        fname = "%s_%04d.vti" % (case, k)
+        write_vti(out_dir / fname, ny, nx, _fields(s))
+        entries.append((float(s.t), fname))
+    pvd = out_dir / ("%s.pvd" % case)
+    write_pvd(pvd, entries)
+    return pvd
+
+
+def main(argv=None) -> int:
+    p = argparse.ArgumentParser(description="Export hyqmom15 campaign cases to ParaView (.vti + .pvd).")
+    p.add_argument("campaign_dir", help="root with one snapshot sub-directory per case")
+    p.add_argument("--case", help="only this case (default: all)")
+    p.add_argument("--out", help="output dir (default: <campaign_dir>/paraview)")
+    args = p.parse_args(argv)
+    root = pathlib.Path(args.campaign_dir)
+    if not root.is_dir():
+        print("campaign dir not found: %s" % root, file=sys.stderr)
+        return 1
+    out_dir = pathlib.Path(args.out) if args.out else root / "paraview"
+    case_dirs = [root / args.case] if args.case else sorted(
+        d for d in root.iterdir() if d.is_dir() and d.name not in ("figures", "h5", "paraview"))
+    n = 0
+    for cd in case_dirs:
+        pvd = export_case(cd, out_dir)
+        if pvd:
+            n += 1
+            print("  %-20s -> %s" % (cd.name, pvd))
+    if args.case and n == 0:
+        print("no snapshots found for case %r under %s" % (args.case, root), file=sys.stderr)
+        return 1
+    print("wrote %d ParaView collections to %s" % (n, out_dir))
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

hyqmom15/plots/diagnostics_plots.py

@@ -190,7 +190,8 @@ def main(argv=None) -> int:
     if args.case:
         case_dirs = [root / args.case]
     else:
-        case_dirs = sorted(d for d in root.iterdir() if d.is_dir() and d.name != "figures")
+        case_dirs = sorted(d for d in root.iterdir()
+                           if d.is_dir() and d.name not in ("figures", "h5", "paraview"))
 
     total = 0
     for cd in case_dirs:

hyqmom15/plots/plot_rie_mom2d_case.py

@@ -196,7 +196,8 @@ def main(argv=None) -> int:
     if args.case:
         case_dirs = [root / args.case]
     else:
-        case_dirs = sorted(d for d in root.iterdir() if d.is_dir() and d.name != "figures")
+        case_dirs = sorted(d for d in root.iterdir()
+                           if d.is_dir() and d.name not in ("figures", "h5", "paraview"))
 
     total = 0
     for cd in case_dirs:

hyqmom15/runs/run_diocotron_periodic.py

@@ -72,7 +72,8 @@ def diocotron_ic(n: int, orientation: str = "standard") -> np.ndarray:
     return init_diocotron_field(periodic_case(n), orientation=orientation).M
 
 
-def build_periodic_sim(n: int, rho_bg: float, name: str = "mom") -> "adc.System":
+def build_periodic_sim(n: int, rho_bg: float, name: str = "mom",
+                       projection: bool = False) -> "adc.System":
     """System periodique fidele au nouveau Matlab : E+B sources, HLL exact, Euler.
 
     omega_c=-20 active la source magnetique (D1). exact_speeds=True donne les
@@ -93,9 +94,11 @@ def build_periodic_sim(n: int, rho_bg: float, name: str = "mom") -> "adc.System"
         rho_background=rho_bg,
         omega_p=CASE.omega_p,                 # borne dt source
         exact_speeds=True,
+        projection=projection,                # True -> projecteur natif relaxation15 (ADC-275)
     )
+    so_name = "hyqmom15_vp_periodic%s.so" % ("_proj" if projection else "")
     compiled = m.compile(
-        os.path.join(case_output_dir("hyqmom15"), "hyqmom15_vp_periodic.so"),
+        os.path.join(case_output_dir("hyqmom15"), so_name),
         adc_include(),
         backend="production",                 # requis pour Euler (cf. ADC-356)
     )