diff --git a/experiments/audit_shape_field_weak.py b/experiments/audit_shape_field_weak.py
new file mode 100644
index 0000000..956a756
--- /dev/null
+++ b/experiments/audit_shape_field_weak.py
@@ -0,0 +1,799 @@
+#!/usr/bin/env python3
+"""Audit frozen shape/field artifacts without reevaluating FEM.
+
+Use the NumPy-enabled Python environment. The SVD is taken directly on the tall
+matrix of saved block actions: squaring its condition number through a Gram
+matrix would obscure precisely the weak response this experiment investigates.
+Physical observations use broken scalar gradients, not RT transport or DG jump
+distributions. Neither these observations nor a small restricted singular value
+identify a physical gauge or establish the rank of the complete Jacobian.
+"""
+import argparse
+import csv
+import hashlib
+import json
+import math
+import os
+from pathlib import Path
+import sys
+
+# Avoid competing with a later numerical campaign through a threaded BLAS pool.
+for _name in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
+    os.environ[_name] = "1"
+import numpy as np
+
+P = "weak_shape_"
+PP = P + "physical_"
+EPS = np.finfo(float).eps
+RTOL = 5e-10
+FIELDS = ("density", "specific_enthalpy", "gravity_potential")
+INTEGRALS = ("base_integral", "pullback_integral", "transport_integral", "eulerian_integral",
+             "base_squared_integral", "pullback_squared_integral", "transport_squared_integral",
+             "eulerian_squared_integral", "pullback_transport_cross_integral")
+SQUARES = INTEGRALS[4:8]
+
+
+def digest(path):
+    h = hashlib.sha256()
+    with Path(path).open("rb") as stream:
+        for chunk in iter(lambda: stream.read(1024 * 1024), b""):
+            h.update(chunk)
+    return h.hexdigest()
+
+
+def rows(path):
+    with Path(path).open(newline="") as stream:
+        reader = csv.DictReader(stream)
+        if not reader.fieldnames or len(set(reader.fieldnames)) != len(reader.fieldnames):
+            raise ValueError(f"Missing or duplicate CSV columns: {path}")
+        result = list(reader)
+    if any(None in row or any(v is None for v in row.values()) for row in result):
+        raise ValueError(f"Malformed CSV row: {path}")
+    return result
+
+
+def finite(value):
+    number = float(value)
+    if not math.isfinite(number):
+        raise ValueError(f"Nonfinite available value: {value}")
+    return number
+
+
+def integer(value):
+    number = finite(value)
+    if number != int(number):
+        raise ValueError(f"Nonintegral count/index: {value}")
+    return int(number)
+
+
+def vector(path):
+    with Path(path).open() as stream:
+        size = int(stream.readline().strip())
+        data = np.loadtxt(stream, dtype=float, ndmin=1)
+    if size <= 0 or data.ndim != 1 or data.size != size or not np.isfinite(data).all():
+        raise ValueError(f"Invalid finite vector: {path}")
+    return data
+
+
+def norm(value):
+    return float(np.linalg.norm(value))
+
+
+def checked_einsum(subscripts, *operands):
+    # NumPy 2.2.6 / Accelerate can flag invalid/divide/overflow for a bounded
+    # 257x6 matrix product whose answer agrees with scalar fsum to 4e-15.
+    # Direct contractions avoid that BLAS path; finite and reconstruction
+    # checks remain mandatory, with no warning suppression or relaxed bounds.
+    if any(not np.isfinite(operand).all() for operand in operands):
+        raise ValueError("Nonfinite contraction input")
+    result = np.einsum(subscripts, *operands, optimize=False)
+    if not np.isfinite(result).all():
+        raise ValueError("Nonfinite contraction result")
+    return result
+
+
+def near(a, b, scale=0.):
+    return math.isfinite(a) and math.isfinite(b) and abs(a - b) <= (
+        RTOL * max(abs(a), abs(b)) + 2048 * EPS * abs(scale) + 1e-28)
+
+
+def vector_near(a, b, scale=0.):
+    return a.shape == b.shape and near(norm(a - b), 0., scale + norm(a) + norm(b))
+
+
+def ratio(a, b):
+    return finite(a / b) if b > 0 else None
+
+
+class Audit:
+    def __init__(self):
+        self.count = 0
+        self.failures = []
+
+    def check(self, name, condition):
+        self.count += 1
+        if not bool(condition):
+            self.failures.append(name)
+
+    def close(self, name, actual, expected, scale=0.):
+        self.check(name, near(actual, expected, scale))
+
+
+def table(path, key):
+    result = {}
+    for row in rows(path):
+        index = key(row)
+        if index in result:
+            raise ValueError(f"Duplicate row {index}: {path}")
+        result[index] = row
+    return result
+
+
+def layout(path, size):
+    result = {}
+    all_rows = rows(path)
+    if any(r["kind"] not in ("value", "residual") for r in all_rows):
+        raise ValueError("Unknown layout kind")
+    for kind in ("value", "residual"):
+        entries = [r for r in all_rows if r["kind"] == kind]
+        offset, names, blocks = 0, set(), []
+        for row in entries:
+            start, count = integer(row["offset"]), integer(row["size"])
+            if start != offset or count <= 0 or row["block"] in names:
+                raise ValueError(f"Invalid {kind} layout partition")
+            names.add(row["block"])
+            blocks.append((row["block"], start, count))
+            offset += count
+        if offset != size:
+            raise ValueError(f"Incomplete {kind} layout partition")
+        result[kind] = blocks
+    return result
+
+
+def restricted_svd(direction, actions, blocks, audit):
+    norms = np.array([norm(direction[o:o+n]) for _, o, n in blocks])
+    active = [i for i, x in enumerate(norms) if x > 0.]
+    for i, x in enumerate(norms):
+        if x == 0.:
+            audit.check(f"zero_direction_block_action:{i}", norm(actions[i]) == 0.)
+    if not active:
+        raise ValueError("No nonzero block directions")
+    matrix = np.column_stack([actions[i] / norms[i] for i in active])
+    if not np.isfinite(matrix).all():
+        raise ValueError("Nonfinite scaled block action")
+    u, singular, vh = np.linalg.svd(matrix, full_matrices=False)
+    weights = vh[-1].copy()
+    # Fix sign for reproducible reporting; repeated singular values still permit
+    # different equally valid vectors, which the JSON explicitly acknowledges.
+    if weights[np.argmax(np.abs(weights))] < 0:
+        weights *= -1
+    combination = np.zeros_like(direction)
+    action = np.zeros_like(direction)
+    coefficients = {}
+    for weight, i in zip(weights, active):
+        name, start, count = blocks[i]
+        coefficient = float(weight / norms[i])
+        combination[start:start+count] = coefficient * direction[start:start+count]
+        action += coefficient * actions[i]
+        coefficients[name] = {"orthonormal_basis_weight": float(weight),
+                              "saved_block_direction_multiplier": coefficient}
+    scale = float(singular[0])
+    audit.close("restricted_svd:unit_direction", norm(combination), 1.)
+    audit.close("restricted_svd:smallest_action_norm", norm(action), float(singular[-1]), scale)
+    audit.check("restricted_svd:action_reconstruction", vector_near(action, checked_einsum("ij,j->i",matrix,weights), scale))
+    audit.check("restricted_svd:factorization", vector_near(matrix, checked_einsum("ij,jk->ik",u*singular,vh), norm(matrix)))
+    audit.check("restricted_svd:finite", np.isfinite(singular).all() and np.isfinite(vh).all())
+    return {
+        "scope": "Jacobian restricted to nonzero normalized saved block directions; not full Jacobian rank or nullspace",
+        "basis": "disjoint block directions divided by their normalized coefficient two-norm",
+        "active_blocks": [blocks[i][0] for i in active], "singular_values": singular.tolist(),
+        "condition_ratio_in_restricted_span": ratio(float(singular[0]), float(singular[-1])),
+        "smallest_direction_norm": norm(combination), "smallest_action_norm": norm(action),
+        "coefficients": coefficients,
+        "reconstructed_direction_float64_sha256": hashlib.sha256(combination.astype("<f8").tobytes()).hexdigest(),
+        "reconstructed_action_float64_sha256": hashlib.sha256(action.astype("<f8").tobytes()).hexdigest(),
+        "smallest_singular_value_gap": float(singular[-2] - singular[-1]) if len(singular) > 1 else None,
+        "qualification": "basis metric is coefficient normalization; vector sign fixed, degenerate subspace basis not unique"
+    }, action
+
+
+def surface_fit(surface, physical, physical_direction, blocks, audit):
+    shape = next(b for b in blocks if b[0] == "surface_deformation")
+    _, offset, size = shape
+    audit.check("surface:complete_parameter_indices", sorted(surface) == list(range(size)))
+    ordered = [surface[i] for i in range(size)]
+    audit.check("surface:unique_volume_dofs", len({integer(r["volume_scalar_true_dof"]) for r in ordered}) == size)
+    audit.check("surface:nonnegative_volume_dofs", all(integer(r["volume_scalar_true_dof"]) >= 0 for r in ordered))
+    xyz = np.array([[finite(r["unit_radial_" + axis]) for axis in "xyz"] for r in ordered])
+    data = np.array([finite(r["physical_direction_parameter"]) for r in ordered])
+    for i, row in enumerate(ordered):
+        r0, r = finite(row["reference_radius"]), finite(row["physical_radius"])
+        b = finite(row["base_parameter"])
+        audit.check(f"surface:radii:{i}", r0 > 0 and r > 0)
+        audit.close(f"surface:unit_radial:{i}", norm(xyz[i]), 1.)
+        audit.close(f"surface:radius_identity:{i}", r, r0 + b, abs(r0) + abs(b))
+        audit.close(f"surface:base:{i}", b, float(physical[offset+i]))
+        audit.close(f"surface:direction:{i}", data[i], float(physical_direction[offset+i]))
+    x, y, z = xyz.T
+    columns = [np.ones(size), x, y, z, x*y, x*z, y*z, x*x-y*y, 3*z*z-1]
+    fits = []
+    for degree, count in ((0, 1), (1, 4), (2, 9)):
+        matrix = np.column_stack(columns[:count])
+        coef, _, rank, singular = np.linalg.lstsq(matrix, data, rcond=None)
+        residual = data - checked_einsum("ij,j->i",matrix,coef)
+        audit.check(f"surface:fit_finite:{degree}", np.isfinite(coef).all() and np.isfinite(singular).all())
+        fits.append({"maximum_degree": degree, "rank": int(rank), "basis_coefficients": coef.tolist(),
+                     "relative_nodal_l2_remainder": ratio(norm(residual), norm(data)),
+                     "nodal_l2_remainder": norm(residual), "singular_values": singular.tolist()})
+    return {"scope": "unweighted nodal least squares in real l=0,1,2 polynomial spans; not area-weighted harmonic coefficients or gauge identification",
+            "basis": ["1", "nx", "ny", "nz", "nx*ny", "nx*nz", "ny*nz", "nx^2-ny^2", "3*nz^2-1"],
+            "direction_nodal_l2": norm(data), "direction_maximum_absolute": float(np.max(np.abs(data))), "fits": fits}
+
+
+def physical_analysis(folder, physical, physical_direction, blocks, audit, required):
+    def read(suffix, key):
+        name = PP + suffix + ".csv"
+        required.add(name)
+        return table(folder / name, key)
+    metadata = read("metadata", lambda r: r["key"])
+    meta = {k: r["value"] for k, r in metadata.items()}
+    audit.check("physical:schema", meta["schema"] == "weak_shape_physical_v1")
+    audit.check("physical:declared_scalar_scope", "broken" in meta["dg_derivatives"] and "not evaluated" in meta["gravity_vector"])
+    elem = read("elements", lambda r: (integer(r["quadrature"]), integer(r["element"]), r["field"]))
+    geometry = read("geometry", lambda r: (integer(r["quadrature"]), integer(r["element"])))
+    affine = read("affine", lambda r: (integer(r["quadrature"]), integer(r["element"])))
+    surface = read("surface_nodes", lambda r: integer(r["parameter"]))
+    summary = read("summary", lambda r: (integer(r["quadrature"]), r["field"], r["metric"]))
+    observed = {key: finite(r["value"]) for key, r in summary.items()}
+    expected_summary = set()
+    def compare(q, field, metric, expected, scale=0.):
+        key = (q, field, metric)
+        expected_summary.add(key)
+        audit.close("physical:summary:" + str(key), observed[key], expected, scale)
+    audit.check("physical:quadrature_coverage", {q for q, _ in geometry} == {10, 14})
+    ids = {e for q, e in geometry if q == 10}
+    audit.check("physical:same_elements_at_both_rules", ids and ids == {e for q, e in geometry if q == 14})
+    audit.check("physical:affine_coverage", set(affine) == set(geometry))
+    audit.check("physical:field_coverage", set(elem) == {(q, e, f) for q, e in geometry for f in FIELDS})
+    outputs = []
+    for q in (10, 14):
+        geo = [geometry[(q, e)] for e in sorted(ids)]
+        aff = [affine[(q, e)] for e in sorted(ids)]
+        for e, g, a in zip(sorted(ids), geo, aff):
+            for key in g:
+                if key not in ("quadrature", "element"):
+                    finite(g[key])
+            volume = finite(g["physical_volume"])
+            audit.check(f"physical:geometry_positive:{q}:{e}", volume > 0 and integer(g["points"]) > 0)
+            for prefix in ("reference", "physical"):
+                lo, hi = finite(g[prefix+"_radius_min"]), finite(g[prefix+"_radius_max"])
+                audit.check(f"physical:radius_range:{q}:{e}:{prefix}", 0 <= lo <= hi)
+            for key in ("base_displacement_squared_integral", "direction_displacement_squared_integral"):
+                audit.check(f"physical:geometry_square:{q}:{e}:{key}", finite(g[key]) >= 0)
+            for key, family in (("density_element_order", "density_family_order"),
+                                ("enthalpy_element_order", "enthalpy_family_order"),
+                                ("potential_element_order", "potential_family_order"),
+                                ("displacement_element_order", "displacement_family_order")):
+                audit.check(f"physical:order:{q}:{e}:{key}", integer(g[key]) == integer(meta[family]))
+            for key in ("reference_volume", "value_error_squared_integral", "reference_value_squared_integral",
+                        "gradient_error_squared_integral", "reference_gradient_squared_integral"):
+                audit.check(f"physical:affine_nonnegative:{q}:{e}:{key}", finite(a[key]) >= 0)
+            audit.check(f"physical:affine_volume:{q}:{e}", finite(a["reference_volume"]) > 0)
+            audit.check(f"physical:affine_points:{q}:{e}", integer(a["points"]) == integer(g["points"]))
+            audit.close(f"physical:affine_identity:{q}:{e}", finite(a["reference_gradient_squared_integral"]), 3*finite(a["reference_volume"]))
+        totals_g = {k: math.fsum(finite(g[k]) for g in geo) for k in
+                    ("physical_volume", "base_displacement_squared_integral", "direction_displacement_squared_integral")}
+        for key, value in totals_g.items():
+            compare(q, "geometry", key, value)
+        compare(q, "geometry", "stellar_elements", len(ids))
+        totals_a = {k: math.fsum(finite(a[k]) for a in aff) for k in
+                    ("reference_volume", "value_error_squared_integral", "reference_value_squared_integral",
+                     "gradient_error_squared_integral", "reference_gradient_squared_integral")}
+        for key, value in totals_a.items():
+            compare(q, "affine", key, value)
+        for label, num, den in (("relative_value_l2_error", "value_error_squared_integral", "reference_value_squared_integral"),
+                                ("relative_gradient_l2_error", "gradient_error_squared_integral", "reference_gradient_squared_integral")):
+            compare(q, "affine", label, math.sqrt(totals_a[num] / totals_a[den]))
+        field_outputs = {}
+        for field in FIELDS:
+            records = [elem[(q, e, field)] for e in sorted(ids)]
+            for e, row, g in zip(sorted(ids), records, geo):
+                values = {k: finite(row[k]) for k in INTEGRALS}
+                audit.check(f"physical:gradient_kind:{q}:{e}:{field}", row["gradient_kind"] == ("continuous_H1" if field == "specific_enthalpy" else "broken_DG"))
+                for key in ("attribute", "geometry", "points"):
+                    audit.check(f"physical:element_binding:{q}:{e}:{field}:{key}", integer(row[key]) == integer(g[key]))
+                audit.close(f"physical:volume_binding:{q}:{e}:{field}", finite(row["physical_volume"]), finite(g["physical_volume"]))
+                audit.check(f"physical:squares:{q}:{e}:{field}", all(values[k] >= 0 for k in SQUARES))
+                for prefix in ("base","pullback","transport","eulerian"):
+                    first,second=values[prefix+"_integral"],values[prefix+"_squared_integral"]
+                    product=finite(row["physical_volume"])*second
+                    audit.check(f"physical:moment_cauchy:{q}:{e}:{field}:{prefix}", first*first <= product+2048*EPS*(first*first+product)+1e-28)
+                # Signed element moments may themselves cancel. Cauchy bounds
+                # their absolute quadrature sums using the saved positive moments.
+                scale = math.sqrt(finite(row["physical_volume"])) * (math.sqrt(values["pullback_squared_integral"])+math.sqrt(values["transport_squared_integral"]))
+                audit.close(f"physical:signed_identity:{q}:{e}:{field}", values["eulerian_integral"], values["pullback_integral"]-values["transport_integral"], scale)
+                d2, t2, e2, cross = (values[k] for k in ("pullback_squared_integral", "transport_squared_integral", "eulerian_squared_integral", "pullback_transport_cross_integral"))
+                audit.close(f"physical:squared_identity:{q}:{e}:{field}", e2, d2+t2-2*cross, d2+t2+2*abs(cross))
+                audit.check(f"physical:cauchy:{q}:{e}:{field}", abs(cross) <= math.sqrt(d2)*math.sqrt(t2) + 2048*EPS*(d2+t2) + 1e-28)
+            totals = {k: math.fsum(finite(r[k]) for r in records) for k in INTEGRALS}
+            for key in (*SQUARES, "pullback_transport_cross_integral"):
+                compare(q, field, key, totals[key])
+            for prefix in ("pullback", "transport", "eulerian"):
+                compare(q, field, prefix+"_l2_norm", math.sqrt(totals[prefix+"_squared_integral"]))
+            d2, t2, e2, cross = (totals[k] for k in ("pullback_squared_integral", "transport_squared_integral", "eulerian_squared_integral", "pullback_transport_cross_integral"))
+            compare(q, field, "squared_norm_identity_error", e2-d2-t2+2*cross, d2+t2+2*abs(cross))
+            audit.close(f"physical:total_squared_identity:{q}:{field}", e2, d2+t2-2*cross, d2+t2+2*abs(cross))
+            # Each element spans a radius interval. Midpoint bins are explicitly
+            # element-group summaries, not integrals over exact spherical shells.
+            bins = []
+            outer = max(finite(g["physical_radius_max"]) for g in geo)
+            for low, high in ((0., .25), (.25, .5), (.5, .75), (.75, .9), (.9, 1.0000001)):
+                chosen = [i for i, g in enumerate(geo) if low <= .5*(finite(g["physical_radius_min"])+finite(g["physical_radius_max"]))/outer < high]
+                bins.append({"normalized_element_midradius_interval": [low, min(high, 1.)], "elements": len(chosen),
+                             "eulerian_squared_integral": math.fsum(finite(records[i]["eulerian_squared_integral"]) for i in chosen),
+                             "pullback_squared_integral": math.fsum(finite(records[i]["pullback_squared_integral"]) for i in chosen)})
+            audit.check(f"physical:radial_bin_coverage:{q}:{field}", sum(b["elements"] for b in bins) == len(ids))
+            top = sorted(zip(sorted(ids), records), key=lambda pair: finite(pair[1]["eulerian_squared_integral"]), reverse=True)[:10]
+            field_outputs[field] = {"integrals": totals, "eulerian_over_pullback_l2": ratio(math.sqrt(e2), math.sqrt(d2)),
+                                    "eulerian_over_pullback_plus_transport_l2": ratio(math.sqrt(e2), math.sqrt(d2)+math.sqrt(t2)),
+                                    "pullback_l2_per_unit_displacement_l2": ratio(math.sqrt(d2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
+                                    "transport_l2_per_unit_displacement_l2": ratio(math.sqrt(t2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
+                                    "eulerian_l2_per_unit_displacement_l2": ratio(math.sqrt(e2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
+                                    "eulerian_rms": math.sqrt(e2/totals_g["physical_volume"]),
+                                    "pullback_transport_cosine": ratio(cross, math.sqrt(d2)*math.sqrt(t2)),
+                                    "radial_element_groups": bins,
+                                    "largest_eulerian_elements": [{"element": e, "squared_integral": finite(r["eulerian_squared_integral"]),
+                                       "fraction": ratio(finite(r["eulerian_squared_integral"]), e2)} for e, r in top]}
+        totals_g["direction_displacement_l2"] = math.sqrt(totals_g["direction_displacement_squared_integral"])
+        totals_g["direction_displacement_rms"] = math.sqrt(totals_g["direction_displacement_squared_integral"]/totals_g["physical_volume"])
+        outputs.append({"quadrature": q, "geometry": totals_g, "affine": totals_a, "fields": field_outputs})
+    for metric, expected in (("accepted_trial_prepared_normalization_unchanged", 1.), ("production_mapping_cache_unchanged", 1.),
+                             ("trial_preparations", 0.), ("volume_direction_actions", 1.)):
+        compare(0, "invariant", metric, expected)
+    for metric in ("maximum_scatter_roundtrip_error", "maximum_surface_trace_error"):
+        key = (0, "reconstruction", metric)
+        expected_summary.add(key)
+        audit.check("physical:"+metric, 0 <= observed[key] <= 512*EPS*max(1., float(np.max(np.abs(physical))), float(np.max(np.abs(physical_direction)))))
+    audit.check("physical:summary_coverage", set(observed) == expected_summary)
+    sensitivity = {}
+    for field in FIELDS:
+        sensitivity[field] = {}
+        for key in (*SQUARES, "pullback_transport_cross_integral"):
+            a, b = (out["fields"][field]["integrals"][key] for out in outputs)
+            sensitivity[field][key] = {"q10": a, "q14": b, "difference": b-a,
+                                       "relative_to_q14": ratio(abs(b-a), abs(b))}
+    return {"metadata": meta, "quadratures": outputs, "quadrature_sensitivity": sensitivity,
+            "surface": surface_fit(surface, physical, physical_direction, blocks, audit),
+            "scope": "recorded stellar broken-scalar integrals; no interface jumps, RT transport, exact shell integrals or continuum gauge test"}
+
+
+def producer_binding(direction_path, folder, audit):
+    producer = direction_path.parent
+    mapping = {"nonlinear_step_physical_base.txt":"physical_base", "nonlinear_step_normalized_base.txt":"base",
+               "nonlinear_step_state_factors.txt":"state_factors", "nonlinear_step_residual_factors.txt":"residual_factors",
+               "nonlinear_step_residual.txt":"residual", "nonlinear_step_action_full.txt":"action"}
+    present = [name for name in mapping if (producer/name).exists()]
+    if not present:
+        return {"status":"producer_neighbors_unavailable", "bound":False}
+    if len(present) != len(mapping) or not (producer/"nonlinear_step_layout.csv").exists():
+        raise ValueError("Incomplete nonlinear-step direction producer artifacts")
+    saved_hashes, equal = {}, {}
+    for name, current in mapping.items():
+        path=producer/name
+        saved_hashes[name]=digest(path)
+        a,b=vector(path),vector(folder/(P+current+".txt"))
+        exact=np.array_equal(a,b)
+        equal[name]=exact
+        audit.check("producer:vector:"+name, exact if "factors" in name else vector_near(a,b))
+    path=producer/"nonlinear_step_layout.csv"
+    saved_hashes[path.name]=digest(path)
+    audit.check("producer:layout", rows(path)==rows(folder/(P+"layout.csv")))
+    run_path=producer.parent/(producer.name+"_run.json")
+    run_bound=False
+    if run_path.exists():
+        run=json.loads(run_path.read_text())
+        for name, actual in saved_hashes.items():
+            audit.check("producer:recorded_output_hash:"+name, run["output_hashes"].get(name)==actual)
+        audit.check("producer:direction_output_hash", run["output_hashes"].get(direction_path.name)==digest(direction_path))
+        saved_hashes[run_path.name]=digest(run_path)
+        run_bound=True
+    return {"status":"producer_neighbors_compared", "bound":run_bound, "producer":str(producer),
+            "hashes":saved_hashes,"vectors_bitwise_equal":equal,
+            "qualification":"neighbor hashes checked against producer run when present; no FEM rerun"}
+
+
+def provenance(folder, required, audit):
+    root = folder.parent
+    source_file, run_file = root / "provenance.json", root / (folder.name + "_run.json")
+    if not source_file.exists() and not run_file.exists():
+        return {"status": "unavailable", "bound": False}
+    if not source_file.exists() or not run_file.exists():
+        raise ValueError("Partially present campaign provenance")
+    source, run = (json.loads(p.read_text()) for p in (source_file, run_file))
+    audit.check("provenance:case_id", run["id"] == folder.name)
+    command = run["command"]
+    final_mode = "--weak-shape-final" in command
+    audit.check("provenance:completed_probe_exit", run["exit_code"] == 0 or (final_mode and run["exit_code"] == 3))
+    audit.check("provenance:required_output_hashes", required <= set(run["output_hashes"]))
+    for name, expected in run["output_hashes"].items():
+        path = folder / name
+        audit.check("provenance:contained_output:"+name, path.resolve().is_relative_to(folder.resolve()))
+        audit.check("provenance:output_hash:"+name, digest(path) == expected)
+    for name, expected in run.get("input_state_hashes", {}).items():
+        audit.check("provenance:input_hash:"+name, digest(Path(name)) == expected)
+    case = next(c for c in source["cases"] if c["id"] == folder.name)
+    expected = [str(root / "plateau_diagnostic"), *case["arguments"], "--mesh", source["mesh"], "--output", str(folder)]
+    audit.check("provenance:command", command == expected)
+    audit.check("provenance:binary_hash", digest(Path(command[0])) == source["binary_sha256"])
+    audit.check("provenance:mesh_hash", digest(Path(source["mesh"])) == source["mesh_sha256"])
+    manifest = root / "cases.json"
+    if manifest.exists():
+        audit.check("provenance:case_manifest", json.loads(manifest.read_text()) == source["cases"])
+    meta_path=folder/"metadata.txt"
+    if meta_path.exists():
+        meta=dict(line.split("=",1) for line in meta_path.read_text().splitlines() if "=" in line)
+        audit.check("provenance:mode", ("--weak-shape-probe" in command) != final_mode)
+        audit.check("provenance:metadata_mode", meta["weak_shape_probe"]==str(int(not final_mode)) and meta["weak_shape_final"]==str(int(final_mode)))
+        audit.check("provenance:metadata_scope", meta["weak_shape_direction_space"]=="normalized_current_factors" and
+                    meta["weak_shape_operator"]=="full_monolithic" and meta["weak_shape_trial_policy"]=="raw_prepare_only")
+        for option,key in (("--linear-tol","linear_tol"),("--restart","restart"),("--max-linear","max_linear"),
+                           ("--max-newton","max_newton"),("--index","n"),("--angular-momentum","J")):
+            if option in command:
+                audit.check("provenance:metadata_argument:"+option, finite(command[command.index(option)+1])==finite(meta[key]))
+        if final_mode:
+            audit.check("provenance:reported_convergence", integer(meta["converged"])==int(run["exit_code"]==0))
+            audit.check("provenance:final_accepted_state", np.array_equal(vector(folder/"accepted_state.txt"),vector(folder/(P+"physical_base.txt"))))
+    producer = None
+    for option, saved in (("--probe-state", "physical_base"), ("--probe-direction", "direction")):
+        if option in command:
+            input_path = Path(command[command.index(option)+1])
+            if not input_path.is_absolute():
+                input_path = Path(source["source"]) / input_path
+            input_path = input_path.resolve()
+            audit.check("provenance:input_bound:"+option, str(input_path) in run.get("input_state_hashes", {}))
+            audit.check("provenance:input_vector:"+option, np.array_equal(vector(input_path), vector(folder/(P+saved+".txt"))))
+            if option == "--probe-direction":
+                producer = producer_binding(input_path,folder,audit)
+    return {"status": "runner_hashes_checked", "bound": True, "source_commit": source["commit"],
+            "provenance_sha256": digest(source_file), "run_sha256": digest(run_file), "binary_sha256": source["binary_sha256"],
+            "direction_producer":producer}
+
+
+def inspect(folder):
+    folder = Path(folder).resolve()
+    audit, required, observations = Audit(), set(), {}
+    try:
+        def get(name):
+            filename = P + name + ".txt"
+            required.add(filename)
+            return vector(folder / filename)
+        required.update(P+n+".csv" for n in ("summary", "layout", "blocks", "trials"))
+        summary = {k: finite(r["value"]) for k, r in table(folder/(P+"summary.csv"), lambda r:r["metric"]).items()}
+        for key, value in {"schema_version": 1, "accepted_commits": 0, "probe_exception": 0,
+                           "accepted_buffers_bitwise_identical": 1, "normalization_factors_bitwise_identical": 1,
+                           "direction_buffers_bitwise_identical": 1, "maximum_restore_error": 0,
+                           "restored_action_error": 0}.items():
+            audit.check("summary:"+key, summary[key] == value)
+        x, xp, f, p, pp, jp, defect, factors, rf = [get(n) for n in
+            ("base", "physical_base", "residual", "direction", "physical_direction", "action", "defect", "state_factors", "residual_factors")]
+        n = x.size
+        audit.check("vector_sizes", all(v.size == n for v in (xp,f,p,pp,jp,defect,factors,rf)))
+        audit.check("summary:state_size", summary["state_size"] == n)
+        audit.check("positive_factors", np.all(factors > 0) and np.all(rf > 0))
+        audit.check("denormalized_base", vector_near(xp, x/factors))
+        audit.check("denormalized_direction", vector_near(pp, p/factors))
+        audit.check("defect_identity", vector_near(defect, f+jp, norm(f)+norm(jp)))
+        audit.check("restored_residual_exact", np.array_equal(f, get("restored_residual")))
+        audit.check("restored_action_exact", np.array_equal(jp, get("restored_action")))
+        fn, pn, jn, dn = map(norm, (f,p,jp,defect))
+        if fn <= 0 or pn <= 0:
+            raise ValueError("Probe requires nonzero base residual and direction")
+        threshold = max(summary["outer_absolute_tolerance"], summary["outer_tolerance"]*fn)
+        audit.check("valid_linear_controls", 0 < summary["outer_tolerance"] < 1 and summary["outer_absolute_tolerance"] >= 0 and
+                    integer(summary["maximum_iterations"]) > 0 and integer(summary["restart"]) > 0)
+        for key, value in {"base_norm":fn, "direction_norm":pn, "action_norm":jn, "verified_defect_norm":dn,
+                           "verified_relative_defect":dn/fn, "verified_threshold":threshold}.items():
+            audit.close("summary:"+key, summary[key], value)
+        audit.check("verified_gate_truth", summary["verified_gate"] in (0,1) and bool(summary["verified_gate"]) == (dn <= threshold))
+        audit.check("supplied_direction_flag", summary["supplied_direction"] in (0,1))
+        if not summary["supplied_direction"]:
+            audit.close("linear_reported_true_norm", summary["linear_reported_true_norm"], dn, fn+jn)
+            audit.check("linear_iterations", 0 <= integer(summary["linear_iterations"]) <= integer(summary["maximum_iterations"]))
+        blocks = layout(folder/(P+"layout.csv"), n)
+        actions = [get("action_b"+str(i)) for i in range(len(blocks["value"]))]
+        audit.check("split_column_count", summary["split_columns"] == len(actions))
+        audit.check("direct_action_count", summary["direct_jacobian_applications"] == len(actions)+2)
+        audit.check("preparation_count", summary["preparations"] == 6 + int(not summary["supplied_direction"]))
+        audit.check("finite_nonnegative_probe_time", summary["total_probe_seconds"] >= 0)
+        total = np.sum(actions, axis=0)
+        action_scale = math.fsum(norm(a) for a in actions)
+        audit.check("split_action_sum", vector_near(total, jp, action_scale))
+        audit.close("reported_split_error", summary["split_action_sum_error"], norm(total-jp), action_scale)
+        reported = table(folder/(P+"blocks.csv"), lambda r:(integer(r["column"]),r["residual_block"]))
+        audit.check("block_table_coverage", set(reported) == {(i,r[0]) for i in range(len(actions)) for r in blocks["residual"]})
+        for i,(name,start,count) in enumerate(blocks["value"]):
+            for rname,offset,size in blocks["residual"]:
+                row = reported[(i,rname)]
+                audit.check(f"block_metadata:{i}:{rname}", row["value_block"]==name and integer(row["value_offset"])==start and
+                            integer(row["value_size"])==count and integer(row["residual_offset"])==offset and integer(row["residual_size"])==size)
+                audit.close(f"block_direction_norm:{i}:{rname}", finite(row["direction_norm"]), norm(p[start:start+count]))
+                audit.close(f"block_action_norm:{i}:{rname}", finite(row["action_norm"]), norm(actions[i][offset:offset+size]))
+        svd, smallest_action = restricted_svd(p, actions, blocks["value"], audit)
+        svd["smallest_action_residual_blocks"] = {name:norm(smallest_action[o:o+s]) for name,o,s in blocks["residual"]}
+        trial_rows = table(folder/(P+"trials.csv"), lambda r:integer(r["trial"]))
+        audit.check("trial_coverage", set(trial_rows)=={0,1,2})
+        trials = []
+        for i,alpha in enumerate((1/128,1/512,1/2048)):
+            row = trial_rows[i]
+            audit.check(f"trial_alpha:{i}", finite(row["alpha"])==alpha)
+            audit.check(f"trial_prepared_flag:{i}", row["prepared"] in ("0","1"))
+            audit.check(f"trial_candidate:{i}", vector_near(get("candidate_a"+str(i)), x+alpha*p, norm(x)+alpha*pn))
+            prediction = f+alpha*jp
+            audit.close(f"trial_prediction_norm:{i}", finite(row["linear_prediction_norm"]), norm(prediction))
+            prepared = row["prepared"]=="1"
+            entry = {"trial":i, "alpha":alpha, "prepared":prepared, "reason":integer(row["reason"]), "stage":integer(row["stage"])}
+            if prepared:
+                audit.check(f"trial_prepared_reasons:{i}", entry["reason"] == entry["stage"] == -1)
+                ft = get("residual_a"+str(i)); remainder = ft-prediction
+                values = {"residual_norm":norm(ft), "remainder_norm":norm(remainder), "remainder_over_alpha":norm(remainder)/alpha}
+                for key,value in values.items():
+                    audit.close(f"trial:{i}:{key}", finite(row[key]), value, (fn+alpha*jn+norm(ft)) / (alpha if key=="remainder_over_alpha" else 1))
+                entry.update(values)
+                entry.update(residual_over_base=norm(ft)/fn, remainder_over_alpha_squared=norm(remainder)/alpha**2,
+                             linear_dot_remainder=float(checked_einsum("i,i->",prediction,remainder)))
+            else:
+                audit.check(f"trial_failed_reasons:{i}", entry["reason"]>=0 and entry["stage"]>=0)
+                audit.check(f"trial_failed_vector_absent:{i}", not (folder/(P+"residual_a"+str(i)+".txt")).exists())
+                audit.check(f"trial_failed_metrics:{i}", all(math.isnan(float(row[k])) for k in ("residual_norm","remainder_norm","remainder_over_alpha")))
+            trials.append(entry)
+        physical_data = physical_analysis(folder,xp,pp,blocks["value"],audit,required)
+        if (folder/"metadata.txt").exists():
+            required.add("metadata.txt")
+            meta=dict(line.split("=",1) for line in (folder/"metadata.txt").read_text().splitlines() if "=" in line)
+            for key,source in (("outer_tolerance","linear_tol"),("restart","restart"),("maximum_iterations","max_linear")):
+                audit.check("metadata:linear_control:"+key,summary[key]==finite(meta[source]))
+            if meta.get("weak_shape_final")=="1":
+                required.add("accepted_state.txt")
+        binding = provenance(folder,required,audit)
+        observations = {"summary":summary, "provenance":binding,
+                        "linear_classification":"verified_newton_direction" if dn<=threshold else "finite_direction_failed_true_gate",
+                        "verified_linear_gate":dn<=threshold, "restricted_block_span_svd":svd, "trials":trials,
+                        "block_action_cancellation_factor":ratio(action_scale,jn), "physical":physical_data,
+                        "input_vector_hashes":{name:digest(folder/(P+name+".txt")) for name in ("physical_base","base","direction","state_factors","residual_factors")},
+                        "layout_sha256":digest(folder/(P+"layout.csv"))}
+    except (OSError, ValueError, KeyError, TypeError, IndexError, StopIteration, OverflowError, ZeroDivisionError, np.linalg.LinAlgError) as error:
+        audit.check("complete_consistent_artifacts: " + str(error),False)
+    passed = not audit.failures
+    return {"directory":str(folder), "artifact_pass":passed, "observations_valid":passed,
+            "checks_count":audit.count, "failures":audit.failures,
+            "required_artifacts":sorted(required), "observations":observations if passed else {}}
+
+
+def comparisons(cases):
+    output=[]
+    for i,left in enumerate(cases):
+        for right in cases[i+1:]:
+            if not left["artifact_pass"] or not right["artifact_pass"]:
+                continue
+            a,b=left["observations"],right["observations"]
+            same_layout=a["layout_sha256"]==b["layout_sha256"]
+            same={k:a["input_vector_hashes"][k]==b["input_vector_hashes"][k] for k in a["input_vector_hashes"]}
+            output.append({"left":left["directory"],"right":right["directory"],"same_layout":same_layout,"same_saved_vectors":same,
+                           "same_state_metric_comparison":same_layout and all(same[k] for k in ("physical_base","base","state_factors","residual_factors")),
+                           "qualification":"different discretizations or factors require physical comparisons; raw coefficient norms are not comparable"})
+    return output
+
+
+def tiny_fixture(folder, failed_gate=False):
+    """Construct a complete eight-coordinate synthetic producer, with no FEM."""
+    folder.mkdir(parents=True,exist_ok=True)
+    def save(name,v):
+        (folder/(P+name+".txt")).write_text(str(len(v))+"\n"+"\n".join(format(float(x),".17g") for x in v)+"\n")
+    def write(name,records):
+        with (folder/name).open("w",newline="") as stream:
+            out=csv.DictWriter(stream,fieldnames=list(records[0]));out.writeheader();out.writerows(records)
+    sf=np.arange(1.,9.);xp=np.array([2.,0.,0.,0.,0.,0.,0.,1.]);x=xp*sf
+    p=np.arange(1.,9.)*1e-4;f=-p;jp=p*(.9 if failed_gate else 1.)
+    for name,v in {"base":x,"physical_base":xp,"residual":f,"direction":p,"physical_direction":p/sf,
+                   "action":jp,"defect":f+jp,"state_factors":sf,"residual_factors":np.ones(8),
+                   "restored_residual":f,"restored_action":jp}.items():save(name,v)
+    blocks=[("density",0,1),("surface_deformation",1,6),("constant",7,1)]
+    write(P+"layout.csv",[{"kind":kind,"block":name,"offset":o,"size":n} for kind in ("value","residual") for name,o,n in blocks])
+    br=[]
+    for i,(name,o,n) in enumerate(blocks):
+        a=np.zeros(8);a[o:o+n]=jp[o:o+n];save("action_b"+str(i),a)
+        for rn,ro,rs in blocks:
+            br.append(dict(column=i,value_block=name,value_offset=o,value_size=n,residual_block=rn,residual_offset=ro,
+                           residual_size=rs,direction_norm=norm(p[o:o+n]),action_norm=norm(a[ro:ro+rs])))
+    write(P+"blocks.csv",br)
+    summary=dict(schema_version=1,state_size=8,supplied_direction=1,accepted_commits=0,outer_tolerance=.03,
+                 outer_absolute_tolerance=0,maximum_iterations=200,restart=160,base_norm=norm(f),direction_norm=norm(p),
+                 action_norm=norm(jp),verified_defect_norm=norm(f+jp),verified_relative_defect=norm(f+jp)/norm(f),
+                 verified_threshold=.03*norm(f),verified_gate=int(not failed_gate),split_columns=3,split_action_sum_error=0,
+                 restored_action_error=0,probe_exception=0,accepted_buffers_bitwise_identical=1,normalization_factors_bitwise_identical=1,
+                 direction_buffers_bitwise_identical=1,maximum_restore_error=0,preparations=6,direct_jacobian_applications=5,total_probe_seconds=1.)
+    write(P+"summary.csv",[dict(metric=k,value=v) for k,v in summary.items()])
+    tr=[]
+    for i,a in enumerate((1/128,1/512,1/2048)):
+        save("candidate_a"+str(i),x+a*p);pred=f+a*jp;ft=pred+a*a*np.ones(8)*1e-5;save("residual_a"+str(i),ft)
+        tr.append(dict(trial=i,alpha=a,prepared=1,reason=-1,stage=-1,residual_norm=norm(ft),remainder_norm=norm(ft-pred),
+                       remainder_over_alpha=norm(ft-pred)/a,linear_prediction_norm=norm(pred)))
+    write(P+"trials.csv",tr)
+    meta=dict(schema="weak_shape_physical_v1",dg_derivatives="broken DG derivatives",gravity_vector="not evaluated",
+              density_family_order="2",enthalpy_family_order="3",potential_family_order="2",displacement_family_order="3")
+    write(PP+"metadata.csv",[dict(key=k,value=v) for k,v in meta.items()])
+    elements=[];geometry=[];affine=[];ps=[]
+    for q in (10,14):
+        g=dict(quadrature=q,element=0,attribute=1,geometry=5,points=8,reference_mesh_geometry_order=4,density_element_order=2,
+               enthalpy_element_order=3,potential_element_order=2,displacement_element_order=3,reference_radius_min=.1,
+               reference_radius_max=.9,physical_radius_min=.1,physical_radius_max=.9,physical_volume=1.,
+               base_displacement_squared_integral=0.,direction_displacement_squared_integral=.2)
+        geometry.append(g)
+        ar=dict(quadrature=q,element=0,points=8,reference_volume=1.,value_error_squared_integral=.01,
+                reference_value_squared_integral=1.,gradient_error_squared_integral=.03,reference_gradient_squared_integral=3.)
+        affine.append(ar)
+        for field in FIELDS:
+            vals=dict(zip(INTEGRALS,(2.,3.,2.,1.,4.,9.,4.,1.,6.)))
+            elements.append(dict(quadrature=q,element=0,attribute=1,geometry=5,field=field,
+                                 gradient_kind="continuous_H1" if field=="specific_enthalpy" else "broken_DG",points=8,physical_volume=1.,**vals))
+            measured={k:vals[k] for k in (*SQUARES,"pullback_transport_cross_integral")}
+            measured.update(pullback_l2_norm=3.,transport_l2_norm=2.,eulerian_l2_norm=1.,squared_norm_identity_error=0.)
+            ps.extend(dict(quadrature=q,field=field,metric=k,value=v) for k,v in measured.items())
+        ps.extend(dict(quadrature=q,field="geometry",metric=k,value=v) for k,v in
+                  dict(stellar_elements=1,physical_volume=1.,base_displacement_squared_integral=0.,direction_displacement_squared_integral=.2).items())
+        am={k:v for k,v in ar.items() if k not in ("quadrature","element","points")}
+        am.update(relative_value_l2_error=.1,relative_gradient_l2_error=.1)
+        ps.extend(dict(quadrature=q,field="affine",metric=k,value=v) for k,v in am.items())
+    ps.extend(dict(quadrature=0,field="invariant",metric=k,value=v) for k,v in
+              dict(accepted_trial_prepared_normalization_unchanged=1,production_mapping_cache_unchanged=1,trial_preparations=0,volume_direction_actions=1).items())
+    ps.extend(dict(quadrature=0,field="reconstruction",metric=k,value=0.) for k in
+              ("maximum_scatter_roundtrip_error","maximum_surface_trace_error"))
+    for name,data in (("elements",elements),("geometry",geometry),("affine",affine),("summary",ps)):write(PP+name+".csv",data)
+    xyz=[[1.,0.,0.],[-1.,0.,0.],[0.,1.,0.],[0.,-1.,0.],[0.,0.,1.],[0.,0.,-1.]]
+    write(PP+"surface_nodes.csv",[dict(parameter=i,volume_scalar_true_dof=i,unit_radial_x=v[0],unit_radial_y=v[1],
+          unit_radial_z=v[2],reference_radius=1.,physical_radius=1.,base_parameter=0.,physical_direction_parameter=p[i+1]/sf[i+1]) for i,v in enumerate(xyz)])
+
+
+def self_test():
+    """Tiny algebraic controls only; never opens campaign artifacts."""
+    import tempfile
+    count=0
+    def require(condition):
+        nonlocal count
+        assert condition
+        count+=1
+    audit=Audit()
+    p=np.array([2.,3.,0.])
+    actions=[np.array([2.,0.,0.]),np.array([-3.,3e-8,0.]),np.zeros(3)]
+    out,action=restricted_svd(p,actions,[("a",0,1),("b",1,1),("zero",2,1)],audit)
+    require(not audit.failures)
+    require(near(out["singular_values"][-1],1e-8/math.sqrt(2),1e-8))
+    require(out["active_blocks"]==["a","b"])
+    require(norm(action)<1e-7)
+    bad=Audit();restricted_svd(p,actions[:2]+[np.ones(3)],[("a",0,1),("b",1,1),("zero",2,1)],bad)
+    require(bool(bad.failures))
+    require(not vector_near(np.array([1.,2.]),np.array([1.,2.001])))
+    require(near(1e-20,0.,1.))
+    # This size triggers spurious Accelerate matmul flags on the investigation
+    # host. Use independent scalar sums, not another BLAS product, as the oracle.
+    a=np.linspace(-1.,1.,257*6).reshape(257,6);v=np.arange(1.,7.)
+    expected=np.array([math.fsum(float(x*y) for x,y in zip(row,v)) for row in a])
+    with np.errstate(divide="raise",over="raise",invalid="raise"):
+        require(vector_near(checked_einsum("ij,j->i",a,v),expected,norm(expected)))
+        b=np.arange(18.).reshape(6,3)/18.
+        expected_matrix=np.array([[math.fsum(float(a[i,k]*b[k,j]) for k in range(6)) for j in range(3)] for i in range(257)])
+        require(vector_near(checked_einsum("ij,jk->ik",a,b),expected_matrix,norm(expected_matrix)))
+    try:checked_einsum("i,i->",np.array([math.inf]),np.ones(1));raise AssertionError("nonfinite contraction accepted")
+    except ValueError:require(True)
+    with tempfile.TemporaryDirectory() as directory:
+        d=Path(directory)
+        (d/"v").write_text("3\n1\n2\n3\n")
+        require(np.array_equal(vector(d/"v"),np.array([1.,2.,3.])))
+        for text in ("2\n1\n", "2\n1\nnan\n", "2\n1\ninf\n"):
+            (d/"v").write_text(text)
+            try:vector(d/"v");raise AssertionError("invalid vector accepted")
+            except ValueError:require(True)
+        (d/"layout").write_text("kind,block,offset,size\nvalue,a,0,2\nresidual,b,0,2\n")
+        require(layout(d/"layout",2)["value"]==[("a",0,2)])
+        (d/"layout").write_text("kind,block,offset,size\nvalue,a,1,2\nresidual,b,0,2\n")
+        try:layout(d/"layout",2);raise AssertionError("bad layout accepted")
+        except ValueError:require(True)
+    xyz=np.array([[1.,0.,0.],[-1.,0.,0.],[0.,1.,0.],[0.,-1.,0.],[0.,0.,1.],[0.,0.,-1.]])
+    values=2.+3.*xyz[:,0]
+    surf={i:{"volume_scalar_true_dof":str(i),"unit_radial_x":str(v[0]),"unit_radial_y":str(v[1]),"unit_radial_z":str(v[2]),
+             "reference_radius":"1","physical_radius":"1","base_parameter":"0","physical_direction_parameter":str(values[i])} for i,v in enumerate(xyz)}
+    audit=Audit();fits=surface_fit(surf,np.zeros(6),values,[("surface_deformation",0,6)],audit)
+    require(not audit.failures)
+    require(fits["fits"][1]["relative_nodal_l2_remainder"]<1e-14)
+    require(fits["fits"][0]["relative_nodal_l2_remainder"]>.1)
+    with tempfile.TemporaryDirectory() as directory:
+        d=Path(directory).resolve()/"case"
+        tiny_fixture(d)
+        result=inspect(d)
+        assert result["artifact_pass"],result["failures"]
+        require(result["observations"]["verified_linear_gate"])
+        require(result["observations"]["provenance"]["status"]=="unavailable")
+        require(result["observations"]["physical"]["quadratures"][0]["fields"]["density"]["eulerian_over_pullback_plus_transport_l2"]==.2)
+        tiny_fixture(d,failed_gate=True)
+        result=inspect(d)
+        require(result["artifact_pass"] and result["observations"]["linear_classification"]=="finite_direction_failed_true_gate")
+        for name in ("defect.txt","restored_action.txt","candidate_a1.txt","physical_base.txt","action_b1.txt"):
+            tiny_fixture(d)
+            path=d/(P+name);lines=path.read_text().splitlines();lines[1]=str(float(lines[1])+.001);path.write_text("\n".join(lines)+"\n")
+            require(not inspect(d)["artifact_pass"])
+        for filename,key in (("elements","pullback_transport_cross_integral"),("summary","value"),("geometry","physical_volume")):
+            tiny_fixture(d);path=d/(PP+filename+".csv");data=rows(path);data[0][key]=str(float(data[0][key])+.1)
+            with path.open("w",newline="") as stream:
+                writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
+            require(not inspect(d)["artifact_pass"])
+        for missing in (False,True):
+            tiny_fixture(d);path=d/(PP+"summary.csv");data=rows(path)
+            if missing:data=[r for r in data if r["metric"]!="production_mapping_cache_unchanged"]
+            else:
+                next(r for r in data if r["metric"]=="production_mapping_cache_unchanged")["value"]="0"
+            with path.open("w",newline="") as stream:
+                writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
+            require(not inspect(d)["artifact_pass"])
+        tiny_fixture(d)
+        path=d/(P+"trials.csv");data=rows(path);data[0].update(prepared="0",reason="2",stage="3",
+            residual_norm="nan",remainder_norm="nan",remainder_over_alpha="nan")
+        with path.open("w",newline="") as stream:
+            writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
+        (d/(P+"residual_a0.txt")).unlink()
+        require(inspect(d)["artifact_pass"])
+        tiny_fixture(d)
+        root=d.parent
+        executable=root/"plateau_diagnostic";executable.write_text("synthetic binary fixture")
+        mesh=root/"mesh";mesh.write_text("synthetic mesh fixture")
+        case={"id":"case","arguments":["--weak-shape-probe"]}
+        prov={"commit":"0"*40,"source":str(root),"cases":[case],"mesh":str(mesh),"mesh_sha256":digest(mesh),"binary_sha256":digest(executable)}
+        (root/"provenance.json").write_text(json.dumps(prov))
+        run={"id":"case","exit_code":0,"input_state_hashes":{},"command":[str(executable),"--weak-shape-probe","--mesh",str(mesh),"--output",str(d)],
+             "output_hashes":{p.name:digest(p) for p in d.iterdir()}}
+        (root/"case_run.json").write_text(json.dumps(run))
+        result=inspect(d)
+        assert result["artifact_pass"],result["failures"]
+        require(result["observations"]["provenance"]["bound"])
+        run["output_hashes"][P+"action.txt"]="0"*64
+        (root/"case_run.json").write_text(json.dumps(run))
+        require(not inspect(d)["artifact_pass"])
+        producer=root/"producer";producer.mkdir()
+        for old,new in {"physical_base":"physical_base","normalized_base":"base","state_factors":"state_factors",
+                        "residual_factors":"residual_factors","residual":"residual","action_full":"action","direction_full":"direction"}.items():
+            (producer/("nonlinear_step_"+old+".txt")).write_bytes((d/(P+new+".txt")).read_bytes())
+        (producer/"nonlinear_step_layout.csv").write_bytes((d/(P+"layout.csv")).read_bytes())
+        (root/"producer_run.json").write_text(json.dumps({"output_hashes":{p.name:digest(p) for p in producer.iterdir()}}))
+        audit=Audit();bound=producer_binding(producer/"nonlinear_step_direction_full.txt",d,audit)
+        require(bound["bound"] and not audit.failures)
+        (producer/"nonlinear_step_state_factors.txt").write_text("8\n"+"2\n"*8)
+        audit=Audit();producer_binding(producer/"nonlinear_step_direction_full.txt",d,audit)
+        require(bool(audit.failures))
+    print(f"{count} tiny synthetic checks passed")
+
+
+def main():
+    parser=argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("directories",nargs="*",type=Path)
+    parser.add_argument("--output",type=Path)
+    parser.add_argument("--self-test",action="store_true")
+    args=parser.parse_args()
+    if args.self_test:
+        self_test();return 0
+    if not args.directories or args.output is None:
+        parser.error("Provide one or more probe directories and --output")
+    destination=args.output.resolve()
+    if destination.exists() or any(destination.is_relative_to(d.resolve()) for d in args.directories):
+        parser.error("Use a fresh output JSON outside the immutable probe directories")
+    cases=[inspect(d) for d in args.directories]
+    result={"schema_version":1,"passed":all(c["artifact_pass"] for c in cases),"checker_sha256":digest(Path(__file__)),
+            "numpy_version":np.__version__,"tolerances":{"relative":RTOL,"cancellation_roundoff_multiplier":2048,"absolute_floor":1e-28},
+            "scope":"Independent recorded-vector algebra/hash and physical-integral identity audit; no FEM rerun, quadrature rerun, full-rank, gauge, convergence or physical-model certificate",
+            "cases":cases,"comparisons":comparisons(cases)}
+    destination.parent.mkdir(parents=True,exist_ok=True)
+    destination.write_text(json.dumps(result,indent=2,allow_nan=False)+"\n")
+    print(json.dumps({"passed":result["passed"],"cases":len(cases),"output":str(destination)}))
+    return 0 if result["passed"] else 1
+
+
+if __name__=="__main__":
+    sys.exit(main())
diff --git a/tests/fem/reference_tables.cpp b/tests/fem/reference_tables.cpp
index 2b9de9b..7d7a458 100644
--- a/tests/fem/reference_tables.cpp
+++ b/tests/fem/reference_tables.cpp
@@ -556,3 +556,120 @@ TEST_CASE("Tensor Gradient Transposes Match Dense Assembly And Adjoint Products"
               mfem::IntRules.Get(mfem::Geometry::CUBE, 2*order+3));
     }
 }
+
+TEST_CASE(
+    "Quartic Curved Hex Geometry Distinguishes Q3 And Q4 Affine Physical Displacements",
+    tags::unit &tags::geometry &tags::quadrature
+) {
+    constexpr int dimension = 3;
+    constexpr double warpAmplitude = 0.2;
+    constexpr double dilation = 0.2;
+    constexpr double exactDeterminant = (1.0 + dilation) * (1.0 + dilation) * (1.0 + dilation);
+    int ranks = 1;
+    MPI_Comm_size(MPI_COMM_WORLD, &ranks);
+
+    // Subdivide only y so every rank can retain the full quartic x variation.
+    // Its coordinate Jacobian is 1 + 2*a*x*(1-x)*(1-2*x) >= 1-a/2 > 0.
+    auto serial = mfem::Mesh::MakeCartesian3D(1, ranks, 1, mfem::Element::HEXAHEDRON);
+    serial.SetCurvature(4, false, dimension, mfem::Ordering::byNODES);
+    mfem::VectorFunctionCoefficient warp(dimension, [](const mfem::Vector &x, mfem::Vector &y) {
+        y = x;
+        const double bubble = x(0) * (1.0 - x(0));
+        y(0) += warpAmplitude * bubble * bubble;
+    });
+    serial.Transform(warp);
+    mfem::ParMesh mesh(MPI_COMM_WORLD, serial);
+    REQUIRE(mesh.GetNodes() != nullptr);
+    CHECK(mesh.GetNodes()->FESpace()->GetMaxElementOrder() == 4);
+
+    // u(x)=s*x pulls back to the quartic coordinate polynomial. This tests the
+    // ambient vector space, not the production surface prescription's r^2 lift.
+    mfem::VectorFunctionCoefficient affineDisplacement(
+        dimension, [](const mfem::Vector &x, mfem::Vector &u) {
+            u = x;
+            u *= dilation;
+        }
+    );
+    const auto &rule = mfem::IntRules.Get(mfem::Geometry::CUBE, 24);
+    for (const int order : {3, 4}) {
+        CAPTURE(order);
+        mfem::H1_FECollection collection(order, dimension, mfem::BasisType::GaussLobatto);
+        mfem::ParFiniteElementSpace space(&mesh, &collection, dimension, mfem::Ordering::byNODES);
+        CHECK(collection.FiniteElementForGeometry(mfem::Geometry::CUBE)->GetDof() ==
+              (order + 1) * (order + 1) * (order + 1));
+        mfem::ParGridFunction field(&space);
+        field.ProjectCoefficient(affineDisplacement);
+        mfem::Vector state;
+        field.GetTrueDofs(state);
+        field.SetFromTrueDofs(state);
+
+        // Accumulate rather than assert per quadrature point, and reject
+        // non-finite samples explicitly so a max reduction cannot hide them.
+        enum Error { coordinates, coordinateJacobian, values, gradients, determinants, count };
+        std::array<double, count> localErrors{};
+        int localFinite = 1;
+        double localGeometryMinimum = std::numeric_limits<double>::infinity();
+        auto recordError = [&](const Error error, const double value) {
+            if (!std::isfinite(value)) localFinite = 0;
+            else localErrors[error] = std::max(localErrors[error], std::abs(value));
+        };
+        mfem::Vector position(dimension), value(dimension);
+        mfem::DenseMatrix gradient(dimension);
+        for (int element = 0; element < mesh.GetNE(); ++element) {
+            auto &transformation = *mesh.GetElementTransformation(element);
+            for (int point = 0; point < rule.GetNPoints(); ++point) {
+                const auto &ip = rule.IntPoint(point);
+                transformation.SetIntPoint(&ip);
+                transformation.Transform(ip, position);
+                const double bubble = ip.x * (1.0 - ip.x);
+                recordError(coordinates, position(0) - (ip.x + warpAmplitude * bubble * bubble));
+                const double expectedJacobian = 1.0 + 2.0 * warpAmplitude * bubble * (1.0 - 2.0 * ip.x);
+                const double coordinateDeterminant = ranks * transformation.Jacobian().Det();
+                recordError(coordinateJacobian, coordinateDeterminant - expectedJacobian);
+                if (std::isfinite(coordinateDeterminant))
+                    localGeometryMinimum = std::min(localGeometryMinimum, coordinateDeterminant);
+
+                field.GetVectorValue(transformation, ip, value);
+                field.GetVectorGradient(transformation, gradient);
+                for (int component = 0; component < dimension; ++component) {
+                    recordError(values, value(component) - dilation * position(component));
+                    for (int axis = 0; axis < dimension; ++axis)
+                        recordError(gradients, gradient(component, axis) -
+                            (component == axis ? dilation : 0.0));
+                }
+                for (int component = 0; component < dimension; ++component) gradient(component, component) += 1.0;
+                recordError(determinants, gradient.Det() - exactDeterminant);
+            }
+        }
+        std::array<double, count> errors{};
+        int finite = 0;
+        double geometryMinimum = 0.0;
+        MPI_Allreduce(localErrors.data(), errors.data(), count, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+        MPI_Allreduce(&localFinite, &finite, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
+        MPI_Allreduce(&localGeometryMinimum, &geometryMinimum, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+        REQUIRE(finite == 1);
+        CHECK(errors[coordinates] < 2.0e-13);
+        CHECK(errors[coordinateJacobian] < 2.0e-12);
+        CHECK(geometryMinimum >= 1.0 - warpAmplitude / 2.0 - 2.0e-12);
+
+        mean_field::fem::PreparedGeometryInspector inspector;
+        const double fastMinimum = inspector.Inspect(space, state);
+        const double referenceMinimum = inspector.Inspect(space, state, true);
+        REQUIRE(std::isfinite(fastMinimum));
+        REQUIRE(std::isfinite(referenceMinimum));
+        CHECK(fastMinimum > 0.0);
+        CHECK(std::abs(fastMinimum - referenceMinimum) < 2.0e-12);
+        CHECK(inspector.GetTensorElementCount() == mesh.GetNE());
+        if (order == 4) {
+            CHECK(errors[values] < 2.0e-13);
+            CHECK(errors[gradients] < 2.0e-12);
+            CHECK(errors[determinants] < 5.0e-12);
+            CHECK(std::abs(fastMinimum - exactDeterminant) < 5.0e-12);
+        } else {
+            CHECK(errors[values] > 1.0e-5);
+            CHECK(errors[gradients] > 1.0e-4);
+            CHECK(errors[determinants] > 1.0e-4);
+            CHECK(std::abs(fastMinimum - exactDeterminant) > 1.0e-4);
+        }
+    }
+}