Files
FESADev/tests/test_compare_hht_dynamic.py

437 lines
22 KiB
Python

"""HHT reference identity, time, projection and common-tolerance regressions."""
import csv
import json
import math
from pathlib import Path
import shutil
import subprocess
import sys
import pytest
h5py = pytest.importorskip("h5py", reason="Run with uv --with h5py for HHT comparison tests.")
np = pytest.importorskip("numpy")
from scripts.verification import compare_hht_dynamic as comparator
@pytest.fixture(scope="module")
def artifacts(tmp_path_factory):
directory = tmp_path_factory.mktemp("hht-comparison")
input_path = comparator.REFERENCE_DIRECTORY / f"{comparator.CASE}.inp"
reference_path = comparator.REFERENCE_DIRECTORY / f"{comparator.CASE} displacements.csv"
nodes, content_identity = comparator.input_inventory(input_path)
reference = comparator.read_reference(reference_path)
result_path = directory / "matching_results.h5"
string_type = h5py.string_dtype("utf-8")
node_type = np.dtype([("internal_node_id", "<u8"), ("instance_name", string_type),
("source_label", string_type), ("coordinates", "<f8", (3,))])
node_rows = np.empty(729, dtype=node_type)
labels = list(nodes)
for index, label in enumerate(labels):
node_rows[index] = (index, "PART-1-1", str(label), nodes[label])
frame_type = np.dtype([("frame_index", "<u8"), ("increment_index", "<u8"),
("time", "<f8"), ("iteration_count", "<u8"), ("is_initial", "u1")])
frames = np.empty(101, dtype=frame_type)
with h5py.File(result_path, "w") as file:
metadata = file.create_group("metadata")
metadata.attrs.update({
"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-1", "step_name": "Step-1", "unit_system_label": "consistent-input-units",
"solver_version": "test-fixture", "source_input_identity": f"path={input_path};content_identity={content_identity}",
"schema_version": np.uint64(1), "increment_count": np.uint64(100), "maximum_increments": np.uint64(1000),
"alpha": -0.05, "beta": 0.275625, "gamma": 0.55, "initial_increment": 0.1, "time_period": 10.0,
})
file.create_dataset("model/nodes", data=node_rows)
for increment in range(101):
time = float(increment) * 0.1
frames[increment] = (increment, increment, time, 0 if increment == 0 else 2, int(increment == 0))
values = np.zeros((729, 6), dtype="<f8")
for label in (2, 21, 22, 69, 70, 165, 309, 442, 586):
values[labels.index(label)] = reference[(increment, label)]
dataset = file.create_dataset(f"steps/Step-1/frames/{increment}/nodal/displacement", data=values)
dataset.attrs.update({
"component_names": "UX,UY,UZ,URX,URY,URZ", "component_unit_dimensions": "length,length,length,1,1,1",
"coordinate_system": "global-cartesian", "location": "nodal", "step_name": "Step-1", "source_step_name": "Step-1",
"rotation_representation": "total-rotation-vector-chart",
"quantity_semantics": "global-translation-and-total-rotation-vector",
"schema_version": np.uint64(1), "frame_index": np.uint64(increment), "increment_index": np.uint64(increment),
"time": time, "is_initial": np.uint8(increment == 0),
})
file.create_dataset("steps/Step-1/accepted_frames", data=frames)
return input_path, reference_path, result_path
def changed_result(tmp_path, artifacts):
target = tmp_path / "changed.h5"
shutil.copyfile(artifacts[2], target)
return target
def reduced_result(tmp_path, artifacts):
path = changed_result(tmp_path, artifacts)
with h5py.File(path, "r+") as file:
file["metadata"].attrs.modify("schema_version", np.uint64(2))
file["metadata"].attrs["time_integration_formulation"] = (
"reduced-coordinate-HHT-with-physical-rate-reconstruction")
for frame in range(101):
root = f"steps/Step-1/frames/{frame}/nodal"
displacement = file[f"{root}/displacement"]
displacement.attrs.modify("schema_version", np.uint64(2))
auxiliary = file.create_dataset(f"{root}/algorithmic_acceleration",
data=np.zeros((729, 6), dtype="<f8"))
for name, value in displacement.attrs.items():
auxiliary.attrs[name] = value
auxiliary.attrs["component_names"] = "AUX1,AUX2,AUX3,AUX4,AUX5,AUX6"
auxiliary.attrs["coordinate_system"] = "fixed-initial-global-chart"
auxiliary.attrs["component_unit_dimensions"] = (
"length/time^2,length/time^2,length/time^2,1/time^2,1/time^2,1/time^2")
auxiliary.attrs["quantity_semantics"] = (
"full-chart-HHT-filter-combination-not-predictor-driver")
return path
def test_reduced_schema_preserves_displacement_comparison(tmp_path, artifacts):
result = comparator.compare(artifacts[0], artifacts[1], reduced_result(tmp_path, artifacts))
assert result["precheck"]["passed"], result["precheck"]
assert result["verdict"] == "passed"
@pytest.mark.parametrize("mutation", ["missing_A", "A_semantics", "A_version", "A_nan", "method"])
def test_reduced_schema_rejects_inconsistent_auxiliary_history(tmp_path, artifacts, mutation):
path = reduced_result(tmp_path, artifacts)
with h5py.File(path, "r+") as file:
auxiliary = file["steps/Step-1/frames/99/nodal/algorithmic_acceleration"]
if mutation == "missing_A":
del file[auxiliary.name]
elif mutation == "A_semantics":
auxiliary.attrs["quantity_semantics"] = "auxiliary-HHT-chart-acceleration"
elif mutation == "A_version":
auxiliary.attrs.modify("schema_version", np.uint64(1))
elif mutation == "A_nan":
auxiliary[0, 0] = float("nan")
else:
file["metadata"].attrs["time_integration_formulation"] = "full-coordinate-HHT"
result = comparator.compare(artifacts[0], artifacts[1], path)
assert not result["precheck"]["passed"]
assert result["families"] == []
def test_original_csv_exact_matching_history_and_deterministic_report(artifacts):
first = comparator.compare(*artifacts)
second = comparator.compare(*artifacts)
assert first == second
assert first["precheck"]["passed"], first["precheck"]
assert first["precheck"]["reference_row_count"] == 909
assert first["precheck"]["checked_result_scalars"] == 441774
assert first["verdict"] == "passed"
assert first["warnings"] == []
assert len(first["families"]) == 202
assert sum(family["row_count"] for family in first["families"]) == 5454
assert [family["quantity"] for family in first["families"][:4]] == ["U", "UR", "U", "UR"]
assert all(family["passed"] for family in first["families"])
assert first["families"][0]["reference_scale"] == 0.0
assert first["families"][6]["hdf5_time"] == 3 * 0.1
assert first["families"][6]["reference_decimal_time"] == "0.3"
assert 3 * 0.1 != float("0.3") # Both legitimate representations were validated independently.
json.dumps(first, allow_nan=False)
@pytest.mark.parametrize("mutation,code", [
("procedure", "wrong-hdf5-attribute-value"),
("mass", "wrong-hdf5-attribute-value"),
("coefficient", "wrong-hdf5-attribute-value"),
("initial_convention", "wrong-hdf5-attribute-value"),
("input_content", "wrong-input-content-identity"),
("component", "wrong-hdf5-attribute-value"),
("rotation", "wrong-hdf5-attribute-value"),
("missing_component_attribute", "missing-hdf5-attribute"),
("nonprojected_nan", "nonfinite-result-value"),
("rotation_inf", "nonfinite-result-value"),
("initial_motion", "nonzero-initial-displacement"),
("duplicate_node", "duplicate-result-node"),
("extra_node_label", "extra-result-node"),
("wrong_instance", "wrong-result-instance"),
("coordinate", "wrong-result-coordinate"),
("wrong_shape", "wrong-hdf5-dataset-layout"),
("float32", "wrong-displacement-type"),
("missing_frame", "result-frame-inventory"),
("extra_frame", "result-frame-inventory"),
("leading_zero_frame", "result-frame-inventory"),
("duplicate_increment", "wrong-result-frame-time"),
("decimal_binary_time", "wrong-result-frame-time"),
("frame_attribute_time", "wrong-hdf5-attribute-value"),
("trial_iteration", "wrong-result-frame-time"),
("initial_iteration", "wrong-result-frame-time"),
])
def test_hdf5_precheck_rejects_before_tolerance(tmp_path, artifacts, mutation, code):
path = changed_result(tmp_path, artifacts)
with h5py.File(path, "r+") as file:
metadata = file["metadata"]
displacement = file["steps/Step-1/frames/3/nodal/displacement"]
nodes = file["model/nodes"]
table = file["steps/Step-1/accepted_frames"]
if mutation in {"procedure", "mass", "coefficient", "initial_convention", "input_content"}:
key, value = {
"procedure": ("analysis_procedure", "linear-static"),
"mass": ("mass_formulation", "consistent-kinetic"),
"coefficient": ("alpha", 0.0),
"initial_convention": ("initial_state", "unloaded"),
"input_content": ("source_input_identity", "wrong-content"),
}[mutation]
metadata.attrs.modify(key, value)
elif mutation in {"component", "rotation", "frame_attribute_time"}:
key, value = {"component": ("component_names", "U1,U2,U3,UR1,UR2,UR3"),
"rotation": ("rotation_representation", "incremental-rotation"),
"frame_attribute_time": ("time", 0.3)}[mutation]
displacement.attrs.modify(key, value)
elif mutation == "missing_component_attribute":
del displacement.attrs["component_names"]
elif mutation in {"nonprojected_nan", "rotation_inf"}:
displacement[0, 0 if mutation == "nonprojected_nan" else 3] = float("nan") if mutation == "nonprojected_nan" else float("inf")
elif mutation == "initial_motion":
file["steps/Step-1/frames/0/nodal/displacement"][0, 0] = 1.0e-100
elif mutation in {"duplicate_node", "extra_node_label", "wrong_instance", "coordinate"}:
rows = nodes[:]
if mutation == "duplicate_node":
rows[1]["source_label"] = rows[0]["source_label"]
elif mutation == "extra_node_label":
rows[0]["source_label"] = "99999"
elif mutation == "wrong_instance":
rows[0]["instance_name"] = "Other"
else:
rows[0]["coordinates"] = [1.0, 2.0, 3.0]
nodes[:] = rows
elif mutation in {"wrong_shape", "float32"}:
name = displacement.name
values = displacement[:]
del file[name]
file.create_dataset(name, data=values[:-1] if mutation == "wrong_shape" else values.astype("<f4"))
elif mutation == "missing_frame":
del file["steps/Step-1/frames/100"]
elif mutation == "extra_frame":
file.create_group("steps/Step-1/frames/101")
elif mutation == "leading_zero_frame":
file.move("steps/Step-1/frames/1", "steps/Step-1/frames/01")
else:
rows = table[:]
if mutation == "duplicate_increment":
rows[3]["increment_index"] = 2
elif mutation == "decimal_binary_time":
rows[3]["time"] = 0.3
elif mutation == "trial_iteration":
rows[3]["iteration_count"] = 101
elif mutation == "initial_iteration":
rows[0]["iteration_count"] = 1
table[:] = rows
report = comparator.compare(artifacts[0], artifacts[1], path)
assert report["precheck"]["passed"] is False
assert report["precheck"]["code"] == code, report["precheck"]
assert report["families"] == []
assert report["verdict"] == "failed"
@pytest.mark.parametrize("mutation,code", [
("duplicate_header", "duplicate-reference-column"),
("missing_header", "missing-reference-column"),
("duplicate_row", "duplicate-reference-row"),
("missing_row", "reference-row-inventory"),
("extra_node", "reference-row-inventory"),
("extra_frame", "reference-row-inventory"),
("time", "wrong-reference-time"),
("instance", "wrong-reference-instance"),
("label", "invalid-node-label"),
("nonfinite_rotation", "nonfinite-reference-value"),
])
def test_reference_precheck_rejects_before_tolerance(tmp_path, artifacts, mutation, code):
with artifacts[1].open(newline="", encoding="utf-8") as stream:
rows = list(csv.reader(stream))
if mutation == "duplicate_header":
rows[0][-1] = rows[0][-2]
elif mutation == "missing_header":
rows[0][-1] = "Not a required component"
elif mutation == "duplicate_row":
rows.append(rows[1])
elif mutation == "missing_row":
rows.pop()
elif mutation in {"extra_node", "extra_frame"}:
row = rows[1].copy()
if mutation == "extra_node":
row[2] = "1"
else:
row[0] = "Increment 101: Step Time = 10.1"
rows.append(row)
elif mutation == "time":
rows[1][0] = "Increment 0: Step Time = 0.0000000001"
elif mutation == "instance":
rows[1][1] = "WRONG-INSTANCE"
elif mutation == "label":
rows[1][2] = "2.0"
else:
rows[1][-1] = "nan"
reference = tmp_path / "changed.csv"
with reference.open("w", newline="", encoding="utf-8") as stream:
csv.writer(stream).writerows(rows)
report = comparator.compare(artifacts[0], reference, artifacts[2])
assert not report["precheck"]["passed"]
assert report["precheck"]["code"] == code, report["precheck"]
assert report["families"] == []
def test_projection_joins_source_identity_and_csv_order_is_irrelevant(tmp_path, artifacts):
path = changed_result(tmp_path, artifacts)
with h5py.File(path, "r+") as file:
nodes = file["model/nodes"][:][::-1]
nodes["internal_node_id"] = np.arange(729, dtype="<u8")
file["model/nodes"][:] = nodes
for frame in range(101):
values = file[f"steps/Step-1/frames/{frame}/nodal/displacement"]
values[:] = values[:][::-1]
with artifacts[1].open(newline="", encoding="utf-8") as stream:
rows = list(csv.reader(stream))
reference = tmp_path / "reordered.csv"
with reference.open("w", newline="", encoding="utf-8") as stream:
csv.writer(stream).writerows([rows[0], *reversed(rows[1:])])
report = comparator.compare(artifacts[0], reference, path)
assert report["precheck"]["passed"], report["precheck"]
assert report["verdict"] == "passed"
@pytest.mark.parametrize("component,passed", [(0, False), (3, True)])
def test_translation_blocks_rotation_warns(tmp_path, artifacts, component, passed):
path = changed_result(tmp_path, artifacts)
with h5py.File(path, "r+") as file:
values = file["steps/Step-1/frames/1/nodal/displacement"]
rows = values[:]
rows[:, component:component + 3] *= 2.0
values[:] = rows
report = comparator.compare(artifacts[0], artifacts[1], path)
assert report["precheck"]["passed"]
assert report["blocking_passed"] is passed
assert len(report["warnings"]) == (2 if passed else 0)
if passed:
assert [warning["code"] for warning in report["warnings"]] == ["row-tolerance-failure", "relative-rms-failure"]
else:
assert report["blocking_failures"] == ["U-global:Step-1:increment-1"]
def test_common_policy_boundaries_and_zero_scale():
near = comparator.evaluate_family(1, "U", [100.0, 2.0, 1.0], [100.0, 1.0, 0.0])
assert near["passed"]
assert [row["branch"] for row in near["rows"]] == ["relative", "near-zero", "near-zero"]
relative = comparator.evaluate_family(1, "U", [105.0], [100.0])
assert relative["rows_passed"] and not relative["rms_passed"]
boundary = comparator.evaluate_family(1, "U", [101.0], [100.0])
assert boundary["passed"]
assert boundary["relative_rms"] == 0.01
outside = comparator.evaluate_family(1, "U", [math.nextafter(101.0, math.inf)], [100.0])
assert outside["rows_passed"] and not outside["rms_passed"]
zero = comparator.evaluate_family(0, "U", [0.0], [0.0])
assert zero["passed"] and zero["relative_rms"] is None
nonzero = comparator.evaluate_family(0, "U", [1.0e-300], [0.0])
assert not nonzero["passed"]
assert nonzero["diagnostic_code"] == "zero-reference-scale-nonzero-error"
json.dumps(nonzero, allow_nan=False)
with pytest.raises(comparator.ComparisonError, match="U"):
comparator.evaluate_family(1, "U", [1.0e308], [-1.0e308])
def test_cli_writes_json_and_exit_two_for_missing_results(tmp_path, artifacts):
report = tmp_path / "evidence.json"
completed = subprocess.run([sys.executable, "-B", str(Path(comparator.__file__)),
"--input", str(artifacts[0]), "--reference", str(artifacts[1]),
"--results", str(tmp_path / "absent.h5"), "--report", str(report)],
check=False, capture_output=True, text=True)
assert completed.returncode == 2, completed.stderr
evidence = json.loads(report.read_text())
assert not evidence["precheck"]["passed"]
assert evidence["families"] == []
DRILLING_MODES = ("direct-initial-director-spring", "legacy-corotated-twist")
def drilling_result(tmp_path, artifacts, mode=DRILLING_MODES[0]):
path = changed_result(tmp_path, artifacts)
with h5py.File(path, "r+") as file:
file["metadata"].attrs["mitc4_drilling_formulation"] = mode
file["metadata"].attrs["drilling_diagnostics_version"] = np.uint64(1)
for frame in range(101):
root = file.require_group(f"steps/Step-1/frames/{frame}/global")
raw = np.array([1., 2., 3., 4., 5., 6.], dtype="<f8")
moment = np.array([4., 5., 6.], dtype="<f8")
corrected = raw.copy()
corrected[3:] -= moment
root.create_dataset("equilibrium", data=raw)
root.create_dataset("drilling_moment", data=moment)
root.create_dataset("equilibrium_with_drilling", data=corrected)
return path
@pytest.mark.parametrize("mode", DRILLING_MODES)
def test_drilling_diagnostics_are_precheck_and_report_all_frames(tmp_path, artifacts, mode):
result = comparator.compare(artifacts[0], artifacts[1], drilling_result(tmp_path, artifacts, mode))
assert result["precheck"]["passed"], result["precheck"]
assert result["verdict"] == "passed"
assert len(result["families"]) == 202
diagnostics = result["drilling_diagnostics"]
assert diagnostics["formulation"] == mode
assert diagnostics["version"] == 1
assert [frame["frame_index"] for frame in diagnostics["frames"]] == list(range(101))
assert diagnostics["frames"][0] == {
"storage_step": "Step-1", "frame_index": 0,
"raw_equilibrium": [1., 2., 3., 4., 5., 6.],
"drilling_moment": [4., 5., 6.],
"equilibrium_with_drilling": [1., 2., 3., 0., 0., 0.]}
def test_markerless_drilling_report_is_null(artifacts):
assert comparator.compare(*artifacts)["drilling_diagnostics"] is None
@pytest.mark.parametrize("mutation", ["mode_missing", "version_missing", "unknown_mode", "unknown_version",
"version_type", "mode_type", "initial_missing", "intermediate_missing", "both_missing",
"wrong_shape", "wrong_type", "big_endian", "nan", "inf", "wrong_algebra", "wrong_force"])
def test_drilling_diagnostics_rejected_before_displacement_tolerance(tmp_path, artifacts, mutation):
path = drilling_result(tmp_path, artifacts)
with h5py.File(path, "r+") as file:
metadata = file["metadata"]
root = file["steps/Step-1/frames/50/global"]
if mutation == "mode_missing":
del metadata.attrs["mitc4_drilling_formulation"]
elif mutation == "version_missing":
del metadata.attrs["drilling_diagnostics_version"]
elif mutation == "unknown_mode":
metadata.attrs["mitc4_drilling_formulation"] = "unknown"
elif mutation == "unknown_version":
metadata.attrs["drilling_diagnostics_version"] = np.uint64(2)
elif mutation == "version_type":
del metadata.attrs["drilling_diagnostics_version"]
metadata.attrs["drilling_diagnostics_version"] = np.uint32(1)
elif mutation == "mode_type":
del metadata.attrs["mitc4_drilling_formulation"]
metadata.attrs["mitc4_drilling_formulation"] = np.bytes_(DRILLING_MODES[0])
elif mutation == "initial_missing":
del file["steps/Step-1/frames/0/global/drilling_moment"]
elif mutation == "intermediate_missing":
del root["equilibrium_with_drilling"]
elif mutation == "both_missing":
del root["drilling_moment"]
del root["equilibrium_with_drilling"]
elif mutation in {"wrong_shape", "wrong_type", "big_endian"}:
del root["drilling_moment"]
root.create_dataset("drilling_moment", data=np.zeros(2 if mutation == "wrong_shape" else 3,
dtype=">f8" if mutation == "big_endian" else ("<f4" if mutation == "wrong_type" else "<f8")))
elif mutation == "nan":
root["drilling_moment"][0] = np.nan
elif mutation == "inf":
root["equilibrium"][0] = np.inf
else:
root["equilibrium_with_drilling"][0 if mutation == "wrong_force" else 3] = 1e-100
result = comparator.compare(artifacts[0], artifacts[1], path)
assert not result["precheck"]["passed"]
assert result["families"] == []