Add frozen smooth-surface restriction study with physical dual maps

experiments/REDUCED_SURFACE_STUDY_2026-09-18.md

@@ -0,0 +1,64 @@
+# One-hour reduced smooth surface study — 2026-09-18
+
+Authorized wall-clock window: 19:30:15–20:30:15 UTC. Branch:
+research/reduced-smooth-surface-20260918, parent 6ceda61.
+
+Keep Q3 displacement/surface, Q4 reference geometry, Q2 DG density/potential,
+Q3 enthalpy and RT2 gravity. No central-density constraint. n=1.5,J=0.
+No production API, operator or extension change. Single MPI rank.
+
+Frozen affine embedding q=q0+B a. The base q0 is never projected. Construct
+real spherical harmonics through degree L at actual surface radial nodes;
+orthonormalize the Q3 interpolants using q14 prepared physical surface area.
+Thin QR of S B and R^-1 B supplies separate normalized state and dual test
+maps. Solve Qf^T Jnormalized Qx p=-Qf^T Fnormalized. The projected full inverse
+Qx^T Mnormalized^-1 Qf is an approximation, not an exact reduced inverse.
+All other unknowns/equations remain active. Cold FGMRES, eta=.03, max200,
+restart160, same gravity-Schur preconditioner (one sweep, relaxation.5).
+
+Priority: fresh full-space direction at the prior Q3 six-update endpoint,
+then L4 there, then L2 and L8 (L6 optional); analytic-seed L4 only if time
+permits. Stop launching FEM cases by20:20; reserve remaining time for independent
+artifact audit and reporting. No convergence claim from frozen steps. Cases,
+builds, and full numerical audits run serially. A library-only Q3 warm rebuild
+may overlap editing experiment-only headers; no production source changes.
+
+Raw independent trial alphas 1,1/4,1/16,1/1024. Preserve original and reduced
+residuals, separate linear gates, full/reduced slopes, excluded shape forces,
+and exact restoration checks. Observational merit criterion is
+norm(Ftrial)<=(1-1e-4 alpha)*norm(Fbase), evaluated separately in both spaces.
+No trial accepted, no SOC, no positivity repair. q14/q18 physical measurements
+at base and the first/largest prepared trial only; all trial residuals saved.
+Direction physical observer q10/q14 reused from the preceding audited study.
+
+Acceptance of an observation requires algebra/transpose/full-rank control tests,
+true reduced gate classification, finite vectors, source/input/output hashes,
+and exact saved-state restoration. Reduced success is not full success.
+Virial1e-5, interior physical errors, density sign, and omitted force remain
+separate measurements. No selected norm reduction can certify a physical root.
+
+## Implementation
+
+Experiment-only headers:
+- `reduced_surface_basis.hpp`: scalar radial-amplitude FE trace basis and
+  quadrature-defined physical surface mass metric, explicit rank/orthogonality
+  checks, and raw data for independent verification.
+- `reduced_surface_operators.hpp`: two separate QR maps, reduced operator,
+  projected flexible preconditioner and deterministic algebra checks.
+- `reduced_surface_diagnostics.hpp`: frozen solve, both original/reduced gates,
+  excluded-force decomposition, raw trials and restoration invariants.
+- `audit_reduced_surface.py`: independent artifact algebra and corruption tests.
+
+The diagnostic executable adds experiment-only flags --reduced-surface-probe
+and --surface-degree (0 means full-space control); no production-facing API is
+modified. Degrees are runtime experimental choices while operator/inverse
+compatibility and the inherited preconditioner application contract remain
+compile-time constrained. Actual mesh dimensions and normalization factors are
+necessarily checked at runtime. Existing physical observers are reused.
+
+Evidence root: `../../diagnostics/reduced_surface_2026-09-18` relative to this
+checkout. The preceding study remains on its two committed research branches.
+
+## Results
+
+Pending the bounded serial campaign and its independent artifact audit.

experiments/audit_reduced_surface.py

@@ -0,0 +1,517 @@
+#!/usr/bin/env python3
+"""Audit saved reduced-surface algebra; a reduced gate is not a full solve.
+
+No FEM action, preparation, or physical quadrature is independently repeated.
+The dual residual restriction is checked against the physical basis and saved
+normalization, rather than assumed to be the state injection's transpose.
+"""
+import argparse
+import json
+import math
+from pathlib import Path
+import sys
+
+from audit_shape_field_weak import (Audit, EPS, RTOL, checked_einsum as einsum,
+                                  digest, finite, integer, layout, near, norm,
+                                  np, rows, table, vector, vector_near)
+from analyze_shape_surface_modes import harmonics
+
+P = "reduced_surface_"
+MAPS = ("basis_physical", "state_basis", "residual_basis", "state_transform", "residual_transform")
+VECTORS = ("base", "physical_base", "residual", "direction", "physical_direction", "action", "defect",
+           "state_factors", "residual_factors", "reduced_residual", "reduced_direction", "reduced_action",
+           "reduced_defect", "excluded_residual", "excluded_defect", "restored_residual", "restored_action")
+ALPHAS = (1., .25, .0625, 1 / 1024)
+
+
+def matrix(path):
+    with Path(path).open() as stream:
+        dimensions = [int(x) for x in stream.readline().split()]
+        values = np.loadtxt(stream, ndmin=1)
+    if len(dimensions) != 2 or min(dimensions) <= 0 or values.ndim != 1 or values.size != math.prod(dimensions) or not np.isfinite(values).all():
+        raise ValueError(f"Invalid finite row-major matrix: {path}")
+    return values.reshape(dimensions)
+
+
+def close_array(audit, label, actual, expected, scale=0.):
+    audit.check(label, vector_near(np.asarray(actual).ravel(), np.asarray(expected).ravel(), scale))
+
+
+class Coordinates:
+    def __init__(self, summary, maps):
+        self.n = integer(summary["full_size"])
+        self.s = integer(summary["shape_size"])
+        self.k = integer(summary["mode_count"])
+        self.vo = integer(summary["shape_value_offset"])
+        self.ro = integer(summary["shape_residual_offset"])
+        self.identity = summary["identity_control"] == 1
+        self.maps = maps
+        if not (0 < self.k <= self.s <= self.n and 0 <= self.vo <= self.n - self.s and 0 <= self.ro <= self.n - self.s):
+            raise ValueError("Invalid coordinate dimensions")
+
+    def restrict(self, values, residual=True):
+        if self.identity:
+            return values.copy()
+        offset = self.ro if residual else self.vo
+        basis = self.maps["residual_basis" if residual else "state_basis"]
+        middle = einsum("ij,i->j", basis, values[offset:offset + self.s])
+        return np.concatenate((values[:offset], middle, values[offset + self.s:]))
+
+    def inject(self, values, residual=False):
+        if self.identity:
+            return values.copy()
+        offset = self.ro if residual else self.vo
+        basis = self.maps["residual_basis" if residual else "state_basis"]
+        middle = einsum("ij,j->i", basis, values[offset:offset + self.k])
+        return np.concatenate((values[:offset], middle, values[offset + self.k:]))
+
+    def excluded(self, values):
+        return values - self.inject(self.restrict(values), residual=True)
+
+    def validate(self, factors, residual_factors, audit):
+        if self.identity:
+            audit.check("identity:full_mode_count", self.k == self.s)
+            return
+        b, px, qf, tx, tf = (self.maps[k] for k in MAPS)
+        audit.check("maps:dimensions", b.shape == px.shape == qf.shape == (self.s, self.k) and tx.shape == tf.shape == (self.k, self.k))
+        sx = factors[self.vo:self.vo + self.s]
+        rf = residual_factors[self.ro:self.ro + self.s]
+        close_array(audit, "maps:state_QR", sx[:, None] * b, einsum("ij,jk->ik", px, tx))
+        close_array(audit, "maps:dual_QR", b / rf[:, None], einsum("ij,jk->ik", qf, tf))
+        for name, basis, transform in (("state", px, tx), ("residual", qf, tf)):
+            close_array(audit, "maps:" + name + "_orthonormal", einsum("ij,ik->jk", basis, basis), np.eye(self.k))
+            close_array(audit, "maps:" + name + "_triangular", np.tril(transform, -1), np.zeros((self.k, self.k)), norm(transform.ravel()))
+            audit.check("maps:" + name + "_nonsingular", np.all(np.abs(np.diag(transform)) > 64 * EPS * max(norm(transform.ravel()), 1e-300)))
+
+
+def basis_package(folder, coords, physical_base, required, audit):
+    first = folder / (P + "basis_raw.txt")
+    if not first.exists():
+        raise ValueError("Missing physical harmonic basis package")
+    def get(name):
+        path = folder / (P + name + ".txt")
+        required.add(path.name)
+        return matrix(path)
+    raw, gram, orthogram, transform, samples = (get(k) for k in
+        ("basis_raw", "gram_raw", "gram_orthonormal", "coefficient_transform", "basis_samples_raw"))
+    weights_path = folder / (P + "basis_weights.txt")
+    required.add(weights_path.name)
+    weights = vector(weights_path)
+    audit.check("basis:positive_quadrature", len(weights) == len(samples) and np.all(weights > 0))
+    close_array(audit, "basis:quadrature_gram", gram, einsum("ni,n,nj->ij", samples, weights, samples))
+    close_array(audit, "basis:physical_from_raw", coords.maps["basis_physical"], einsum("ij,jk->ik", raw, transform))
+    actual = einsum("ij,ik,kl->jl", transform, gram, transform)
+    close_array(audit, "basis:transformed_gram", orthogram, actual)
+    close_array(audit, "basis:mass_orthonormal", actual, np.eye(coords.k))
+    node_path, mode_path, meta_path = (folder / (P + n + ".csv") for n in ("basis_nodes", "basis_modes", "basis_metadata"))
+    required.update(p.name for p in (node_path, mode_path, meta_path))
+    nodes, modes = rows(node_path), rows(mode_path)
+    audit.check("basis:node_indices", [integer(r["parameter"]) for r in nodes] == list(range(coords.s)))
+    node_base = np.array([finite(r["base_parameter"]) for r in nodes])
+    reference_radius = np.array([finite(r["reference_radius"]) for r in nodes])
+    physical_radius = np.array([finite(r["physical_radius"]) for r in nodes])
+    close_array(audit, "basis:frozen_base_parameters", node_base, physical_base[coords.vo:coords.vo + coords.s])
+    close_array(audit, "basis:physical_radius", physical_radius, reference_radius + node_base)
+    audit.check("basis:positive_radii", np.all(reference_radius > 0) and np.all(physical_radius > 0))
+    points = np.array([[finite(r["unit_radial_" + axis]) for axis in "xyz"] for r in nodes])
+    lengths = np.sqrt(einsum("ij,ij->i", points, points))
+    close_array(audit, "basis:unit_directions", lengths, np.ones(coords.s))
+    maximum_degree = max(integer(r["degree"]) for r in modes)
+    expected, labels = harmonics(points / lengths[:, None], maximum_degree)
+    audit.check("basis:mode_indices", [integer(r["column"]) for r in modes] == list(range(coords.k)))
+    audit.check("basis:mode_order", [(integer(r["degree"]), integer(r["m"])) for r in modes] == [(l["degree"], l["m"]) for l in labels])
+    close_array(audit, "basis:independent_harmonics", raw, expected)
+    text_keys = {"metric_scope", "scalar_finite_element_collection", "displacement_finite_element_collection"}
+    metadata = {k: r["value"] if k in text_keys else finite(r["value"])
+                for k, r in table(meta_path, lambda r: r["metric"]).items()}
+    for key, value in {"maximum_degree": maximum_degree, "modes": coords.k, "surface_parameters": coords.s,
+                       "quadrature": 14, "surface_quadrature_points": len(weights), "accepted_buffers_unchanged": 1,
+                       "prepared_state_unchanged": 1, "normalization_factors_unchanged": 1,
+                       "vector_displacement_metric_claimed": 0, "base_projection_performed": 0,
+                       "trial_preparations": 0, "radial_trace_validation_actions": 1}.items():
+        audit.check("basis:metadata:" + key, metadata[key] == value)
+    audit.check("basis:metric_scope", metadata["metric_scope"] == "scalar_radial_amplitude_FE_trace_mass_on_current_physical_surface")
+    audit.check("basis:reported_error_bounds", 0 <= metadata["maximum_gram_orthogonality_error"] <= 1e-10 and
+                0 <= metadata["maximum_sample_reconstruction_error"] <= 1e-10 and 0 <= metadata["adjoint_relative_error"] <= 1e-12)
+    for j, row in enumerate(modes):
+        pivot, relative = finite(row["pivot_norm"]), finite(row["relative_pivot"])
+        audit.check(f"basis:independent_mode:{j}", pivot > 0 and relative > metadata["rank_relative_pivot_threshold"] > 0)
+        audit.close(f"basis:pivot_transform:{j}", pivot * transform[j, j], 1.)
+        audit.close(f"basis:relative_pivot:{j}", relative, pivot / math.sqrt(gram[j, j]))
+    audit.close("basis:surface_area", metadata["physical_surface_area"], math.fsum(weights))
+    return {"status": "nodal_basis_and_saved_quadrature_gram_checked", "metadata": metadata,
+            "scope": "Saved FE trace samples and weights are checked algebraically; physical quadrature is not rerun."}
+
+
+def provenance(folder, required, audit):
+    run_path = folder.parent / (folder.name + "_run.json")
+    source_path = folder.parent / "provenance.json"
+    if not run_path.exists() and not source_path.exists():
+        return {"status": "unavailable", "bound": False}
+    if not run_path.exists() or not source_path.exists():
+        raise ValueError("Incomplete campaign provenance")
+    run, source = json.loads(run_path.read_text()), json.loads(source_path.read_text())
+    audit.check("provenance:required_hashes", required <= set(run["output_hashes"]))
+    for name, expected in run["output_hashes"].items():
+        audit.check("provenance:output:" + name, digest(folder / name) == expected)
+    for name, expected in run.get("input_state_hashes", {}).items():
+        audit.check("provenance:input:" + name, digest(name) == expected)
+    command = run.get("command", run.get("cmd"))
+    if command is None:
+        raise ValueError("Missing recorded command")
+    case = next(c for c in source["cases"] if c["id"] == folder.name)
+    expected = [str(folder.parent / "plateau_diagnostic"), *case["arguments"], "--mesh", source["mesh"], "--output", str(folder)]
+    audit.check("provenance:command", command == expected)
+    audit.check("provenance:binary", digest(command[0]) == source["binary_sha256"])
+    audit.check("provenance:mesh", digest(source["mesh"]) == source["mesh_sha256"])
+    if "--probe-state" in command:
+        path = Path(command[command.index("--probe-state") + 1])
+        if not path.is_absolute():
+            path = Path(source["source"]) / path
+        audit.check("provenance:base_input_bound", str(path.resolve()) in run.get("input_state_hashes", {}))
+        close_array(audit, "provenance:physical_base_input", vector(path), vector(folder / (P + "physical_base.txt")))
+    audit.check("provenance:successful_probe_exit", run["exit_code"] == 0)
+    return {"status": "runner_hashes_checked", "bound": True, "source_commit": source["commit"],
+            "run_sha256": digest(run_path), "provenance_sha256": digest(source_path)}
+
+
+def inspect(folder):
+    folder = Path(folder).resolve()
+    audit, required, observations = Audit(), set(), {}
+    try:
+        def get(name):
+            path = folder / (P + name + ".txt")
+            required.add(path.name)
+            return vector(path)
+        summary_path = folder / (P + "summary.csv")
+        trials_path = folder / (P + "trials.csv")
+        layout_path = folder / (P + "layout.csv")
+        required.update(p.name for p in (summary_path, trials_path, layout_path))
+        summary = {k: finite(r["value"]) for k, r in table(summary_path, lambda r: r["metric"]).items()}
+        audit.check("identity:boolean", summary["identity_control"] in (0, 1))
+        maps = {}
+        for name in MAPS:
+            path = folder / (P + name + ".txt")
+            if summary["identity_control"]:
+                audit.check("identity:matrix_absent:" + name, not path.exists())
+            else:
+                required.add(path.name)
+                maps[name] = matrix(path)
+        coordinates = Coordinates(summary, maps)
+        n, s, k, vo, ro = coordinates.n, coordinates.s, coordinates.k, coordinates.vo, coordinates.ro
+        degree = integer(summary["maximum_degree"])
+        audit.check("prescribed_basis_size", (degree == 0 and k == s) if coordinates.identity else (degree in (2, 4, 6, 8) and k == (degree + 1) ** 2))
+        audit.check("reduced_size", summary["reduced_size"] == n - s + k)
+        audit.check("schema_version", summary["schema_version"] == 1)
+        audit.check("reported_algebra_selfcheck", summary["algebra_self_check_passed"] == 1)
+        blocks = layout(layout_path, n)
+        audit.check("layout:value_shape", ("surface_deformation", vo, s) in blocks["value"])
+        audit.check("layout:residual_shape", ("surface_shape_balance", ro, s) in blocks["residual"])
+        data = {name: get(name) for name in VECTORS}
+        for name, values in data.items():
+            audit.check("size:" + name, len(values) == (n - s + k if name.startswith("reduced_") else n))
+        sx, rf = data["state_factors"], data["residual_factors"]
+        audit.check("normalization:positive", np.all(sx > 0) and np.all(rf > 0))
+        coordinates.validate(sx, rf, audit)
+        for key, value in {"accepted_buffers_bitwise_identical": 1, "normalization_factors_bitwise_identical": 1,
+                           "direction_buffers_bitwise_identical": 1, "restored_residual_error": 0,
+                           "restored_action_error": 0, "probe_exception": 0, "accepted_commits": 0}.items():
+            audit.check("invariant:" + key, summary[key] == value)
+        audit.check("restoration:residual_exact", np.array_equal(data["residual"], data["restored_residual"]))
+        audit.check("restoration:action_exact", np.array_equal(data["action"], data["restored_action"]))
+        close_array(audit, "normalization:base", data["physical_base"], data["base"] / sx)
+        close_array(audit, "normalization:direction", data["physical_direction"], data["direction"] / sx)
+        close_array(audit, "injection:direction", data["direction"], coordinates.inject(data["reduced_direction"]))
+        for name in ("residual", "action", "defect"):
+            # Forming R(F+Jp) and RF+RJp rounds differently when F and Jp
+            # cancel. The allowance must scale with the operands, not only
+            # with the potentially zero defect being checked.
+            scale = norm(data["residual"]) + norm(data["action"]) if name == "defect" else 0.
+            close_array(audit, "restriction:" + name, data["reduced_" + name], coordinates.restrict(data[name]), scale)
+        close_array(audit, "defect:full", data["defect"], data["residual"] + data["action"], norm(data["residual"]) + norm(data["action"]))
+        close_array(audit, "defect:reduced", data["reduced_defect"], data["reduced_residual"] + data["reduced_action"],
+                    norm(data["reduced_residual"]) + norm(data["reduced_action"]))
+        for name in ("residual", "defect"):
+            excluded = coordinates.excluded(data[name])
+            close_array(audit, "excluded:" + name, data["excluded_" + name], excluded, norm(data[name]))
+            audit.close("pythagoras:" + name, norm(data[name]) ** 2,
+                        norm(data["reduced_" + name]) ** 2 + norm(excluded) ** 2, norm(data[name]) ** 2)
+        eta = summary["linear_relative_tolerance"]
+        absolute = summary["linear_absolute_tolerance"]
+        audit.check("linear:controls", 0 < eta < 1 and absolute >= 0 and integer(summary["linear_maximum_iterations"]) > 0)
+        gates, slopes = {}, {}
+        for label, prefix in (("full", ""), ("reduced", "reduced_")):
+            fn, dn = norm(data[prefix + "residual"]), norm(data[prefix + "defect"])
+            audit.check(label + ":nonzero_rhs", fn > 0)
+            passed = dn <= max(absolute, eta * fn)
+            audit.check(label + ":reported_gate", summary[label + "_verified_gate"] in (0, 1) and bool(summary[label + "_verified_gate"]) == passed)
+            audit.close(label + ":relative_defect", summary[label + "_relative_defect"], dn / fn)
+            gates[label] = {"verified_gate": passed, "residual_norm": fn, "defect_norm": dn, "relative_defect": dn / fn}
+            slopes[label] = float(einsum("i,i->", data[prefix + "residual"], data[prefix + "action"]))
+            audit.close(label + ":reported_norm", summary[label + "_base_norm"], fn)
+            audit.close(label + ":reported_slope", summary[label + "_merit_slope"], slopes[label], fn * norm(data[prefix + "action"]))
+        audit.close("reported_direction_norm", summary["direction_norm"], norm(data["direction"]))
+        audit.close("reported_excluded_base", summary["excluded_base_norm"], norm(data["excluded_residual"]), norm(data["residual"]))
+        audit.close("linear:reported_true_norm", summary["linear_reported_true_norm"], norm(data["reduced_defect"]), norm(data["reduced_residual"]))
+        audit.check("linear:iteration_budget", 0 <= integer(summary["linear_iterations"]) <= integer(summary["linear_maximum_iterations"]))
+        trials = rows(trials_path)
+        audit.check("trial:complete_grid", len(trials) == len(ALPHAS))
+        trial_observations = []
+        for index, row in enumerate(trials):
+            alpha = finite(row["alpha"])
+            audit.check(f"trial:{index}:alpha", index < len(ALPHAS) and alpha == ALPHAS[index])
+            audit.check(f"trial:{index}:sample", integer(row["sample"]) == index)
+            prepared = integer(row["prepared"])
+            audit.check(f"trial:{index}:prepared_flag", prepared in (0, 1))
+            result = {"sample": index, "alpha": alpha, "prepared": bool(prepared), "reason": integer(row["reason"]), "stage": integer(row["stage"])}
+            candidate_path = folder / (P + f"candidate_a{index}.txt")
+            if prepared or candidate_path.exists():
+                candidate = get(f"candidate_a{index}")
+                close_array(audit, f"trial:{index}:candidate", candidate, data["base"] + alpha * data["direction"])
+            if prepared:
+                audit.check(f"trial:{index}:prepared_status", result["reason"] == result["stage"] == -1)
+                physical_candidate = get(f"physical_candidate_a{index}")
+                close_array(audit, f"trial:{index}:physical_candidate", physical_candidate, candidate / sx)
+                full, reduced = get(f"residual_a{index}"), get(f"reduced_residual_a{index}")
+                close_array(audit, f"trial:{index}:restriction", reduced, coordinates.restrict(full))
+                excluded = coordinates.excluded(full)
+                full_norm, reduced_norm, excluded_norm = norm(full), norm(reduced), norm(excluded)
+                remainder = full - data["residual"] - alpha * data["action"]
+                reduced_remainder = reduced - data["reduced_residual"] - alpha * data["reduced_action"]
+                for label, actual in (("full_norm", full_norm), ("reduced_norm", reduced_norm), ("excluded_shape_norm", excluded_norm),
+                                      ("full_remainder_norm", norm(remainder)), ("reduced_remainder_norm", norm(reduced_remainder))):
+                    audit.close(f"trial:{index}:" + label, finite(row[label]), actual, norm(data["residual"]))
+                audit.close(f"trial:{index}:pythagoras", full_norm ** 2, reduced_norm ** 2 + excluded_norm ** 2, full_norm ** 2)
+                for label, actual in (("full", full_norm), ("reduced", reduced_norm)):
+                    accepted = actual <= (1 - 1e-4 * alpha) * gates[label]["residual_norm"]
+                    audit.check(f"trial:{index}:" + label + "_accept", integer(row[label + "_merit_pass"]) in (0, 1) and bool(integer(row[label + "_merit_pass"])) == accepted)
+                    result[label + "_merit_pass"] = accepted
+                result.update(full_norm=full_norm, reduced_norm=reduced_norm, excluded_shape_norm=excluded_norm,
+                              full_remainder_norm=norm(remainder), reduced_remainder_norm=norm(reduced_remainder))
+            else:
+                audit.check(f"trial:{index}:failure_status", result["reason"] >= 0 and result["stage"] >= 0)
+                for name in (f"residual_a{index}", f"reduced_residual_a{index}"):
+                    audit.check(f"trial:{index}:unavailable:" + name, not (folder / (P + name + ".txt")).exists())
+                for label in ("full_norm", "reduced_norm", "excluded_shape_norm", "full_remainder_norm", "reduced_remainder_norm"):
+                    audit.check(f"trial:{index}:unavailable:" + label, math.isnan(float(row[label])))
+                audit.check(f"trial:{index}:failure_not_accepted", integer(row["full_merit_pass"]) == integer(row["reduced_merit_pass"]) == 0)
+            trial_observations.append(result)
+        package = basis_package(folder, coordinates, data["physical_base"], required, audit) if not coordinates.identity else {"status": "identity_control"}
+        if not coordinates.identity:
+            audit.check("basis:degree_matches_summary", package["metadata"]["maximum_degree"] == degree)
+        binding = provenance(folder, required, audit)
+        observations = {"summary": summary, "linear": gates, "merit_directional_derivatives": slopes,
+                        "classification": "full_and_reduced_gates_pass" if gates["full"]["verified_gate"] and gates["reduced"]["verified_gate"] else
+                        "reduced_gate_only_not_full_Newton" if gates["reduced"]["verified_gate"] else "finite_reduced_solve_failed_gate",
+                        "base_excluded_shape_norm": norm(data["excluded_residual"]), "trials": trial_observations,
+                        "basis": package, "provenance": binding}
+    except (ValueError, OSError, KeyError, TypeError, IndexError, StopIteration, np.linalg.LinAlgError) as error:
+        audit.check("complete_readable_consistent_artifacts:" + str(error), False)
+    return {"directory": str(folder), "artifact_pass": not audit.failures, "checks_count": audit.count,
+            "failures": audit.failures, "required_artifacts": sorted(required), "observations": observations if not audit.failures else {}}
+
+
+def self_test():
+    count = 0
+    def checked(condition):
+        nonlocal count
+        assert condition
+        count += 1
+    # Deliberately different value/residual factors and shape offsets catch the
+    # tempting but invalid reuse of the state basis as a residual test basis.
+    b = np.array([[1., .4], [1., -.1], [1., -.3]])
+    sx = np.array([2., 3., 5., 7., 11., 13., 17.])
+    rf = np.array([.2, .3, .7, 1.1, 1.3, 1.7, 1.9])
+    summary = dict(full_size=7, shape_size=3, mode_count=2, shape_value_offset=1, shape_residual_offset=3, identity_control=0)
+    px, tx = np.linalg.qr(sx[1:4, None] * b)
+    qf, tf = np.linalg.qr(b / rf[3:6, None])
+    maps = dict(zip(MAPS, (b, px, qf, tx, tf)))
+    c = Coordinates(summary, maps)
+    a = Audit(); c.validate(sx, rf, a); checked(not a.failures)
+    z = np.arange(1., 7.)
+    f = np.arange(-3., 4.)
+    checked(vector_near(c.restrict(c.inject(z), False), z))
+    checked(vector_near(c.restrict(c.inject(z, True)), z))
+    checked(near(norm(f) ** 2, norm(c.restrict(f)) ** 2 + norm(c.excluded(f)) ** 2))
+    physical_dual = einsum("ij,i->j", b, f[3:6] / rf[3:6])
+    checked(vector_near(einsum("ji,j->i", tf, c.restrict(f)[3:5]), physical_dual))
+    bad = Coordinates(summary, {**maps, "residual_basis": px})
+    a = Audit(); bad.validate(sx, rf, a); checked("maps:dual_QR" in a.failures)
+    a = Audit(); c.validate(sx, rf[::-1], a); checked(bool(a.failures))
+    checked(norm(c.restrict(c.excluded(f))) < 1e-14)
+    # Full-rank bases may rotate the coordinates. Verify the recovered full
+    # direction, not equality of the rotated coefficient arrays.
+    bfull = np.array([[1., .4, .2], [1., -.1, .7], [1., -.3, -.6]])
+    px, tx = np.linalg.qr(sx[1:4, None] * bfull)
+    qf, tf = np.linalg.qr(bfull / rf[3:6, None])
+    full = Coordinates({**summary, "mode_count": 3}, dict(zip(MAPS, (bfull, px, qf, tx, tf))))
+    op = np.diag(np.arange(2., 9.)) + .03 * np.ones((7, 7))
+    columns = []
+    for i in range(7):
+        columns.append(full.restrict(einsum("ij,j->i", op, full.inject(np.eye(7)[:, i]))))
+    reduced_op = np.column_stack(columns)
+    result = full.inject(np.linalg.solve(reduced_op, -full.restrict(f)))
+    checked(vector_near(result, np.linalg.solve(op, -f)))
+    identity = Coordinates({**summary, "mode_count": 3, "identity_control": 1}, {})
+    checked(np.array_equal(identity.inject(f), f) and np.array_equal(identity.restrict(f), f))
+    checked(np.array_equal(identity.excluded(f), np.zeros(7)))
+    # A reduced zero defect can coexist with a substantial discarded full one.
+    discarded = c.excluded(f)
+    checked(norm(discarded) > .1 and norm(c.restrict(discarded)) < 1e-14)
+    import csv
+    import tempfile
+    def fixture(folder, identity=False):
+        folder.mkdir()
+        def save(name, array):
+            array = np.asarray(array)
+            with (folder / (P + name + ".txt")).open("w") as stream:
+                stream.write((f"{len(array)}" if array.ndim == 1 else f"{array.shape[0]} {array.shape[1]}") + "\n")
+                np.savetxt(stream, array.ravel(), fmt="%.17g")
+        def csv_save(name, columns, records):
+            with (folder / (P + name + ".csv")).open("w", newline="") as stream:
+                writer = csv.writer(stream); writer.writerow(columns); writer.writerows(records)
+        n, s, k, vo, ro = 18, 12, 9, 1, 4
+        points = []
+        for i in range(s):
+            z = 1 - 2 * (i + .5) / s
+            phi = i * math.pi * (3 - math.sqrt(5))
+            points.append((math.sqrt(1 - z * z) * math.cos(phi), math.sqrt(1 - z * z) * math.sin(phi), z))
+        points = np.asarray(points)
+        raw, modes = harmonics(points, 2)
+        weights = np.ones(s)
+        gram = einsum("ni,n,nj->ij", raw, weights, raw)
+        transform = np.linalg.inv(np.linalg.cholesky(gram).T)
+        b = einsum("ij,jk->ik", raw, transform)
+        sx, rf = np.linspace(1., 3., n), np.linspace(.4, 2., n)
+        px, tx = np.linalg.qr(sx[vo:vo + s, None] * b)
+        qf, tf = np.linalg.qr(b / rf[ro:ro + s, None])
+        maps = dict(zip(MAPS, (b, px, qf, tx, tf)))
+        summary = dict(schema_version=1, full_size=n, shape_size=s, mode_count=s if identity else k,
+                       shape_value_offset=vo, shape_residual_offset=ro, identity_control=int(identity),
+                       maximum_degree=0 if identity else 2, reduced_size=n if identity else n - s + k,
+                       linear_relative_tolerance=.03, linear_absolute_tolerance=0., linear_maximum_iterations=200,
+                       linear_restart=160, linear_status=0, linear_iterations=7, linear_seconds=.01,
+                       algebra_self_check_passed=1, accepted_commits=0, accepted_buffers_bitwise_identical=1,
+                       normalization_factors_bitwise_identical=1, direction_buffers_bitwise_identical=1,
+                       restored_residual_error=0, restored_action_error=0, preparations=8, probe_exception=0)
+        c = Coordinates(summary, {} if identity else maps)
+        reduced_direction = np.sin(np.arange(summary["reduced_size"]) + 1.)
+        direction = c.inject(reduced_direction)
+        action = direction.copy()
+        residual = -action + c.excluded(np.cos(np.arange(n)))
+        base = np.linspace(.1, .3, n)
+        data = dict(base=base, physical_base=base / sx, residual=residual, direction=direction,
+                    physical_direction=direction / sx, action=action, defect=residual + action,
+                    state_factors=sx, residual_factors=rf, reduced_direction=reduced_direction,
+                    restored_residual=residual, restored_action=action)
+        for name in ("residual", "action", "defect"):
+            data["reduced_" + name] = c.restrict(data[name])
+        for name in ("residual", "defect"):
+            data["excluded_" + name] = c.excluded(data[name])
+        for label, prefix in (("full", ""), ("reduced", "reduced_")):
+            fn, dn = norm(data[prefix + "residual"]), norm(data[prefix + "defect"])
+            summary[label + "_base_norm"] = fn
+            summary[label + "_relative_defect"] = dn / fn
+            summary[label + "_verified_gate"] = int(dn <= .03 * fn)
+            summary[label + "_merit_slope"] = float(einsum("i,i->", data[prefix + "residual"], data[prefix + "action"]))
+        summary["linear_reported_true_norm"] = norm(data["reduced_defect"])
+        summary["direction_norm"] = norm(direction)
+        summary["excluded_base_norm"] = norm(data["excluded_residual"])
+        for name, values in data.items():
+            save(name, values)
+        csv_save("summary", ["metric", "value"], summary.items())
+        csv_save("layout", ["kind", "block", "offset", "size"],
+                 [("value", "density", 0, vo), ("value", "surface_deformation", vo, s), ("value", "tail", vo + s, n - vo - s),
+                  ("residual", "head", 0, ro), ("residual", "surface_shape_balance", ro, s), ("residual", "tail", ro + s, n - ro - s)])
+        trials = []
+        for index, alpha in enumerate(ALPHAS):
+            save(f"candidate_a{index}", base + alpha * direction)
+            if index == 0:
+                trials.append((index, alpha, 0, 0, 0, "nan", "nan", "nan", 0, 0, "nan", "nan"))
+                continue
+            f = residual + alpha * action + .025 * alpha ** 2 * residual
+            fr = c.restrict(f)
+            remainder = f - residual - alpha * action
+            save(f"physical_candidate_a{index}", (base + alpha * direction) / sx)
+            save(f"residual_a{index}", f); save(f"reduced_residual_a{index}", fr)
+            trials.append((index, alpha, 1, -1, -1, norm(f), norm(fr), norm(c.excluded(f)),
+                           int(norm(f) <= (1 - 1e-4 * alpha) * norm(residual)),
+                           int(norm(fr) <= (1 - 1e-4 * alpha) * norm(data["reduced_residual"])),
+                           norm(remainder), norm(c.restrict(remainder))))
+        csv_save("trials", ["sample", "alpha", "prepared", "reason", "stage", "full_norm", "reduced_norm", "excluded_shape_norm",
+                            "full_merit_pass", "reduced_merit_pass", "full_remainder_norm", "reduced_remainder_norm"], trials)
+        if not identity:
+            for name, values in maps.items():
+                save(name, values)
+            for name, values in {"basis_raw": raw, "gram_raw": gram, "coefficient_transform": transform,
+                                 "gram_orthonormal": einsum("ij,ik,kl->jl", transform, gram, transform),
+                                 "basis_samples_raw": raw, "basis_weights": weights}.items():
+                save(name, values)
+            csv_save("basis_nodes", ["parameter", "unit_radial_x", "unit_radial_y", "unit_radial_z", "reference_radius", "physical_radius", "base_parameter"],
+                     [(i, *point, 1., 1 + data["physical_base"][vo + i], data["physical_base"][vo + i]) for i, point in enumerate(points)])
+            csv_save("basis_modes", ["column", "degree", "m", "pivot_norm", "relative_pivot"],
+                     [(i, m["degree"], m["m"], 1 / transform[i, i], 1 / (transform[i, i] * math.sqrt(gram[i, i]))) for i, m in enumerate(modes)])
+            metadata = dict(maximum_degree=2, modes=k, surface_parameters=s, quadrature=14, surface_quadrature_points=s,
+                            physical_surface_area=s, accepted_buffers_unchanged=1, prepared_state_unchanged=1, normalization_factors_unchanged=1,
+                            vector_displacement_metric_claimed=0, base_projection_performed=0, trial_preparations=0,
+                            radial_trace_validation_actions=1, metric_scope="scalar_radial_amplitude_FE_trace_mass_on_current_physical_surface",
+                            maximum_gram_orthogonality_error=1e-15, maximum_sample_reconstruction_error=0.,
+                            adjoint_relative_error=0., rank_relative_pivot_threshold=1e-12)
+            csv_save("basis_metadata", ["metric", "value"], metadata.items())
+        return save
+    with tempfile.TemporaryDirectory() as tmp:
+        folder = Path(tmp) / "reduced"
+        save = fixture(folder)
+        result = inspect(folder)
+        assert result["artifact_pass"], result["failures"]
+        checked(result["artifact_pass"])
+        checked(result["observations"]["classification"] == "reduced_gate_only_not_full_Newton")
+        for name in ("residual_basis", "candidate_a1", "excluded_residual", "gram_raw", "state_factors"):
+            path = folder / (P + name + ".txt")
+            original = path.read_bytes()
+            values = matrix(path) if name in ("residual_basis", "gram_raw") else vector(path)
+            values.flat[0] += .02
+            save(name, values)
+            checked(not inspect(folder)["artifact_pass"])
+            path.write_bytes(original)
+        summary_path = folder / (P + "summary.csv")
+        original = summary_path.read_text()
+        summary_path.write_text(original.replace("full_verified_gate,0", "full_verified_gate,1"))
+        checked(not inspect(folder)["artifact_pass"])
+        summary_path.write_text(original)
+        save("residual_a0", np.ones(18))
+        checked(not inspect(folder)["artifact_pass"])
+        (folder / (P + "residual_a0.txt")).unlink()
+        identity_folder = Path(tmp) / "identity"
+        fixture(identity_folder, True)
+        result = inspect(identity_folder)
+        checked(result["artifact_pass"] and result["observations"]["classification"] == "full_and_reduced_gates_pass")
+    print(f"{count} tiny algebra, end-to-end artifact and corruption checks passed")
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("directories", nargs="*")
+    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 probe directories and a fresh output JSON")
+    output = args.output.resolve()
+    if output.exists() or any(output.is_relative_to(Path(d).resolve()) for d in args.directories):
+        parser.error("Use a fresh output outside 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(__file__),
+              "tolerances": {"relative": RTOL, "cancellation_epsilon_multiplier": 2048}, "cases": cases,
+              "scope": "Independent saved-vector and map algebra only. Reduced success does not certify the unrestricted equations, FEM evaluation, or physical validity."}
+    output.parent.mkdir(parents=True, exist_ok=True)
+    output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n")
+    print(f"{'PASS' if result['passed'] else 'FAIL'}: {len(cases)} cases; {output}")
+    return 0 if result["passed"] else 1
+
+
+if __name__ == "__main__":
+    sys.exit(main())

experiments/plateau_diagnostic.cpp

@@ -31,6 +31,7 @@ import mean_field;
 #include "plateau_positivity_strategy_diagnostics.hpp"
 #include "plateau_polynomial_certificate.hpp"
 #include "shape_field_weak_diagnostics.hpp"
+#include "reduced_surface_diagnostics.hpp"
 
 using namespace mean_field;
 namespace {
@@ -42,6 +43,8 @@ struct Options {
     bool strategyProbe = false, positivityProbe = false, positivityCertificate = false;
     bool preconditionerProbe = false, nonlinearStepProbe = false;
     bool weakShapeProbe = false, weakShapeFinal = false;
+    bool reducedSurfaceProbe = false, reducedSurfaceAlgebraCheck = false;
+    int surfaceDegree = 4;
     bool coupledTrialProbe = false, coupledSeedControl = false;
     bool dampedCoupledProbe = false;
     std::string coupledSearch = "near";
@@ -115,6 +118,12 @@ Options Parse(int argc, char **argv) {
             o.preconditionerProbe = true;
         else if (a == "--nonlinear-step-probe")
             o.nonlinearStepProbe = true;
+        else if (a == "--reduced-surface-probe")
+            o.reducedSurfaceProbe = true;
+        else if (a == "--reduced-surface-algebra-check")
+            o.reducedSurfaceAlgebraCheck = true;
+        else if (a == "--surface-degree")
+            o.surfaceDegree = std::stoi(v());
         else if (a == "--weak-shape-probe")
             o.weakShapeProbe = true;
         else if (a == "--weak-shape-final")
@@ -177,6 +186,18 @@ Options Parse(int argc, char **argv) {
     }
     if (o.output.empty() || o.n < 1 || o.n >= 5 || o.maxNewton < 0 || o.maxLinear < 1 || o.restart < 1)
         throw std::invalid_argument("Invalid options");
+    if (o.reducedSurfaceProbe || o.reducedSurfaceAlgebraCheck) {
+        if ((o.reducedSurfaceProbe && o.reducedSurfaceAlgebraCheck) || o.maxNewton != 0 ||
+            o.strategy != "baseline" || o.probes || o.pressureProbe || o.shapeCurvature ||
+            o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.positivityCertificate ||
+            o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe ||
+            o.weakShapeProbe || o.weakShapeFinal || !o.probeDirection.empty() ||
+            !o.replayCandidate.empty() || !o.captureIterations.empty() || !o.noCentral ||
+            (o.surfaceDegree != 0 && o.surfaceDegree != 2 && o.surfaceDegree != 4 &&
+             o.surfaceDegree != 6 && o.surfaceDegree != 8) ||
+            !std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.)
+            throw std::invalid_argument("Reduced-surface study requires a standalone unbordered frozen probe; degree 0 is the full control");
+    }
     if (o.weakShapeProbe || o.weakShapeFinal) {
         if ((o.weakShapeProbe && o.weakShapeFinal) || o.probes || o.pressureProbe || o.shapeCurvature ||
             o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.positivityCertificate ||
@@ -234,8 +255,8 @@ Options Parse(int argc, char **argv) {
             throw std::invalid_argument("Newton preconditioner variants are independent, gravity-schur, stellar-ldu, or stellar-ldu-gravity-schur");
         if (o.preconditionerVariant != "independent") {
             if (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe ||
-                o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe) ||
+                o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe && !o.reducedSurfaceProbe) ||
-                (!o.probeDirection.empty() && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe) ||
+                (!o.probeDirection.empty() && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe && !o.reducedSurfaceProbe) ||
                 !o.replayCandidate.empty())
                 throw std::invalid_argument("A custom Newton preconditioner cannot be combined with frozen or certificate probes");
             if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" &&
@@ -529,7 +550,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
              << "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart
              << "\nfrozen_probe_backend=" << (o.preconditionerProbe ? "direct_mfem_fgmres_preconditioner_benchmark" :
                  (o.coupledTrialProbe || o.dampedCoupledProbe) ? "supplied_outer_direction_verified_with_full_J_and_auxiliary_fgmres" :
-                 (o.nonlinearStepProbe || ((o.weakShapeProbe || o.weakShapeFinal) && o.probeDirection.empty()))
+                 (o.reducedSurfaceProbe || o.nonlinearStepProbe || ((o.weakShapeProbe || o.weakShapeFinal) && o.probeDirection.empty()))
                      ? "prepared_fgmres_cold_with_independent_true_gate" :
                  o.positivityCertificate ? "polynomial_certificate_no_linear_solve" :
                  o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve")
@@ -539,6 +560,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
              << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart
              << "\npreconditioner_probe=" << o.preconditionerProbe
              << "\nnonlinear_step_probe=" << o.nonlinearStepProbe
+             << "\nreduced_surface_probe=" << o.reducedSurfaceProbe << "\nsurface_degree=" << o.surfaceDegree
              << "\nweak_shape_probe=" << o.weakShapeProbe << "\nweak_shape_final=" << o.weakShapeFinal
              << "\ndamped_coupled_probe=" << o.dampedCoupledProbe
              << "\ncoupled_search=" << o.coupledSearch
@@ -564,6 +586,23 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
                  << "weak_shape_operator=full_monolithic\n"
                  << "weak_shape_trial_policy=raw_prepare_only\n";
     metadata.flush();
+    if (o.reducedSurfaceProbe) {
+        metadata << "reduced_surface_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n'
+                 << "reduced_surface_trial_policy=raw_prepare_only_no_commit\n"
+                 << "probe_state=" << o.probeState << '\n';
+        metadata.flush();
+        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
+            mfem::Vector base(s.acceptedNormalizedState);
+            if (!o.probeState.empty()) {
+                const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
+                s.normalizedOperator->NormalizeState(physical, base);
+            }
+            experiment::plateau_diagnostics::ProbeReducedSurface(s, f, base, o.output, o.surfaceDegree, o.n, scales.K, {
+                .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0.,
+                           .maximumIterations = o.maxLinear}, .restartLength = o.restart});
+        });
+        return 0;
+    }
     if (o.weakShapeProbe) {
         metadata << "weak_shape_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n'
                  << "probe_state=" << o.probeState << "\nprobe_direction=" << o.probeDirection << '\n';
@@ -1088,7 +1127,14 @@ int main(int argc, char **argv) {
             throw std::runtime_error("Single rank only");
         mfem::Device device("cpu");
         std::cout << std::setprecision(12);
-        result = Run(Parse(argc, argv));
+        const auto options = Parse(argc, argv);
+        if (options.reducedSurfaceAlgebraCheck) {
+            const auto checks = experiment::plateau_diagnostics::RunReducedSurfaceAlgebraSelfCheck(MPI_COMM_WORLD);
+            std::cout << "Reduced surface algebra self-check passed: " << checks.checks
+                      << " checks, maximum error " << checks.maximumError << std::endl;
+        } else {
+            result = Run(options);
+        }
     } catch (const std::exception &e) {
         std::cerr << "Diagnostic failure: " << e.what() << std::endl;
         result = 2;

experiments/reduced_surface_basis.hpp

@@ -0,0 +1,339 @@
+#pragma once
+
+// Include after import mean_field. This is an experiment-local coordinate
+// construction, not a change to the stellar model or its deformation policy.
+#include <algorithm>
+#include <array>
+#include <chrono>
+#include <cmath>
+#include <cstring>
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <limits>
+#include <numbers>
+#include <stdexcept>
+#include <string>
+#include <string_view>
+#include <vector>
+
+#include <mfem.hpp>
+#include <mpi.h>
+#include "plateau_coupled_coarse.hpp"
+
+namespace experiment::plateau_diagnostics {
+
+struct ReducedSurfaceBasis {
+    mfem::DenseMatrix physicalBasis;
+    mfem::DenseMatrix rawPhysicalBasis;
+    mfem::DenseMatrix rawGram;
+    mfem::DenseMatrix orthonormalGram;
+    mfem::DenseMatrix coefficientTransform;
+    int maximumDegree{};
+    double physicalSurfaceArea{};
+    double maximumGramOrthogonalityError{};
+};
+
+namespace reduced_surface_basis_detail {
+inline void Require(bool condition, const std::string &message) {
+    if (!condition) throw std::runtime_error("Reduced surface basis: " + message);
+}
+
+inline bool Identical(const mfem::Vector &a, const mfem::Vector &b) {
+    return a.Size() == b.Size() && (a.Size() == 0 ||
+        std::memcmp(a.GetData(), b.GetData(), sizeof(mfem::real_t) * a.Size()) == 0);
+}
+
+inline std::ofstream Output(const std::filesystem::path &directory, std::string_view name,
+                            std::string_view extension = ".txt") {
+    const auto path = directory / ("reduced_surface_" + std::string(name) + std::string(extension));
+    Require(!std::filesystem::exists(path), "refusing to overwrite " + path.string());
+    std::ofstream stream(path);
+    stream.exceptions(std::ios::badbit | std::ios::failbit);
+    stream << std::setprecision(std::numeric_limits<double>::max_digits10);
+    return stream;
+}
+
+inline void Save(const std::filesystem::path &directory, std::string_view name,
+                 const mfem::DenseMatrix &matrix) {
+    auto stream = Output(directory, name);
+    stream << matrix.Height() << ' ' << matrix.Width() << '\n';
+    for (int i = 0; i < matrix.Height(); ++i)
+        for (int j = 0; j < matrix.Width(); ++j) {
+            Require(std::isfinite(matrix(i, j)), "nonfinite matrix output");
+            stream << matrix(i, j) << '\n';
+        }
+}
+
+// The normalized recurrence and signed-m convention match the existing
+// analyze_shape_surface_modes.py analysis, including Condon--Shortley phase.
+inline std::vector<double> Harmonics(int degree, double x, double y, double z) {
+    std::array<std::array<double, 9>, 9> p{};
+    const double phi = std::atan2(y, x), rxy = std::hypot(x, y);
+    p[0][0] = 1. / std::sqrt(4. * std::numbers::pi);
+    for (int m = 0; m <= degree; ++m) {
+        if (m > 0) p[m][m] = -std::sqrt((2. * m + 1.) / (2. * m)) * rxy * p[m - 1][m - 1];
+        if (m < degree) p[m + 1][m] = std::sqrt(2. * m + 3.) * z * p[m][m];
+        for (int l = m + 2; l <= degree; ++l) {
+            const double a = std::sqrt((4. * l * l - 1.) / (l * l - m * m));
+            const double b = std::sqrt(((2. * l + 1.) * ((l - 1.) * (l - 1.) - m * m)) /
+                                       ((2. * l - 3.) * (l * l - m * m)));
+            p[l][m] = a * z * p[l - 1][m] - b * p[l - 2][m];
+        }
+    }
+    std::vector<double> values;
+    values.reserve((degree + 1) * (degree + 1));
+    for (int l = 0; l <= degree; ++l)
+        for (int m = -l; m <= l; ++m) {
+            double value = p[l][std::abs(m)];
+            if (m < 0) value *= std::sqrt(2.) * std::sin(-m * phi);
+            if (m > 0) value *= std::sqrt(2.) * std::cos(m * phi);
+            Require(std::isfinite(value), "nonfinite harmonic value");
+            values.push_back(value);
+        }
+    return values;
+}
+
+inline long double ColumnDot(const mfem::DenseMatrix &a, int i,
+                             const mfem::DenseMatrix &b, int j) {
+    long double value = 0.;
+    for (int row = 0; row < a.Height(); ++row)
+        value += static_cast<long double>(a(row, i)) * b(row, j);
+    return value;
+}
+
+struct Face {
+    int boundary{};
+    mfem::Array<int> scalarDofs;
+    mfem::Vector displacement;
+    mfem::DenseMatrix coefficients;
+};
+} // namespace reduced_surface_basis_detail
+
+template <class State>
+ReducedSurfaceBasis MakeReducedSurfaceBasis(State &state, const mean_field::fem::FEM &fem,
+                                             int maximumDegree, const std::filesystem::path &output) {
+    namespace d = reduced_surface_basis_detail;
+    using Schema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema;
+    const auto start = std::chrono::steady_clock::now();
+    d::Require(maximumDegree == 2 || maximumDegree == 4 || maximumDegree == 6 || maximumDegree == 8,
+               "maximum degree must be one of 2, 4, 6, 8");
+    const auto &problem = state.Problem();
+    int ranks = 0;
+    MPI_Comm_size(problem.GetCommunicator(), &ranks);
+    d::Require(ranks == 1, "this bounded experiment supports exactly one MPI rank");
+    d::Require(state.normalizedOperator->IsPrepared() && problem.GetPhysicalOperator().IsPrepared(),
+               "a prepared frozen base is required");
+    const auto generation = problem.GetPreparationGeneration();
+    const mfem::Vector accepted(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual),
+        acceptedPhysical(state.AcceptedPhysicalState()), trial(state.trialNormalizedState),
+        trialResidual(state.trialNormalizedResidual), candidatePhysical(state.candidatePhysicalState),
+        preparedPhysical(state.normalizedOperator->GetPhysicalState()),
+        preparedResidual(state.normalizedOperator->GetPhysicalResidual()),
+        stateFactors(state.normalizedOperator->GetNormalization().StateFactors()),
+        residualFactors(state.normalizedOperator->GetNormalization().ResidualFactors());
+    const auto &physical = problem.GetPhysicalOperator();
+    const mfem::Vector baseDisplacement(physical.GetGeneratedVolumeDisplacement());
+    // This verifies the actual production radial trace, rather than relying
+    // only on its erased runtime descriptor. It does not prepare any state.
+    const auto radial = coupled_shape_detail::MakeDefaultRadialSurface(problem);
+    const int n = radial.parameterCount(), k = (maximumDegree + 1) * (maximumDegree + 1);
+    d::Require(n >= k && fem.mesh && fem.surfaceDeformationFes && fem.displacementFes,
+               "incomplete spaces or fewer surface parameters than requested modes");
+    d::Require(!fem.surfaceDeformationFes->Nonconforming() &&
+               !fem.surfaceDeformationFes->IsVariableOrder(), "conforming fixed-order scalar traces are required");
+    const auto &parameters = physical.GetSurfaceDeformationParameters();
+    d::Require(parameters.Size() == n, "surface parameter ordering mismatch");
+    ReducedSurfaceBasis result;
+    result.maximumDegree = maximumDegree;
+    result.rawPhysicalBasis.SetSize(n, k);
+    for (int i = 0; i < n; ++i) {
+        const auto values = d::Harmonics(maximumDegree, radial.radialDirection(i, 0),
+                                        radial.radialDirection(i, 1), radial.radialDirection(i, 2));
+        for (int j = 0; j < k; ++j) result.rawPhysicalBasis(i, j) = values[j];
+        d::Require(std::isfinite(parameters(i)) && radial.referenceRadius(i) + parameters(i) > 0.,
+                   "nonpositive or nonfinite nodal physical radius");
+    }
+
+    mfem::ParGridFunction displacement(fem.displacementFes.get()), scalar(fem.surfaceDeformationFes.get());
+    displacement.SetFromTrueDofs(baseDisplacement);
+    std::vector<d::Face> faces;
+    constexpr int attribute = Schema::template boundary_attribute<mean_field::utils::domain::StellarSurface>();
+    constexpr int quadrature = 14;
+    int sampleCount = 0;
+    // Tagged internal interfaces are included by GetNBE; a Boundary-only face
+    // restriction would silently omit the stellar/vacuum interface here.
+    for (int b = 0; b < fem.surfaceDeformationFes->GetNBE(); ++b) {
+        if (fem.surfaceDeformationFes->GetBdrAttribute(b) != attribute) continue;
+        d::Face face;
+        face.boundary = b;
+        const auto &element = *fem.surfaceDeformationFes->GetBE(b);
+        mfem::Array<int> vectorDofs;
+        d::Require(fem.surfaceDeformationFes->GetBdrElementVDofs(b, face.scalarDofs) == nullptr &&
+                   fem.displacementFes->GetBdrElementVDofs(b, vectorDofs) == nullptr,
+                   "nontrivial boundary DOF transformations are unsupported");
+        d::Require(face.scalarDofs.Size() == element.GetDof() && vectorDofs.Size() == 3 * element.GetDof(),
+                   "the scalar and vector trace elements disagree");
+        displacement.GetSubVector(vectorDofs, face.displacement);
+        face.coefficients.SetSize(element.GetDof(), k);
+        sampleCount += mfem::IntRules.Get(element.GetGeomType(), quadrature).GetNPoints();
+        faces.push_back(std::move(face));
+    }
+    d::Require(sampleCount > 0, "empty stellar surface quadrature");
+    mfem::Vector mode(n), ambient(fem.surfaceDeformationFes->GetTrueVSize()), local;
+    for (int j = 0; j < k; ++j) {
+        for (int i = 0; i < n; ++i) mode(i) = result.rawPhysicalBasis(i, j);
+        radial.surfaceDofMap().scatter(mode, ambient);
+        scalar.SetFromTrueDofs(ambient);
+        for (auto &face : faces) {
+            scalar.GetSubVector(face.scalarDofs, local);
+            for (int i = 0; i < local.Size(); ++i) face.coefficients(i, j) = local(i);
+        }
+    }
+    mfem::DenseMatrix samples(sampleCount, k), weighted(sampleCount, k);
+    mfem::Vector weights(sampleCount), shape;
+    mfem::DenseMatrix derivative, currentJacobian(3, 2);
+    int sample = 0;
+    long double surfaceArea = 0.;
+    for (const auto &face : faces) {
+        const auto &element = *fem.surfaceDeformationFes->GetBE(face.boundary);
+        auto &transform = *fem.surfaceDeformationFes->GetBdrElementTransformation(face.boundary);
+        const auto &rule = mfem::IntRules.Get(element.GetGeomType(), quadrature);
+        const int ndof = element.GetDof();
+        shape.SetSize(ndof);
+        derivative.SetSize(ndof, 2);
+        for (int q = 0; q < rule.GetNPoints(); ++q, ++sample) {
+            const auto &point = rule.IntPoint(q);
+            transform.SetIntPoint(&point);
+            element.CalcShape(point, shape);
+            element.CalcDShape(point, derivative);
+            currentJacobian = transform.Jacobian();
+            // Different polynomial orders in the mesh and displacement are
+            // kept distinct: this is the exact trace Jacobian of X + u_h.
+            for (int c = 0; c < 3; ++c)
+                for (int a = 0; a < 2; ++a)
+                    for (int i = 0; i < ndof; ++i)
+                        currentJacobian(c, a) += face.displacement(i + c * ndof) * derivative(i, a);
+            const double weight = point.weight * currentJacobian.Weight();
+            d::Require(std::isfinite(weight) && weight > 0., "invalid physical surface area weight");
+            weights(sample) = weight;
+            surfaceArea += weight;
+            for (int j = 0; j < k; ++j) {
+                long double value = 0.;
+                for (int i = 0; i < ndof; ++i) value += static_cast<long double>(shape(i)) * face.coefficients(i, j);
+                samples(sample, j) = static_cast<double>(value);
+                weighted(sample, j) = std::sqrt(weight) * samples(sample, j);
+            }
+        }
+    }
+    result.physicalSurfaceArea = static_cast<double>(surfaceArea);
+    result.rawGram.SetSize(k, k);
+    for (int i = 0; i < k; ++i)
+        for (int j = 0; j < k; ++j)
+            result.rawGram(i, j) = static_cast<double>(d::ColumnDot(weighted, i, weighted, j));
+
+    // Ordered, twice-reorthogonalized columns preserve nested harmonic spans.
+    // Rank loss rejects the requested space; it never silently changes K.
+    mfem::DenseMatrix orthogonal(weighted);
+    result.coefficientTransform.SetSize(k, k);
+    result.coefficientTransform = 0.;
+    std::vector<double> pivots(k), relativePivots(k);
+    for (int j = 0; j < k; ++j) {
+        result.coefficientTransform(j, j) = 1.;
+        const double originalNorm = std::sqrt(static_cast<double>(d::ColumnDot(weighted, j, weighted, j)));
+        for (int pass = 0; pass < 2; ++pass)
+            for (int i = 0; i < j; ++i) {
+                const double projection = static_cast<double>(d::ColumnDot(orthogonal, i, orthogonal, j));
+                for (int row = 0; row < sampleCount; ++row)
+                    orthogonal(row, j) -= projection * orthogonal(row, i);
+                for (int row = 0; row < k; ++row)
+                    result.coefficientTransform(row, j) -= projection * result.coefficientTransform(row, i);
+            }
+        pivots[j] = std::sqrt(static_cast<double>(d::ColumnDot(orthogonal, j, orthogonal, j)));
+        relativePivots[j] = pivots[j] / originalNorm;
+        d::Require(std::isfinite(relativePivots[j]) && originalNorm > 0. && relativePivots[j] > 1e-12,
+                   "the requested harmonic trace space is numerically rank deficient");
+        for (int row = 0; row < sampleCount; ++row) orthogonal(row, j) /= pivots[j];
+        for (int row = 0; row < k; ++row) result.coefficientTransform(row, j) /= pivots[j];
+    }
+    result.physicalBasis.SetSize(n, k);
+    mfem::Mult(result.rawPhysicalBasis, result.coefficientTransform, result.physicalBasis);
+    mfem::DenseMatrix reconstructed(sampleCount, k);
+    mfem::Mult(weighted, result.coefficientTransform, reconstructed);
+    double reconstructionError = 0.;
+    for (int i = 0; i < sampleCount; ++i)
+        for (int j = 0; j < k; ++j)
+            reconstructionError = std::max(reconstructionError, std::abs(reconstructed(i, j) - orthogonal(i, j)));
+    result.orthonormalGram.SetSize(k, k);
+    for (int i = 0; i < k; ++i)
+        for (int j = 0; j < k; ++j) {
+            const double value = static_cast<double>(d::ColumnDot(reconstructed, i, reconstructed, j));
+            d::Require(std::isfinite(value), "nonfinite reconstructed mass Gram matrix");
+            result.orthonormalGram(i, j) = value;
+            result.maximumGramOrthogonalityError = std::max(result.maximumGramOrthogonalityError,
+                                                             std::abs(value - (i == j ? 1. : 0.)));
+        }
+    d::Require(result.maximumGramOrthogonalityError <= 1e-10 && reconstructionError <= 1e-10,
+               "physical mass orthonormalization failed its reconstruction check");
+    mfem::Vector a(k), r(n), ba(n), btr(k);
+    for (int i = 0; i < k; ++i) a(i) = std::cos(.21 * (i + 1));
+    for (int i = 0; i < n; ++i) r(i) = std::sin(.37 * (i + 1));
+    result.physicalBasis.Mult(a, ba);
+    result.physicalBasis.MultTranspose(r, btr);
+    const double adjointError = std::abs((ba * r) - (a * btr)) /
+        std::max({1., ba.Norml2() * r.Norml2(), a.Norml2() * btr.Norml2()});
+    d::Require(std::isfinite(adjointError) && adjointError <= 1e-12, "basis transpose identity failed");
+
+    const bool acceptedUnchanged = d::Identical(accepted, state.acceptedNormalizedState) &&
+        d::Identical(acceptedResidual, state.acceptedNormalizedResidual) &&
+        d::Identical(acceptedPhysical, state.AcceptedPhysicalState());
+    const bool preparedUnchanged = generation == problem.GetPreparationGeneration() &&
+        state.normalizedOperator->IsPrepared() && d::Identical(preparedPhysical, state.normalizedOperator->GetPhysicalState()) &&
+        d::Identical(preparedResidual, state.normalizedOperator->GetPhysicalResidual()) &&
+        d::Identical(baseDisplacement, physical.GetGeneratedVolumeDisplacement()) &&
+        d::Identical(trial, state.trialNormalizedState) && d::Identical(trialResidual, state.trialNormalizedResidual) &&
+        d::Identical(candidatePhysical, state.candidatePhysicalState);
+    const bool factorsUnchanged = d::Identical(stateFactors, state.normalizedOperator->GetNormalization().StateFactors()) &&
+        d::Identical(residualFactors, state.normalizedOperator->GetNormalization().ResidualFactors());
+    d::Require(acceptedUnchanged && preparedUnchanged && factorsUnchanged, "a solver state buffer changed");
+    std::filesystem::create_directories(output);
+    d::Save(output, "basis_raw", result.rawPhysicalBasis);
+    d::Save(output, "basis_physical", result.physicalBasis);
+    d::Save(output, "gram_raw", result.rawGram);
+    d::Save(output, "gram_orthonormal", result.orthonormalGram);
+    d::Save(output, "coefficient_transform", result.coefficientTransform);
+    d::Save(output, "basis_samples_raw", samples);
+    auto weightOutput = d::Output(output, "basis_weights");
+    weightOutput << sampleCount << '\n';
+    for (int i = 0; i < sampleCount; ++i) weightOutput << weights(i) << '\n';
+    auto nodes = d::Output(output, "basis_nodes", ".csv");
+    nodes << "parameter,unit_radial_x,unit_radial_y,unit_radial_z,reference_radius,physical_radius,base_parameter\n";
+    for (int i = 0; i < n; ++i)
+        nodes << i << ',' << radial.radialDirection(i, 0) << ',' << radial.radialDirection(i, 1) << ','
+              << radial.radialDirection(i, 2) << ',' << radial.referenceRadius(i) << ','
+              << radial.referenceRadius(i) + parameters(i) << ',' << parameters(i) << '\n';
+    auto modes = d::Output(output, "basis_modes", ".csv");
+    modes << "column,degree,m,pivot_norm,relative_pivot\n";
+    int column = 0;
+    for (int l = 0; l <= maximumDegree; ++l)
+        for (int m = -l; m <= l; ++m, ++column)
+            modes << column << ',' << l << ',' << m << ',' << pivots[column] << ',' << relativePivots[column] << '\n';
+    auto metadata = d::Output(output, "basis_metadata", ".csv");
+    metadata << "metric,value\nmaximum_degree," << maximumDegree << "\nmodes," << k
+             << "\nsurface_parameters," << n << "\nquadrature," << quadrature
+             << "\nsurface_quadrature_points," << sampleCount << "\nphysical_surface_area," << result.physicalSurfaceArea
+             << "\nmaximum_gram_orthogonality_error," << result.maximumGramOrthogonalityError
+             << "\nmaximum_sample_reconstruction_error," << reconstructionError << "\nadjoint_relative_error," << adjointError
+             << "\naccepted_buffers_unchanged," << acceptedUnchanged << "\nprepared_state_unchanged," << preparedUnchanged
+             << "\nnormalization_factors_unchanged," << factorsUnchanged
+             << "\nmetric_scope,scalar_radial_amplitude_FE_trace_mass_on_current_physical_surface"
+             << "\nscalar_trace_order," << fem.surfaceDeformationFes->GetBE(faces.front().boundary)->GetOrder()
+             << "\nscalar_finite_element_collection," << fem.surfaceDeformationFes->FEColl()->Name()
+             << "\ndisplacement_finite_element_collection," << fem.displacementFes->FEColl()->Name()
+             << "\nstellar_surface_faces," << faces.size()
+             << "\nvector_displacement_metric_claimed,0\nbase_projection_performed,0\ntrial_preparations,0"
+             << "\nradial_trace_validation_actions,1\nrank_relative_pivot_threshold,1e-12"
+             << "\nconstruction_seconds," << std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count() << '\n';
+    return result;
+}
+} // namespace experiment::plateau_diagnostics

experiments/reduced_surface_diagnostics.hpp

@@ -0,0 +1,270 @@
+#pragma once
+
+// A frozen, affine surface restriction tests approximation-space compatibility
+// without changing the production problem or committing a filtered model.
+#include "plateau_nonlinear_step_diagnostics.hpp"
+#include "plateau_coupled_trial_physical.hpp"
+#include "shape_field_physical_diagnostics.hpp"
+#include "reduced_surface_basis.hpp"
+#include "reduced_surface_operators.hpp"
+#include <memory>
+
+namespace experiment::plateau_diagnostics {
+namespace reduced_surface_probe_detail {
+inline void SaveMatrix(const std::filesystem::path &path, const mfem::DenseMatrix &matrix) {
+    auto file = nonlinear_step_detail::Open(path);
+    file << matrix.Height() << ' ' << matrix.Width() << '\n';
+    for (int i = 0; i < matrix.Height(); ++i)
+        for (int j = 0; j < matrix.Width(); ++j) {
+            if (!std::isfinite(matrix(i, j))) throw std::runtime_error("Nonfinite reduced surface matrix.");
+            file << matrix(i, j) << '\n';
+        }
+}
+} // namespace reduced_surface_probe_detail
+
+template <class State>
+void ProbeReducedSurface(State &state, const mean_field::fem::FEM &fem,
+                         const mfem::Vector &normalizedBase, const std::filesystem::path &output,
+                         int maximumDegree, double n, double K, NonlinearStepProbeOptions options = {}) {
+    using namespace mean_field;
+    namespace d = nonlinear_step_detail;
+    namespace pd = reduced_surface_probe_detail;
+    const auto start = d::Clock::now();
+    options.linear.Validate();
+    int ranks = 0;
+    MPI_Comm_size(state.Problem().GetCommunicator(), &ranks);
+    const int size = state.Problem().StateSize();
+    if (ranks != 1 || size != state.Problem().EquationSize() || normalizedBase.Size() != size ||
+        !d::Finite(normalizedBase) || options.restartLength < 1)
+        throw std::invalid_argument("Reduced surface probes require one rank and finite square root data.");
+    const auto algebra = RunReducedSurfaceAlgebraSelfCheck(state.Problem().GetCommunicator());
+    const bool identity = maximumDegree == 0;
+    const mfem::Vector base(normalizedBase);
+    const auto &normalization = state.normalizedOperator->GetNormalization();
+    const mfem::Vector stateFactors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors());
+    const mfem::Vector savedAccepted(state.acceptedNormalizedState), savedResidual(state.acceptedNormalizedResidual),
+        savedPhysical(state.AcceptedPhysicalState()), savedCorrection(state.normalizedCorrection),
+        savedPhysicalCorrection(state.physicalCorrection), savedVolumeDirection(state.volumeDisplacementDirection);
+    const double savedMinimumJacobian = state.acceptedMinimumJacobianDeterminant;
+    const auto values = state.Problem().GetManifest().valueBlocks();
+    const auto rows = state.Problem().GetManifest().residualBlocks();
+    const auto shapeValue = shape_field_physical_detail::RequireBlock(values, "surface_deformation", size);
+    const auto shapeRow = shape_field_physical_detail::RequireBlock(rows, "surface_shape_balance", size);
+    if (shapeValue.size != shapeRow.size) throw std::invalid_argument("Shape dimensions disagree.");
+    std::filesystem::create_directories(output);
+    const auto file = [&](const std::string &name) { return output / ("reduced_surface_" + name); };
+    auto summary = d::Open(file("summary.csv"));
+    summary << "metric,value\n";
+    const auto metric = [&](const char *key, auto value) { summary << key << ',' << value << '\n'; };
+    auto layout = d::Open(file("layout.csv"));
+    layout << "kind,block,offset,size\n";
+    for (const auto &b : values) layout << "value," << b.stableId << ',' << b.offset << ',' << b.size << '\n';
+    for (const auto &b : rows) layout << "residual," << b.stableId << ',' << b.offset << ',' << b.size << '\n';
+    layout.flush();
+    mfem::Vector residual(size), direction(size), action(size), physicalBase, physicalDirection;
+    bool haveBase = false, finalized = false;
+    int preparations = 0;
+    double restoreError = 0.;
+    const auto prepare = [&](const mfem::Vector &candidate) {
+        state.trialNormalizedState = candidate;
+        ++preparations;
+        return state.PrepareTrial();
+    };
+    const auto restore = [&]() {
+        if (!prepare(base)) throw std::runtime_error("Reduced surface probe cannot restore base.");
+        if (!haveBase) { residual = state.trialNormalizedResidual; haveBase = true; }
+        if (!d::Finite(state.trialNormalizedResidual)) throw std::runtime_error("Nonfinite base residual.");
+        restoreError = std::max(restoreError, d::Difference(residual, state.trialNormalizedResidual));
+    };
+    const auto finish = [&]() {
+        finalized = true;
+        restore();
+        state.physicalInverse->Refresh();
+        state.normalizedCorrection = savedCorrection;
+        state.physicalCorrection = savedPhysicalCorrection;
+        state.volumeDisplacementDirection = savedVolumeDirection;
+        const bool acceptedSame = d::Identical(savedAccepted, state.acceptedNormalizedState) &&
+            d::Identical(savedResidual, state.acceptedNormalizedResidual) &&
+            d::Identical(savedPhysical, state.AcceptedPhysicalState()) &&
+            savedMinimumJacobian == state.acceptedMinimumJacobianDeterminant;
+        const bool factorsSame = d::Identical(stateFactors, normalization.StateFactors()) &&
+            d::Identical(residualFactors, normalization.ResidualFactors());
+        metric("accepted_buffers_bitwise_identical", acceptedSame);
+        metric("normalization_factors_bitwise_identical", factorsSame);
+        metric("direction_buffers_bitwise_identical", d::Identical(savedCorrection, state.normalizedCorrection) &&
+            d::Identical(savedPhysicalCorrection, state.physicalCorrection) &&
+            d::Identical(savedVolumeDirection, state.volumeDisplacementDirection));
+        metric("restored_residual_error", restoreError);
+        d::Save(file("restored_residual.txt"), state.trialNormalizedResidual);
+        metric("preparations", preparations);
+        metric("total_probe_seconds", d::Seconds(start));
+        summary.flush();
+        if (!acceptedSame || !factorsSame || restoreError != 0.)
+            throw std::runtime_error("Reduced surface probe changed frozen invariants.");
+    };
+    try {
+        restore();
+        state.physicalInverse->Refresh();
+        physicalBase = state.normalizedOperator->GetPhysicalState();
+        metric("schema_version", 1);
+        metric("maximum_degree", maximumDegree);
+        metric("identity_control", identity);
+        metric("full_size", size);
+        metric("shape_size", shapeValue.size);
+        metric("shape_value_offset", shapeValue.offset);
+        metric("shape_residual_offset", shapeRow.offset);
+        metric("linear_relative_tolerance", options.linear.relativeTolerance);
+        metric("linear_absolute_tolerance", options.linear.absoluteTolerance);
+        metric("linear_maximum_iterations", options.linear.maximumIterations);
+        metric("linear_restart", options.restartLength);
+        metric("accepted_commits", 0);
+        metric("algebra_self_check_passed", 1);
+        metric("algebra_self_check_count", algebra.checks);
+        metric("algebra_self_check_maximum_error", algebra.maximumError);
+
+        std::unique_ptr<ReducedSurfaceCoordinates> coordinates;
+        if (!identity) {
+            auto basis = MakeReducedSurfaceBasis(state, fem, maximumDegree, output);
+            coordinates = std::make_unique<ReducedSurfaceCoordinates>(basis.physicalBasis,
+                stateFactors, residualFactors, shapeValue.offset, shapeRow.offset, state.Problem().GetCommunicator());
+            pd::SaveMatrix(file("state_basis.txt"), coordinates->StateBasis());
+            pd::SaveMatrix(file("residual_basis.txt"), coordinates->ResidualBasis());
+            pd::SaveMatrix(file("state_transform.txt"), coordinates->StateTransform());
+            pd::SaveMatrix(file("residual_transform.txt"), coordinates->ResidualTransform());
+        }
+        const int reducedSize = identity ? size : coordinates->ReducedSize();
+        metric("mode_count", identity ? shapeValue.size : coordinates->ModeCount());
+        metric("reduced_size", reducedSize);
+        const auto restrictResidual = [&](const mfem::Vector &full, mfem::Vector &reduced) {
+            if (identity) reduced = full; else coordinates->RestrictResidual(full, reduced);
+        };
+        const auto excluded = [&](const mfem::Vector &full, mfem::Vector &result) {
+            mfem::Vector reduced, retained;
+            restrictResidual(full, reduced);
+            if (identity) retained = reduced; else coordinates->InjectResidual(reduced, retained);
+            result = full; result -= retained;
+        };
+        mfem::Vector reducedResidual, reducedDirection(reducedSize), reducedAction, reducedDefect;
+        restrictResidual(residual, reducedResidual);
+        reducedDirection = 0.;
+        mfem::Vector rhs(reducedResidual); rhs *= -1.;
+        const auto solve = [&](auto &operation, auto &inverse) {
+            auto backend = solver::linear::prepareLinearBackend(
+                solver::linear::FGMRES({.restartLength = options.restartLength, .printLevel = -1}),
+                operation, inverse, state.Problem().GetCommunicator());
+            const auto result = backend.Solve(rhs, reducedDirection, options.linear);
+            metric("linear_status", int(result.status));
+            metric("linear_iterations", result.iterations);
+            metric("linear_seconds", result.solveSeconds);
+            metric("linear_reported_true_norm", result.trueResidualNorm);
+        };
+        if (identity) solve(*state.normalizedOperator, *state.normalizedInverse);
+        else {
+            ReducedSurfaceOperator operation(*state.normalizedOperator, *coordinates);
+            ReducedSurfacePreconditioner inverse(*state.normalizedInverse, *coordinates);
+            solve(operation, inverse);
+        }
+        restore();
+        if (identity) direction = reducedDirection; else coordinates->InjectState(reducedDirection, direction);
+        if (!d::Finite(direction) || d::Norm(direction) == 0.) throw std::runtime_error("Invalid reduced direction.");
+        state.normalizedOperator->Mult(direction, action);
+        if (!d::Finite(action)) throw std::runtime_error("Nonfinite full Jacobian action.");
+        restrictResidual(action, reducedAction);
+        mfem::Vector defect(action); defect += residual;
+        reducedDefect = reducedAction; reducedDefect += reducedResidual;
+        mfem::Vector excludedResidual, excludedDefect;
+        excluded(residual, excludedResidual); excluded(defect, excludedDefect);
+        const double fullBaseNorm = d::Norm(residual), reducedBaseNorm = d::Norm(reducedResidual);
+        for (double norm : {fullBaseNorm, reducedBaseNorm, d::Norm(defect), d::Norm(reducedDefect),
+                            d::Norm(direction), d::Norm(excludedResidual), d::Norm(excludedDefect)})
+            if (!std::isfinite(norm)) throw std::runtime_error("Nonfinite reduced diagnostic norm.");
+        if (!(fullBaseNorm > 0.) || !(reducedBaseNorm > 0.))
+            throw std::runtime_error("This comparison requires nonzero full and retained residuals.");
+        metric("full_base_norm", fullBaseNorm);
+        metric("reduced_base_norm", reducedBaseNorm);
+        metric("excluded_base_norm", d::Norm(excludedResidual));
+        metric("full_relative_defect", d::Ratio(d::Norm(defect), fullBaseNorm));
+        metric("reduced_relative_defect", d::Ratio(d::Norm(reducedDefect), reducedBaseNorm));
+        metric("full_verified_gate", d::Norm(defect) <= options.linear.ConvergenceThreshold(fullBaseNorm));
+        metric("reduced_verified_gate", d::Norm(reducedDefect) <= options.linear.ConvergenceThreshold(reducedBaseNorm));
+        metric("full_merit_slope", d::Dot(residual, action, 0, size));
+        metric("reduced_merit_slope", d::Dot(reducedResidual, reducedAction, 0, reducedSize));
+        metric("direction_norm", d::Norm(direction));
+        state.normalizedOperator->DenormalizeState(direction, physicalDirection);
+        for (const auto &[name, vector] : std::initializer_list<std::pair<const char *, const mfem::Vector *>>{
+                 {"base", &base}, {"physical_base", &physicalBase}, {"residual", &residual},
+                 {"direction", &direction}, {"physical_direction", &physicalDirection}, {"action", &action},
+                 {"defect", &defect}, {"state_factors", &stateFactors}, {"residual_factors", &residualFactors},
+                 {"reduced_residual", &reducedResidual}, {"reduced_direction", &reducedDirection},
+                 {"reduced_action", &reducedAction}, {"reduced_defect", &reducedDefect},
+                 {"excluded_residual", &excludedResidual}, {"excluded_defect", &excludedDefect}})
+            d::Save(file(std::string(name) + ".txt"), *vector);
+        summary.flush();
+        std::cout << "Reduced surface degree " << maximumDegree << ": linear iterations complete; full relative defect "
+                  << d::Ratio(d::Norm(defect), fullBaseNorm) << ", reduced "
+                  << d::Ratio(d::Norm(reducedDefect), reducedBaseNorm) << std::endl;
+        MeasureShapeFieldPhysical(state, fem, physicalBase, physicalDirection, output);
+        auto physical = d::Open(file("physical.csv"));
+        physical << "sample,quadrature,metric,value\n";
+        const auto measurePhysical = [&](int sample) {
+            coupled_physical_detail::PreparedSnapshotView<State> view(state);
+            experiment::polytrope_validation::PhysicalState snapshot(view, fem);
+            for (int q : {14, 18})
+                for (const auto &[key, value] : MeasurePhysical(snapshot, n, K, q))
+                    physical << sample << ',' << q << ',' << key << ',' << value << '\n';
+            physical.flush();
+        };
+        measurePhysical(-1);
+        auto trials = d::Open(file("trials.csv"));
+        trials << "sample,alpha,prepared,reason,stage,full_norm,reduced_norm,excluded_shape_norm,"
+                  "full_merit_pass,reduced_merit_pass,full_remainder_norm,reduced_remainder_norm\n";
+        bool measuredTrial = false;
+        int sample = 0;
+        for (double alpha : {1., .25, .0625, 1./1024.}) {
+            // The affine base is retained verbatim: no projection of q0 and no
+            // enthalpy repair, SOC, or accepted-state update belongs in this test.
+            mfem::Vector candidate(base); candidate.Add(alpha, direction);
+            const auto prepared = prepare(candidate);
+            d::Save(file("candidate_a" + std::to_string(sample) + ".txt"), candidate);
+            if (prepared) {
+                const mfem::Vector trial(state.trialNormalizedResidual);
+                if (!d::Finite(trial)) throw std::runtime_error("Nonfinite reduced trial residual.");
+                mfem::Vector reducedTrial, omitted, remainder(trial), reducedRemainder;
+                restrictResidual(trial, reducedTrial); excluded(trial, omitted);
+                remainder -= residual; remainder.Add(-alpha, action);
+                restrictResidual(remainder, reducedRemainder);
+                const double fullNorm = d::Norm(trial), reducedNorm = d::Norm(reducedTrial);
+                trials << sample << ',' << alpha << ",1,-1,-1," << fullNorm << ',' << reducedNorm << ','
+                       << d::Norm(omitted) << ',' << (fullNorm <= (1.-1e-4*alpha)*fullBaseNorm) << ','
+                       << (reducedNorm <= (1.-1e-4*alpha)*reducedBaseNorm) << ','
+                       << d::Norm(remainder) << ',' << d::Norm(reducedRemainder) << '\n';
+                d::Save(file("residual_a" + std::to_string(sample) + ".txt"), trial);
+                d::Save(file("reduced_residual_a" + std::to_string(sample) + ".txt"), reducedTrial);
+                mfem::Vector physicalCandidate;
+                state.normalizedOperator->DenormalizeState(candidate, physicalCandidate);
+                d::Save(file("physical_candidate_a" + std::to_string(sample) + ".txt"), physicalCandidate);
+                if (!measuredTrial) { measurePhysical(sample); measuredTrial = true; }
+            } else {
+                trials << sample << ',' << alpha << ",0," << int(prepared.error().reason) << ','
+                       << int(prepared.error().stage) << ",nan,nan,nan,0,0,nan,nan\n";
+            }
+            trials.flush();
+            ++sample;
+        }
+        restore();
+        mfem::Vector restoredAction(size);
+        state.normalizedOperator->Mult(direction, restoredAction);
+        d::Save(file("restored_action.txt"), restoredAction);
+        metric("restored_action_error", d::Difference(action, restoredAction));
+        if (!d::Identical(action, restoredAction)) throw std::runtime_error("Restored Jacobian changed.");
+        finish();
+        metric("probe_exception", 0);
+        summary.flush();
+    } catch (...) {
+        const auto failure = std::current_exception();
+        metric("probe_exception", 1);
+        if (!finalized) finish();
+        std::rethrow_exception(failure);
+    }
+}
+} // namespace experiment::plateau_diagnostics

experiments/reduced_surface_operators.hpp

@@ -0,0 +1,527 @@
+#pragma once
+
+// Include after importing mean_field. These maps change the discrete shape
+// space; they are not an invertible preconditioner of the original equations.
+#include <algorithm>
+#include <cmath>
+#include <concepts>
+#include <cstdint>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <type_traits>
+#include <vector>
+
+#include <mfem.hpp>
+#include <mpi.h>
+
+namespace experiment::plateau_diagnostics {
+    struct ReducedSurfaceCoordinateOptions final {
+        double relativeRankTolerance{1.e-12};
+        double orthogonalityTolerance{5.e-12};
+    };
+
+    struct ReducedSurfaceCoordinateReport final {
+        double stateOrthogonalityError{0.0};
+        double residualOrthogonalityError{0.0};
+        double stateReconstructionError{0.0};
+        double residualReconstructionError{0.0};
+        double stateMinimumIndependentFraction{1.0};
+        double residualMinimumIndependentFraction{1.0};
+    };
+
+    namespace reduced_surface_detail {
+        inline void RequireFinite(const mfem::Vector &value, const char *description) {
+            for (int i = 0; i < value.Size(); ++i)
+                if (!std::isfinite(value[i])) throw std::invalid_argument(description);
+        }
+
+        inline double StableNorm(const mfem::Vector &value) {
+            double norm = 0.0;
+            for (int i = 0; i < value.Size(); ++i) norm = std::hypot(norm, value[i]);
+            return norm;
+        }
+
+        inline double DotColumn(const mfem::DenseMatrix &matrix, int column, const mfem::Vector &value) {
+            long double result = 0.0L;
+            for (int i = 0; i < matrix.Height(); ++i)
+                result += static_cast<long double>(matrix(i, column)) * value[i];
+            return static_cast<double>(result);
+        }
+
+        struct QRReport final {
+            double orthogonalityError{0.0};
+            double reconstructionError{0.0};
+            double minimumIndependentFraction{1.0};
+        };
+
+        // Reorthogonalization preserves the advertised coordinate identities
+        // without a new dense linear-algebra dependency. Rejecting a dependent
+        // mode keeps the scientific degree/mode ordering explicit in artifacts.
+        inline QRReport ThinQR(
+            const mfem::DenseMatrix &input,
+            mfem::DenseMatrix &basis,
+            mfem::DenseMatrix &transform,
+            ReducedSurfaceCoordinateOptions options
+        ) {
+            const int rows = input.Height(), columns = input.Width();
+            basis.SetSize(rows, columns);
+            transform.SetSize(columns, columns);
+            basis = 0.0;
+            transform = 0.0;
+            QRReport report;
+            mfem::Vector column(rows);
+            for (int j = 0; j < columns; ++j) {
+                for (int i = 0; i < rows; ++i) column[i] = input(i, j);
+                RequireFinite(column, "The reduced surface basis contains nonfinite entries.");
+                const double originalNorm = StableNorm(column);
+                if (!std::isfinite(originalNorm) || originalNorm <= 0.0)
+                    throw std::invalid_argument("The reduced surface basis has a zero or overflowing column.");
+                for (int pass = 0; pass < 2; ++pass) {
+                    for (int previous = 0; previous < j; ++previous) {
+                        const double projection = DotColumn(basis, previous, column);
+                        if (!std::isfinite(projection))
+                            throw std::invalid_argument("The reduced surface QR projection is nonfinite.");
+                        transform(previous, j) += projection;
+                        for (int i = 0; i < rows; ++i) column[i] -= projection * basis(i, previous);
+                    }
+                }
+                const double independentNorm = StableNorm(column);
+                const double fraction = independentNorm / originalNorm;
+                if (!std::isfinite(fraction) || fraction <= options.relativeRankTolerance)
+                    throw std::invalid_argument("The reduced surface basis is numerically rank deficient.");
+                report.minimumIndependentFraction = std::min(report.minimumIndependentFraction, fraction);
+                transform(j, j) = independentNorm;
+                for (int i = 0; i < rows; ++i) basis(i, j) = column[i] / independentNorm;
+            }
+            double inputNorm = 0.0, errorNorm = 0.0;
+            for (int j = 0; j < columns; ++j) {
+                for (int l = 0; l < columns; ++l) {
+                    long double product = 0.0L;
+                    for (int i = 0; i < rows; ++i)
+                        product += static_cast<long double>(basis(i, j)) * basis(i, l);
+                    report.orthogonalityError = std::max(
+                        report.orthogonalityError, std::abs(static_cast<double>(product) - (j == l ? 1.0 : 0.0))
+                    );
+                }
+                for (int i = 0; i < rows; ++i) {
+                    long double reconstructed = 0.0L;
+                    for (int l = 0; l <= j; ++l)
+                        reconstructed += static_cast<long double>(basis(i, l)) * transform(l, j);
+                    inputNorm = std::hypot(inputNorm, input(i, j));
+                    errorNorm = std::hypot(errorNorm, static_cast<double>(reconstructed) - input(i, j));
+                }
+            }
+            report.reconstructionError = errorNorm / inputNorm;
+            if (!std::isfinite(report.reconstructionError) ||
+                !std::isfinite(report.orthogonalityError) ||
+                report.orthogonalityError > options.orthogonalityTolerance ||
+                report.reconstructionError > options.orthogonalityTolerance)
+                throw std::invalid_argument("The reduced surface QR failed its coordinate consistency check.");
+            return report;
+        }
+    } // namespace reduced_surface_detail
+
+    class ReducedSurfaceCoordinates final {
+    public:
+        ReducedSurfaceCoordinates(
+            const mfem::DenseMatrix &physicalBasis,
+            const mfem::Vector &stateFactors,
+            const mfem::Vector &residualFactors,
+            int shapeValueOffset,
+            int shapeResidualOffset,
+            MPI_Comm communicator,
+            ReducedSurfaceCoordinateOptions options = {}
+        )
+            : m_physicalBasis(physicalBasis), m_stateFactors(stateFactors), m_residualFactors(residualFactors),
+              m_shapeValueOffset(shapeValueOffset), m_shapeResidualOffset(shapeResidualOffset) {
+            int ranks = 0;
+            if (MPI_Comm_size(communicator, &ranks) != MPI_SUCCESS || ranks != 1)
+                throw std::invalid_argument("The reduced surface experiment currently requires one MPI rank.");
+            if (!std::isfinite(options.relativeRankTolerance) || options.relativeRankTolerance <= 0.0 ||
+                options.relativeRankTolerance >= 1.0 || !std::isfinite(options.orthogonalityTolerance) ||
+                options.orthogonalityTolerance <= 0.0 || options.orthogonalityTolerance >= 1.0)
+                throw std::invalid_argument("Invalid reduced surface coordinate tolerances.");
+            const int size = stateFactors.Size(), shape = physicalBasis.Height(), modes = physicalBasis.Width();
+            if (size != residualFactors.Size() || shape <= 0 || modes <= 0 || modes > shape || shape > size ||
+                shapeValueOffset < 0 || shapeValueOffset > size - shape ||
+                shapeResidualOffset < 0 || shapeResidualOffset > size - shape)
+                throw std::invalid_argument("Invalid reduced surface dimensions or distinct block offsets.");
+            for (int i = 0; i < size; ++i)
+                if (!std::isfinite(stateFactors[i]) || stateFactors[i] <= 0.0 ||
+                    !std::isfinite(residualFactors[i]) || residualFactors[i] <= 0.0)
+                    throw std::invalid_argument("Reduced surface normalization factors must be finite and positive.");
+
+            mfem::DenseMatrix stateImage(shape, modes), residualDualImage(shape, modes);
+            for (int j = 0; j < modes; ++j) {
+                for (int i = 0; i < shape; ++i) {
+                    stateImage(i, j) = stateFactors[shapeValueOffset + i] * physicalBasis(i, j);
+                    residualDualImage(i, j) = physicalBasis(i, j) / residualFactors[shapeResidualOffset + i];
+                }
+            }
+            // State and residual scales act on different layouts and are not
+            // interchangeable. Q_r^T Fhat = T_r^{-T} B^T Fphysical preserves
+            // precisely the physical Galerkin equations with an invertible
+            // reduced residual scaling, even when Q_x differs from Q_r.
+            const auto state = reduced_surface_detail::ThinQR(stateImage, m_stateBasis, m_stateTransform, options);
+            const auto residual = reduced_surface_detail::ThinQR(
+                residualDualImage, m_residualBasis, m_residualTransform, options
+            );
+            m_report = {
+                state.orthogonalityError, residual.orthogonalityError,
+                state.reconstructionError, residual.reconstructionError,
+                state.minimumIndependentFraction, residual.minimumIndependentFraction
+            };
+        }
+
+        [[nodiscard]] int FullSize() const noexcept { return m_stateFactors.Size(); }
+        [[nodiscard]] int ReducedSize() const noexcept { return FullSize() - ShapeSize() + ModeCount(); }
+        [[nodiscard]] int ShapeSize() const noexcept { return m_physicalBasis.Height(); }
+        [[nodiscard]] int ModeCount() const noexcept { return m_physicalBasis.Width(); }
+        [[nodiscard]] int ShapeValueOffset() const noexcept { return m_shapeValueOffset; }
+        [[nodiscard]] int ShapeResidualOffset() const noexcept { return m_shapeResidualOffset; }
+        [[nodiscard]] const mfem::DenseMatrix &PhysicalBasis() const noexcept { return m_physicalBasis; }
+        [[nodiscard]] const mfem::DenseMatrix &StateBasis() const noexcept { return m_stateBasis; }
+        [[nodiscard]] const mfem::DenseMatrix &ResidualBasis() const noexcept { return m_residualBasis; }
+        [[nodiscard]] const mfem::DenseMatrix &StateTransform() const noexcept { return m_stateTransform; }
+        [[nodiscard]] const mfem::DenseMatrix &ResidualTransform() const noexcept { return m_residualTransform; }
+        [[nodiscard]] const ReducedSurfaceCoordinateReport &Report() const noexcept { return m_report; }
+
+        [[nodiscard]] bool MatchesNormalization(
+            const mfem::Vector &stateFactors, const mfem::Vector &residualFactors
+        ) const noexcept {
+            if (stateFactors.Size() != FullSize() || residualFactors.Size() != FullSize()) return false;
+            for (int i = 0; i < FullSize(); ++i)
+                if (stateFactors[i] != m_stateFactors[i] || residualFactors[i] != m_residualFactors[i]) return false;
+            return true;
+        }
+
+        void InjectState(const mfem::Vector &reduced, mfem::Vector &full) const {
+            Inject(reduced, full, m_shapeValueOffset, m_stateBasis);
+        }
+        void RestrictState(const mfem::Vector &full, mfem::Vector &reduced) const {
+            Restrict(full, reduced, m_shapeValueOffset, m_stateBasis);
+        }
+        void RestrictResidual(const mfem::Vector &full, mfem::Vector &reduced) const {
+            Restrict(full, reduced, m_shapeResidualOffset, m_residualBasis);
+        }
+        void InjectResidual(const mfem::Vector &reduced, mfem::Vector &full) const {
+            Inject(reduced, full, m_shapeResidualOffset, m_residualBasis);
+        }
+
+    private:
+        void Inject(
+            const mfem::Vector &reduced, mfem::Vector &full, int offset, const mfem::DenseMatrix &basis
+        ) const {
+            if (reduced.Size() != ReducedSize() || (&reduced == &full))
+                throw std::invalid_argument("Reduced surface injection requires distinct correctly sized vectors.");
+            reduced_surface_detail::RequireFinite(reduced, "Cannot inject a nonfinite reduced surface vector.");
+            full.SetSize(FullSize());
+            for (int i = 0; i < offset; ++i) full[i] = reduced[i];
+            for (int i = 0; i < ShapeSize(); ++i) {
+                long double value = 0.0L;
+                for (int j = 0; j < ModeCount(); ++j)
+                    value += static_cast<long double>(basis(i, j)) * reduced[offset + j];
+                full[offset + i] = static_cast<double>(value);
+            }
+            for (int i = offset + ShapeSize(); i < FullSize(); ++i)
+                full[i] = reduced[i - ShapeSize() + ModeCount()];
+            reduced_surface_detail::RequireFinite(full, "Reduced surface injection overflowed.");
+        }
+
+        void Restrict(
+            const mfem::Vector &full, mfem::Vector &reduced, int offset, const mfem::DenseMatrix &basis
+        ) const {
+            if (full.Size() != FullSize() || (&full == &reduced))
+                throw std::invalid_argument("Reduced surface restriction requires distinct correctly sized vectors.");
+            reduced_surface_detail::RequireFinite(full, "Cannot restrict a nonfinite full surface vector.");
+            reduced.SetSize(ReducedSize());
+            for (int i = 0; i < offset; ++i) reduced[i] = full[i];
+            for (int j = 0; j < ModeCount(); ++j) {
+                long double value = 0.0L;
+                for (int i = 0; i < ShapeSize(); ++i)
+                    value += static_cast<long double>(basis(i, j)) * full[offset + i];
+                reduced[offset + j] = static_cast<double>(value);
+            }
+            for (int i = offset + ShapeSize(); i < FullSize(); ++i)
+                reduced[i - ShapeSize() + ModeCount()] = full[i];
+            reduced_surface_detail::RequireFinite(reduced, "Reduced surface restriction overflowed.");
+        }
+
+        mfem::DenseMatrix m_physicalBasis, m_stateBasis, m_residualBasis, m_stateTransform, m_residualTransform;
+        mfem::Vector m_stateFactors, m_residualFactors;
+        int m_shapeValueOffset, m_shapeResidualOffset;
+        ReducedSurfaceCoordinateReport m_report;
+    };
+
+    template <std::derived_from<mfem::Operator> FullOperator>
+    class ReducedSurfaceOperator final : public mfem::Operator {
+    public:
+        ReducedSurfaceOperator(const FullOperator &operation, const ReducedSurfaceCoordinates &coordinates)
+            : mfem::Operator(coordinates.ReducedSize()), m_operation(&operation), m_coordinates(&coordinates),
+              m_fullDirection(coordinates.FullSize()), m_fullAction(coordinates.FullSize()) {
+            Refresh();
+        }
+
+        // Refresh does not prepare the physical problem. It adopts an already
+        // prepared generation while requiring the frozen coordinate scales.
+        void Refresh() {
+            m_current = false;
+            if (m_operation->Height() != m_coordinates->FullSize() ||
+                m_operation->Width() != m_coordinates->FullSize())
+                throw std::invalid_argument("Reduced surface operator needs the matching square full Jacobian.");
+            if constexpr (requires { m_operation->IsPrepared(); })
+                if (!m_operation->IsPrepared()) throw std::invalid_argument("The full Jacobian is not prepared.");
+            if constexpr (requires { m_operation->GetNormalization(); }) {
+                const auto &normalization = m_operation->GetNormalization();
+                if (!m_coordinates->MatchesNormalization(normalization.StateFactors(), normalization.ResidualFactors()))
+                    throw std::invalid_argument("The frozen reduced surface normalization has changed.");
+            }
+            if constexpr (requires { m_operation->GetProblem().GetPreparationGeneration(); })
+                m_generation = m_operation->GetProblem().GetPreparationGeneration();
+            m_current = true;
+        }
+
+        [[nodiscard]] bool IsCurrent() const {
+            if (!m_current || m_operation->Height() != m_coordinates->FullSize() ||
+                m_operation->Width() != m_coordinates->FullSize()) return false;
+            if constexpr (requires { m_operation->IsPrepared(); })
+                if (!m_operation->IsPrepared()) return false;
+            if constexpr (requires { m_operation->GetProblem().GetPreparationGeneration(); })
+                if (m_generation != m_operation->GetProblem().GetPreparationGeneration()) return false;
+            if constexpr (requires { m_operation->GetNormalization(); }) {
+                const auto &normalization = m_operation->GetNormalization();
+                if (!m_coordinates->MatchesNormalization(normalization.StateFactors(), normalization.ResidualFactors()))
+                    return false;
+            }
+            return true;
+        }
+
+        void Mult(const mfem::Vector &direction, mfem::Vector &action) const override {
+            if (!IsCurrent()) throw std::invalid_argument("The reduced surface Jacobian is stale.");
+            m_coordinates->InjectState(direction, m_fullDirection);
+            m_operation->Mult(m_fullDirection, m_fullAction);
+            m_coordinates->RestrictResidual(m_fullAction, action);
+        }
+
+        [[nodiscard]] const FullOperator &GetFullOperator() const noexcept { return *m_operation; }
+
+    private:
+        const FullOperator *m_operation;
+        const ReducedSurfaceCoordinates *m_coordinates;
+        mutable mfem::Vector m_fullDirection, m_fullAction;
+        std::uint64_t m_generation{0};
+        bool m_current{false};
+    };
+
+    template <class FullInverse>
+        requires std::derived_from<FullInverse, mfem::Solver> &&
+                 mean_field::solver::LinearPreconditionerApplicationContractAvailable<FullInverse>
+    class ReducedSurfacePreconditioner final : public mfem::Solver {
+    public:
+        // Projecting an inverse is an approximation, not an inverse identity.
+        // In particular a flexible underlying application remains flexible.
+        static constexpr auto applicationContract =
+            mean_field::solver::linearPreconditionerApplicationContract<FullInverse>;
+
+        ReducedSurfacePreconditioner(FullInverse &inverse, const ReducedSurfaceCoordinates &coordinates)
+            : mfem::Solver(coordinates.ReducedSize()), m_inverse(&inverse), m_coordinates(&coordinates),
+              m_fullResidual(coordinates.FullSize()), m_fullCorrection(coordinates.FullSize()) {
+            if (!IsCurrent()) throw std::invalid_argument("The full reduced-surface inverse is not current.");
+            iterative_mode = false;
+        }
+
+        void SetOperator(const mfem::Operator &operation) override {
+            if (operation.Height() != Height() || operation.Width() != Width())
+                throw std::invalid_argument("The projected inverse received an incompatible reduced operator.");
+            // The physical inverse is already bound to the authoritative full
+            // operator. Passing this reduced operator into it would be invalid.
+            if (!IsCurrent()) throw std::invalid_argument("The projected inverse is stale.");
+        }
+
+        [[nodiscard]] bool IsCurrent() const {
+            if (m_inverse->Height() != m_coordinates->FullSize() ||
+                m_inverse->Width() != m_coordinates->FullSize() || m_inverse->iterative_mode) return false;
+            if constexpr (requires { m_inverse->IsCurrent(); })
+                if (!m_inverse->IsCurrent()) return false;
+            return true;
+        }
+
+        void Mult(const mfem::Vector &residual, mfem::Vector &correction) const override {
+            if (!IsCurrent()) throw std::invalid_argument("The projected inverse is stale.");
+            m_coordinates->InjectResidual(residual, m_fullResidual);
+            m_fullCorrection = 0.0;
+            m_inverse->Mult(m_fullResidual, m_fullCorrection);
+            m_coordinates->RestrictState(m_fullCorrection, correction);
+        }
+
+    private:
+        FullInverse *m_inverse;
+        const ReducedSurfaceCoordinates *m_coordinates;
+        mutable mfem::Vector m_fullResidual, m_fullCorrection;
+    };
+
+    struct ReducedSurfaceAlgebraSelfCheckReport final {
+        int checks{0};
+        double maximumError{0.0};
+        bool rejectedDependentBasis{false};
+        bool rejectedNonfiniteBasis{false};
+        bool rejectedStaleOperator{false};
+    };
+
+    namespace reduced_surface_detail {
+        class SelfCheckOperator final : public mfem::Operator {
+        public:
+            explicit SelfCheckOperator(int size) : mfem::Operator(size), matrix(size) {
+                matrix = 0.0;
+                for (int i = 0; i < size; ++i) {
+                    matrix(i, i) = 2.0 + i;
+                    for (int j = 0; j < size; ++j)
+                        if (i != j) matrix(i, j) = 0.03 * (1 + i - 2 * j);
+                }
+            }
+            void Mult(const mfem::Vector &x, mfem::Vector &y) const override {
+                ++applications;
+                y.SetSize(Height());
+                matrix.Mult(x, y);
+            }
+            [[nodiscard]] bool IsPrepared() const noexcept { return true; }
+            [[nodiscard]] const SelfCheckOperator &GetProblem() const noexcept { return *this; }
+            [[nodiscard]] std::uint64_t GetPreparationGeneration() const noexcept { return generation; }
+            mfem::DenseMatrix matrix;
+            std::uint64_t generation{1};
+            mutable int applications{0};
+        };
+
+        class SelfCheckInverse final : public mfem::Solver {
+        public:
+            static constexpr auto applicationContract =
+                mean_field::preconditioning::ApplicationContract::stationary_linear;
+            explicit SelfCheckInverse(int size) : mfem::Solver(size) { iterative_mode = false; }
+            void SetOperator(const mfem::Operator &) override {}
+            [[nodiscard]] bool IsCurrent() const noexcept { return true; }
+            void Mult(const mfem::Vector &x, mfem::Vector &y) const override {
+                ++applications;
+                y.SetSize(Height());
+                for (int i = 0; i < Height(); ++i) y[i] = x[i] / (2.0 + i);
+            }
+            mutable int applications{0};
+        };
+    } // namespace reduced_surface_detail
+
+    // Tiny deterministic algebra only: this exercises different value/row
+    // offsets and non-proportional scales without constructing any FEM state.
+    [[nodiscard]] inline ReducedSurfaceAlgebraSelfCheckReport RunReducedSurfaceAlgebraSelfCheck(MPI_Comm comm) {
+        ReducedSurfaceAlgebraSelfCheckReport report;
+        const auto check = [&report](double error) {
+            ++report.checks;
+            if (!std::isfinite(error) || error > 2.e-12)
+                throw std::runtime_error("Reduced surface algebra self-check failed.");
+            report.maximumError = std::max(report.maximumError, std::abs(error));
+        };
+        const auto vectorError = [&check](const mfem::Vector &a, const mfem::Vector &b) {
+            if (a.Size() != b.Size()) throw std::runtime_error("Reduced surface self-check size mismatch.");
+            double error = 0.0;
+            for (int i = 0; i < a.Size(); ++i) error = std::max(error, std::abs(a[i] - b[i]));
+            check(error);
+        };
+        mfem::Vector stateFactors(7), residualFactors(7), full(7), reduced(6);
+        for (int i = 0; i < 7; ++i) {
+            stateFactors[i] = 0.7 + 0.3 * i;
+            residualFactors[i] = 1.9 - 0.13 * i;
+            full[i] = 0.2 - 0.17 * i;
+        }
+        for (int i = 0; i < 6; ++i) reduced[i] = 0.11 * (i + 1);
+        mfem::DenseMatrix physical(3, 2);
+        physical(0, 0) = 1.0; physical(1, 0) = 0.5; physical(2, 0) = 0.0;
+        physical(0, 1) = 0.2; physical(1, 1) = -0.3; physical(2, 1) = 1.0;
+        ReducedSurfaceCoordinates coordinates(physical, stateFactors, residualFactors, 1, 4, comm);
+        check(coordinates.Report().stateOrthogonalityError);
+        check(coordinates.Report().residualOrthogonalityError);
+        mfem::Vector injected, restricted, roundtrip, other;
+        coordinates.InjectState(reduced, injected);
+        coordinates.RestrictState(injected, roundtrip);
+        vectorError(reduced, roundtrip);
+        coordinates.RestrictState(full, restricted);
+        check(std::abs(injected * full - reduced * restricted));
+        coordinates.InjectResidual(reduced, injected);
+        coordinates.RestrictResidual(injected, roundtrip);
+        vectorError(reduced, roundtrip);
+        coordinates.RestrictResidual(full, restricted);
+        check(std::abs(injected * full - reduced * restricted));
+
+        // Verify the physical dual identity independently of InjectResidual:
+        // T_r^T rhat_shape = B^T (D_r^{-1} Fhat)_shape.
+        for (int j = 0; j < 2; ++j) {
+            double lhs = 0.0, rhs = 0.0;
+            for (int l = 0; l < 2; ++l) lhs += coordinates.ResidualTransform()(l, j) * restricted[4 + l];
+            for (int i = 0; i < 3; ++i) rhs += physical(i, j) * full[4 + i] / residualFactors[4 + i];
+            check(std::abs(lhs - rhs));
+        }
+        // The two QR bases really differ in this fixture; transpose reuse
+        // would therefore fail the physical dual identity above.
+        double mapDifference = 0.0;
+        for (int j = 0; j < 2; ++j)
+            for (int i = 0; i < 3; ++i)
+                mapDifference += std::abs(coordinates.StateBasis()(i, j) - coordinates.ResidualBasis()(i, j));
+        if (mapDifference < 1.e-3) throw std::runtime_error("Reduced surface dual-map fixture is ineffective.");
+        ++report.checks;
+
+        reduced_surface_detail::SelfCheckOperator fullOperator(7);
+        ReducedSurfaceOperator operation(fullOperator, coordinates);
+        operation.Mult(reduced, roundtrip);
+        coordinates.InjectState(reduced, injected);
+        fullOperator.Mult(injected, other);
+        coordinates.RestrictResidual(other, restricted);
+        vectorError(roundtrip, restricted);
+        reduced_surface_detail::SelfCheckInverse fullInverse(7);
+        ReducedSurfacePreconditioner inverse(fullInverse, coordinates);
+        inverse.SetOperator(operation);
+        const int oldApplications = fullOperator.applications;
+        inverse.Mult(reduced, roundtrip);
+        coordinates.InjectResidual(reduced, injected);
+        fullInverse.Mult(injected, other);
+        coordinates.RestrictState(other, restricted);
+        vectorError(roundtrip, restricted);
+        check(fullOperator.applications - oldApplications);
+
+        ++fullOperator.generation;
+        try { operation.Mult(reduced, roundtrip); }
+        catch (const std::invalid_argument &) { report.rejectedStaleOperator = true; }
+        if (!report.rejectedStaleOperator || operation.IsCurrent())
+            throw std::runtime_error("Reduced surface stale-generation guard failed.");
+        ++report.checks;
+        operation.Refresh();
+        operation.Mult(reduced, roundtrip);
+
+        mfem::DenseMatrix identity(3);
+        identity = 0.0;
+        for (int i = 0; i < 3; ++i) identity(i, i) = 1.0;
+        ReducedSurfaceCoordinates fullCoordinates(identity, stateFactors, residualFactors, 1, 4, comm);
+        fullCoordinates.InjectState(full, injected);
+        vectorError(full, injected);
+        fullCoordinates.RestrictResidual(full, restricted);
+        vectorError(full, restricted);
+        ReducedSurfaceOperator fullRecovery(fullOperator, fullCoordinates);
+        fullRecovery.Mult(full, roundtrip);
+        fullOperator.Mult(full, other);
+        vectorError(roundtrip, other);
+        ReducedSurfacePreconditioner inverseRecovery(fullInverse, fullCoordinates);
+        inverseRecovery.Mult(full, roundtrip);
+        fullInverse.Mult(full, other);
+        vectorError(roundtrip, other);
+
+        mfem::DenseMatrix invalid(physical);
+        for (int i = 0; i < 3; ++i) invalid(i, 1) = invalid(i, 0);
+        try { ReducedSurfaceCoordinates rejected(invalid, stateFactors, residualFactors, 1, 4, comm); }
+        catch (const std::invalid_argument &) { report.rejectedDependentBasis = true; }
+        invalid = physical;
+        invalid(0, 0) = std::numeric_limits<double>::quiet_NaN();
+        try { ReducedSurfaceCoordinates rejected(invalid, stateFactors, residualFactors, 1, 4, comm); }
+        catch (const std::invalid_argument &) { report.rejectedNonfiniteBasis = true; }
+        if (!report.rejectedDependentBasis || !report.rejectedNonfiniteBasis)
+            throw std::runtime_error("Reduced surface basis validity guards failed.");
+        report.checks += 2;
+        return report;
+    }
+} // namespace experiment::plateau_diagnostics