"""Compare the approved 101-frame HHT reference without modifying artifacts. Run with ``uv run --with h5py python -B scripts/verification/compare_hht_dynamic.py``. Only the documented displacement/rotation reference quantities are compared; this is not an independent validation of inertia, energy, or convergence. """ from __future__ import annotations import argparse import csv from decimal import Decimal, InvalidOperation import hashlib import json import math from pathlib import Path import re from typing import Any import h5py import numpy as np ROOT = Path(__file__).resolve().parents[2] CASE = "shellNL-dyn-meshfine" REFERENCE_DIRECTORY = ROOT / "reference" / CASE NODE_LABELS = (2, 21, 22, 69, 70, 165, 309, 442, 586) INSTANCE = "PART-1-1" STEP = "Step-1" NODE_COUNT = 729 FRAME_COUNT = 101 INCREMENT = 0.1 COMPONENTS = ("U-U1", "U-U2", "U-U3", "UR-UR1", "UR-UR2", "UR-UR3") CSV_COLUMNS = ("Frame", "Part Instance Name", "Node Label", *COMPONENTS) FRAME_PATTERN = re.compile( r"Increment\s+(\d+):\s*Step Time\s*=\s*" r"([+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?)" ) POLICY = { "policy_id": "common-reference-tolerance", "near_zero_ratio": 0.01, "relative_tolerance": 0.05, "relative_rms_tolerance": 0.01, } class ComparisonError(ValueError): """An identity/schema/finite failure that must precede tolerance evaluation.""" def __init__(self, code: str, message: str): super().__init__(message) self.code = code def require(condition: bool, code: str, message: str) -> None: if not condition: raise ComparisonError(code, message) def positive_label(text: str) -> int: require(bool(re.fullmatch(r"\d+", text)), "invalid-node-label", repr(text)) label = int(text) require(0 < label <= 2**63 - 1, "invalid-node-label", repr(text)) return label def canonical_instance(text: str) -> str: return text.translate(str.maketrans("abcdefghijklmnopqrstuvwxyz", "ABCDEFGHIJKLMNOPQRSTUVWXYZ")) def input_inventory(path: Path) -> tuple[dict[int, tuple[float, ...]], str]: """Read identity and schedule authority from this specific approved INP.""" raw = path.read_bytes() nodes: dict[int, tuple[float, ...]] = {} parts: list[str] = [] instances: list[tuple[str, str]] = [] steps: list[dict[str, str]] = [] dynamic_rows: list[tuple[float, float]] = [] keyword = "" in_part = False for line_number, line in enumerate(raw.decode("utf-8-sig").splitlines(), 1): text = line.strip() if not text or text.startswith("**"): continue fields = [field.strip() for field in text.split(",")] if text.startswith("*"): keyword = fields[0].upper() parameters = dict( (name.upper(), value.strip()) for field in fields[1:] if "=" in field for name, value in [field.split("=", 1)] ) if keyword == "*PART": in_part = True parts.append(canonical_instance(parameters.get("NAME", ""))) elif keyword == "*END PART": in_part = False elif keyword == "*INSTANCE": instances.append((canonical_instance(parameters.get("NAME", "")), canonical_instance(parameters.get("PART", "")))) elif keyword == "*STEP": steps.append(parameters) elif keyword == "*DYNAMIC": require([field.upper() for field in fields] == ["*DYNAMIC", "DIRECT"], "wrong-input-procedure", f"line {line_number}") continue if keyword == "*NODE": require(in_part and len(fields) == 4, "invalid-input-node", f"line {line_number}") label = positive_label(fields[0]) require(label not in nodes, "duplicate-input-node", str(label)) coordinates = tuple(float(value) for value in fields[1:]) require(all(math.isfinite(value) for value in coordinates), "nonfinite-input-node", str(label)) nodes[label] = coordinates elif keyword == "*DYNAMIC": if fields[-1] == "": fields.pop() require(len(fields) == 2, "wrong-input-schedule", f"line {line_number}") dynamic_rows.append((float(fields[0]), float(fields[1]))) require(parts == ["PART-1"] and instances == [(INSTANCE, "PART-1")], "wrong-input-instance", "Expected the original single identity instance.") require(len(nodes) == NODE_COUNT and set(NODE_LABELS) <= nodes.keys(), "wrong-input-node-inventory", f"Expected {NODE_COUNT} source nodes.") require(len(steps) == 1 and steps[0].get("NAME") == STEP and steps[0].get("NLGEOM", "").upper() == "YES" and steps[0].get("INC") == "1000", "wrong-input-step", str(steps)) require(dynamic_rows == [(INCREMENT, 10.0)], "wrong-input-schedule", str(dynamic_rows)) value = 14695981039346656037 for byte in raw: value = ((value ^ byte) * 1099511628211) & (2**64 - 1) return nodes, f"fnv1a64:{value:016x}" def read_reference(path: Path) -> dict[tuple[int, int], tuple[float, ...]]: rows: dict[tuple[int, int], tuple[float, ...]] = {} with path.open(encoding="utf-8-sig", newline="") as stream: reader = csv.reader(stream) header = [field.strip() for field in next(reader, [])] require(len(header) == len(set(header)), "duplicate-reference-column", str(header)) require(set(CSV_COLUMNS) <= set(header), "missing-reference-column", str(header)) for line_number, fields in enumerate(reader, 2): if not fields or all(not field.strip() for field in fields): continue require(len(fields) == len(header), "invalid-reference-row", f"line {line_number}") row = dict(zip(header, (field.strip() for field in fields))) match = FRAME_PATTERN.fullmatch(row["Frame"]) require(match is not None, "invalid-reference-frame", f"line {line_number}") increment = int(match.group(1)) time = Decimal(match.group(2)) require(time.is_finite() and time == Decimal(increment) / Decimal(10), "wrong-reference-time", f"line {line_number}: {row['Frame']}") require(canonical_instance(row["Part Instance Name"]) == INSTANCE, "wrong-reference-instance", f"line {line_number}") label = positive_label(row["Node Label"]) key = (increment, label) require(key not in rows, "duplicate-reference-row", str(key)) values = tuple(float(row[component]) for component in COMPONENTS) require(all(math.isfinite(value) for value in values), "nonfinite-reference-value", str(key)) rows[key] = values expected = {(frame, label) for frame in range(FRAME_COUNT) for label in NODE_LABELS} require(rows.keys() == expected, "reference-row-inventory", f"missing={sorted(expected - rows.keys())}; extra={sorted(rows.keys() - expected)}") return rows def text_value(value: Any) -> str: return value.decode("utf-8") if isinstance(value, bytes) else str(value) def string_attribute(obj: Any, name: str, expected: str | None = None) -> str: require(name in obj.attrs, "missing-hdf5-attribute", f"{obj.name}@{name}") attribute = obj.attrs.get_id(name) string = h5py.check_string_dtype(attribute.dtype) require(attribute.shape == () and string is not None and string.length is None and string.encoding == "utf-8", "wrong-hdf5-attribute-type", f"{obj.name}@{name}") value = text_value(obj.attrs[name]) require(bool(value) and (expected is None or value == expected), "wrong-hdf5-attribute-value", f"{obj.name}@{name}: {value!r}") return value def numeric_attribute(obj: Any, name: str, dtype: str, expected: int | float) -> None: require(name in obj.attrs, "missing-hdf5-attribute", f"{obj.name}@{name}") attribute = obj.attrs.get_id(name) require(attribute.shape == () and attribute.dtype == np.dtype(dtype), "wrong-hdf5-attribute-type", f"{obj.name}@{name}") value = obj.attrs[name] require(math.isfinite(value) and value == expected, "wrong-hdf5-attribute-value", f"{obj.name}@{name}: {value!r}") def dataset(file: h5py.File, path: str, shape: tuple[int, ...]) -> h5py.Dataset: require(path in file and isinstance(file[path], h5py.Dataset), "missing-hdf5-dataset", path) result = file[path] require(result.shape == shape and result.chunks is None and result.compression is None, "wrong-hdf5-dataset-layout", path) return result def read_results(path: Path, nodes: dict[int, tuple[float, ...]], content_identity: str) -> tuple[list[np.ndarray], dict[str, Any]]: """Precheck every result node/frame/component, then project by source key.""" projected: list[np.ndarray] = [] with h5py.File(path, "r") as file: require("metadata" in file and isinstance(file["metadata"], h5py.Group), "missing-hdf5-metadata", "/metadata") metadata = file["metadata"] for name, expected in { "feature_id": "hht-implicit-dynamics", "analysis_procedure": "hht-implicit-dynamics", "mass_formulation": "nodal-lumped-kinetic", "internal_formulation": "FESA-MITC4-total-lagrangian", "initial_state": "zero-q-v-step-load-0+", "acceleration_convention": "physical-endpoint-and-auxiliary-separated", "source_step_name": STEP, "step_name": STEP, "unit_system_label": "consistent-input-units", }.items(): string_attribute(metadata, name, expected) string_attribute(metadata, "solver_version") source_identity = string_attribute(metadata, "source_input_identity") require(source_identity.endswith(";content_identity=" + content_identity), "wrong-input-content-identity", source_identity) require("schema_version" in metadata.attrs, "missing-hdf5-attribute", "/metadata@schema_version") version_attribute = metadata.attrs.get_id("schema_version") require(version_attribute.shape == () and version_attribute.dtype == np.dtype(" dict[str, Any]: """Apply the common policy exactly, including zero-scale and family RMS.""" require(len(result) == len(reference) and len(result) > 0, "scalar-inventory", quantity) require(all(math.isfinite(value) for value in (*result, *reference)), "nonfinite-comparison-value", quantity) scale = max(abs(value) for value in reference) band = POLICY["near_zero_ratio"] * scale rows = [] error_norm = 0.0 for index, (actual, expected) in enumerate(zip(result, reference)): error = abs(actual - expected) require(math.isfinite(error), "nonfinite-comparison-value", quantity) error_norm = math.hypot(error_norm, error) require(math.isfinite(error_norm), "nonfinite-comparison-value", quantity) relative_error = None if scale == 0.0: branch, threshold, passed = "zero-scale-exact", 0.0, error == 0.0 elif abs(expected) <= band: branch, threshold, passed = "near-zero", band, error <= band else: branch = "relative" threshold = POLICY["relative_tolerance"] * abs(expected) relative_error = error / abs(expected) require(math.isfinite(relative_error), "nonfinite-comparison-value", quantity) passed = relative_error <= POLICY["relative_tolerance"] rows.append({"instance": INSTANCE, "node_label": NODE_LABELS[index // 3], "component": f"{quantity}{index % 3 + 1}", "fesa_value": actual, "reference_value": expected, "absolute_error": error, "relative_error": relative_error, "branch": branch, "threshold": threshold, "passed": passed}) rows_passed = all(row["passed"] for row in rows) rms = error_norm / math.sqrt(len(rows)) / scale if scale else None require(rms is None or math.isfinite(rms), "nonfinite-comparison-value", quantity) rms_passed = rows_passed if scale == 0.0 else rms <= POLICY["relative_rms_tolerance"] return {"identity": f"{quantity}-global:{STEP}:increment-{increment}", "case": CASE, "step": STEP, "increment": increment, "hdf5_time": float(increment) * INCREMENT, "reference_decimal_time": str(Decimal(increment) / Decimal(10)), "quantity": quantity, "components": [f"{quantity}{i}" for i in (1, 2, 3)], "unit_dimension": "length" if quantity == "U" else "rad", "coordinate_system": "global-cartesian" if quantity == "U" else "fixed-initial-global-total-rotation-vector", "impact": "blocking" if quantity == "U" else "warning-only", "row_count": len(rows), "reference_scale": scale, "near_zero_band": band, "rows": rows, "rows_passed": rows_passed, "relative_rms": rms, "rms_passed": rms_passed, "diagnostic_code": "zero-reference-scale-nonzero-error" if scale == 0.0 and not rows_passed else None, "passed": rows_passed and rms_passed} def compare(input_path: Path, reference_path: Path, results_path: Path) -> dict[str, Any]: report: dict[str, Any] = { "schema": "hht-reference-comparison-v1", "case": CASE, "policy": dict(POLICY), "identities": {}, "precheck": {"passed": False}, "families": [], "blocking_failures": [], "warnings": [], "blocking_passed": False, "verdict": "failed", } try: for name, path in (("input", input_path), ("reference", reference_path), ("results", results_path)): report["identities"][name] = {"path": str(path.absolute())} with path.open("rb") as stream: report["identities"][name]["sha256"] = hashlib.file_digest(stream, "sha256").hexdigest() nodes, content_identity = input_inventory(input_path) reference = read_reference(reference_path) result, evidence = read_results(results_path, nodes, content_identity) families = [] for frame in range(FRAME_COUNT): for quantity, first in (("U", 0), ("UR", 3)): actual = [float(value) for row in result[frame] for value in row[first:first + 3]] expected = [value for label in NODE_LABELS for value in reference[(frame, label)][first:first + 3]] families.append(evaluate_family(frame, quantity, actual, expected)) report["precheck"] = {"passed": True, "reference_row_count": len(reference), "reference_frame_count": FRAME_COUNT, **evidence} report["families"] = families for family in families: if family["passed"]: continue if family["impact"] == "blocking": report["blocking_failures"].append(family["identity"]) else: if not family["rows_passed"]: report["warnings"].append({"family": family["identity"], "code": "row-tolerance-failure"}) if not family["rms_passed"]: report["warnings"].append({"family": family["identity"], "code": "relative-rms-failure"}) report["blocking_passed"] = not report["blocking_failures"] report["verdict"] = "passed" if report["blocking_passed"] else "failed" except (ComparisonError, OSError, ValueError, InvalidOperation, UnicodeError) as error: report["precheck"] = {"passed": False, "code": getattr(error, "code", "invalid-comparison-input"), "message": str(error)} return report def main() -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--input", type=Path, default=REFERENCE_DIRECTORY / f"{CASE}.inp") parser.add_argument("--reference", type=Path, default=REFERENCE_DIRECTORY / f"{CASE} displacements.csv") parser.add_argument("--results", type=Path) parser.add_argument("--report", type=Path, required=True) args = parser.parse_args() results = args.results or args.input.with_name(args.input.stem + "_results.h5") if args.report.resolve() in {args.input.resolve(), args.reference.resolve(), results.resolve()}: parser.error("The report path must not overwrite an input, reference, or solver result.") report = compare(args.input, args.reference, results) args.report.parent.mkdir(parents=True, exist_ok=True) args.report.write_text(json.dumps(report, indent=2, sort_keys=True, allow_nan=False) + "\n", encoding="utf-8") print(json.dumps({"verdict": report["verdict"], "precheck": report["precheck"], "blocking_failure_count": len(report["blocking_failures"]), "warning_count": len(report["warnings"]), "report": str(args.report)}, allow_nan=False)) return 2 if not report["precheck"]["passed"] else (0 if report["blocking_passed"] else 1) if __name__ == "__main__": raise SystemExit(main())