diff --git a/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md b/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md index 6e9ec1e..cb5918e 100644 --- a/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md +++ b/experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md @@ -62,4 +62,99 @@ comparison, not evidence of high-n/high-rotation convergence. ## Results -Pending recorded runs. No production default has been changed on the baseline. +Baseline observations are recorded below. No production default has been changed on this baseline branch. + +### Reproducible baseline + +Evidence root: `../../diagnostics/shape_field_weak_response_2026-09-18` relative +to this checkout. Frozen probes `baseline_probes_v1` use source `d79db94`; the +six-update control `q3_live_v1` uses `2c4a5e0`. Both include immutable binaries, +mesh/input hashes, commands, clean-source provenance, and independent elapsed +and CPU clocks. All full-J linear gates passed. Independent audits +`baseline_probes_audit_v2.json` and `q3_live_audit_v1.json` pass; accepted buffers, +factors and restored actions agree exactly. Q3 curved/tensor/geometry regression: +5 cases, 21,368 assertions passed. + +The first checker version omitted the new cache-preservation invariant from its +schema. Version 2 requires it. A small finite synthetic matrix reproduced +NumPy/Accelerate matmul warnings; direct finite-checked contractions avoid that +backend path, without changing tolerances. Raw FEM artifacts remain unchanged. +The artifact checker passes 38 tiny algebra/corruption fixtures. The angular +analyzer independently checks its degree-8 recurrence, orthogonality, cubic +symmetries, rank, and file binding on analytic fixtures. + +| Observation | Seed | Old late state (12 updates) | Fresh control (6 updates) | Coupled endpoint (6 updates) | +|---|---:|---:|---:|---:| +| Normalized residual | 2.73659e-4 | 1.39670e-5 | 1.49720e-5 | 1.24583e-5 | +| Sampled block-span minimum singular response | 5.80574e-2 | 1.43543e-4 | 8.16558e-4 | 6.94894e-4 | +| Sampled surface radius span / nodal mean | 3.42403e-5 | 1.24122e-2 | 1.09849e-2 | 1.11125e-2 | + +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, +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. + +Late q14 Eulerian/pullback change norm ratios are density 0.05769, enthalpy +0.18870, potential 0.01288, versus seed 0.3243, 0.3080, 0.1993. This is substantial +transport compensation but not physically invisible motion. Broken DG gradients +omit jumps. The maximum q10/q14 relative difference in the three original cases' +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; +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. + +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, +but adding surface modes could also worsen weak response. + +### Structural interpretation and a rejected shortcut + +At J=0, bulk Bernoulli rows are unweighted moments. Their direct shape derivative +at fixed reference-field coefficients is `integral(q * B * div(v))`; the EOS +shape column similarly contains its local closure imbalance. Thus pointwise +near-equilibrium shape information increasingly relies on gravity and mechanical +coupling. Small projected residuals alone do not prove small pointwise imbalance. +Current shape rows are the prescribed lift's transpose applied to volume force. +For regular surface h~s, density~s^n and pressure~s^(n+1), weakening near-surface +coupling as n increases is plausible. DG jumps and the independent RT field +prevent concluding exact discrete redundancy from the continuum identities. + +Relevant assembly: `libmeanfield/impl/operators/prepared_hydrostatic_equilibrium.cpp` +(lines 864–911, 1354–1397), `operators/kernels/barotropic_closure_kernels.cpp` +(line 655), and `operators/prepared_stellar_equilibrium.cpp` (lines 806–811). + +A naive replacement of mechanical shape rows by surface Bernoulli moments is +rank deficient at the desired pointwise B=0, J=0 root with the current trace +spaces. The actual closed stellar surface has 96 quads, 192 edges, 98 vertices. +Continuous surface Qp has 96*p*p+2 degrees of freedom; broken Q(p-1) potential +traces supply only 96*p*p. Q3 surface therefore has 866 rows versus 864 DG Q2 +stellar trace coefficients, after eliminating enforced h=0 trace variations. +C adds no independent trace dimension; the direct shape measure term vanishes +with B. Surface-only Q4 widens this hypothetical deficit to 1538-864=674. +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. + +### 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. diff --git a/experiments/analyze_shape_surface_modes.py b/experiments/analyze_shape_surface_modes.py new file mode 100644 index 0000000..26ef777 --- /dev/null +++ b/experiments/analyze_shape_surface_modes.py @@ -0,0 +1,366 @@ +#!/usr/bin/env python3 +"""Describe saved radial surface directions by unweighted nodal angular fits. + +These are projections of the sampled coefficient values, not surface-integral +harmonic coefficients. Cubic symmetry can suggest a mesh-related pattern, but +does not establish its cause, physical accuracy, or a nullspace of the solver. +No FEM operations or full-state vector reads are performed. +""" +import argparse +import csv +import hashlib +import itertools +import json +import math +import os +from pathlib import Path +import sys + +for _name in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"): + os.environ[_name] = "1" +import numpy as np + +SURFACE_FILE = "weak_shape_physical_surface_nodes.csv" +DEGREES = (0, 2, 4, 6, 8) +EPS = np.finfo(float).eps +NODE_TOLERANCE = 2e-10 + + +def digest(path): + return hashlib.sha256(Path(path).read_bytes()).hexdigest() + + +def contract(subscripts, *arrays): + # Direct contraction avoids the Accelerate matmul flag defect documented + # in audit_shape_field_weak.py; finite results remain an explicit check. + if any(not np.isfinite(a).all() for a in arrays): + raise ValueError("Nonfinite contraction input") + result = np.einsum(subscripts, *arrays, optimize=False) + if not np.isfinite(result).all(): + raise ValueError("Nonfinite contraction output") + return result + + +def norm(array): + result = math.sqrt(float(contract("i,i->", array, array))) + if not math.isfinite(result): + raise ValueError("Nonfinite norm") + return result + + +def harmonics(points, maximum_degree=8): + """Continuous-L2-normalized real Y_lm, including Condon--Shortley phase. + + Recurring on normalized associated Legendre functions avoids separately + generating large factorials and then cancelling them in normalization. + Negative m denotes sine, positive m cosine, with sqrt(2) for nonzero m. + """ + if maximum_degree < 0 or maximum_degree > 8: + raise ValueError("This bounded diagnostic supports degrees zero through eight") + z, phi = points[:, 2], np.arctan2(points[:, 1], points[:, 0]) + radial = np.hypot(points[:, 0], points[:, 1]) + legendre = {(0, 0): np.full(len(points), 1 / math.sqrt(4 * math.pi))} + for m in range(maximum_degree + 1): + if m: + legendre[m, m] = -math.sqrt((2 * m + 1) / (2 * m)) * radial * legendre[m - 1, m - 1] + if m < maximum_degree: + legendre[m + 1, m] = math.sqrt(2 * m + 3) * z * legendre[m, m] + for degree in range(m + 2, maximum_degree + 1): + a = math.sqrt((4 * degree * degree - 1) / (degree * degree - m * m)) + b = math.sqrt(((2 * degree + 1) * ((degree - 1) ** 2 - m * m)) / + ((2 * degree - 3) * (degree * degree - m * m))) + legendre[degree, m] = a * z * legendre[degree - 1, m] - b * legendre[degree - 2, m] + labels, columns = [], [] + for degree in range(maximum_degree + 1): + for m in range(-degree, degree + 1): + value = legendre[degree, abs(m)] + if m: + value = math.sqrt(2) * value * (np.sin(-m * phi) if m < 0 else np.cos(m * phi)) + labels.append({"degree": degree, "m": m}) + columns.append(value) + matrix = np.column_stack(columns) + if not np.isfinite(matrix).all(): + raise ValueError("Nonfinite harmonic basis") + return matrix, labels + + +def fit(matrix, values, column_scale=False): + """Tall-matrix SVD with explicit reconstruction/orthogonality checks.""" + if not np.isfinite(matrix).all() or not np.isfinite(values).all(): + raise ValueError("Nonfinite fit input") + scales = np.sqrt(contract("ij,ij->j", matrix, matrix)) if column_scale else np.ones(matrix.shape[1]) + scales = np.where(scales > 0, scales, 1.) + scaled = matrix / scales + u, singular, vh = np.linalg.svd(scaled, full_matrices=False) + if not all(np.isfinite(a).all() for a in (u, singular, vh)): + raise ValueError("Nonfinite SVD") + cutoff = EPS * max(matrix.shape) * singular[0] + rank = int(np.count_nonzero(singular > cutoff)) + scores = contract("ij,i->j", u[:, :rank], values) + coefficients = contract("ji,j->i", vh[:rank], scores / singular[:rank]) / scales + projected = contract("ij,j->i", matrix, coefficients) + remainder = values - projected + reconstruction = contract("ij,jk->ik", u * singular, vh) + factorization_error = norm((reconstruction - scaled).ravel()) + factorization_limit = 4096 * EPS * max(norm(scaled.ravel()), 1.) + if factorization_error > factorization_limit: + raise ValueError("SVD reconstruction failed") + # The retained singular directions define the actual numerical projection; + # omitted singular directions need not be resolved by the node set. + orthogonality = norm(contract("ij,i->j", u[:, :rank], remainder)) + magnitude = norm(values) + limit = 8192 * EPS * max(magnitude, norm(projected), 1e-300) + if orthogonality > limit: + raise ValueError("Projection orthogonality failed") + error = norm(remainder) + return {"sample_count": len(values), "column_count": matrix.shape[1], "rank": rank, + "full_column_rank": rank == matrix.shape[1], "rank_cutoff": float(cutoff), + "singular_values": singular.tolist(), "column_scaling": "unit nodal L2" if column_scale else "none", + "column_scales": scales.tolist(), + "retained_condition_ratio": float(singular[0] / singular[rank - 1]) if rank else None, + "basis_coefficients": coefficients.tolist(), "nodal_l2_remainder": error, + "relative_nodal_l2_remainder": error / magnitude if magnitude else None, + "explained_nodal_squared_norm_fraction": 1 - (error / magnitude) ** 2 if magnitude else None, + "svd_reconstruction_error": factorization_error, "retained_projection_orthogonality_error": orthogonality} + + +def cubic_columns(points): + squared = points * points + c4 = np.sum(squared * squared, axis=1) - 3 / 5 + # On the unit sphere, the reflection/permutation-invariant polynomial ring + # through degree eight is spanned by 1, e2, e3, e2^2. C4 is an affine + # replacement for e2; e3's spherical mean is 1/105. These columns are not + # individually harmonic or orthogonal, except the pure l=4 C4 column. + e3 = np.prod(squared, axis=1) - 1 / 105 + return np.column_stack((np.ones(len(points)), c4, e3, c4 * c4)) + + +def symmetry_projection(points, values): + """Project onto all signed-permutation invariants if nodes form orbits. + + A permutation of the entire node set preserves the unweighted metric. + Without that closure, a nodal group average is not a justified orthogonal + projection; invariant polynomial fits remain separately available. + """ + buckets = {} + for i, point in enumerate(points): + key = tuple(np.rint(point / NODE_TOLERANCE).astype(np.int64)) + buckets.setdefault(key, []).append(i) + shifts = tuple(itertools.product((-1, 0, 1), repeat=3)) + transforms = [] + maximum_error = 0. + for permutation in itertools.permutations(range(3)): + for signs in itertools.product((-1, 1), repeat=3): + image_indices = [] + for point in points[:, permutation] * np.asarray(signs): + key = tuple(np.rint(point / NODE_TOLERANCE).astype(np.int64)) + candidates = [] + for shift in shifts: + for index in buckets.get(tuple(key[k] + shift[k] for k in range(3)), ()): + distance = float(np.max(np.abs(point - points[index]))) + if distance <= NODE_TOLERANCE: + candidates.append((distance, index)) + if len(candidates) != 1: + return {"status": "unavailable_node_cloud_not_uniquely_closed", "node_tolerance": NODE_TOLERANCE, + "matching_node_count_at_failure": len(candidates), + "scope": "No full cubic projection claimed; polynomial invariant fits remain valid."} + distance, index = candidates[0] + maximum_error = max(maximum_error, distance) + image_indices.append(index) + if len(set(image_indices)) != len(points): + raise ValueError("Cubic transformation did not define a nodal permutation") + transforms.append(image_indices) + indices = np.asarray(transforms, dtype=int) + projected = np.mean(values[indices], axis=0) + error = norm(values - projected) + magnitude = norm(values) + idempotency = norm(np.mean(projected[indices], axis=0) - projected) + cross = abs(float(contract("i,i->", projected, values - projected))) + if idempotency > 4096 * EPS * max(magnitude, 1e-300) or cross > 4096 * EPS * max(magnitude ** 2, 1e-300): + raise ValueError("Cubic nodal projection identities failed") + return {"status": "closed_48_element_group", "node_tolerance": NODE_TOLERANCE, + "maximum_node_matching_error": maximum_error, "nodal_l2_remainder": error, + "relative_nodal_l2_remainder": error / magnitude if magnitude else None, + "invariant_nodal_squared_norm_fraction": 1 - (error / magnitude) ** 2 if magnitude else None, + "projection_idempotency_error": idempotency, "projection_cross_term": cross, + "scope": "All signed coordinate permutations in the saved Cartesian axes; unweighted nodal group average."} + + +def read_surface(path): + with path.open(newline="") as stream: + reader = csv.DictReader(stream) + required = {"parameter", "volume_scalar_true_dof", "unit_radial_x", "unit_radial_y", "unit_radial_z", + "reference_radius", "physical_radius", "base_parameter", "physical_direction_parameter"} + if not reader.fieldnames or len(set(reader.fieldnames)) != len(reader.fieldnames) or not required <= set(reader.fieldnames): + raise ValueError("Missing or duplicate surface CSV columns") + rows = list(reader) + if not rows or any(None in r or any(v is None for v in r.values()) for r in rows): + raise ValueError("Empty or malformed surface CSV") + data = {k: np.asarray([float(r[k]) for r in rows]) for k in required} + if any(not np.isfinite(a).all() for a in data.values()): + raise ValueError("Nonfinite surface data") + if not np.array_equal(data["parameter"], np.arange(len(rows))): + raise ValueError("Noncontiguous surface parameter indices") + identifiers = data["volume_scalar_true_dof"] + if np.any(identifiers < 0) or np.any(identifiers != np.floor(identifiers)) or len(set(identifiers)) != len(rows): + raise ValueError("Invalid or repeated surface degree of freedom") + if np.any(data["reference_radius"] <= 0) or np.any(data["physical_radius"] <= 0): + raise ValueError("Nonpositive surface radius") + discrepancy = np.abs(data["physical_radius"] - data["reference_radius"] - data["base_parameter"]) + if np.any(discrepancy > 2048 * EPS * (data["physical_radius"] + data["reference_radius"] + np.abs(data["base_parameter"]))): + raise ValueError("Physical radius identity failed") + points = np.column_stack([data["unit_radial_" + axis] for axis in "xyz"]) + lengths = np.sqrt(contract("ij,ij->i", points, points)) + unit_error = float(np.max(np.abs(lengths - 1))) + if unit_error > 5e-10: + raise ValueError("Surface radial directions are not unit length") + # Remove decimal-rounding length error before the spherical recurrence. + return points / lengths[:, None], data["physical_direction_parameter"], unit_error, data + + +def analyze(path): + path = Path(path).resolve() + if path.is_dir(): + path /= SURFACE_FILE + points, values, unit_error, surface_data = read_surface(path) + source_hash = digest(path) + binding = {"status": "input_hash_only"} + run_path = path.parent.parent / (path.parent.name + "_run.json") + if run_path.exists(): + run = json.loads(run_path.read_text()) + if run.get("output_hashes", {}).get(path.name) != source_hash: + raise ValueError("Surface file does not match runner output hash") + binding = {"status": "runner_surface_hash_checked", "run_path": str(run_path), "run_sha256": digest(run_path), + "scope": "Only this surface artifact binding; no independent audit of the complete FEM campaign."} + matrix, labels = harmonics(points) + harmonic_fits = [{"maximum_degree": degree, **fit(matrix[:, :(degree + 1) ** 2], values)} for degree in DEGREES] + cubic = cubic_columns(points) + invariant_fits = [{"maximum_degree": degree, **fit(cubic[:, :count], values, True)} + for degree, count in ((0, 1), (2, 1), (4, 2), (6, 3), (8, 4))] + return {"surface_file": str(path), "surface_sha256": source_hash, "binding": binding, + "sample_count": len(points), "maximum_original_unit_length_error": unit_error, + "base_surface_radius": { + "scope": "unweighted nodal statistics; extrema bound observed asphericity, not full-surface extrema", + "reference_minimum": float(np.min(surface_data["reference_radius"])), + "reference_maximum": float(np.max(surface_data["reference_radius"])), + "physical_minimum": float(np.min(surface_data["physical_radius"])), + "physical_maximum": float(np.max(surface_data["physical_radius"])), + "physical_mean": float(np.mean(surface_data["physical_radius"])), + "physical_span_over_mean": float(np.ptp(surface_data["physical_radius"]) / np.mean(surface_data["physical_radius"])), + "maximum_absolute_base_parameter": float(np.max(np.abs(surface_data["base_parameter"])))}, + "direction_nodal_l2": norm(values), "direction_maximum_absolute": float(np.max(np.abs(values))), + "harmonic_basis": {"normalization": "continuous unit-sphere L2 orthonormal real Y_lm with Condon--Shortley phase", + "ordering": labels, "negative_m": "sine", "positive_m": "cosine"}, + "harmonic_fits": harmonic_fits, + "cubic_invariant_basis": ["1", "x^4+y^4+z^4-3/5", "x^2*y^2*z^2-1/105", "(x^4+y^4+z^4-3/5)^2"], + "cubic_invariant_fits": invariant_fits, "full_cubic_nodal_projection": symmetry_projection(points, values)} + + +def self_test(): + # Product quadrature integrates products through l=8 exactly, independently + # checking the recurrence's normalization, phase and degree organization. + z, weights = np.polynomial.legendre.leggauss(9) + phi = np.arange(18) * (2 * math.pi / 18) + points = np.asarray([(math.sqrt(1 - q * q) * math.cos(t), math.sqrt(1 - q * q) * math.sin(t), q) for q in z for t in phi]) + quadrature = np.repeat(weights, len(phi)) * (2 * math.pi / len(phi)) + matrix, labels = harmonics(points) + gram = contract("ki,k,kj->ij", matrix, quadrature, matrix) + assert np.max(np.abs(gram - np.eye(81))) < 8e-14 + assert np.max(np.abs(matrix[:, 2] - math.sqrt(3 / (4 * math.pi)) * points[:, 2])) < 1e-15 + assert np.max(np.abs(matrix[:, 3] + math.sqrt(3 / (4 * math.pi)) * points[:, 0])) < 1e-15 + for degree in range(9): + columns = matrix[:, degree * degree:(degree + 1) ** 2] + assert np.max(np.abs(np.sum(columns * columns, axis=1) - (2 * degree + 1) / (4 * math.pi))) < 2e-14 + cubic = cubic_columns(points) + assert fit(matrix[:, :9], cubic[:, 1])["relative_nodal_l2_remainder"] > .8 + assert fit(matrix[:, :25], cubic[:, 1])["relative_nodal_l2_remainder"] < 1e-13 + for degree in DEGREES: + values = matrix[:, degree * degree] + result = fit(matrix[:, :(degree + 1) ** 2], values) + assert result["rank"] == (degree + 1) ** 2 and result["relative_nodal_l2_remainder"] < 1e-13 + assert fit(matrix[:, :1], np.zeros(len(points)))["relative_nodal_l2_remainder"] is None + assert fit(np.column_stack((matrix[:, 0], matrix[:, 0])), matrix[:, 0])["rank"] == 1 + cloud = [] + for seed in ((1, 2, 3), (1, 3, 5), (2, 3, 7), (1, 4, 9)): + for permutation in itertools.permutations(seed): + for signs in itertools.product((-1, 1), repeat=3): + point = np.asarray(permutation, dtype=float) * signs + cloud.append(point / np.linalg.norm(point)) + cloud = np.asarray(cloud) + cubic = cubic_columns(cloud) + invariant = contract("ij,j->i", cubic, np.asarray((.4, 2., -7., .6))) + assert fit(cubic, invariant, True)["relative_nodal_l2_remainder"] < 1e-13 + assert symmetry_projection(cloud, invariant)["relative_nodal_l2_remainder"] < 1e-13 + antisymmetric = cloud[:, 0] + result = symmetry_projection(cloud, invariant + antisymmetric) + assert abs(result["nodal_l2_remainder"] - norm(antisymmetric)) < 1e-12 + assert symmetry_projection(cloud[:-1], invariant[:-1])["status"] == "unavailable_node_cloud_not_uniquely_closed" + # Bind the analysis to a tiny real CSV and reject corruption in both parser + # and runner certificate, without reading any campaign artifacts. + import tempfile + with tempfile.TemporaryDirectory() as tmp: + root = Path(tmp) + folder = root / "case" + folder.mkdir() + path = folder / SURFACE_FILE + columns = ["parameter", "volume_scalar_true_dof", "unit_radial_x", "unit_radial_y", "unit_radial_z", + "reference_radius", "physical_radius", "base_parameter", "physical_direction_parameter"] + with path.open("w", newline="") as stream: + writer = csv.writer(stream) + writer.writerow(columns) + for i, (point, value) in enumerate(zip(cloud, invariant)): + writer.writerow([i, i, *point, 1., 1., 0., value]) + run = root / "case_run.json" + run.write_text(json.dumps({"output_hashes": {SURFACE_FILE: digest(path)}})) + assert analyze(folder)["binding"]["status"] == "runner_surface_hash_checked" + run.write_text(json.dumps({"output_hashes": {SURFACE_FILE: "0" * 64}})) + try: + analyze(folder) + except ValueError as exc: + assert "runner output hash" in str(exc) + else: + raise AssertionError("Corrupted source hash passed") + path.write_text(path.read_text().replace(",1.0,1.0,0.0,", ",nan,1.0,0.0,", 1)) + try: + read_surface(path) + except ValueError as exc: + assert "Nonfinite" in str(exc) + else: + raise AssertionError("Nonfinite source passed") + print("Tiny harmonic, invariant, symmetry, rank and artifact-binding fixtures passed") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("inputs", nargs="*", help="Probe directories or saved surface-node CSV files") + parser.add_argument("--output", type=Path) + parser.add_argument("--self-test", action="store_true") + args = parser.parse_args() + if args.self_test: + if args.inputs or args.output: + parser.error("Self-test does not consume campaign inputs or create an output") + self_test() + return 0 + if not args.inputs or not args.output: + parser.error("Provide surface inputs and a fresh output JSON") + output = args.output.resolve() + if output.exists() or any(output.is_relative_to(Path(p).resolve() if Path(p).is_dir() else Path(p).resolve().parent) for p in args.inputs): + parser.error("Use a fresh output outside the immutable input directories") + cases = [] + for path in args.inputs: + try: + cases.append({"analysis_pass": True, **analyze(path)}) + except (ValueError, OSError, KeyError, np.linalg.LinAlgError) as exc: + cases.append({"input": str(Path(path).resolve()), "analysis_pass": False, "error": str(exc)}) + result = {"schema_version": 1, "passed": all(c["analysis_pass"] for c in cases), + "analyzer_sha256": digest(__file__), "numpy_version": np.__version__, + "scope": "Unweighted nodal projections of the physical radial direction. Not area-weighted harmonic power, not a continuum symmetry or solver-nullspace certificate.", + "cases": cases} + 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)} surface analyses; {output}") + return 0 if result["passed"] else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/experiments/plot_shape_surface_comparison.py b/experiments/plot_shape_surface_comparison.py new file mode 100644 index 0000000..c7128d7 --- /dev/null +++ b/experiments/plot_shape_surface_comparison.py @@ -0,0 +1,58 @@ +#!/usr/bin/env python3 +"""Plot observed nonrotating surface asphericity at saved radial nodes.""" +import argparse +import csv +import os +from pathlib import Path +for name in ('OPENBLAS_NUM_THREADS', 'OMP_NUM_THREADS', 'VECLIB_MAXIMUM_THREADS'): + os.environ[name] = '1' +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('directories', type=Path, nargs=2) + parser.add_argument('--output', type=Path, required=True) + args = parser.parse_args() + if args.output.exists(): + parser.error('Use a fresh output file.') + data = [] + for directory, expected_order in zip(args.directories, (3, 4)): + metadata = dict(line.split("=", 1) for line in (directory / "metadata.txt").read_text().splitlines() if "=" in line) + with (directory / "weak_shape_physical_metadata.csv").open() as stream: + physical_metadata = {r["key"]: r["value"] for r in csv.DictReader(stream)} + if float(metadata["n"]) != 1.5 or float(metadata["J"]) != 0 or int(physical_metadata["displacement_family_order"]) != expected_order: + 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)) + 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): + raise ValueError('Invalid physical surface samples.') + 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))) + 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')): + 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.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) + args.output.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(args.output, dpi=200) + plt.close(fig) + + +if __name__ == '__main__': + main()