Preserve completed weak-response study while retaining Q3 baseline

experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md

@@ -9,7 +9,7 @@ that the Q3 displacement/Q4 reference geometry mismatch contributes to weak
 response. Cancellation by itself is not evidence of a gauge or bad element orders.
 
 The baseline branch is `research/shape-field-weak-response-20260918`, based on
-`db854c0f7cdf5137e377273499f0ac5b6053b752`. The isolated order variant will be
+`db854c0f7cdf5137e377273499f0ac5b6053b752`. The isolated order variant is
 `research/shape-field-displacement-q4-20260918`. Preserve original checkout,
 immutable executables, source revisions, mesh/input hashes and serial run records.
 
@@ -62,7 +62,8 @@ comparison, not evidence of high-n/high-rotation convergence.
 
 ## Results
 
-Baseline observations are recorded below. No production default has been changed on this baseline branch.
+The baseline branch retains Q3 displacement. The separate Q4 branch is an
+isolated numerical experiment, not a recommendation to change the default.
 
 ### Reproducible baseline
 
@@ -91,7 +92,9 @@ symmetries, rank, and file binding on analytic fixtures.
 
 The largest sampled singular response remains approximately 22.1–22.6. The SVD
 uses the six nonzero disjoint value-block pieces of each p, normalized by their
-coefficient norms. It is not a full-J spectrum or rank test. In particular,
+coefficient norms. The direction-dependent subspace changes between observations;
+the comparison combines changes in the operator and the probe subspace. It is
+not a full-J spectrum or rank test. In particular,
 comparing the 12-step late state with six-step endpoints does not demonstrate a
 causal benefit of the coupled correction: the matched six-step control has a
 slightly stronger minimum response than the coupled endpoint.
@@ -105,16 +108,20 @@ Eulerian squared integrals is 1.69e-5, an integration sensitivity observation.
 The radius observations expose a validation gap: the old late state has nodal
 radii 0.997674–1.010093, while its q18 P(h)-based virial error is 4.06909e-6.
 The fresh six-step control has q18 virial error 4.22590e-6 and reproduces the old
-six-step residual trajectory. Neither is converged or physically certified;
+six-step residual trajectory. Its sampled hydrostatic-gradient RMS increases
+from 0.00152810 to 0.0106811 despite the lower residual and virial error.
+Neither is converged or physically certified;
 negative density samples remain. Virial success can coexist with appreciable
 nonrotating surface asphericity at these unfinished states.
 
 `q3_surface_modes_v1.json` contains unweighted fits, not surface-integral power.
 The 12-step late direction is 99.9995107% cubic-invariant in sampled squared
 norm, but degree <=8 spherical harmonics leave 0.893679 of its nodal norm.
-Degree <=2 leaves 0.981646. Thus this is predominantly finer cubic-symmetric
+Degree <=2 leaves 0.981646. Thus most of the nodal direction is unresolved by
-structure, not simple expansion, translation, or quadrupolar deformation.
+degree <=8 fits and is nearly cubic-invariant. This is consistent with
-This is consistent with mesh/discretization coupling, without proving its cause.
+mesh/discretization coupling, without identifying its source: cubic symmetry
+also exists in the bulk and gravity discretizations. These are not surface-area
+spectra or evidence that the surface discretization alone causes the problem.
 
 The affine control on the actual stellar mesh gives Q3 relative value error
 9.01277e-5 and gradient error 1.23472e-3. This warrants the isolated Q4 comparison,
@@ -148,13 +155,198 @@ This concerns that proposed replacement, not a proven deficiency of the current
 mechanical Jacobian. Any boundary formulation requires a compatible trace/test
 space and explicit rank/stability analysis.
 
-An invertible transport coordinate/preconditioner change preserves roots and may
+An invertible linear/tangent coordinate or preconditioner change preserves the
-improve linear algebra, but cannot itself change the exact Newton direction or
+exact mapped physical Newton correction at a fixed state. It can improve
-remove a true weak physical response. Shape representation, step control, and
+coordinate-dependent singular values and inexact linear solves. A nonlinear
-operator formulation remain distinct hypotheses.
+coordinate map or retraction can also change finite trial curvature; none of
+these observations alone removes weak physical response. Shape representation,
+step control, and operator formulation remain distinct hypotheses.
 
 ### Q4 intervention
 
-Pending on its isolated branch. Its rationale is to remove measured ambient
+Applied 2026-09-18 at approximately 17:45 UTC on
-affine representation error and test the effect on the observed weak direction;
+`research/shape-field-displacement-q4-20260918`, branched from baseline `211be17`.
-the additional surface modes and quadrature changes are explicit confounders.
+Changed `libmeanfield/interface/field/field_registry.cppm`,
+`Displacement::vectorOrder`, from `3 + uniformPolynomialOrderIncrement` to
+`4 + uniformPolynomialOrderIncrement`; configured uniform increment stays zero.
+The dated source comment records the rationale. Eleven expectations and their
+comments in six test files change to enforce the new displacement degree and
+derived quadrature, while assertions for all other field orders stay fixed.
+
+The purpose is to remove measured ambient affine representation error and test
+its effect on the weak direction. The shared scalar surface increases from 866
+to 1538 parameters, so this also adds shape modes. Derived pressure-force order
+increases by one, gravity-force order by two, and geometry-inspection request
+from 18 to 24, including the geometry-admissibility sampling; common physical
+diagnostic quadratures remain 10/14 (weak-direction
+observations) and 14/18 (state metrics). The surface metric, normalization and
+preconditioner also change with their FE space. Equal relative linear tolerance
+does not imply equal physical error or an identical line-search norm. The default
+r-squared radial lift stays
+unchanged and is not made globally affine by this ambient enrichment.
+
+Use the separate `cmake-build-release-displacement-q4` directory with the same
+LLVM 22.1.8, Release/IPO, MPI and library configuration. Rebuild all affected
+modules and registered-space tests together. Q4 cases start from the analytic
+seed; no coefficient snapshot is transferred between orders. The fixed
+200-iteration linear budget and eta=0.03 remain unchanged. A 3600-second process
+limit bounds each Q4 numerical case; no silent control relaxation is permitted.
+
+### Completed Q3/Q4 comparison
+
+Q4 source is `4be1a23c3ec2650af6a2a49528a8f5248cb001c3`. Its only production
+change is the compile-time displacement degree at
+`libmeanfield/interface/field/field_registry.cppm:163`; all instrumentation and
+analysis remain in `experiments/`. No physics-facing API changed. Build and
+dependency configuration values match the Q3 build. The Q4 regression campaign
+passed 22 cases and 380,856 assertions: registry/quadrature contracts, field
+spaces, curved geometry, surface/radial/domain Jacobians, Lane–Emden projection,
+and prepared pressure/gravity force Jacobians. Cases, builds and full artifact
+audits were serialized.
+
+`q4_v1` contains the completed seed and live campaigns. The independent audit
+`q4_audit_v1.json` passes 101,282 and 101,310 checks respectively, including raw
+vector algebra, runner hashes, physical-integral identities, and exact buffer,
+normalization, residual and Jacobian-action restoration. An artifact audit pass
+does not turn a failed solve into a successful solve. The endpoint direction is
+explicitly classified `finite_direction_failed_true_gate`.
+
+The analytic seed has identical physical volume diagnostics in both variants.
+Q4 reduces the ambient affine relative value error from 9.01277e-5 to
+2.11629e-16 and gradient error from 1.23472e-3 to 3.24701e-15. Nevertheless,
+the seed Newton direction's stellar displacement L2 grows from 0.00121967 to
+0.0385749, a factor of 31.63. The seed linear solve grows from 26 to 67 iterations;
+both meet the same independently checked eta=0.03 criterion. All three raw small
+Q4 seed trials prepare and exhibit a quadratic residual remainder. This is not
+evidence of a missing first-order derivative in the tested direction.
+
+The fixed six-iteration-budget comparison gives:
+
+| Measurement | Q3 displacement / surface | Q4 displacement / surface |
+|---|---:|---:|
+| Surface parameters | 866 | 1538 |
+| Total root unknowns | 178073 | 178745 |
+| Accepted updates | 6 | 3 |
+| Termination | Nonlinear iteration limit | Fourth linear solve exhausted 200 iterations |
+| First three accepted-step time | 269.31 s | 505.65 s |
+| Residual after three updates / own initial residual | 0.056755 | 0.906188 |
+| Final residual / own initial residual | 0.054710 | 0.906188 |
+| All Newton-iteration time, including failed attempt | 547.34 s | 638.42 s |
+| Outer Krylov iterations, including failed attempt | 354 | 492 |
+| Endpoint probe time, excluded above | 93.07 s | 222.45 s |
+| q18 P(h) potential-energy virial error | 4.22590e-6 | 8.89666e-5 |
+| q18 P(h) force virial error | 2.71367e-6 | 8.27467e-6 |
+| q18 mass error | 3.41685e-6 | 1.96178e-5 |
+| q18 scaled EOS density mismatch RMS | 1.27948e-4 | 1.61156e-4 |
+| q18 hydrostatic-gradient RMS | 0.0106811 | 0.325413 |
+| q18 gravity versus broken potential-gradient RMS | 0.0121494 | 0.326881 |
+| Negative density mass fraction | 2.89683e-5 | 4.19384e-5 |
+
+These are differing stopped endpoints, not equal-accuracy or converged solutions.
+The three-update comparison isolates the accepted-count difference in the table,
+but all the previously listed order/space/quadrature/scaling confounders remain.
+Iteration time excludes initialization, final physical measurements and the
+diagnostic probe; total live-process times were 705.76 s and 944.98 s.
+The stdlib helper `compare_shape_order_runs.py` binds the small comparison inputs
+to runner hashes and reproduces `q3_q4_comparison_v1.json`.
+
+Q3 accepts steps 1, 1/16, 1/16, 1/16, 1/32, 1/32. Q4 accepts approximately
+0.1770, 0.04425, 0.04425. Its fourth linear solve and independently repeated cold
+endpoint solve both stop after 200 iterations with relative full-J defect
+0.0868564, above 0.03. The finite failed direction remains a directional
+observation only. It is not used as an equally accurate Newton comparison.
+
+Q4 also changes the observed physical character of the seed direction: the q14
+Eulerian/pullback norm ratios for density, enthalpy and potential become
+13.0825, 13.4289 and 13.2393. Thus the added deformation is not largely cancelled
+by scalar transport, unlike the late Q3 direction. These are broken scalar
+observations, not a complete distributional or RT transport analysis. The
+largest q10/q14 relative change in Q4 endpoint Eulerian squared integrals is
+4.35e-5; this is an integration sensitivity check, not discretization convergence.
+
+`q4_surface_modes_v1.json` passes its independent angular analysis. At the stopped
+Q4 endpoint, sampled radii span 0.927235–1.049839, or 12.2559% of their nodal
+mean. The Q3 endpoint's observed span is 1.09849%. These are different nodal
+sample sets, not common-point or full-surface extrema. In particular, zero
+initial radius variation at Q4 nodes does not make the curved reference surface
+an exact sphere: the same physical mesh sampled at Q3 nodes has a small nonzero
+span. The Q4 seed direction's degree <=8 nodal fit leaves 0.982810 of its norm;
+its sampled squared norm is 99.9999989% cubic-invariant. The corresponding failed
+endpoint direction leaves 0.988475 of its norm. This remains fine angular
+structure, without identifying which discretization component causes it.
+
+The saved figure `surface_order_comparison_v2.png` uses a common color scale and
+labels the different sample and accepted-update counts. Version 1 is retained
+as an unused layout draft; version 2 fixes the caption/colorbar overlap.
+
+### Conclusions and the next isolated hypotheses
+
+**Do not promote the Q4-only change.** It repairs ambient affine representation
+but worsens useful nonlinear progress, linear cost, interior gradient checks,
+and potential-energy virial balance in this case. This rules out order matching
+alone as a sufficient repair under the tested configuration. It does not show
+that the Q3 representation error contributes nothing, nor that a compatible
+enrichment of the other fields could not help.
+
+The leading hypothesis is an unfavorable combination of surface freedom and
+the resolving power of the coupled field equations. The observed directions
+have strong fine angular structure, and adding surface freedom produced much
+larger physical motion without useful residual reduction. Normalization and
+preconditioning remain possible contributors; full-J residual accuracy does
+not bound physical correction error for a weak operator. More Krylov work might
+recover a linear solution, but the first three gate-passing Q4 directions already
+gave poor accepted progress. No exact nullspace or discrete stability theorem
+has been established for the present mechanical operator.
+
+The next direct experiment should separate independent surface parameters from
+ambient displacement interpolation. A fixed map `q=q0+B*a` with small smooth
+bases of increasing degree can test this hypothesis around the same frozen
+state. Inject physical variations and restrict the existing dual shape force
+with `B^T`; the reduced blocks are the corresponding Galerkin restrictions.
+Pull back the surface metric and preconditioning surrogate consistently, then
+apply state/residual normalization. Do not transpose an already normalized value
+basis and assume it is the physical dual restriction. A surface-mass-orthonormal
+basis is a controlled choice. The deformation contract already separates
+parameter count from boundary-displacement size; the default problem constructor,
+PhysicalRiesz surface metric, and material-surface preconditioners still contain
+nodal assumptions that must be changed together.
+
+That reduced solve changes the approximation/test space. Its excluded original
+shape forces and full original residual must remain visible. A frozen `q0`
+retains its existing surface distortion; success there would justify a later
+analytic-seed sequence with increasing surface resolution, not certify a smooth
+physical root. Suggested first frozen bases use spherical degrees 2/4/6/8
+(9/25/49/81 parameters). High-rotation work must demonstrate adequate shape
+resolution rather than enforce a low-mode answer. This experiment has not been
+implemented here.
+
+Existing `experiments/README.md`, `gravity_completed_rigid_motion.cpp`, and the
+surface-mode/reachability/frequency probes provide relevant mode construction
+and gravity-completion infrastructure to reuse. The proposed addition is the
+fully coupled restricted Newton solve with a consistent dual restriction and
+excluded-force checks. Historical CSVs in the original checkout were not rerun
+or incorporated into this controlled order comparison.
+
+A separate extension test addresses the virial dilation mode. Q4 ambient fields
+can represent the reference position X, but the production logical-radius-squared
+extension E0 does not generate X from its surface trace b. At a retained logical
+radius 1/8 vertex, lifting b gives only one eighth of the required affine value;
+this remains true with Q4 interpolation. A fixed rank-one correction
+`E(z)=E0(z)+(X-E0(b))*ell(z)`, with `ell(b)=1`, preserves the supplied trace and
+gives `E(b)=X`. Its exact transpose is
+`E0^T*f + ell^T*((X-E0(b))^T*f)`; its Jacobian is constant. With an interior
+state displacement `u=E(z)`, direction `b+z` then generates current physical
+dilation `X+u`. Reference-coordinate access, MPI weights, and adjoint virtual
+work require explicit tests. Keep the fixed-infinity vacuum policy unchanged.
+
+This would place the discrete dilation test in the mechanical force space; it
+does not automatically enforce the potential-energy virial using the separate
+mixed potential. The Q4 endpoint demonstrates why both independent virials are
+needed. The rank-one lift also does not remove the observed fine angular modes,
+so it is a separate structural control, not a presumed cure. It remains
+unimplemented and must not be bundled with surface reduction in a causal test.
+
+All reported states remain nonconverged and include negative density samples;
+P(rho)-based virials are undefined there. The required 1e-5 virial milestone and
+internal pointwise validation have not been achieved by this order intervention.
+No claim about higher-n or rapidly rotating convergence follows from this study.

experiments/compare_shape_order_runs.py

@@ -0,0 +1,312 @@
+#!/usr/bin/env python3
+"""Summarize completed Q3/Q4 six-iteration-budget controls; standard library only.
+
+Usage: compare_shape_order_runs.py Q3_ENDPOINT Q4_ENDPOINT --output FRESH.json
+Each endpoint must have its sibling ENDPOINT_run.json completion/hash record.
+The parent campaign provenance.json supplies the recorded mesh/source identities.
+Exit 0 means comparable complete artifacts, including bounded solver failures;
+exit 2 means missing, inconsistent or unsuitable inputs. No physical pass is implied.
+"""
+
+import argparse
+import csv
+import hashlib
+import io
+import json
+import math
+from pathlib import Path
+import statistics
+import sys
+
+
+PHYSICAL = (
+    "mass", "mass_absolute_error", "volume", "minimum_mapping_determinant",
+    "enthalpy_pressure_integral", "enthalpy_virial_error", "enthalpy_virial_signed",
+    "enthalpy_force_virial_error", "enthalpy_force_virial_signed",
+    "eos_density_mismatch_rms", "eos_density_mismatch_rms_scaled",
+    "eos_density_mismatch_maximum", "eos_density_mismatch_maximum_scaled",
+    "minimum_density", "maximum_density", "negative_density_samples",
+    "negative_density_volume_fraction", "negative_density_mass_fraction",
+    "negative_density_mass_magnitude", "minimum_enthalpy", "negative_enthalpy_samples",
+    "bernoulli_raw_rms", "bernoulli_raw_rms_scaled", "hydrostatic_gradient_rms",
+    "gravity_gradient_vs_broken_potential_gradient_rms",
+    "pressure_from_density_defined_everywhere", "pressure_from_enthalpy_defined_everywhere",
+    "virial_error", "force_virial_error", "complete_volume_measurement",
+    "stellar_samples", "valid_stellar_samples", "invalid_stellar_samples",
+    "invalid_density_eos_samples", "invalid_enthalpy_eos_samples",
+    "normalized_unbordered_residual", "central_border", "polytropic_index", "polytropic_constant",
+)
+CONTROLS = (
+    "n", "K", "rho_c", "J", "M", "R_reference", "central_phase", "max_newton",
+    "restart", "max_linear", "linear_tol", "atol", "rtol", "strategy",
+    "gravity_tolerance", "gravity_iterations", "gravity_restart", "weak_shape_final",
+    "context_preconditioner_variant", "context_gravity_schur_sweeps",
+    "context_gravity_schur_relaxation",
+)
+WARNINGS = [
+    "Across-order normalized residual norms and coefficient spectra are not physical accuracy measures. "
+    "Only within-variant residual reductions are compared.",
+    "Q3 and Q4 have different surface sample/parameter spaces; raw nodal norms and unweighted harmonic "
+    "fits are not directly comparable surface integrals. Surface-mode audits are not consumed here.",
+    "q18 denotes the physical diagnostic integration-rule order, not a field order or converged quadrature limit. "
+    "P(h) energy/force virials are reported separately from density-based virials, which may be undefined.",
+    "q14 direction norms describe each endpoint's independently computed Newton direction, not the same "
+    "vector across orders. Density/potential gradients are broken element-local derivatives; jumps are excluded.",
+    "A failed endpoint verified_gate remains a failed linear observation; its direction norms do not describe "
+    "a direction satisfying the requested Newton forcing tolerance.",
+    "The affine representation check projects reference physical coordinates into the displacement space; "
+    "it does not certify that the production r^2 radial lift is affine.",
+    "Six configured iterations are a bounded comparison, not a convergence certificate. "
+    "Use check_plateau_physics.py for the existing solver/virial milestone gate.",
+]
+
+
+def require(condition, message):
+    if not condition:
+        raise ValueError(message)
+
+
+def number(value, name, finite=True):
+    result = float(value)
+    require(not finite or math.isfinite(result), f"Nonfinite {name}")
+    return result
+
+
+def integer(value, name):
+    result = number(value, name)
+    require(result.is_integer(), f"Noninteger {name}")
+    return int(result)
+
+
+def close(a, b, name):
+    require(math.isclose(a, b, rel_tol=1e-10, abs_tol=1e-18), f"Inconsistent {name}: {a} vs {b}")
+
+
+def pairs(rows, key="metric"):
+    result = {}
+    for row in rows:
+        name = row[key]
+        require(name and name not in result, f"Empty or duplicate key: {name!r}")
+        result[name] = row["value"]
+    return result
+
+
+class Inputs:
+    def __init__(self, directory):
+        self.directory = directory.resolve()
+        self.hashes = {}
+
+    def text(self, path):
+        path = path.resolve()
+        data = path.read_bytes()
+        self.hashes[str(path)] = {"sha256": hashlib.sha256(data).hexdigest(), "bytes": len(data)}
+        return data.decode("utf-8")
+
+    def csv(self, name, columns):
+        reader = csv.DictReader(io.StringIO(self.text(self.directory / name)))
+        require(reader.fieldnames is not None and set(columns).issubset(reader.fieldnames),
+                f"Missing columns in {name}: {columns}")
+        require(len(reader.fieldnames) == len(set(reader.fieldnames)), f"Duplicate columns in {name}")
+        rows = list(reader)
+        require(all(None not in row and all(v is not None for v in row.values()) for row in rows),
+                f"Truncated or malformed row in {name}")
+        return rows
+
+
+def physical(source, name, metadata):
+    rows = source.csv(name, ("quadrature", "metric", "value"))
+    selected = pairs([r for r in rows if integer(r["quadrature"], name) == 18])
+    missing = set(PHYSICAL) - selected.keys()
+    require(not missing, f"Missing q18 metrics in {name}: {sorted(missing)}")
+    values = {k: number(selected[k], k, finite=False) for k in PHYSICAL}
+    close(values["polytropic_index"], 1.5, name + " polytropic index")
+    close(values["polytropic_constant"], number(metadata["K"], "K"), name + " polytropic constant")
+    require(values["central_border"] == 0, f"Unexpected central border in {name}")
+    require(math.isfinite(values["normalized_unbordered_residual"]) and
+            values["normalized_unbordered_residual"] >= 0, f"Invalid residual in {name}")
+    # Undefined density-based quantities are observations, not silently zeroed errors.
+    unavailable = {k: str(v) for k, v in values.items() if not math.isfinite(v)}
+    return {"quadrature_order": 18, "metrics": {k: v if math.isfinite(v) else None for k, v in values.items()},
+            "unavailable_metrics": unavailable}
+
+
+def summarize(directory, displacement_order):
+    source = Inputs(directory)
+    metadata = pairs([dict(zip(("metric", "value"), line.split("=", 1)))
+                      for line in source.text(directory / "metadata.txt").splitlines() if line.strip()])
+    require(set(CONTROLS).issubset(metadata), "Missing run configuration metadata")
+    require(number(metadata["n"], "n") == 1.5 and number(metadata["J"], "J") == 0,
+            "Expected n=1.5, J=0")
+    require(integer(metadata["central_phase"], "central_phase") == 0, "Central phase must be off")
+    require(integer(metadata["max_newton"], "max_newton") == 6, "Expected a six-iteration budget")
+    require(metadata.get("converged") in ("0", "1"), "Missing completed solver outcome")
+    accepted = integer(metadata["accepted_steps"], "accepted_steps")
+    require(0 <= accepted <= 6, "Accepted-step count outside configured budget")
+    require(metadata["converged"] == "1" or metadata.get("failure"), "Missing solver-failure explanation")
+    runner_path = directory.parent / (directory.name + "_run.json")
+    runner = json.loads(source.text(runner_path))
+    require(runner.get("finished_utc") and runner.get("exit_code") in (0, 3),
+            "Runner incomplete, interrupted, or exited without a normal solver outcome")
+    require((runner["exit_code"] == 0) == (metadata["converged"] == "1"), "Runner/solver outcome mismatch")
+    provenance = json.loads(source.text(directory.parent / "provenance.json"))
+    provenance = {k: provenance[k] for k in ("commit", "branch", "source", "mesh", "mesh_sha256",
+                                            "binary_sha256", "runtime_environment")}
+    require(all(isinstance(provenance[k], str) and len(provenance[k]) == 64
+                for k in ("mesh_sha256", "binary_sha256")), "Missing recorded mesh/binary hashes")
+
+    fe = pairs(source.csv("weak_shape_physical_metadata.csv", ("key", "value")), key="key")
+    expected = {"density": 2, "enthalpy": 3, "potential": 2, "gravity_flux": 2,
+                "displacement": displacement_order}
+    for field, order in expected.items():
+        require(integer(fe[field + "_family_order"], field) == order, f"Unexpected {field} order")
+    for key, collection in {"reference_mesh_geometry_collection": "H1_3D_P4",
+                            "displacement_collection": f"H1_3D_P{displacement_order}",
+                            "surface_deformation_collection": f"H1_3D_P{displacement_order}"}.items():
+        require(fe[key] == collection, f"Unexpected {key}: {fe[key]}")
+    weak = {k: number(v, k, finite=False) for k, v in
+            pairs(source.csv("weak_shape_summary.csv", ("metric", "value"))).items()}
+    for key in ("accepted_buffers_bitwise_identical", "normalization_factors_bitwise_identical",
+                "direction_buffers_bitwise_identical"):
+        require(weak[key] == 1, f"Endpoint probe changed {key}")
+    require(weak["probe_exception"] == 0 and weak["accepted_commits"] == 0, "Incomplete or mutating endpoint probe")
+
+    layout = {kind: {} for kind in ("value", "residual")}
+    for row in source.csv("weak_shape_layout.csv", ("kind", "block", "offset", "size")):
+        kind, block = row["kind"], row["block"]
+        require(kind in layout and block not in layout[kind], "Invalid/duplicate layout block")
+        size = integer(row["size"], "block size")
+        require(size >= 0 and integer(row["offset"], "block offset") == sum(layout[kind].values()),
+                "Noncontiguous or negative block layout")
+        layout[kind][block] = size
+    require(all(sum(blocks.values()) == weak["state_size"] for blocks in layout.values()), "Layout size mismatch")
+    require(layout["value"]["surface_deformation"] == layout["residual"]["surface_shape_balance"],
+            "Surface parameter/equation count mismatch")
+
+    initial = physical(source, "initial_physical.csv", metadata)
+    final = physical(source, "final_physical.csv", metadata)
+    initial_f = initial["metrics"]["normalized_unbordered_residual"]
+    final_f = final["metrics"]["normalized_unbordered_residual"]
+    require(initial_f > 0, "Initial residual must be positive for this comparison")
+    close(final_f, weak["base_norm"], "final physical/probe residual")
+    history = source.csv("newton.csv", ("iteration", "accepted", "F", "alpha", "trials", "linear_status",
+                         "linear_iterations", "linear_true", "seconds", "linear_seconds",
+                         "prepare_seconds", "refresh_seconds", "linear_rhs_norm"))
+    steps = []
+    for i, row in enumerate(history):
+        require(integer(row["iteration"], "iteration") == i and i < 6, "Invalid iteration sequence")
+        step = {k: integer(row[k], k) for k in ("iteration", "accepted", "trials", "linear_status", "linear_iterations")}
+        require(step["accepted"] in (0, 1) and step["trials"] >= 0 and step["linear_iterations"] >= 0,
+                "Invalid iteration count/flag")
+        step.update({k: number(row[k], k) for k in ("F", "alpha", "seconds", "linear_seconds",
+                    "prepare_seconds", "refresh_seconds")})
+        require(all(step[k] >= 0 for k in ("F", "seconds", "linear_seconds", "prepare_seconds", "refresh_seconds")),
+                "Negative residual/timing")
+        relative = number(row["linear_true"], "linear_true", finite=False)
+        step["linear_relative_true_defect"] = relative if math.isfinite(relative) else None
+        step["linear_relative_true_defect_unavailable"] = None if math.isfinite(relative) else str(relative)
+        if step["accepted"]:
+            require(0 < step["alpha"] <= 1 and step["linear_status"] == 0 and math.isfinite(relative),
+                    "Accepted step lacks a successful finite linear report")
+            require(relative <= number(metadata["linear_tol"], "linear_tol") * (1 + 1e-10) + 1e-15,
+                    "Accepted step exceeds the recorded linear tolerance")
+        else:
+            require(i == len(history) - 1, "Nonfinal rejected iteration")
+        steps.append(step)
+    require(sum(s["accepted"] for s in steps) == accepted, "History/metadata accepted count mismatch")
+    if steps:
+        close(steps[-1]["F"], final_f, "last recorded/final residual")
+        rhs = number(history[0]["linear_rhs_norm"], "first linear RHS", finite=False)
+        if math.isfinite(rhs):
+            close(rhs, initial_f, "initial residual/first linear RHS")
+    else:
+        close(initial_f, final_f, "zero-iteration endpoint residual")
+
+    direction = {}
+    physical_rows = source.csv("weak_shape_physical_summary.csv", ("quadrature", "field", "metric", "value"))
+    for row in physical_rows:
+        if integer(row["quadrature"], "direction quadrature") == 14:
+            field = direction.setdefault(row["field"], {})
+            require(row["metric"] not in field, "Duplicate q14 direction metric")
+            field[row["metric"]] = number(row["value"], "q14 direction metric")
+    for field in ("density", "specific_enthalpy", "gravity_potential"):
+        require(all(k in direction.get(field, {}) for k in
+                    ("pullback_l2_norm", "transport_l2_norm", "eulerian_l2_norm")), f"Missing q14 {field} norms")
+    require(all(k in direction.get("affine", {}) for k in ("relative_value_l2_error", "relative_gradient_l2_error")),
+            "Missing q14 affine errors")
+    squared = direction["geometry"]["direction_displacement_squared_integral"]
+    require(squared >= 0, "Negative displacement squared integral")
+    direction["geometry"]["direction_displacement_l2_norm"] = math.sqrt(squared)
+    for path, digest in source.hashes.items():
+        if Path(path).parent == source.directory:
+            require(runner["output_hashes"].get(Path(path).name) == digest["sha256"],
+                    f"Runner hash missing/mismatched: {path}")
+
+    times = [s["seconds"] for s in steps]
+    endpoint = {k: weak[k] for k in ("state_size", "linear_status", "linear_iterations", "linear_seconds",
+                                   "verified_gate", "verified_relative_defect", "total_probe_seconds")}
+    endpoint_unavailable = {k: str(v) for k, v in endpoint.items() if not math.isfinite(v)}
+    return {"directory": str(source.directory), "input_files": source.hashes,
+            "recorded_campaign_provenance": provenance,
+            "controls": {k: metadata[k] for k in CONTROLS}, "orders": expected,
+            "finite_element_metadata": fe, "block_sizes": layout,
+            "outcome": {"runner_completed": True, "runner_exit_code": runner["exit_code"],
+                        "converged": metadata["converged"] == "1", "failure": metadata.get("failure"),
+                        "accepted_updates": accepted, "completed_six_updates": accepted == 6},
+            "newton": {"recorded_iterations": len(steps), "steps": steps,
+                       "accepted_alphas": [s["alpha"] for s in steps if s["accepted"]],
+                       "krylov_iterations_total": sum(s["linear_iterations"] for s in steps),
+                       "iteration_seconds_total": sum(times),
+                       "iteration_seconds_mean": statistics.mean(times) if times else None,
+                       "iteration_seconds_median": statistics.median(times) if times else None,
+                       "linear_seconds_total": sum(s["linear_seconds"] for s in steps),
+                       "preparation_seconds_total": sum(s["prepare_seconds"] for s in steps),
+                       "refresh_seconds_total": sum(s["refresh_seconds"] for s in steps),
+                       "timing_scope": "AfterIteration records only; excludes initialization, final physical measurement "
+                                       "and endpoint weak-shape probe. Components are not asserted to partition wall time.",
+                       "initial_normalized_residual": initial_f, "final_normalized_residual": final_f,
+                       "final_over_initial_normalized_residual": final_f / initial_f},
+            "initial_physical_q18": initial, "final_physical_q18": final,
+            "endpoint_probe": {k: v if math.isfinite(v) else None for k, v in endpoint.items()},
+            "endpoint_probe_unavailable_metrics": endpoint_unavailable,
+            "endpoint_direction_q14": direction}
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("q3", type=Path)
+    parser.add_argument("q4", type=Path)
+    parser.add_argument("--output", required=True, type=Path)
+    args = parser.parse_args()
+    try:
+        require(args.q3.resolve() != args.q4.resolve(), "Q3 and Q4 must be distinct endpoint directories")
+        require(not args.output.exists(), "Output must be fresh")
+        require(not any(args.output.resolve().is_relative_to(p.resolve()) for p in (args.q3, args.q4)),
+                "Write the comparison outside the immutable endpoint directories")
+        q3, q4 = summarize(args.q3.resolve(), 3), summarize(args.q4.resolve(), 4)
+        require(q3["controls"] == q4["controls"], "Physics/solver controls differ beyond displacement order")
+        require(q3["recorded_campaign_provenance"]["mesh_sha256"] ==
+                q4["recorded_campaign_provenance"]["mesh_sha256"], "Recorded reference meshes differ")
+        require(q3["recorded_campaign_provenance"]["runtime_environment"] ==
+                q4["recorded_campaign_provenance"]["runtime_environment"], "Recorded runtime environments differ")
+        result = {"schema_version": 1, "artifact_comparison_complete": True,
+                  "scope": "Recorded endpoint comparison; no FEM rerun or physical acceptance certificate",
+                  "comparison": {"matched_maximum_newton_iterations": 6,
+                                 "equal_accepted_counts": q3["outcome"]["accepted_updates"] ==
+                                                          q4["outcome"]["accepted_updates"],
+                                 "both_completed_six_updates": q3["outcome"]["completed_six_updates"] and
+                                                               q4["outcome"]["completed_six_updates"]},
+                  "warnings": WARNINGS, "q3": q3, "q4": q4}
+        # Deliberately omit an across-order normalized-residual ratio.
+        text = json.dumps(result, indent=2, allow_nan=False) + "\n"
+        with args.output.open("x") as stream:
+            stream.write(text)
+        print(f"Wrote {args.output}; accepted updates: Q3={q3['outcome']['accepted_updates']}, "
+              f"Q4={q4['outcome']['accepted_updates']}. Not a physical-validation certificate.")
+        return 0
+    except (OSError, ValueError, KeyError, TypeError) as error:
+        print(f"Comparison incomplete or unsuitable: {error}", file=sys.stderr)
+        return 2
+
+
+if __name__ == "__main__":
+    sys.exit(main())

experiments/plot_shape_surface_comparison.py

@@ -28,6 +28,8 @@ def main():
             raise ValueError("This comparison requires n=1.5, J=0, and Q3 followed by Q4 displacement.")
         with (directory / 'weak_shape_physical_surface_nodes.csv').open() as stream:
             rows = list(csv.DictReader(stream))
+        with (directory / 'newton.csv').open() as stream:
+            accepted = sum(int(row['accepted']) for row in csv.DictReader(stream))
         points = np.array([[float(row['unit_radial_'+a]) for a in 'xyz'] for row in rows])
         radii = np.array([float(row['physical_radius']) for row in rows])
         if not np.isfinite(points).all() or not np.isfinite(radii).all() or np.any(radii <= 0):
@@ -35,20 +37,21 @@ def main():
         latitude = np.arcsin(np.clip(points[:, 2], -1, 1))
         longitude = np.arctan2(points[:, 1], points[:, 0])
         deviation = 100 * (radii / np.mean(radii) - 1)
-        data.append((longitude, latitude, deviation, len(rows)))
+        data.append((longitude, latitude, deviation, len(rows), accepted))
     limit = max(float(np.max(np.abs(row[2]))) for row in data)
     if limit == 0:
         limit = 1e-12
     fig, axes = plt.subplots(1, 2, figsize=(12, 5), subplot_kw={'projection': 'mollweide'}, layout='constrained')
-    for axis, (longitude, latitude, deviation, count), title in zip(axes, data, ('Q3 displacement / surface', 'Q4 displacement / surface')):
+    fig.get_layout_engine().set(rect=(0, .07, 1, .93))
+    for axis, (longitude, latitude, deviation, count, accepted), title in zip(axes, data, ('Q3 displacement / surface', 'Q4 displacement / surface')):
         picture = axis.scatter(longitude, latitude, c=deviation, s=14, cmap='RdBu_r', vmin=-limit, vmax=limit, linewidths=0)
-        axis.set_title(f'{title}\n{count} radial nodes')
+        axis.set_title(f'{title}\n{accepted} accepted updates; {count} radial nodes')
         axis.grid(True, alpha=.25)
         axis.tick_params(labelsize=8)
     fig.colorbar(picture, ax=list(axes), orientation='horizontal', shrink=.65, pad=.08,
                  label='Radius relative to its unweighted nodal mean (%)')
     fig.suptitle('Nonrotating n = 1.5: surface samples after the bounded order trial', fontsize=14)
-    fig.text(.5, .01, 'Common color scale. Nodal observations only; neither endpoint is certified by this plot.', ha='center', fontsize=9)
+    fig.text(.5, .01, 'Common color scale; different nodal sampling. Neither endpoint is certified by this plot.', ha='center', fontsize=9)
     args.output.parent.mkdir(parents=True, exist_ok=True)
     fig.savefig(args.output, dpi=200)
     plt.close(fig)