tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Record audited shape-field compensation and cubic surface structure before order trial

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-18 13:45:27 -04:00
parent 2c4a5e05a6
commit 211be17190
3 changed files with 520 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

97
experiments/SHAPE_FIELD_WEAK_RESPONSE_2026-09-18.md
View File

@@ -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.

366
experiments/analyze_shape_surface_modes.py Normal file
View File

@@ -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())

58
experiments/plot_shape_surface_comparison.py Normal file
View File

@@ -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()