From 6ceda61bddc0d4881d9f7ac98d170d07980a963d Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 18 Sep 2026 14:27:15 -0400 Subject: [PATCH] Preserve completed weak-response study while retaining Q3 baseline --- .../SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md | 220 +++++++++++- experiments/compare_shape_order_runs.py | 312 ++++++++++++++++++ experiments/plot_shape_surface_comparison.py | 11 +- 3 files changed, 525 insertions(+), 18 deletions(-) create mode 100644 experiments/compare_shape_order_runs.py diff --git a/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md b/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md index cb5918e..0f7a523 100644 --- a/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md +++ b/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 -structure, not simple expansion, translation, or quadrupolar deformation. -This is consistent with mesh/discretization coupling, without proving its cause. +Degree <=2 leaves 0.981646. Thus most of the nodal direction is unresolved by +degree <=8 fits and is nearly cubic-invariant. This is consistent with +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 -improve linear algebra, but cannot itself change the exact Newton direction or -remove a true weak physical response. Shape representation, step control, and -operator formulation remain distinct hypotheses. +An invertible linear/tangent coordinate or preconditioner change preserves the +exact mapped physical Newton correction at a fixed state. It can improve +coordinate-dependent singular values and inexact linear solves. A nonlinear +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 -affine representation error and test the effect on the observed weak direction; -the additional surface modes and quadrature changes are explicit confounders. +Applied 2026-09-18 at approximately 17:45 UTC on +`research/shape-field-displacement-q4-20260918`, branched from baseline `211be17`. +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. diff --git a/experiments/compare_shape_order_runs.py b/experiments/compare_shape_order_runs.py new file mode 100644 index 0000000..cddcd05 --- /dev/null +++ b/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()) diff --git a/experiments/plot_shape_surface_comparison.py b/experiments/plot_shape_surface_comparison.py index c7128d7..c75d218 100644 --- a/experiments/plot_shape_surface_comparison.py +++ b/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)