KOKO\Mimi
750 lines
44 KiB
Python
750 lines
44 KiB
Python
"""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", "<u8", 1)
|
|
for name, expected in {
|
|
"feature_id": "hex8-structural-analysis", "analysis_procedure": CASES[case],
|
|
"element_formulation": "FESA-Hex8", "integration_rule": "gauss-2x2x2",
|
|
"mass_formulation": "consistent-reference", "source_step_name": STEP,
|
|
"unit_system_label": "user-consistent-unspecified",
|
|
"coordinate_convention": "global-cartesian; reference-strain; current-cauchy-stress",
|
|
"component_order": "node-major; UX,UY,UZ",
|
|
"source_input_identity": f"path={model['source_path']};content_identity={model['content_identity']}",
|
|
}.items():
|
|
string_attribute(metadata, name, expected)
|
|
string_attribute(metadata, "solver_version")
|
|
nodes = table(file, "model/nodes", NODE_COUNT, [
|
|
("internal_node_id", "<u8"), ("instance_name", "string"),
|
|
("source_label", "string"), ("coordinates", ("<f8", (3,)))])
|
|
node_map = {}
|
|
for i, row in enumerate(nodes):
|
|
node = label(text(row["source_label"]))
|
|
require(int(row["internal_node_id"]) == i and canonical(text(row["instance_name"])) == INSTANCE,
|
|
"wrong-result-node-identity", str(i))
|
|
require(node not in node_map and node in model["nodes"], "result-node-inventory", str(node))
|
|
node_map[node] = i
|
|
require(node_map.keys() == model["nodes"].keys(), "result-node-inventory", str(path))
|
|
elements = table(file, "model/elements", ELEMENT_COUNT, [
|
|
("internal_element_id", "<u8"), ("instance_name", "string"), ("source_label", "string"),
|
|
("source_element_type", "string"), ("internal_formulation", "string"),
|
|
("node_internal_ids", ("<u8", (8,))), ("solid_section_internal_id", "<u8"),
|
|
("material_internal_id", "<u8")])
|
|
element_map = {}
|
|
for i, row in enumerate(elements):
|
|
element = label(text(row["source_label"]))
|
|
require(int(row["internal_element_id"]) == i and canonical(text(row["instance_name"])) == INSTANCE,
|
|
"wrong-result-element-identity", str(i))
|
|
require(element not in element_map and element in model["elements"],
|
|
"result-element-inventory", str(element))
|
|
require(text(row["source_element_type"]) == "C3D8"
|
|
and text(row["internal_formulation"]) == "FESA-Hex8"
|
|
and tuple(row["node_internal_ids"]) == tuple(node_map[n] for n in model["elements"][element])
|
|
and row["solid_section_internal_id"] == 0 and row["material_internal_id"] == 0,
|
|
"wrong-result-element-connectivity", str(element))
|
|
element_map[element] = i
|
|
materials = table(file, "model/solid/materials", 1, [
|
|
("internal_material_id", "<u8"), ("name", "string"), ("youngs_modulus", "<f8"),
|
|
("poisson_ratio", "<f8"), ("density_present", "u1"), ("density", "<f8")])
|
|
material = materials[0]
|
|
require(material["internal_material_id"] == 0 and canonical(text(material["name"])) == "MATERIAL-1"
|
|
and material["density_present"] == 1, "wrong-result-material-identity", str(path))
|
|
sections = table(file, "model/solid/sections", 1, [
|
|
("internal_section_id", "<u8"), ("source_file", "string"), ("source_line", "<u8"),
|
|
("source_elset", "string"), ("material_internal_id", "<u8")])
|
|
section = sections[0]
|
|
require(section["internal_section_id"] == 0 and section["material_internal_id"] == 0
|
|
and text(section["source_file"]) == model["source_path"]
|
|
and section["source_line"] == model["section_line"]
|
|
and canonical(text(section["source_elset"])) == canonical(model["section_elset"]),
|
|
"wrong-result-section-identity", str(path))
|
|
points = table(file, "model/solid_integration_points", ELEMENT_COUNT * 8, [
|
|
("internal_element_id", "<u8"), ("point_id", "u1"),
|
|
("natural_coordinates", ("<f8", (3,))), ("reference_coordinates", ("<f8", (3,)))])
|
|
for i, row in enumerate(points):
|
|
require(row["internal_element_id"] == i // 8 and row["point_id"] == i % 8 + 1,
|
|
"wrong-result-point-identity", str(i))
|
|
mask = dataset(file, "model/nodal_constraint_mask", (NODE_COUNT, 3))
|
|
require(mask.dtype == np.dtype("u1"), "wrong-result-mask-type", str(path))
|
|
string_attribute(mask, "component_names", "UX,UY,UZ")
|
|
mask = mask[:]
|
|
expected_mask = np.array([[(node, c) in model["constrained"] for c in (1, 2, 3)]
|
|
for node in node_map], dtype="u1")
|
|
require(np.array_equal(mask, expected_mask), "wrong-result-constraint-identity", str(path))
|
|
require("steps" in file and set(file["steps"]) == {STEP}, "wrong-result-step", str(path))
|
|
step = file[f"steps/{STEP}"]
|
|
expected_step_groups = {"eigen"} if frequency else ({"frames", "accepted_frames"} if dynamic or nonlinear else {"frames"})
|
|
require(set(step) == expected_step_groups, "wrong-result-step-inventory", str(path))
|
|
if frequency:
|
|
modal = read_modal_tables(file, arrays, int((mask == 0).sum()))
|
|
else:
|
|
count = 101 if dynamic else (10 if nonlinear else 1)
|
|
accepted = None
|
|
if dynamic or nonlinear:
|
|
coordinate = "time" if dynamic else "load_factor"
|
|
flag = "is_initial" if dynamic else "is_factor_one"
|
|
accepted = table(file, f"steps/{STEP}/accepted_frames", count, [
|
|
("frame_index", "<u8"), ("increment_index", "<u8"), (coordinate, "<f8"),
|
|
("iteration_count", "<u8"), (flag, "u1")])
|
|
for frame, row in enumerate(accepted):
|
|
increment = frame if dynamic else frame + 1
|
|
value = increment * 0.1
|
|
expected_flag = int(frame == (0 if dynamic else count - 1))
|
|
require(row["frame_index"] == frame and row["increment_index"] == increment
|
|
and row[coordinate] == value and row[flag] == expected_flag
|
|
and (row["iteration_count"] == 0 if dynamic and frame == 0
|
|
else 1 <= row["iteration_count"] <= 100),
|
|
"wrong-result-frame-schedule", str(frame))
|
|
frame_root = f"steps/{STEP}/frames"
|
|
require(set(file[frame_root]) == {str(i) for i in range(count)},
|
|
"result-frame-inventory", str(path))
|
|
for frame in range(count):
|
|
increment = frame if dynamic else (frame + 1 if nonlinear else 10)
|
|
value = increment * 0.1 if dynamic or nonlinear else 1.0
|
|
frames[increment] = read_frame(file, frame, increment, value, dynamic, nonlinear, arrays)
|
|
# Identity and schema checks above deliberately finish before physical checks.
|
|
arrays.extend([("model/nodes/coordinates", nodes["coordinates"]),
|
|
("model/points/natural", points["natural_coordinates"]),
|
|
("model/points/reference", points["reference_coordinates"]),
|
|
("model/material", np.array([material[n] for n in ("youngs_modulus", "poisson_ratio", "density")]))])
|
|
evidence = {"source_input_identity": string_attribute(metadata, "source_input_identity"),
|
|
"source_node_count": NODE_COUNT, "source_element_count": ELEMENT_COUNT,
|
|
"accepted_frame_count": len(frames), "constrained_scalar_count": int(mask.sum()),
|
|
"result_node_map": node_map, "result_element_map": element_map}
|
|
return {"frames": frames, "modal": modal, "nodes": nodes, "points": points,
|
|
"material": material, "arrays": arrays, "mask": mask, "evidence": evidence}
|
|
|
|
|
|
def physical_dataset(file, path, shape, names, units, location, frame=None,
|
|
increment=None, coordinate=None, dynamic=False, extra=None):
|
|
value = dataset(file, path, shape)
|
|
require(value.dtype == np.dtype("<f8"), "wrong-result-value-type", path)
|
|
attributes = {"component_names": names, "unit_dimensions": units,
|
|
"coordinate_system": "global-cartesian", "location": location}
|
|
if frame is not None:
|
|
attributes["step_name"] = STEP
|
|
attributes.update(extra or {})
|
|
for name, expected in attributes.items():
|
|
string_attribute(value, name, expected)
|
|
if frame is not None:
|
|
numeric_attribute(value, "frame_index", "<u8", frame)
|
|
numeric_attribute(value, "increment_index", "<u8", increment)
|
|
numeric_attribute(value, "time" if dynamic else "load_factor", "<f8", coordinate)
|
|
return value[()]
|
|
|
|
|
|
def read_frame(file, frame, increment, coordinate, dynamic, nonlinear, arrays):
|
|
root = f"steps/{STEP}/frames/{frame}"
|
|
require(set(file[root]) == {"nodal", "solid", "global"}, "wrong-result-frame-inventory", root)
|
|
definitions = {
|
|
"nodal/displacement": ((NODE_COUNT, 3), "U1,U2,U3", "length", "nodal", {}),
|
|
"nodal/reaction": ((NODE_COUNT, 3), "RF1,RF2,RF3", "force", "nodal", {}),
|
|
"nodal/residual": ((NODE_COUNT, 3), "R1,R2,R3", "force", "nodal", {}),
|
|
"solid/stress": ((ELEMENT_COUNT, 8, 6), "S11,S22,S33,S12,S13,S23", "force/length^2",
|
|
"integration-point", {"stress_measure": "cauchy"}),
|
|
"solid/strain": ((ELEMENT_COUNT, 8, 6), "E11,E22,E33,E12,E13,E23", "dimensionless",
|
|
"integration-point", {"strain_measure": "green-lagrange" if nonlinear else "infinitesimal",
|
|
"shear_convention": "tensor"}),
|
|
"global/strain_energy": ((), "scalar", "force*length", "global",
|
|
{"definition": "sum(W*detJ0)" if nonlinear else "0.5*u^T*K0*u"}),
|
|
}
|
|
if dynamic:
|
|
definitions.update({
|
|
"nodal/velocity": ((NODE_COUNT, 3), "V1,V2,V3", "length/time", "nodal", {}),
|
|
"nodal/acceleration": ((NODE_COUNT, 3), "A1,A2,A3", "length/time^2", "nodal", {"acceleration_kind": "physical"}),
|
|
"nodal/algorithmic_acceleration": ((NODE_COUNT, 3), "A1,A2,A3", "length/time^2", "nodal", {"acceleration_kind": "algorithmic"}),
|
|
"global/kinetic_energy": ((), "scalar", "force*length", "global", {"definition": "0.5*v^T*M*v"}),
|
|
"global/external_work": ((), "scalar", "force*length", "global", {"definition": "sum(0.5*(F_old+F_new)*delta_u)"}),
|
|
"global/energy_balance_error": ((), "scalar", "force*length", "global", {"definition": "T+U-T0-U0-Wext"}),
|
|
"global/linear_momentum": ((3,), "P1,P2,P3", "mass*length/time", "global", {}),
|
|
"global/angular_momentum": ((3,), "L1,L2,L3", "mass*length^2/time", "global", {"origin": "global0"}),
|
|
})
|
|
for group in ("nodal", "solid", "global"):
|
|
expected = {key.split("/", 1)[1] for key in definitions if key.startswith(group + "/")}
|
|
require(set(file[f"{root}/{group}"]) == expected, "wrong-result-quantity-inventory", root + "/" + group)
|
|
values = {}
|
|
for quantity, (shape, names, units, location, extra) in definitions.items():
|
|
path = root + "/" + quantity
|
|
values[quantity] = physical_dataset(file, path, shape, names, units, location,
|
|
frame, increment, coordinate, dynamic, extra)
|
|
arrays.append((path, values[quantity]))
|
|
return values
|
|
|
|
|
|
def read_modal_tables(file, arrays, free_count):
|
|
root = f"steps/{STEP}/eigen"
|
|
require(set(file[root]) == {"modes", "frequencies", "vectors", "verification"},
|
|
"wrong-modal-result-inventory", root)
|
|
modes = table(file, root + "/modes", 10, [("mode_id", "<u8"), ("eigenvalue", "<f8"),
|
|
("classification", "u1"), ("cluster_id", "<u8")])
|
|
for name, component, units in (("modes", "eigenvalue", "1/time^2"),
|
|
("frequencies", "frequency", "1/time")):
|
|
obj = file[root + "/" + name]
|
|
for attr, expected in {"component_names": component, "unit_dimensions": units,
|
|
"coordinate_system": "global-cartesian", "location": "global"}.items():
|
|
string_attribute(obj, attr, expected)
|
|
frequencies = table(file, root + "/frequencies", sum(row["eigenvalue"] >= 0 for row in modes),
|
|
[("mode_id", "<u8"), ("frequency", "<f8")])
|
|
require(list(modes["mode_id"]) == list(range(1, 11)) and (modes["classification"] <= 1).all()
|
|
and modes["cluster_id"][0] == 1
|
|
and all(b in (a, a + 1) for a, b in zip(modes["cluster_id"], modes["cluster_id"][1:])),
|
|
"wrong-modal-mode", root)
|
|
require(list(frequencies["mode_id"]) == [row["mode_id"] for row in modes if row["eigenvalue"] >= 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", "<u8"), ("solver_column", "<u8"), ("eigenvalue", "<f8"),
|
|
("modal_mass", "<f8"), ("rayleigh_quotient", "<f8"), ("eta_r", "<f8"),
|
|
("eta_norm", "<f8"), ("eta_energy", "<f8")])
|
|
require(list(pairs["sorted_rank"]) == list(range(1, pair_count + 1))
|
|
and len(set(pairs["solver_column"])) == pair_count
|
|
and (pairs["solver_column"] < pair_count).all(), "wrong-modal-pair-identity", verification)
|
|
for name in ("eta_gram", "max_eta_orth"):
|
|
value = dataset(file, verification + "/" + name, ())
|
|
require(value.dtype == np.dtype("<f8"), "wrong-result-value-type", name)
|
|
arrays.append((verification + "/" + name, value[()]))
|
|
arrays.extend([(root + "/vectors", vectors), (root + "/lambda", modes["eigenvalue"]),
|
|
(root + "/frequency", frequencies["frequency"])])
|
|
arrays.extend([(verification + "/" + name, pairs[name]) for name in pairs.dtype.names[2:]])
|
|
return {"modes": modes, "frequencies": frequencies, "pairs": pairs}
|
|
|
|
|
|
def check_result_values(result, model):
|
|
for path, values in result["arrays"]:
|
|
require(np.isfinite(values).all(), "nonfinite-result-value", path)
|
|
for row in result["nodes"]:
|
|
node = label(text(row["source_label"]))
|
|
require(tuple(row["coordinates"]) == model["nodes"][node], "wrong-result-coordinate", str(node))
|
|
material = result["material"]
|
|
require(tuple(material[name] for name in ("youngs_modulus", "poisson_ratio")) == model["elastic"]
|
|
and material["density"] == model["density"], "wrong-result-material-value", model["case"])
|
|
for i, row in enumerate(result["points"]):
|
|
element = list(result["evidence"]["result_element_map"])[i // 8]
|
|
expected, scale = point_coordinates(model, element, i % 8 + 1)
|
|
require(np.array_equal(row["natural_coordinates"], NATURAL_POINTS[i % 8])
|
|
and (np.abs(row["reference_coordinates"] - expected) <= 16 * np.finfo(float).eps * scale).all(),
|
|
"result-coordinate-convention", str(i))
|
|
for increment, values in result["frames"].items():
|
|
require((values["nodal/reaction"][result["mask"] == 0] == 0.0).all(),
|
|
"nonzero-free-result-reaction", str(increment))
|
|
if model["case"].startswith("hexa_dyn") and increment == 0:
|
|
require((values["nodal/displacement"] == 0.0).all() and (values["nodal/velocity"] == 0.0).all(),
|
|
"nonzero-initial-result-state", str(increment))
|
|
if result["modal"]:
|
|
modes, frequencies, pairs = (result["modal"][key] for key in ("modes", "frequencies", "pairs"))
|
|
require((np.diff(modes["eigenvalue"]) >= 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())
|