diff --git a/docs/hex8-structural-analysis/reference-model.md b/docs/hex8-structural-analysis/reference-model.md index d0ca90e..71cb15b 100644 --- a/docs/hex8-structural-analysis/reference-model.md +++ b/docs/hex8-structural-analysis/reference-model.md @@ -32,8 +32,9 @@ Frequency는 original unstressed 상태의10 modes 요청이며 결과 Value는 ## Quantity와 identity -Static U1/U2/U3 global, RF1/RF2/RF3의 constrained projection, integration-point Cauchy -S11/S22/S33/S12/S13/S23 각각 blocking family다. U/RF는442-row 전체 inventory를 먼저 확인하고, +Static은 U1/U2/U3의 U-global, RF1/RF2/RF3 constrained projection의 RF-global, +integration-point Cauchy S11/S22/S33/S12/S13/S23의 S-global을 각각 독립 blocking family로 비교한다. +U/RF는442-row 전체 inventory를 먼저 확인하고, RF만 실제 constrained components로 투영한다. Free Abaqus RF와 FESA free residual을 비교하지 않는다. Stress는300×8×6 전부 포함하며 elemental-nodal stress·절점 외삽·평균으로 변환하지 않는다. Dynamic U는101 frame 각각442×3인 독립 blocking family다. Frequency Value는 mode별 singleton이며 diff --git a/docs/hex8-structural-analysis/verification.md b/docs/hex8-structural-analysis/verification.md new file mode 100644 index 0000000..366bc4e --- /dev/null +++ b/docs/hex8-structural-analysis/verification.md @@ -0,0 +1,137 @@ +# Hex8 original comparison evidence + +This page records T8's comparison-only Python implementation and its observed +results. It does not replace the feature PRD, numerical controls or completion +checklist. The comparator reads official rounded Float64 quantities; it does not +recompute the internal high+low solver acceptance gates from those observables. + +## Entry point and identities + +`uv run --with h5py --with pytest python scripts/verification/compare_hex8.py` +requires `--case`, the actual scratch `--input`, official `--results`, and +`--report`. `--references` defaults to `reference/hexa`. The optional +`--execution-record` consumes a closed native CLI record; it never executes its +command. A supplied record must have native exit0, unchanged original/copy byte +identities, and the exact input and result hashes. Without that record the report +explicitly certifies file comparison only, not successful solver execution. + +For example, after a real original run has naturally completed: + +```powershell +uv run --with h5py --with pytest python scripts/verification/compare_hex8.py ` + --case hexa_linearstatic --input /hexa_linearstatic.inp ` + --results /hexa_linearstatic_results.h5 ` + --execution-record /native.json --report /comparison.json +``` + +The report path may not overwrite an input, reference, execution record or result. +Authoritative repository originals remain protected even with alternate copied +references; existing symlink/junction and hardlink aliases are also checked. +Exit0 means the declared comparison passes; exit1 means a numerical family fails; +exit2 means an identity/schema/finite/precondition error or missing modal reference. +The tool is not a solver launcher, result generator, reference repair or release +readiness check. + +The fixed source inventory is the [reference-model](reference-model.md) inventory: +442 nodes, 300 C3D8 elements and 13 frozen original files. Input and CSV raw SHA256 +must match the documented Windows working-tree bytes. The HDF5 metadata must +identify the actual lexically normalized absolute scratch path and that file's FNV1a64 content +identity; a different path with identical bytes is not silently accepted. + +Complete CSV identities are checked before any physical values are converted. +Complete required Solid model/frame/physical schemas, UTF-8 attributes and explicit +LE F64/U64/U8 widths are then checked. Result arrays join the source node or +source element/point labels through their internal ID maps; CSV row order is not +an identity. Physical values, including nonprojected RF rows and unused exported +CSV derived quantities, must be finite before tolerance evaluation. + +Static reference increment10 maps to the completed load-factor1 frame, independent +of its HDF5 storage frame index. RF validates all442 source rows before projecting +the actual51 constrained scalar components. It never compares free RF with free +residual. Stress compares all300×8×6 unaveraged global Cauchy components. Fixed +IntPt1–8 conventions are checked without nearest-coordinate remapping. Point +coordinates are auxiliary: CSV coordinate errors use the approved +`16*float32-eps*sum(abs(Na*Xa_axis))` source-resolution budget; the scalar stress +policy is unchanged. + +Dynamic comparisons require exactly101 accepted frames. Integer increment0–100 +is matched first, source decimal time must equal n/10, and each stored time must +be binary64 `n*0.1`. Each frame's1326 U scalars form a separate blocking family. +There is no interpolation, nearest-time matching, scale sharing between frames, +rotation family or acceleration/energy reference claim. + +## Common numerical policy + +Every static case has independent U-global, constrained RF-global and S-global +blocking families. Every dynamic frame has its own U-global blocking family. +For each family, only reference values define S=max(abs(reference)) and B=.01S. +Rows with abs(reference)<=B require absolute error<=B; other rows require relative +error<=.05. All-row RMS/S must also be<=.01. Boundaries are inclusive, and S=0 +requires every actual value to be exactly0. No absolute floor, denominator clamp, +row deletion or cancellation of signed errors is applied. Reports preserve each +scalar identity/value/error/branch/pass and family scale/band/max-absolute-error, +absolute/relative RMS, source SHA and policy. + +The independent typed modal comparison behavior requires explicit source step, +mode1–10, eigenvalue or frequency quantity, and matching `1/time^2` or `1/time` +units. Each mode is a singleton using the same RMS1% policy. This is not an +assumption about the absent original Value file's eventual CSV layout. Actual +modal original comparison remains **NOT RUN** until a genuine reference is +provided and its semantics and hashes are confirmed. + +## Current observations + +The initial seven behavior assertions failed before implementation (pytest native +exit1). The first expanded suite passed63/63 with native exit0. Those tests include +original-sized independently declared HDF5 schema fixtures and reference values; +they are comparator behavior fixtures, not actual FESA source validation evidence. +The affected Hex8 plus legacy HHT comparator suite passed119/119 (native exit0, +62.51s) before independent review fixes. Those bounded fixes were reproduced by +safe behavior RED tests and the five targeted regressions then passed5/5. Final affected checks after all three independent findings were corrected passed +124/124 with native exit0 (81.41s), with raw stdout and JUnit retained under +`.harness/hex8-evidence/t8/final-affected/`. The preliminary whole `uv ... pytest -q` attempt and the separate exact +configured `uv run --with h5py --with pytest python -B -m pytest -v -rs` +(`docs/operations.md:87`, adding h5py) both stopped during collection with native +exit2 and five legacy operational import errors (`msvc_harness` and removed `scripts.hooks`). It is not a whole-suite PASS; +no tests were excluded and no retired orchestration was restored. + +Actual official linear-static and nonlinear-static results passed source, +model, schema, point, finite and frame prechecks and failed all three numerical +families. Their closed CLI evidence records both have native exit0 and matching +original/copy/result hashes; their comparator commands have native exit1. +These are genuine comparison failures under the unchanged policy, not reader +precheck failures or solver native execution failures. Frequency official readback +passes its prechecks; its comparator exits2 with `not-run`, not PASS. + +| Case / blocking quantity | Scalars | Failing rows | S | Maximum absolute error | RMS/S | +| --- | ---: | ---: | ---: | ---: | ---: | +| linear static U | 1326 | 79 | .00186131289 | .00007588175038730229 | .016739519723527385 | +| linear static RF | 51 | 34 | 513497.094 | 78630.61640144575 | .06818892528188494 | +| linear static S | 14400 | 7783 | 3430145 | 982855.3871926877 | .04270903417463437 | +| nonlinear static U | 1326 | 79 | .00186110672 | .00007586510202476195 | .016737790679783753 | +| nonlinear static RF | 51 | 34 | 513445.594 | 78623.61011241788 | .0681899286018377 | +| nonlinear static S | 14400 | 7782 | 3429817.75 | 982753.151651503 | .04270850513087158 | + +Official result SHA256 values are +`3ef92c35ea83ddb4d11437ff87ef15eca7dac179d378eda19570abae579b09d0` (linear), +`0095ebf2d360a8bd8a71212af0dd500ed611a56540a173a39d9c8d88d8e05bbf` (nonlinear), +and `0ec4d2b09250bab17a8532228cd34a312518c04023076f971e929b387031c6b3` (frequency). +Scratch input paths are +`C:/git/FESADev/.harness/hex8-evidence/t8-original-cli//.inp`. +Each exact command, closed native exit and full row report is retained under +`C:/git/FESADev/.dryforge/worktrees/hex8-t8-reference/.harness/hex8-evidence/t8/official-comparisons//`. +The command records identify the actual scratch input, original reference directory, +HDF5 and closed native CLI record consumed. Comparator implementation review and +remaining naturally completed dynamic originals are tracked separately. + +| Actual original | Current comparison state | +| --- | --- | +| `hexa_linearstatic` | Full precheck PASS; U/RF/S numerical FAIL (comparator native1) | +| `hexa_nonlinearstatic` | Full precheck PASS; U/RF/S numerical FAIL (comparator native1) | +| `hexa_dynlinear` | NOT RUN by comparator; root-owned actual CLI run pending | +| `hexa_dynnonlinear` | NOT RUN by comparator; root-owned actual CLI run pending | +| `hexa_freq` | Official HDF5 precheck PASS; modal Value absent, actual comparison NOT RUN (native2) | + +Full MSVC build/CTest, original CLI orchestration, independent numerical/physics +investigation of the reference conflict and whole-feature completion are root-owned +and are not inferred from this Python implementation or its fixture tests. diff --git a/scripts/verification/compare_hex8.py b/scripts/verification/compare_hex8.py new file mode 100644 index 0000000..5b178fa --- /dev/null +++ b/scripts/verification/compare_hex8.py @@ -0,0 +1,749 @@ +"""Compare original Hex8 CSVs with official Solid HDF5; never launch a solver. + +This entry point is intentionally restricted to reference/hexa's fixed original +models. Identity/schema checks precede finite checks and the common tolerance. +Rounded HDF5 observables are compared; internal high+low solver gates are not +recomputed from them. Missing original modal Value data remain NOT RUN. +""" + +import argparse +import csv +from decimal import Decimal, InvalidOperation +import hashlib +import json +import math +import os +from pathlib import Path +import re +from typing import Any + +import h5py +import numpy as np + + +ROOT = Path(__file__).resolve().parents[2] +REFERENCE_DIRECTORY = ROOT / "reference/hexa" +STEP = "Step-1" +INSTANCE = "PART-1-1" +NODE_COUNT = 442 +ELEMENT_COUNT = 300 +POLICY = {"relative_tolerance": 0.05, "relative_rms_tolerance": 0.01, + "near_zero_ratio": 0.01, "zero_scale": "exact-zero"} +CASES = { + "hexa_linearstatic": "linear-static", + "hexa_nonlinearstatic": "nonlinear-static-full-newton", + "hexa_dynlinear": "hht-linear-dynamics", + "hexa_dynnonlinear": "hht-nonlinear-dynamics", + "hexa_freq": "linear-frequency", +} +# Working-tree bytes (Windows CRLF), not Git's normalized LF blobs. +SOURCE_SHA256 = { + "hexa_dynlinear.inp": "f8db81b8ea56b689ab0ad892e6a84b7c6a0162c9d3ff95328ea99722f3de078e", + "hexa_dynlinear_U.csv": "623c20489e8cc18d23fa1728563ce5a26059f321861514ef8569e3d204ef6cb5", + "hexa_dynnonlinear.inp": "1685f86f844e8020bf46e1747b949699ac9aff0d7570bd21ec95109399365ff9", + "hexa_dynnonlinear_U.csv": "f5832a24577f54838bb49dec86a0dc9668cf43e78953f718a9d3635e0f574226", + "hexa_freq.inp": "42cbcec2a8b738f8b6d2c1e13114b8840dccdde6fc3e6f8a43a3272e85d82347", + "hexa_linearstatic.inp": "67183909c5dc962a4d4b65ffb81f7bd65b72e89a5d015c8369a8e0226347a647", + "hexa_linearstatic_RF.csv": "3e04e71ec075ffef098cf4870d77f9c400557fccb375369da8bf67cf67c1e36b", + "hexa_linearstatic_S.csv": "8840a502b9d28a91e5685349301a56fed37d2ac62327b71304716b714db577e8", + "hexa_linearstatic_U.csv": "1ff49aea3ec601cfd796d6a2fb6a4f91be3fc3de9426cdf3744d981216d4bd5f", + "hexa_nonlinearstatic.inp": "7858254025b96c6bfe3fb6498605ecf9c473425c2bc3593f5fbdc6902fd7c6cb", + "hexa_nonlinearstatic_RF.csv": "3ba67a4fb75b2ce3edee13d35b2ea66ac7d52fc16e612c086ed1122a834fee7b", + "hexa_nonlinearstatic_S.csv": "a8866a730f0d65b1e3921d504d255b72158b4a2b9ee30fc691712d94982b07ba", + "hexa_nonlinearstatic_U.csv": "bd6c2fe2b2bb2ffb0bc2ea52f67513c3967a396b5707bbebc9533fde7be25ac6", +} +BASE_COLUMNS = ["ODB Name", "Step", "Frame", "Part Instance Name"] +TAIL_COLUMNS = ["X", "Y", "Z", "Section Name", "Material Name", "Section Point"] +COMPONENTS = {"U": ["U1", "U2", "U3"], "RF": ["RF1", "RF2", "RF3"], + "S": ["S11", "S22", "S33", "S12", "S13", "S23"]} +STRESS_DERIVED = ["Mises", "Max. Principal", "Max. Principal (Abs)", + "Mid. Principal", "Min. Principal", "Tresca", "Pressure", "Third Invariant"] +NATURAL_NODES = np.array([[-1, -1, -1], [1, -1, -1], [1, 1, -1], [-1, 1, -1], + [-1, -1, 1], [1, -1, 1], [1, 1, 1], [-1, 1, 1]]) +NATURAL_POINTS = np.array([[-1, -1, -1], [1, -1, -1], [-1, 1, -1], [1, 1, -1], + [-1, -1, 1], [1, -1, 1], [-1, 1, 1], [1, 1, 1]]) / math.sqrt(3.0) +FRAME_PATTERN = re.compile( + r"Increment\s+(\d+):\s*Step Time\s*=\s*([+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?)" +) + + +class ComparisonError(ValueError): + def __init__(self, code: str, message: str): + super().__init__(message) + self.code = code + + +def require(condition: bool, code: str, message: str): + if not condition: + raise ComparisonError(code, message) + + +def canonical(value: str) -> str: + """Abaqus ASCII case folding only; distinct instances never coalesce.""" + return value.translate(str.maketrans("abcdefghijklmnopqrstuvwxyz", "ABCDEFGHIJKLMNOPQRSTUVWXYZ")) + + +def label(value: str) -> int: + require(re.fullmatch(r"[0-9]+", value) is not None, "invalid-source-label", value) + number = int(value) + require(0 < number <= 2**63 - 1, "invalid-source-label", value) + return number + + +def lexical_absolute(path: Path) -> str: + """Match C++ absolute.lexically_normal without resolving physical aliases.""" + return Path(os.path.abspath(path)).as_posix() + + +def sha256(path: Path) -> str: + with path.open("rb") as stream: + return hashlib.file_digest(stream, "sha256").hexdigest() + + +def fnv1a64(raw: bytes) -> str: + value = 14695981039346656037 + for byte in raw: + value = ((value ^ byte) * 1099511628211) & (2**64 - 1) + return f"fnv1a64:{value:016x}" + + +def input_inventory(path: Path, case: str) -> dict[str, Any]: + """Read only the frozen single-part identity-instance original subset.""" + require(case in CASES, "unknown-source-case", case) + require(sha256(path) == SOURCE_SHA256[f"{case}.inp"], "source-input-hash-mismatch", str(path)) + nodes, elements, sets, boundaries = {}, {}, {}, [] + keyword, options, in_assembly = "", {}, False + model = {"case": case, "nodes": nodes, "elements": elements, "node_sets": sets, + "section_line": None, "source_path": lexical_absolute(path), + "content_identity": fnv1a64(path.read_bytes())} + for line_number, raw in enumerate(path.read_text(encoding="utf-8-sig").splitlines(), 1): + line = raw.strip() + if not line or line.startswith("**"): + continue + if line.startswith("*"): + tokens = [token.strip() for token in line[1:].split(",")] + keyword = canonical(tokens[0]) + options = {canonical(token.split("=", 1)[0]): token.split("=", 1)[1].strip() + for token in tokens[1:] if "=" in token} + if keyword == "ASSEMBLY": + in_assembly = True + elif keyword == "END ASSEMBLY": + in_assembly = False + elif keyword == "INSTANCE": + require(canonical(options["NAME"]) == INSTANCE and len(tokens) == 3, + "unsupported-source-instance", line) + elif keyword == "SOLID SECTION": + model["section_line"] = line_number + model["section_elset"] = options["ELSET"] + model["material"] = options["MATERIAL"] + elif keyword == "STEP": + require(options["NAME"] == STEP, "wrong-source-step", line) + continue + tokens = [token.strip() for token in line.split(",") if token.strip()] + if keyword == "NODE": + node = label(tokens[0]) + require(node not in nodes and len(tokens) == 4, "source-node-inventory", line) + nodes[node] = tuple(float(value) for value in tokens[1:]) + elif keyword == "ELEMENT": + element = label(tokens[0]) + require(element not in elements and len(tokens) == 9 and options["TYPE"] == "C3D8", + "source-element-inventory", line) + elements[element] = tuple(label(value) for value in tokens[1:]) + elif keyword == "NSET" and in_assembly: + name = canonical(options["NSET"]) + require(canonical(options["INSTANCE"]) == INSTANCE, "wrong-source-instance", line) + sets.setdefault(name, []).extend(label(value) for value in tokens) + elif keyword == "BOUNDARY": + first, last = int(tokens[1]), int(tokens[2]) + require(len(tokens) == 3 or float(tokens[3]) == 0.0, "unsupported-source-boundary", line) + boundaries.append((canonical(tokens[0]), first, last)) + elif keyword == "DENSITY": + model["density"] = float(tokens[0]) + elif keyword == "ELASTIC": + model["elastic"] = tuple(float(value) for value in tokens) + require(set(nodes) == set(range(1, NODE_COUNT + 1)) + and set(elements) == set(range(1, ELEMENT_COUNT + 1)), + "source-model-inventory", "Expected original 442 nodes and 300 elements.") + constrained = set() + for name, first, last in boundaries: + require(name in sets and 1 <= first <= last <= 6, "source-boundary-identity", name) + constrained.update((node, component) for node in sets[name] + for component in range(first, min(last, 3) + 1)) + model["constrained"] = constrained + return model + + +def point_coordinates(model, element: int, point: int): + xi = NATURAL_POINTS[point - 1] + shape = np.prod(1.0 + NATURAL_NODES * xi, axis=1) / 8.0 + coordinates = np.array([model["nodes"][node] for node in model["elements"][element]]) + terms = shape[:, None] * coordinates + return np.sum(terms, axis=0), np.sum(np.abs(terms), axis=0) + + +def read_reference(path: Path, case: str, quantity: str, model): + """Validate all row identities before converting any physical values.""" + require(quantity in COMPONENTS, "unknown-quantity", quantity) + dynamic = case.startswith("hexa_dyn") + require(not dynamic or quantity == "U", "wrong-reference-quantity", quantity) + extra = (["Element Label", "IntPt"] if quantity == "S" else ["Node Label"]) + suffix = ([f"S-{name}" for name in STRESS_DERIVED] if quantity == "S" + else [f"{quantity}-Magnitude"]) + headers = BASE_COLUMNS + extra + TAIL_COLUMNS + suffix + [f"{quantity}-{c}" for c in COMPONENTS[quantity]] + records = {} + with path.open(newline="", encoding="utf-8-sig") as stream: + reader = csv.reader(stream) + fields = [field.strip() for field in next(reader, [])] + require(len(fields) == len(set(fields)), "duplicate-reference-header", str(path)) + require(set(fields) == set(headers), "wrong-reference-columns", str(fields)) + for row_number, values in enumerate(reader, 2): + require(len(values) == len(fields), "wrong-reference-row-width", str(row_number)) + row = dict(zip(fields, [value.strip() for value in values])) + require(row["ODB Name"].replace("\\", "/").rsplit("/", 1)[-1] == f"{case}.odb", + "wrong-reference-case", str(row_number)) + require(row["Step"] == STEP, "wrong-reference-step", str(row_number)) + require(canonical(row["Part Instance Name"]) == INSTANCE, + "wrong-reference-instance", str(row_number)) + match = FRAME_PATTERN.fullmatch(row["Frame"]) + require(match is not None, "invalid-reference-frame", str(row_number)) + increment = int(match[1]) + time = Decimal(match[2]) + require(time.is_finite() and increment in (range(101) if dynamic else (10,)) + and time == (Decimal(increment) / 10 if dynamic else Decimal(1)), + "wrong-reference-frame-time", str(row_number)) + entity = label(row["Element Label" if quantity == "S" else "Node Label"]) + point = label(row["IntPt"]) if quantity == "S" else 0 + require(entity in model["elements" if quantity == "S" else "nodes"], + "extra-reference-entity", str(row_number)) + require(quantity != "S" or 1 <= point <= 8, "wrong-reference-point", str(row_number)) + require(row["Section Point"] == "" + and canonical(row["Material Name"]) == ("MATERIAL-1" if quantity == "S" else "NA") + and canonical(row["Section Name"]) == ("ASSEMBLY_PART-1-1_SET-1" if quantity == "S" else "NA"), + "wrong-reference-section", str(row_number)) + key = (increment, entity, point) + require(key not in records, "duplicate-reference-identity", str(key)) + records[key] = row + entities = model["elements" if quantity == "S" else "nodes"] + expected = {(n, entity, p) for n in (range(101) if dynamic else (10,)) + for entity in entities for p in (range(1, 9) if quantity == "S" else (0,))} + require(records.keys() == expected, "reference-inventory", "Missing reference identities.") + return records + + +def reference_values(records, quantity, model): + result = {} + for key in sorted(records): + row = records[key] + values = [float(row[f"{quantity}-{component}"]) for component in COMPONENTS[quantity]] + coordinates = np.array([float(row[axis]) for axis in ("X", "Y", "Z")]) + # Every exported numeric column is checked, including unused derived values. + derived = STRESS_DERIVED if quantity == "S" else ["Magnitude"] + other = [float(row[f"{quantity}-{name}"]) for name in derived] + require(all(math.isfinite(value) for value in values + other) and np.isfinite(coordinates).all(), + "nonfinite-reference-value", str(key)) + if quantity == "S": + expected, scale = point_coordinates(model, key[1], key[2]) + else: + expected = np.array(model["nodes"][key[1]]) + scale = np.abs(expected) + budget = 16.0 * np.finfo(np.float32).eps * scale + require((np.abs(coordinates - expected) <= budget).all(), + "reference-coordinate-convention", str(key)) + result[key] = values + return result + + +def evaluate_family(identity, result, reference): + require(len(result) == len(reference) == len(identity["rows"]) and len(result) > 0, + "scalar-inventory", identity["quantity"]) + require(all(math.isfinite(value) for value in (*result, *reference)), + "nonfinite-comparison-value", "Comparison values must be finite.") + scale = max(abs(value) for value in reference) + band = POLICY["near_zero_ratio"] * scale + rows, error_norm = [], 0.0 + for key, actual, expected in zip(identity["rows"], result, reference): + error = abs(actual - expected) + require(math.isfinite(error), "nonfinite-comparison-value", str(key)) + error_norm = math.hypot(error_norm, error) + require(math.isfinite(error_norm), "nonfinite-comparison-value", str(key)) + relative = 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, threshold = "relative", POLICY["relative_tolerance"] * abs(expected) + relative = error / abs(expected) + require(math.isfinite(relative), "nonfinite-comparison-value", str(key)) + passed = relative <= POLICY["relative_tolerance"] + rows.append({**key, "fesa_value": actual, "reference_value": expected, + "absolute_error": error, "relative_error": relative, "branch": branch, + "threshold": threshold, "passed": passed}) + rows_passed = all(row["passed"] for row in rows) + absolute_rms = error_norm / math.sqrt(len(rows)) + relative_rms = absolute_rms / scale if scale else None + require(relative_rms is None or math.isfinite(relative_rms), "nonfinite-comparison-value", "RMS") + rms_passed = rows_passed if not scale else relative_rms <= POLICY["relative_rms_tolerance"] + return {**{k: v for k, v in identity.items() if k != "rows"}, "impact": "blocking", + "row_count": len(rows), "reference_scale": scale, "near_zero_band": band, + "maximum_absolute_error": max(row["absolute_error"] for row in rows), + "absolute_rms": absolute_rms, "relative_rms": relative_rms, + "policy": dict(POLICY), "rows": rows, "rows_passed": rows_passed, + "rms_passed": rms_passed, "passed": rows_passed and rms_passed, + "diagnostic_code": "zero-reference-scale-nonzero-error" if not scale and not rows_passed else None} + + +def compare_modal(reference, result, quantity: str, unit: str): + """Typed modal rows, not an invented format for the absent original CSV. + + Each row explicitly supplies step, mode, quantity, unit_dimensions and value. + Callers must obtain those semantics from a genuine supplied artifact first. + Eigenvalues and frequencies are never silently converted into one another. + """ + units = {"eigenvalue": "1/time^2", "frequency": "1/time"} + require(quantity in units and unit == units[quantity], "wrong-modal-unit", unit) + inventories = [] + for rows in (reference, result): + mapping = {} + for row in rows: + require(row["quantity"] == quantity, "wrong-modal-quantity", str(row["mode"])) + require(row["unit_dimensions"] == unit, "wrong-modal-unit", str(row["mode"])) + require(row["step"] == STEP, "wrong-modal-step", str(row["mode"])) + mode = row["mode"] + require(type(mode) is int and mode in range(1, 11), "wrong-modal-mode", str(mode)) + require(mode not in mapping, "duplicate-modal-mode", str(mode)) + mapping[mode] = row["value"] + require(set(mapping) == set(range(1, 11)), "modal-mode-inventory", "Lowest ten modes required.") + inventories.append(mapping) + families = [evaluate_family( + {"case": "hexa_freq", "step": STEP, "mode": mode, "quantity": quantity, + "unit_dimensions": unit, "rows": [{"mode": mode, "component": quantity}]}, + [inventories[1][mode]], [inventories[0][mode]]) for mode in range(1, 11)] + return {"families": families, "passed": all(family["passed"] for family in families)} + + +def text(value) -> str: + return value.decode("utf-8", errors="strict") if isinstance(value, bytes) else str(value) + + +def string_attribute(obj, name, expected=None): + require(name in obj.attrs, "missing-result-attribute", f"{obj.name}:{name}") + attr = obj.attrs.get_id(name) + string = h5py.check_string_dtype(attr.dtype) + require(attr.shape == () and string is not None and string.length is None + and string.encoding == "utf-8", "wrong-result-attribute-type", f"{obj.name}:{name}") + value = text(obj.attrs[name]) + require(bool(value) and (expected is None or value == expected), + "wrong-result-attribute", f"{obj.name}:{name}={value}") + return value + + +def numeric_attribute(obj, name, dtype, expected): + require(name in obj.attrs, "missing-result-attribute", f"{obj.name}:{name}") + attr = obj.attrs.get_id(name) + require(attr.shape == () and attr.dtype == np.dtype(dtype), + "wrong-result-attribute-type", f"{obj.name}:{name}") + require(obj.attrs[name] == expected, "wrong-result-attribute", f"{obj.name}:{name}") + + +def dataset(file, path, shape): + require(path in file and isinstance(file[path], h5py.Dataset), "missing-result-dataset", path) + value = file[path] + require(value.shape == shape and value.chunks is None and value.compression is None, + "wrong-result-shape", path) + return value + + +def table(file, path, count, fields): + value = dataset(file, path, (count,)) + require(value.dtype.names == tuple(name for name, _ in fields), "wrong-result-table-schema", path) + for name, dtype in fields: + actual = value.dtype[name] + if dtype == "string": + info = h5py.check_string_dtype(actual) + valid = info is not None and info.length is None and info.encoding == "utf-8" + else: + valid = actual == np.dtype(dtype) + require(valid, "wrong-result-table-schema", f"{path}:{name}") + return value[:] + + +def read_results(path: Path, model): + """Read complete official schema/identities; defer physical finite checks.""" + case = model["case"] + dynamic, nonlinear = case.startswith("hexa_dyn"), "nonlinear" in case + frequency = case == "hexa_freq" + arrays, frames, modal = [], {}, {} + with h5py.File(path, "r") as file: + require("metadata" in file and isinstance(file["metadata"], h5py.Group), + "missing-result-metadata", str(path)) + metadata = file["metadata"] + numeric_attribute(metadata, "schema_version", "= 0 for row in modes), + [("mode_id", "= 0], + "wrong-modal-frequency-inventory", root) + vectors = physical_dataset(file, root + "/vectors", (10, NODE_COUNT, 3), + "UX,UY,UZ", "mass-normalized", "nodal", + extra={"normalization": "phi^T*M*phi=1"}) + verification = root + "/verification" + require(set(file[verification]) == {"pairs", "eta_gram", "max_eta_orth"}, + "wrong-modal-verification-inventory", verification) + pair_object = file[verification + "/pairs"] + require(isinstance(pair_object, h5py.Dataset) and len(pair_object.shape) == 1, + "wrong-modal-pair-shape", verification + "/pairs") + pair_count = pair_object.shape[0] + require(10 <= pair_count <= free_count, "wrong-modal-pair-count", str(pair_count)) + pairs = table(file, verification + "/pairs", pair_count, [ + ("sorted_rank", "= 0).all() and (np.diff(pairs["eigenvalue"]) >= 0).all(), + "wrong-modal-order", model["case"]) + for row in frequencies: + expected = math.sqrt(modes[row["mode_id"] - 1]["eigenvalue"]) / (2 * math.pi) + require(row["frequency"] == expected, "wrong-modal-frequency", str(row["mode_id"])) + + +def compare(case: str, input_path: Path, results_path: Path, reference_directory: Path = REFERENCE_DIRECTORY, + execution_record: Path | None = None): + report = {"schema": "hex8-reference-comparison-v1", "case": case, "policy": dict(POLICY), + "identities": {}, "precheck": {"passed": False}, "families": [], + "blocking_failures": [], "blocking_passed": False, "verdict": "failed"} + try: + require(case in CASES, "unknown-source-case", case) + quantities = ["U"] if case.startswith("hexa_dyn") else (["U", "RF", "S"] if case != "hexa_freq" else []) + paths = {"input": input_path, "results": results_path} + paths.update({quantity: reference_directory / f"{case}_{quantity}.csv" for quantity in quantities}) + for name, path in paths.items(): + report["identities"][name] = {"path": lexical_absolute(path), "sha256": sha256(path)} + model = input_inventory(input_path, case) + report["input_copy_matches_approved_bytes"] = True + report["approved_input_sha256"] = SOURCE_SHA256[f"{case}.inp"] + report["execution"] = {"verified": False, "scope": "file comparison only; execution record not supplied"} + if execution_record is not None: + execution = json.loads(execution_record.read_text(encoding="utf-8-sig")) + command = execution.get("command") + require(execution.get("case") == case and execution.get("native_exit") == 0 + and type(execution.get("native_exit")) is int + and execution.get("original_sha256") == report["identities"]["input"]["sha256"] + and execution.get("scratch_sha256") == report["identities"]["input"]["sha256"] + and execution.get("result_sha256") == report["identities"]["results"]["sha256"] + and execution.get("source_unchanged") is True + and execution.get("scratch_unchanged") is True + and execution.get("result_exists") is True + and isinstance(command, list) and len(command) == 2 + and all(isinstance(item, str) for item in command) + and Path(lexical_absolute(Path(command[1]))) == Path(lexical_absolute(input_path)), + "invalid-execution-evidence", str(execution_record)) + report["execution"] = {"verified": True, "native_exit": 0, "command": command, + "record_path": lexical_absolute(execution_record), + "record_sha256": sha256(execution_record)} + records = {} + for quantity in quantities: + require(report["identities"][quantity]["sha256"] == SOURCE_SHA256[f"{case}_{quantity}.csv"], + "reference-hash-mismatch", str(paths[quantity])) + records[quantity] = read_reference(paths[quantity], case, quantity, model) + result = read_results(results_path, model) + references = {quantity: reference_values(rows, quantity, model) for quantity, rows in records.items()} + check_result_values(result, model) + report["precheck"] = {"passed": True, **result["evidence"], + "reference_row_counts": {q: len(rows) for q, rows in records.items()}, + "coordinate_check": "auxiliary; fixed identity; 16*float32-eps uncancelled source-coordinate budget"} + # Internal maps are useful to validation but unnecessary in the public report. + del report["precheck"]["result_node_map"], report["precheck"]["result_element_map"] + if case == "hexa_freq": + report["verdict"] = "not-run" + report["reason"] = "Original modal Value reference is absent; no format/quantity/unit assumption made." + return report + families = [] + for increment in sorted(result["frames"]): + if not case.startswith("hexa_dyn") and increment != 10: + continue + values = result["frames"][increment] + for quantity in quantities: + actual, expected, identities = [], [], [] + for key in sorted(references[quantity]): + if key[0] != increment: + continue + _, entity, point = key + index = result["evidence"]["result_element_map" if quantity == "S" else "result_node_map"][entity] + row = values[{"U": "nodal/displacement", "RF": "nodal/reaction", "S": "solid/stress"}[quantity]][index] + if quantity == "S": + row = row[point - 1] + for c, component in enumerate(COMPONENTS[quantity]): + if quantity == "RF" and (entity, c + 1) not in model["constrained"]: + continue + identities.append({"instance": INSTANCE, "component": component, + "element_label" if quantity == "S" else "node_label": entity, + **({"integration_point": point} if quantity == "S" else {})}) + actual.append(float(row[c])) + expected.append(references[quantity][key][c]) + identity = {"identity": f"{case}:{STEP}:increment-{increment}:{quantity}-global", + "case": case, "step": STEP, "increment": increment, "quantity": quantity, + "unit_dimensions": {"U": "length", "RF": "force", "S": "force/length^2"}[quantity], + "coordinate_system": "global-cartesian", "rows": identities, + "reference_sha256": report["identities"][quantity]["sha256"]} + if case.startswith("hexa_dyn"): + identity.update({"hdf5_time": increment * 0.1, + "reference_decimal_time": str(Decimal(increment) / 10)}) + families.append(evaluate_family(identity, actual, expected)) + report["families"] = families + report["blocking_failures"] = [f["identity"] for f in families if not f["passed"]] + report["blocking_passed"] = not report["blocking_failures"] + report["verdict"] = "passed" if report["blocking_passed"] else "failed" + except (ComparisonError, OSError, ValueError, InvalidOperation, UnicodeError, KeyError, TypeError) 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("--case", choices=CASES, required=True) + parser.add_argument("--input", type=Path, required=True, help="Actual byte-identical scratch input read by the solver.") + parser.add_argument("--results", type=Path, required=True) + parser.add_argument("--references", type=Path, default=REFERENCE_DIRECTORY) + parser.add_argument("--execution-record", type=Path, help="Closed native CLI record with input/result hashes; never executed here.") + parser.add_argument("--report", type=Path, required=True) + args = parser.parse_args() + protected = {args.input.resolve(), args.results.resolve()} + if args.execution_record is not None: + protected.add(args.execution_record.resolve()) + protected.update(path.resolve() for path in args.references.glob("*")) + protected.update(path.resolve() for path in REFERENCE_DIRECTORY.glob("*")) + report_path = args.report.resolve() + hardlink_alias = args.report.exists() and any( + path.exists() and args.report.samefile(path) for path in protected) + if report_path in protected or hardlink_alias: + parser.error("The report path must not overwrite an input, reference or result.") + report = compare(args.case, args.input, args.results, args.references, args.execution_record) + args.report.parent.mkdir(parents=True, exist_ok=True) + args.report.write_text(json.dumps(report, sort_keys=True, indent=2, allow_nan=False) + "\n", encoding="utf-8") + print(json.dumps({"case": args.case, "verdict": report["verdict"], "precheck": report["precheck"], + "blocking_failure_count": len(report["blocking_failures"]), "report": str(args.report)}, allow_nan=False)) + return 2 if not report["precheck"]["passed"] or report["verdict"] == "not-run" else (0 if report["blocking_passed"] else 1) + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/test_compare_hex8.py b/tests/test_compare_hex8.py new file mode 100644 index 0000000..6034614 --- /dev/null +++ b/tests/test_compare_hex8.py @@ -0,0 +1,774 @@ +"""Independent comparison contracts; generated result fixtures are not source evidence.""" + +import csv +import importlib.util +import json +import math +import os +from pathlib import Path +import shutil +import subprocess +import sys + +import h5py +import numpy as np +import pytest + + +ROOT = Path(__file__).resolve().parents[1] +SPEC = importlib.util.spec_from_file_location( + "compare_hex8", ROOT / "scripts/verification/compare_hex8.py" +) +COMPARE = importlib.util.module_from_spec(SPEC) +SPEC.loader.exec_module(COMPARE) + + +def family(actual, reference): + identities = [{"node_label": i + 1, "component": "U1"} for i in range(len(actual))] + return COMPARE.evaluate_family( + {"case": "hexa_linearstatic", "step": "Step-1", "increment": 10, + "quantity": "U", "unit_dimensions": "length", "rows": identities}, + actual, reference, + ) + + +def test_zero_reference_scale_requires_exact_zero(): + report = family([0.0, 1e-300], [0.0, 0.0]) + assert not report["passed"] + assert report["diagnostic_code"] == "zero-reference-scale-nonzero-error" + + +def test_row_pass_does_not_override_one_percent_rms(): + report = family([1.02] * 10, [1.0] * 10) + assert report["rows_passed"] + assert not report["rms_passed"] + assert not report["passed"] + + +def test_near_zero_band_uses_reference_scale_only(): + report = family([100.0, 1.0], [100.0, 0.0]) + assert report["passed"] + assert report["rows"][1]["branch"] == "near-zero" + assert not family([100.0, 1.0001], [100.0, 0.0])["passed"] + + +def test_nonfinite_family_fails_before_tolerance(): + with pytest.raises(COMPARE.ComparisonError, match="finite"): + family([float("nan")], [1.0]) + + +def test_header_strip_collision_rejected_before_bad_numeric(tmp_path): + path = tmp_path / "collision.csv" + path.write_text("Node Label, Node Label\n1,NaN\n", encoding="utf-8") + with pytest.raises(COMPARE.ComparisonError) as error: + COMPARE.read_reference(path, "hexa_linearstatic", "U", {}) + assert error.value.code == "duplicate-reference-header" + + +def test_modal_quantity_and_units_are_explicit(): + rows = [{"step": "Step-1", "mode": i, "quantity": "eigenvalue", + "unit_dimensions": "1/time^2", "value": float(i)} for i in range(1, 11)] + with pytest.raises(COMPARE.ComparisonError) as error: + COMPARE.compare_modal(rows, rows, "frequency", "1/time") + assert error.value.code == "wrong-modal-quantity" + + +def test_modal_singletons_use_common_one_percent_rms(): + reference = [{"step": "Step-1", "mode": i, "quantity": "eigenvalue", + "unit_dimensions": "1/time^2", "value": float(i)} for i in range(1, 11)] + actual = [dict(row, value=row["value"] * 1.02) for row in reference] + report = COMPARE.compare_modal(reference, actual, "eigenvalue", "1/time^2") + assert not report["passed"] + assert len(report["families"]) == 10 + assert all(row["rows_passed"] and not row["rms_passed"] for row in report["families"]) + +# These schema declarations are independent of the comparator's declarations. +TEXT = h5py.string_dtype(encoding="utf-8") +REFERENCE = ROOT / "reference/hexa" + + +def write_strings(obj, attributes): + for name, value in attributes.items(): + obj.attrs.create(name, value, dtype=TEXT) + + +def compound(group, name, schema, rows): + data = np.empty(len(rows), dtype=np.dtype(schema)) + for i, row in enumerate(rows): + data[i] = row + return group.create_dataset(name, data=data) + + +@pytest.fixture(scope="module") +def source_models(): + return {case: COMPARE.input_inventory(REFERENCE / f"{case}.inp", case) + for case in ("hexa_linearstatic", "hexa_nonlinearstatic", "hexa_dynlinear", "hexa_dynnonlinear", "hexa_freq")} + + +def source_rows(case, quantity): + with (REFERENCE / f"{case}_{quantity}.csv").open(newline="", encoding="utf-8") as stream: + reader = csv.reader(stream) + headers = [value.strip() for value in next(reader)] + return [dict(zip(headers, [value.strip() for value in row])) for row in reader] + + +def make_result(path, model): + """Build an original-sized schema fixture, not a validated physical solution.""" + case = model["case"] + dynamic = case in ("hexa_dynlinear", "hexa_dynnonlinear") + nonlinear = case in ("hexa_nonlinearstatic", "hexa_dynnonlinear") + frequency = case == "hexa_freq" + reference_values = {} + if not frequency: + for quantity in (["U"] if dynamic else ["U", "RF", "S"]): + table_rows = {} + for row in source_rows(case, quantity): + increment = int(row["Frame"].split(":")[0].split()[-1]) + entity = int(row["Element Label" if quantity == "S" else "Node Label"]) + point = int(row["IntPt"]) if quantity == "S" else 0 + names = ("S11", "S22", "S33", "S12", "S13", "S23") if quantity == "S" else tuple(f"{quantity}{c}" for c in (1, 2, 3)) + table_rows[increment, entity, point] = [float(row[f"{quantity}-{name}"]) for name in names] + reference_values[quantity] = table_rows + with h5py.File(path, "w") as file: + metadata = file.create_group("metadata") + metadata.attrs.create("schema_version", np.uint64(1), dtype="f8", "wrong-result-value-type")]) +def test_wrong_shape_or_dtype_is_not_sliced_into_supported_components(tmp_path, result_fixtures, shape, dtype, code): + path = copy_result(tmp_path, result_fixtures) + with h5py.File(path, "r+") as file: + key = "steps/Step-1/frames/0/nodal/displacement" + del file[key] + file.create_dataset(key, shape=shape, dtype=dtype) + assert compare_fixture(path)["precheck"]["code"] == code + + +def test_signed_ids_and_fixed_ascii_are_not_official_compound_schema(tmp_path, result_fixtures): + path = copy_result(tmp_path, result_fixtures) + with h5py.File(path, "r+") as file: + old = file["model/nodes"][:] + del file["model/nodes"] + compound(file["model"], "nodes", [("internal_node_id", "