Files
MeanField_SERiF/tools/validate_ester.py

336 lines
23 KiB
Python

#!/usr/bin/env python3
"""Run fresh SERiF solves and compare the unchanged physical requests with the archived ESTER study.
Example:
python tools/validate_ester.py --report ../MeanField/reports/ester_validation_2026-09-21 \
--executable build/tests/solver/ester_validation --output /tmp/serif-ester
Omit --executable to analyze already generated SERiF output. No archived solver
field is used as a new measurement. Archive coefficients are warm initial guesses;
solve.json must demonstrate independent SERiF nonlinear convergence.
"""
import argparse
import csv
import hashlib
import importlib.util
import json
import math
from pathlib import Path
import subprocess
import sys
import numpy as np
CASES = ("n1_J0", "n1_J0p05", "n1_J0p1", "n1p5_J0", "n1p625_J0")
SCALARS = {"rho": "rho", "h": "h", "P": "pressure_h", "Phi": "phi", "bernoulli_defect": "bernoulli_defect"}
def read_rows(path):
"""Retain archived point order and names, including failed point locations."""
with path.open() as source:
return list(csv.DictReader(source))
def column(rows, name):
"""Read one named numerical measurement from every requested point."""
return np.asarray([float(row[name]) for row in rows])
def magnitude(values):
"""Use absolute scalar values or Euclidean vector lengths for error norms."""
return np.linalg.norm(values, axis=1) if values.ndim == 2 else np.abs(values)
def norms(values, mask, scale, weights=None):
"""Equal-weight requested-grid norms, as specified by the original protocol."""
selected = magnitude(values[mask]) / scale
if not len(selected) or not np.isfinite(selected).all():
raise ValueError("A required physical comparison contains missing or nonfinite samples")
selected_weights = np.ones(len(selected)) if weights is None else weights[mask]
if not np.isfinite(selected_weights).all() or np.any(selected_weights <= 0):
raise ValueError("A required physical comparison has invalid quadrature weights")
return {"count": len(selected), "rms_scaled": float(np.sqrt(np.average(selected**2, weights=selected_weights))),
"max_scaled": float(np.max(selected))}
def digest(path):
"""Record the exact bytes of a source or independently generated measurement."""
return hashlib.sha256(path.read_bytes()).hexdigest()
def compare_surface(report, case, output, models):
"""Repeat the original actual-surface quadrature and sampled-axis measurements."""
old = read_rows(report / "meanfield_runs" / case / "surface.csv")
new = read_rows(output / "surface.csv")
identity = ("kind", "boundary_element", "sample", "reference_u", "reference_v")
if [[row[key] for key in identity[:3]] for row in old] != [[row[key] for key in identity[:3]] for row in new]:
raise ValueError(f"{case}: actual surface sample identities changed")
for key in identity[3:]:
if not np.array_equal(column(old, key), column(new, key)):
raise ValueError(f"{case}: actual surface reference requests changed")
base = models["base"]
scale = base.globals["R_eq"]
kind = np.asarray([row["kind"] for row in old])
old_radius, radius = column(old, "radius"), column(new, "radius")
xyz = np.stack([column(new, c) for c in ("x", "y", "z")], axis=1)
if not np.isfinite(xyz).all() or not np.isfinite(radius).all() or np.any(radius <= 0):
raise ValueError(f"{case}: a required actual surface sample is nonfinite or nonpositive")
if not np.allclose(np.linalg.norm(xyz, axis=1), radius, rtol=32*np.finfo(float).eps, atol=1e-15):
raise ValueError(f"{case}: actual surface positions and radii disagree")
reference = {level: model.surface_radius(column(new, "theta")) for level, model in models.items()}
old_reference = base.surface_radius(column(old, "theta"))
archived = json.loads((report / "analysis" / case / "metrics.json").read_text())["surface"]
quadrature = {}
for rule in ("q14", "q18"):
mask = kind == rule
fresh = norms(radius - reference["base"], mask, scale, column(new, "area_weight"))
reconstructed = norms(old_radius - old_reference, mask, scale, column(old, "area_weight"))
recorded = {key: archived[rule][key] for key in ("rms_scaled", "max_scaled")}
if fresh["count"] != archived[rule]["eligible"] or any(abs(reconstructed[key] - recorded[key]) > 128*np.finfo(float).eps for key in recorded):
raise ValueError(f"{case}: archived {rule} surface metric did not reproduce")
sensitivity = {level: norms(reference[level] - reference["base"], mask, scale, column(new, "area_weight")) for level in ("low", "high")}
quadrature[rule] = {"fresh": fresh, "archived": recorded, "archived_reconstructed": reconstructed,
"fresh_minus_archived": norms(radius - old_radius, mask, scale, column(old, "area_weight")),
"area": float(column(new, "area_weight")[mask].sum()), "archived_area": archived[rule]["area"],
"ester_grid_change": sensitivity,
"rms_change": fresh["rms_scaled"] - recorded["rms_scaled"], "max_change": fresh["max_scaled"] - recorded["max_scaled"],
"within_archived_plus_high_grid_change": all(fresh[key] <= recorded[key] + sensitivity["high"][key] for key in recorded)}
# These are selected actual nodes, never extrema or ideal harmonic radii.
axis_tolerance = 2e-12 * scale
old_xyz = np.stack([column(old, c) for c in ("x", "y", "z")], axis=1)
axes = {}
# Freeze the original node visits and exact-coordinate deduplication. Repeating that
# deduplication on fresh floating-point positions can change weights through roundoff alone.
for name, select in {"equator": np.abs(old_xyz[:, 2]) <= axis_tolerance,
"north_pole": (np.hypot(old_xyz[:, 0], old_xyz[:, 1]) <= axis_tolerance) & (old_xyz[:, 2] > 0),
"south_pole": (np.hypot(old_xyz[:, 0], old_xyz[:, 1]) <= axis_tolerance) & (old_xyz[:, 2] < 0)}.items():
indices = np.flatnonzero((kind == "node") & select)
if not len(indices):
raise ValueError(f"{case}: missing required {name} actual surface nodes")
_, unique = np.unique(old_xyz[indices], axis=0, return_index=True)
indices = indices[unique]
transverse = np.abs(xyz[indices, 2]) if name == "equator" else np.hypot(xyz[indices, 0], xyz[indices, 1])
if np.any(transverse > axis_tolerance):
raise ValueError(f"{case}: a required {name} surface request moved off its physical axis")
values = radius[indices]
axes[name] = {"mean": float(values.mean()), "minimum": float(values.min()), "maximum": float(values.max()),
"archived_node_visits": int(((kind == "node") & select).sum()), "fixed_archived_unique_coordinate_samples": len(values),
"coordinate_tolerance": axis_tolerance}
return {"quadrature": quadrature, "axis_samples": axes}
def compare_globals(report, case, output, old, new, models, surface):
"""Compare measured center traces, typed global coefficients and actual surface nodes."""
center = [index for index, row in enumerate(old) if row["id"] == "center:00000"]
if len(center) != 1 or any(float(new[center[0]][c]) != 0 for c in ("x", "y", "z")):
raise ValueError(f"{case}: the explicit center trace is missing")
index = center[0]
summary = json.loads((output / "physical_summary.json").read_text())
measured = {name: summary[name] for name in ("Omega", "C", "surface_volume_radius_q14", "surface_volume_radius_q18")}
measured.update(rho_c=float(new[index]["rho"]), h_c=float(new[index]["h"]))
bulk_rows = read_rows(output / "bulk.csv")
bulk = {(row["quadrature"], row["metric"]): float(row["value"]) for row in bulk_rows}
if len(bulk) != len(bulk_rows):
raise ValueError(f"{case}: duplicate physical bulk measurements")
for quadrature in ("14", "18"):
for quantity, metric in {"M": "mass", "I": "moment_of_inertia_z", "J": "angular_momentum_z", "R_volume": "volume_radius"}.items():
measured[quantity + "_q" + quadrature] = bulk[(quadrature, metric)]
for axis, values in surface["axis_samples"].items():
quantity = ("R_eq" if axis == "equator" else "R_p") + "_sampled_" + axis + "_mean"
measured[quantity] = values["mean"]
records = {(row["quantity"], row["quadrature"]): row for row in read_rows(report / "analysis" / "global_comparison.csv") if row["case"] == case}
result = {}
for name, value in measured.items():
if name.startswith("surface_volume_radius_q"):
recorded = records[("R_volume", name.rsplit("q", 1)[1])]
quantity = "R_volume"
elif name.endswith(("_q14", "_q18")):
quantity, quadrature = name.rsplit("_q", 1)
recorded = records[(quantity, quadrature)]
else:
recorded = records[(name, "14" if name in ("Omega", "C") else "")]
quantity = name.split("_sampled_")[0]
scale = float(recorded["normalization"])
if not math.isfinite(value):
raise ValueError(f"{case}: {name} is nonfinite")
reference = {level: model.globals[quantity] for level, model in models.items()}
previous = float(recorded["meanfield"])
if any(abs(reference[level] - float(recorded["ester_" + level])) > 128*np.finfo(float).eps*max(1, abs(reference[level])) for level in reference):
raise ValueError(f"{case}: archived {name} ESTER global reference no longer reproduces")
fresh_error = (value - reference["base"]) / scale
old_error = (previous - reference["base"]) / scale
sensitivity = {level: (reference[level] - reference["base"]) / scale for level in ("low", "high")}
result[name] = {"fresh": value, "archived": previous, "ester": reference, "normalization": scale,
"fresh_scaled_difference": fresh_error, "archived_scaled_difference": old_error,
"fresh_minus_archived_scaled": (value - previous) / scale, "ester_grid_change_scaled": sensitivity,
"within_archived_plus_high_grid_change": bool(abs(fresh_error) <= abs(old_error) + abs(sensitivity["high"]))}
return {"quantities": result, "center_boundary_trace_possible": bool(float(new[index]["boundary_trace_possible"])),
"surface_volume_measurement": "Fresh divergence-theorem boundary quadrature compared with original volume-integral radii.",
"physical_summary": summary}
def compare(report, case, output, reference_class):
"""Keep the original common-material masks fixed so missing samples cannot improve a score."""
solve = json.loads((output / "solve.json").read_text())
if solve["case"] != case or solve["converged"] is not True or not math.isfinite(solve["final_residual_norm"]) or not 0.0 <= solve["final_residual_norm"] <= 1e-10:
raise ValueError(f"{case}: fresh SERiF nonlinear qualification failed")
metadata = dict(line.split("=", 1) for line in (report / "inputs" / case / "metadata.txt").read_text().splitlines() if "=" in line)
if digest(report / "inputs" / case / "mesh.smesh") != metadata["mesh_sha256"]:
raise ValueError(f"{case}: the reference mesh no longer matches its recorded hash")
catalog = json.loads((report / "reference_catalog.json").read_text())[case]
old = read_rows(report / "meanfield_runs" / case / "samples.csv")
new = read_rows(output / "samples.csv")
if [row["id"] for row in old] != [row["id"] for row in new]:
raise ValueError(f"{case}: point identity or order changed")
xyz = np.stack([column(new, c) for c in ("x", "y", "z")], axis=1)
if not np.array_equal(xyz, np.stack([column(old, c) for c in ("x", "y", "z")], axis=1)):
raise ValueError(f"{case}: requested physical coordinates changed")
old_located = column(old, "located") == 1
new_located = column(new, "located") == 1
if np.any(old_located & ~new_located):
raise ValueError(f"{case}: {np.count_nonzero(old_located & ~new_located)} previously located samples were lost")
def physical_coordinates(rows, located):
"""Validate actual located positions against their unchanged physical requests."""
mapped = np.stack([column(rows, "mapped_" + c) for c in ("x", "y", "z")], axis=1)
error = np.linalg.norm(mapped - xyz, axis=1)
tolerance = 2e-11 * np.maximum(1.0, np.linalg.norm(xyz, axis=1))
roundoff = 64 * np.finfo(float).eps * np.maximum(1.0, np.linalg.norm(xyz, axis=1))
if not np.isfinite(mapped[located]).all() or np.any(error[located] > tolerance[located] + roundoff[located]):
raise ValueError(f"{case}: sampled physical location exceeds the original locator tolerance")
if np.any(np.abs(error[located] - column(rows, "location_error")[located]) > roundoff[located]):
raise ValueError(f"{case}: sampled location error is inconsistent with mapped coordinates")
coordinates = xyz.copy()
coordinates[located] = mapped[located]
return coordinates
old_xyz = physical_coordinates(old, old_located)
new_xyz = physical_coordinates(new, new_located)
models = {level: reference_class(report / "ester_runs" / Path(catalog[level]["directory"]).name,
n=float(metadata["n"]), K=float(metadata["K"]), M=float(metadata["M"]), G=float(metadata["G"]))
for level in ("base", "low", "high")}
if not all(model.globals["reference_converged"] for model in models.values()):
raise ValueError(f"{case}: an ESTER reference did not meet its recorded convergence criterion")
match_tolerance = json.loads((report / "matching_plan.json").read_text())["angular_momentum_absolute_tolerance"]
if any(abs(model.globals["J"] - float(metadata["J"])) > match_tolerance + 16 * np.finfo(float).eps for model in models.values()):
raise ValueError(f"{case}: an ESTER reference does not match the physical angular momentum")
reference = {level: model.evaluate(new_xyz) for level, model in models.items()}
old_reference = models["base"].evaluate(old_xyz)
base = models["base"]
radii = np.linalg.norm(old_xyz, axis=1)
theta = np.arccos(np.clip(np.divide(old_xyz[:, 2], radii, out=np.ones(len(radii)), where=radii > 0), -1, 1))
common = old_located & (column(old, "stellar_material") == 1) & old_reference["inside_star"]
strict = common & (column(old, "boundary_trace_possible") == 0) & (radii / base.surface_radius(theta) < 1 - 1e-8)
if np.any(strict & (column(new, "stellar_material") != 1)):
raise ValueError(f"{case}: a required interior physical request changed material support")
if np.any(strict & (column(new, "boundary_trace_possible") != 0)):
raise ValueError(f"{case}: a required interior sample moved onto an ambiguous FE trace")
scales = {"rho": base.rho_c, "h": base.H, "P": base.rho_c * base.H / (base.n + 1),
"Phi": base.G * base.M / base.globals["R_eq"], "gradPhi": base.G * base.M / base.globals["R_eq"]**2,
"gradPhi_broken": base.G * base.M / base.globals["R_eq"]**2, "bernoulli_defect": base.H}
recorded_metrics = {(row["category"], row["region"], row["field"]): row
for row in read_rows(report / "analysis" / "metrics.csv") if row["case"] == case and row["group"] == "all"}
def field_values(rows, name):
"""Read the specified physical scalar or one of the two gravity gradients."""
if name in ("gradPhi", "gradPhi_broken"):
prefix = "mixed_dphi_d" if name == "gradPhi" else "broken_phi_d"
return np.stack([column(rows, prefix + c) for c in ("x", "y", "z")], axis=1)
return column(rows, SCALARS[name])
def field_comparison(name, mask, category, region="all"):
"""Reconstruct the exact archived metric before comparing a new independently sampled field."""
archived_name = "gradPhi_mixed" if name == "gradPhi" else name
reference_name = "gradPhi" if name == "gradPhi_broken" else name
recorded = recorded_metrics[(category, region, archived_name)]
if int(recorded["eligible"]) != int(mask.sum()):
raise ValueError(f"{case}: {category}/{region} support no longer reproduces the original study")
old_field, new_field = field_values(old, name), field_values(new, name)
fresh = norms(new_field - reference["base"][reference_name], mask, scales[name])
preserved = norms(new_field - old_field, mask, scales[name])
sensitivity = {level: norms(reference[level][reference_name] - reference["base"][reference_name], mask, scales[name]) for level in ("low", "high")}
archived = {key: float(recorded[key]) for key in ("rms_scaled", "max_scaled")}
reconstructed = norms(old_field - old_reference[reference_name], mask, scales[name])
if any(abs(reconstructed[key] - archived[key]) > 128 * np.finfo(float).eps * max(1.0, abs(archived[key]))
for key in ("rms_scaled", "max_scaled")):
raise ValueError(f"{case}: archived {category}/{region}/{name} metrics did not reproduce at their original mapped coordinates")
return {"fresh": fresh, "archived": archived, "archived_reconstructed": reconstructed,
"reference_relocation": norms(reference["base"][reference_name] - old_reference[reference_name], mask, scales[name]),
"fresh_minus_archived": preserved, "ester_grid_change": sensitivity,
"rms_change": fresh["rms_scaled"] - archived["rms_scaled"],
"max_change": fresh["max_scaled"] - archived["max_scaled"],
"within_archived_plus_high_grid_change": all(fresh[key] <= archived[key] + sensitivity["high"][key]
for key in ("rms_scaled", "max_scaled"))}
fields = {name: field_comparison(name, strict, "common_material_strict") for name in (*SCALARS, "gradPhi", "gradPhi_broken")}
# Preserve the original radial-region definitions, including gravity in physical vacuum.
fraction = radii / base.surface_radius(theta)
regions = {"all": old_located, "core": old_located & (fraction < .5),
"envelope": old_located & (fraction >= .5) & (fraction < .9),
"surface_layer": old_located & (fraction >= .9) & (fraction <= 1),
"outside_reference": old_located & (fraction > 1)}
gravity = {region: {name: field_comparison(name, mask, "located_all_materials", region) for name in ("Phi", "gradPhi", "gradPhi_broken")}
for region, mask in regions.items() if mask.any()}
def material_diagnostics(rows, located):
"""Expose negative material tails and undefined pressure without clipping them."""
material = located & (column(rows, "stellar_material") == 1)
density, enthalpy, pressure = (column(rows, name) for name in ("rho", "h", "pressure_h"))
if not np.isfinite(density[material]).all() or not np.isfinite(enthalpy[material]).all():
raise ValueError(f"{case}: material density or enthalpy contains a nonfinite value")
return {"material_points": int(material.sum()), "negative_density_points": int((material & (density < 0)).sum()),
"negative_enthalpy_points": int((material & (enthalpy < 0)).sum()),
"nonfinite_pressure_points": int((material & ~np.isfinite(pressure)).sum()),
"minimum_density": float(density[material].min()), "minimum_enthalpy": float(enthalpy[material].min())}
tails = {"archived": material_diagnostics(old, old_located), "fresh": material_diagnostics(new, new_located)}
surface = compare_surface(report, case, output, models)
globals_comparison = compare_globals(report, case, output, old, new, models, surface)
provenance = {str(path): digest(path) for path in [report / "inputs" / case / "accepted_state.txt",
report / "inputs" / case / "mesh.smesh", output / "accepted_state.txt", output / "samples.csv", report / "ester_reference.py",
report / "inputs" / case / "layout.csv", report / "inputs" / case / "metadata.txt",
output / "physical_summary.json", output / "surface.csv", output / "bulk.csv"]}
result = {"case": case, "fresh_solve": solve, "strict_common_material_points": int(strict.sum()),
"located_points": int(new_located.sum()), "fields": fields, "located_gravity_by_region": gravity, "surface": surface, "globals": globals_comparison, "material_tail_diagnostics": tails, "sha256": provenance,
"ester_reference_sha256": {level: model.hashes for level, model in models.items()},
"qualification": "Fresh nonlinear atol=1e-10; fixed archived physical requests, strict material masks and all-located gravity masks; no universal physical percentage tolerance exists in the source protocol.",
"envelope": "Archived error plus the separately measured high-minus-base ESTER grid change is reported, not substituted for the recorded errors."}
qualified = list(fields.values()) + [item for region in gravity.values() for item in region.values()]
qualified += list(surface["quadrature"].values()) + list(globals_comparison["quantities"].values())
result["within_archived_plus_high_grid_change"] = all(item["within_archived_plus_high_grid_change"] for item in qualified)
(output / "comparison.json").write_text(json.dumps(result, indent=2, allow_nan=False) + "\n")
return result
def main():
"""Optionally solve selected stars, then write their independently checked comparisons."""
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--report", type=Path, required=True)
parser.add_argument("--output", type=Path, required=True)
parser.add_argument("--executable", type=Path)
parser.add_argument("--cases", nargs="+", choices=CASES, default=CASES)
args = parser.parse_args()
args.output.mkdir(parents=True, exist_ok=True)
sys.dont_write_bytecode = True
specification = importlib.util.spec_from_file_location("serif_ester_reference", args.report / "ester_reference.py")
reference = importlib.util.module_from_spec(specification)
specification.loader.exec_module(reference)
results = []
for case in args.cases:
target = args.output / case
target.mkdir(parents=True, exist_ok=True)
if args.executable:
with (target / "solve.log").open("w") as log:
subprocess.run([str(args.executable.resolve()), str(args.report.resolve()), case, str(target.resolve())],
stdout=log, stderr=subprocess.STDOUT, check=True)
result = compare(args.report, case, target, reference.Model)
results.append(result)
print(case, "converged", result["fresh_solve"]["final_residual_norm"], "Newton steps", result["fresh_solve"]["accepted_iterations"], flush=True)
for field, data in result["fields"].items():
print(f" {field}: RMS {data['fresh']['rms_scaled']:.9g}; archived {data['archived']['rms_scaled']:.9g}; max physical drift {data['fresh_minus_archived']['max_scaled']:.3g}", flush=True)
for rule, data in result["surface"]["quadrature"].items():
print(f" surface {rule}: RMS {data['fresh']['rms_scaled']:.9g}; archived {data['archived']['rms_scaled']:.9g}; max radius drift {data['fresh_minus_archived']['max_scaled']:.3g}", flush=True)
print(" all reported errors within archived plus measured ESTER high-grid change:",
result["within_archived_plus_high_grid_change"], flush=True)
(args.output / "comparison.json").write_text(json.dumps(results, indent=2, allow_nan=False) + "\n")
if __name__ == "__main__":
main()