#!/usr/bin/env python3 """Audit frozen density-order controls and short coupled campaigns. Run only after the serial FEM workload. Saved polynomial algebra, sums, input binding and reported gates are checked; mapped weights and local M,b are not exported, so this is not an independent FEM projection/integration. No sampled result is promoted to a global positivity or converged-physics certificate. """ from __future__ import annotations import argparse from collections import defaultdict from decimal import Decimal, localcontext import json import math from pathlib import Path import tempfile import unittest import audit_negative_density_regions as common EPS = 2.220446049250313e-16 ORDERS = (14, 22, 30) CONTROLS = ("eos_nodal_interpolation", "eos_l2_q26") BOUNDS = common.NEGATIVE_BOUNDS SUMS = common.SUM_FIELDS + ( "negative_mass_uncertainty_scale", "first_moment_x", "first_moment_y", "first_moment_z", "inertia_z", "negative_enthalpy_samples", "resolved_negative_intervals", "uncertain_negative_intervals", "resolved_witness_negative_volume", "resolved_witness_negative_mass", "uncertain_negative_volume", "uncertain_negative_mass") EOS_SUMS = ("target_mass", "eos_squared_integral", "signed_eos_integral") require, close, rows, metadata, keyed, digest = ( common.require, common.close, common.rows, common.metadata, common.keyed, common.digest) class Campaign(common.Campaign): def run(self, name): if name not in self.runs: r = json.loads(self.record(self.path / f"{name}_run.json").read_text()) require(r["id"] == name and r["exit_code"] in (0, 3) and r["inputs_unchanged"] is True, f"Incomplete/changed run {name}") expected = [str(self.path / "plateau_diagnostic"), *self.cases[name]["arguments"], "--mesh", self.provenance["mesh"], "--output", str(self.path / name)] require(r["command"] == expected, f"Manifest/command mismatch: {name}") for path, value in r["input_state_hashes"].items(): require(digest(self.record(Path(path))) == value, f"Changed input: {path}") self.runs[name] = r return self.runs[name] def vector(path): with path.open() as stream: size = int(next(stream)) result = [float(line) for line in stream] require(size == len(result) and size > 0 and all(map(math.isfinite, result)), f"Invalid vector {path}") return result def numeric(row, omit=("control",)): return {k: float(v) for k, v in row.items() if k not in omit} def enthalpy_observations(records, key): values = [row[key] for row in records] require(all(map(math.isfinite, values)), "Nonfinite saved enthalpy observation") return {"samples": len(values), "minimum_enthalpy": min(values) if values else None, "negative_enthalpy_samples": sum(value < 0 for value in values)} def poly_value(coefficients, x): value = coefficients[-1] for coefficient in reversed(coefficients[:-1]): value = value * x + coefficient return value def decimal_roots(coefficients): """Bracket derivative-separated real roots using 60-digit arithmetic. Unlike the producer's compensated binary arithmetic, the saved decimal coefficients define this independent polynomial. Exact repeated roots may move under serialization; shallow intervals retain their uncertainty label. """ c = list(coefficients) while len(c) > 1 and c[-1] == 0: c.pop() if len(c) == 1: return [] if len(c) == 2: root = -c[0] / c[1] return [root] if 0 <= root <= 1 else [] derivative = [i * c[i] for i in range(1, len(c))] cuts = sorted(set([Decimal(0), Decimal(1), *decimal_roots(derivative)])) result = [point for point in cuts if poly_value(c, point) == 0] for left, right in zip(cuts, cuts[1:]): fl, fr = poly_value(c, left), poly_value(c, right) if not fl or not fr or (fl < 0) == (fr < 0): continue for _ in range(100): middle = (left + right) / 2 fm = poly_value(c, middle) if not fm: left = right = middle break if (fm < 0) == (fl < 0): left, fl = middle, fm else: right = middle result.append((left + right) / 2) return sorted(set(result)) def independent_partition(coefficient_strings): with localcontext() as context: context.prec = 60 c = [Decimal(v) for v in coefficient_strings] cuts = sorted(set([Decimal(0), Decimal(1), *decimal_roots(c)])) # A repeated exact root can acquire a pocket at the last Decimal digit # after the derivative root is rounded. It has no distinct binary64 # endpoints and cannot correspond to an exported interval. negative = [(float(a), float(b)) for a, b in zip(cuts, cuts[1:]) if float(a) < float(b) and poly_value(c, (a + b) / 2) < 0] critical = decimal_roots([i * c[i] for i in range(1, len(c))]) return negative, [float(x) for x in critical] def interval_line(line, gauss_nodes): first = line[0] coefficient_keys = [f"c{i}" for i in range(5)] c = [float(first[k]) for k in coefficient_keys] degree = int(first["degree"]) require(degree in (2, 3, 4) and all(c[i] == 0 for i in range(degree + 1, 5)), "Wrong polynomial degree") require(all(all(row[k] == first[k] for k in coefficient_keys) for row in line), "Inconsistent line polynomial") scale = math.fsum(map(abs, c)) tolerance = 2048 * EPS * scale expected, critical = independent_partition([first[k] for k in coefficient_keys]) parsed = sorted([numeric(row) for row in line], key=lambda r: r["t_left"]) for index, r in enumerate(parsed): common.finite(r) require(r["axis"] in (0, 1, 2) and r["orientation"] in (-1, 1), "Invalid interval axis") for coordinate in (r["u"], r["v"]): require(min(abs(coordinate - node) for node in gauss_nodes) < 2e-14, "Interval is not at its Gauss node") a, b = r["t_left"], r["t_right"] require(0 <= a < b <= 1 and (index == 0 or parsed[index - 1]["t_right"] <= a), "Overlapping/outside interval") close(r["reconstruction_tolerance"], tolerance, "Coefficient-scaled interval tolerance", atol=1e-30) require(r["negative_volume"] > 0 and r["negative_mass"] >= 0, "Invalid negative measure") require(r["resolved_negative_witness"] in (0, 1) and bool(r["resolved_negative_witness"]) == (r["minimum_density"] < -tolerance), "Witness classification changed") samples = [a, b, (a + b) / 2, a + .25 * (b - a), a + .75 * (b - a)] samples += [x for x in critical if a < x < b] # Evaluate the decimal polynomial at the saved endpoints/critical points; # ordinary binary Horner could invent signs in the shallow Q4 tails. with localcontext() as context: context.prec = 60 exact = [Decimal(first[k]) for k in coefficient_keys] values = [float(poly_value(exact, Decimal.from_float(x))) for x in samples] allowance = tolerance + 128 * EPS * scale require(max(values) <= allowance, "Positive interior of alleged negative interval") require(abs(r["minimum_density"] - min(values)) <= allowance, "Polynomial/FE minimum disagrees") for endpoint in (a, b): if endpoint not in (0., 1.): require(abs(float(poly_value(exact, Decimal.from_float(endpoint)))) <= allowance, "Interval endpoint is not a root") require(r["negative_mass"] / r["negative_volume"] <= -min(values) + allowance, "Negative mass exceeds polynomial minimum times volume") # Only lines whose every exported pocket has a resolved FE witness demand # one-to-one root matching. Rounded near-tangent tails are retained, not # silently recast as a sign certificate or discarded from the integrals. resolved = all(r["resolved_negative_witness"] for r in parsed) if resolved: require(len(parsed) == len(expected), "Missing/extra resolved negative interval") for r, (left, right) in zip(parsed, expected): for got, wanted in ((r["t_left"], left), (r["t_right"], right)): if wanted in (0., 1.): require(got == wanted, "Missing cell endpoint") else: derivative = abs(math.fsum(i * c[i] * wanted ** (i - 1) for i in range(1, 5))) # Coefficient-scaled function uncertainty divided by the # root slope. Saturating the slope prevents an unbounded # location allowance at a multiple root. allowance = max(128 * EPS, tolerance / max(derivative, math.sqrt(tolerance * scale))) close(got, wanted, "Independent polynomial root", rtol=0, atol=allowance) return parsed, resolved def measure_check(r): allowed = list(BOUNDS) if r["negative_intervals"] == 0 else [] if r["negative_enthalpy_samples"]: allowed += EOS_SUMS if r["mass"] == 0: allowed += ["negative_mass_fraction"] common.finite(r, allowed) require(r["degree"] in (2, 3, 4) and r["interval_order"] == 14 and r["tangential_order"] in ORDERS, "Unexpected density/quadrature degree") require(r["volume"] > 0 and 0 <= r["negative_volume"] <= r["volume"] and r["negative_mass"] >= 0, "Invalid region measure") close(r["negative_volume_fraction"], r["negative_volume"] / r["volume"], "Negative volume fraction") if r["mass"]: close(r["negative_mass_fraction"], r["negative_mass"] / abs(r["mass"]), "Negative mass fraction") require(r["maximum_reconstruction_error"] <= r["maximum_reconstruction_tolerance"], "Reconstruction gate failed") require(r["eos_defined_everywhere_sampled"] == (r["negative_enthalpy_samples"] == 0), "EOS definition flag") for key in ("negative_intervals", "resolved_negative_intervals", "uncertain_negative_intervals", "negative_enthalpy_samples", "tangent_lines", "integration_points"): require(r[key] >= 0 and r[key].is_integer(), f"Invalid count {key}") require(r["resolved_negative_intervals"] + r["uncertain_negative_intervals"] == r["negative_intervals"], "Negative classification counts") for quantity in ("mass", "volume"): close(r[f"resolved_witness_negative_{quantity}"] + r[f"uncertain_negative_{quantity}"], r[f"negative_{quantity}"], "Resolved/uncertain measure decomposition") if r["negative_intervals"] == 0: require(r["negative_volume"] == r["negative_mass"] == 0 and all(math.isnan(r[k]) for k in BOUNDS), "Empty-negative-region convention") else: require(0 <= r["negative_interval_t_minimum"] < r["negative_interval_t_maximum"] <= 1 and 0 <= r["negative_radius_minimum"] <= r["negative_radius_maximum"], "Negative region bounds") if not r["negative_enthalpy_samples"]: require(r["eos_squared_integral"] >= 0 and r["target_mass"] >= 0, "Invalid EOS integral") close(r["signed_eos_integral"], r["mass"] - r["target_mass"], "Signed EOS integral", rtol=2e-10, atol=256 * EPS * max(abs(r["mass"]), r["target_mass"])) else: require(all(math.isnan(r[k]) for k in EOS_SUMS), "Undefined EOS integral fabricated") def audit_measures(get, expected_controls, n, K): summary = keyed(rows(get("density_order_summary.csv")), lambda r: (r["control"], int(r["degree"]), int(r["tangential_order"]))) require(set(summary) == {(c, d, q) for c, d in expected_controls for q in ORDERS}, "Missing order/control summary") elements = keyed(rows(get("density_order_elements.csv")), lambda r: (r["control"], int(r["degree"]), int(r["tangential_order"]), int(r["element"]))) support = {element for _, _, _, element in elements} require(set(elements) == {(*key, e) for key in summary for e in support}, "Element support differs between controls") totals = {key: numeric(row) for key, row in summary.items()} elemental = {key: numeric(row) for key, row in elements.items()} for key, total in totals.items(): measure_check(total) subset = [elemental[(*key, e)] for e in support] for r in subset: measure_check(r) require(r["tangent_lines"] == (key[2] // 2 + 1) ** 2 and r["witness_element"] == r["element"] and r["axis"] in (0, 1, 2) and r["orientation"] in (-1, 1), "Element axis/count mismatch") require(r["integration_points"] % 8 == 0 and 8 * r["tangent_lines"] <= r["integration_points"] <= 8 * (key[1] + 1) * r["tangent_lines"], "Interval quadrature count mismatch") for field in SUMS: scale = math.fsum(abs(r[field]) for r in subset) close(math.fsum(r[field] for r in subset), total[field], "Element/global " + field, atol=128 * EPS * max(scale, 1e-25)) if total["eos_defined_everywhere_sampled"]: for field in EOS_SUMS: close(math.fsum(r[field] for r in subset), total[field], "Element/global " + field, atol=128 * EPS * math.fsum(abs(r[field]) for r in subset)) for field in ("minimum_density", "minimum_enthalpy"): close(min(r[field] for r in subset), total[field], "Element/global minimum") witness = elemental[(*key, int(total["witness_element"]))] for field in ("minimum_density", "witness_h", "witness_t", "witness_u", "witness_v", "witness_x", "witness_y", "witness_z"): close(witness[field], total[field], "Global witness binding", atol=1e-15) nonempty = [r for r in subset if r["negative_intervals"]] if nonempty: for field, operation in zip(BOUNDS, (min, max, min, max)): close(operation(r[field] for r in nonempty), total[field], "Global negative bound") for field in ("maximum_reconstruction_error", "maximum_reconstruction_tolerance"): close(max(r[field] for r in subset), total[field], "Global reconstruction maximum") lines, pieces = defaultdict(list), defaultdict(list) for row in rows(get("density_order_intervals.csv")): key = row["control"], int(row["degree"]), int(row["element"]) require((*key[:2], 30, key[2]) in elements, "Interval unsupported element/control") lines[(*key, row["u"], row["v"])].append(row) resolved_lines = 0 gauss_nodes = [x for x, _ in common.gauss_rule(16)] for key, line in lines.items(): parsed, resolved = interval_line(line, gauss_nodes) resolved_lines += resolved element = elemental[(*key[:2], 30, key[2])] require(all(r["axis"] == element["axis"] and r["orientation"] == element["orientation"] for r in parsed), "Interval geometry differs from element") pieces[key[:3]].extend(parsed) for control, degree in expected_controls: for element in support: saved = pieces[(control, degree, element)] total = elemental[(control, degree, 30, element)] require(len(saved) == total["negative_intervals"], "Finest interval coverage") for field in ("negative_mass", "negative_volume"): close(math.fsum(r[field] for r in saved), total[field], "Interval/element " + field) for resolved, prefix in ((True, "resolved_witness"), (False, "uncertain")): subset = [r for r in saved if bool(r["resolved_negative_witness"]) == resolved] count_field = "resolved_negative_intervals" if resolved else "uncertain_negative_intervals" require(len(subset) == total[count_field], "Interval resolved/uncertain count") for quantity in ("mass", "volume"): close(math.fsum(r[f"negative_{quantity}"] for r in subset), total[f"{prefix}_negative_{quantity}"], "Interval classification sum") if saved: close(min(r["t_left"] for r in saved), total["negative_interval_t_minimum"], "Finest left bound") close(max(r["t_right"] for r in saved), total["negative_interval_t_maximum"], "Finest right bound") surface = defaultdict(list) for row in rows(get("density_order_surface.csv")): key = row["control"], int(row["degree"]) require(key in expected_controls and int(row["element"]) in support, "Surface unsupported control") r = numeric(row) common.finite(r, ("rho_eos", "error") if r["h"] < 0 else ()) require(r["u"] in (0, .25, .5, .75, 1) and r["v"] in (0, .25, .5, .75, 1), "Surface sampling changed") require(r["eos_defined"] == (r["h"] >= 0) and r["resolved_negative"] == (r["rho"] < -r["reconstruction_tolerance"]), "Surface classification contradiction") if r["h"] >= 0: close(r["rho_eos"], (r["h"] / ((n + 1) * K)) ** n, "Surface EOS", atol=1e-30) close(r["error"], r["rho"] - r["rho_eos"], "Surface EOS difference", atol=1e-30) else: require(math.isnan(r["rho_eos"]) and math.isnan(r["error"]), "Undefined surface EOS fabricated") surface[key].append(r) result = {} for control, degree in sorted(expected_controls): key = control, degree traces = surface[key] require(traces, "Missing surface samples") unique = {(r["face"], r["u"], r["v"]) for r in traces} require(len(unique) == len(traces) == 25 * len({r["face"] for r in traces}), "Surface coverage") witnesses = [row for element in support for row in pieces[(control, degree, element)]] result[f"Q{degree}/{control}"] = { "summary": {str(q): totals[(control, degree, q)] for q in ORDERS}, "summary_enthalpy_scope": "Integration nodes only; EOS availability flags do not cover boundary samples or density-minimum witnesses.", "surface": {"samples": len(traces), "minimum_density": min(r["rho"] for r in traces), "negative_samples": sum(r["rho"] < 0 for r in traces), "resolved_negative_samples": sum(r["resolved_negative"] for r in traces), "enthalpy": {**enthalpy_observations(traces, "h"), "scope": "Saved one-sided 5-by-5 face samples; counts are not volume or area fractions."}}, "negative_interval_witness_enthalpy": {**enthalpy_observations(witnesses, "witness_h"), "scope": "One saved density-minimum witness per exported negative interval on order30 tangential lines; this is not an enthalpy-extremum search."}, "negative_mass_q30_minus_q22": totals[(control, degree, 30)]["negative_mass"] - totals[(control, degree, 22)]["negative_mass"], "negative_volume_q30_minus_q22": totals[(control, degree, 30)]["negative_volume"] - totals[(control, degree, 22)]["negative_volume"]} return result, support, totals, {"negative_lines": len(lines), "fully_resolved_negative_lines": resolved_lines} def audit_selfchecks(table): expected = {"positive_constant": (0., 0.), "negative_constant": (1., 2.), "linear_crossing": (.25, .03125), "quadratic_interior_pocket": (math.sqrt(6) / 15, math.sqrt(6) / 2250), "cubic_interior_pocket": (.5, 1 / 96), "quartic_two_pockets": (.5, 19 / 15360), "quartic_double_roots_nonnegative": (0., 0.), **{f"reconstruct_t_power_{d}": (0., 0.) for d in (2, 3, 4)}} table = keyed(table, lambda row: row["case"]) require(table.keys() == expected.keys(), "Self-check coverage") for name, (length, mass) in expected.items(): r = numeric(table[name], ("case",)) common.finite(r) require(r["passed"] == 1, "Producer root self-check failed") for key, value in (("negative_length", length), ("expected_negative_length", length), ("negative_mass", mass), ("expected_negative_mass", mass)): close(r[key], value, "Independent analytic self-check " + name, atol=2e-18) def audit_projection(get, support, n, K, totals): elements = keyed(rows(get("density_order_projection_elements.csv")), lambda r: (int(r["degree"]), int(r["element"]))) require(set(elements) == {(d, e) for d in (2, 3, 4) for e in support}, "Projection element coverage") by_element = defaultdict(list) vectors = {(d, c): vector(get(f"density_order_Q{d}_{c}_coefficients.txt")) for d in (2, 3, 4) for c in CONTROLS} seen = defaultdict(set) degree_nodes = {degree: [x for x, _ in common.gauss_rule(degree + 1)] for degree in (2, 3, 4)} for row in rows(get("density_order_nodes.csv")): r = numeric(row, ()) common.finite(r) degree, element, local, glob = (int(r[k]) for k in ("degree", "element", "local_dof", "global_dof")) require((degree, element) in elements and r["h"] >= 0, "Invalid nodal projection support/enthalpy") require(0 <= local < (degree + 1) ** 3 and glob not in seen[degree], "Duplicate/invalid DG dof") require(glob == element * (degree + 1) ** 3 + local, "Serial DG element ordering changed") seen[degree].add(glob) nodes = degree_nodes[degree] for coordinate in (r["xi"], r["eta"], r["zeta"]): require(min(abs(coordinate - x) for x in nodes) < 2e-14, "Wrong interpolation nodes") close(r["target_rho"], (r["h"] / ((n + 1) * K)) ** n, "Nodal EOS target", atol=1e-30) require(r["interpolation_coefficient"] == r["target_rho"], "Interpolation is not EOS at the node") for control, field in zip(CONTROLS, ("interpolation_coefficient", "projection_coefficient")): require(0 <= glob < len(vectors[(degree, control)]) and vectors[(degree, control)][glob] == r[field], "Node/coefficient vector disagreement") by_element[(degree, element)].append(r) result = {} for (degree, element), row in elements.items(): r = numeric(row, ("basis_family",)) common.finite(r) node_rows = by_element[(degree, element)] require(len(node_rows) == r["dofs"] == (degree + 1) ** 3 and {int(v["local_dof"]) for v in node_rows} == set(range((degree + 1) ** 3)), "Nodal element coverage") require(len({(v["xi"], v["eta"], v["zeta"]) for v in node_rows}) == len(node_rows), "Duplicate tensor nodes") require(r["quadrature_request"] == 26 and r["quadrature_achieved"] == 27 and r["quadrature_points"] == 14 ** 3 and "L2" in row["basis_family"], "Projection rule/basis changed") require(r["rhs_norm"] > 0 and 0 <= r["normal_defect_norm"] and r["normal_relative_defect"] < 1e-10, "Reported local normal solve gate failed") close(r["normal_relative_defect"], r["normal_defect_norm"] / r["rhs_norm"], "Normal defect ratio") require(0 <= r["constant_coefficient_error"] < 1e-10 and 0 <= r["constant_sample_error"] < 1e-10, "Constant reconstruction failed") close(r["projection_mass"], r["target_mass"], "Local projection mass", rtol=2e-11) require(r["projection_eos_squared"] >= 0 and r["interpolation_eos_squared"] >= 0 and r["projection_eos_squared"] <= r["interpolation_eos_squared"] + 512 * EPS * r["target_mass"] ** 2 / r["volume"], "Reported orthogonal projection exceeds interpolation error") for field, operation, source in (("interpolation_coefficient_minimum", min, "interpolation_coefficient"), ("projection_coefficient_minimum", min, "projection_coefficient"), ("projection_coefficient_maximum", max, "projection_coefficient")): require(r[field] == operation(v[source] for v in node_rows), "Projection coefficient extrema") for degree in (2, 3, 4): data = [numeric(row, ("basis_family",)) for (d, _), row in elements.items() if d == degree] totals26 = {k: math.fsum(r[k] for r in data) for k in ("volume", "target_mass", "interpolation_mass", "projection_mass", "interpolation_eos_squared", "projection_eos_squared")} for control, key in zip(CONTROLS, ("interpolation_mass", "projection_mass")): fine = totals[(control, degree, 30)] close(totals26[key], fine["mass"], "Q26/sliced mass", rtol=2e-9, atol=1e-13) close(totals26["volume"], fine["volume"], "Q26/sliced volume", rtol=2e-9, atol=1e-13) for control in CONTROLS: require(all(value == 0 for i, value in enumerate(vectors[(degree, control)]) if i not in seen[degree]), "Control density outside stellar support is nonzero") result[f"Q{degree}"] = {**totals26, "maximum_reported_relative_normal_defect": max(r["normal_relative_defect"] for r in data), "active_density_dofs": len(seen[degree]), "full_density_storage": len(vectors[(degree, CONTROLS[0])])} return result def audit_error_quadrature(get, controls, support, totals, frozen): table = keyed(rows(get("density_order_error_quadrature.csv")), lambda r: (r["control"], int(r["degree"]), int(r["quadrature_request"]), int(r["element"]))) require(set(table) == {(c, d, q, e) for c, d in controls for q in (26, 30) for e in support}, "Independent error-rule coverage") projection = keyed(rows(get("density_order_projection_elements.csv")), lambda r: (int(r["degree"]), int(r["element"]))) if frozen else {} results = {} for control, degree in sorted(controls): result = {} for order in (26, 30): records = [] for element in support: r = numeric(table[(control, degree, order, element)]) common.finite(r, EOS_SUMS if r["negative_enthalpy_samples"] else ()) require(r["quadrature_achieved"] == order + 1 and r["quadrature_points"] == (order // 2 + 1) ** 3, "Independent error quadrature changed") require(r["volume"] > 0 and r["negative_enthalpy_samples"] >= 0 and r["eos_defined_everywhere_sampled"] == (r["negative_enthalpy_samples"] == 0), "Independent error validity") if r["negative_enthalpy_samples"]: require(all(math.isnan(r[k]) for k in EOS_SUMS), "Undefined independent EOS integrals fabricated") else: require(r["target_mass"] >= 0 and r["eos_squared_integral"] >= 0, "Invalid independent EOS integrals") close(r["signed_eos_integral"], r["mass"] - r["target_mass"], "Independent signed EOS error", atol=256 * EPS * max(abs(r["mass"]), r["target_mass"])) if frozen and order == 26: p = numeric(projection[(degree, element)], ("basis_family",)) label = "interpolation" if control == CONTROLS[0] else "projection" for field, source in (("volume", "volume"), ("target_mass", "target_mass"), ("mass", label + "_mass"), ("eos_squared_integral", label + "_eos_squared")): require(r[field] == p[source], "Q26 error/projection assembly binding") records.append(r) valid = all(r["eos_defined_everywhere_sampled"] for r in records) global_values = {key: math.fsum(r[key] for r in records) for key in ("volume", "mass", "negative_enthalpy_samples")} global_values.update({key: math.fsum(r[key] for r in records) if valid else math.nan for key in EOS_SUMS}) global_values["eos_defined_everywhere_sampled"] = valid global_values["minimum_density"] = min(r["minimum_density"] for r in records) global_values["minimum_enthalpy"] = min(r["minimum_enthalpy"] for r in records) global_values["eos_rms"] = math.sqrt(global_values["eos_squared_integral"] / global_values["volume"]) if valid else math.nan result[str(order)] = global_values for key in ("mass", "volume"): close(global_values[key], totals[(control, degree, 30)][key], "Unsliced/sliced " + key, rtol=2e-9, atol=1e-13) result["eos_squared_q30_minus_q26"] = result["30"]["eos_squared_integral"] - result["26"]["eos_squared_integral"] result["sliced_eos_squared_minus_q30"] = totals[(control, degree, 30)]["eos_squared_integral"] - result["30"]["eos_squared_integral"] results[f"Q{degree}/{control}"] = result return results def bind_density(state, full_density, block, degree, support): """Bind the supported serial, discontinuous tensor element ordering. Frozen nodal exports separately validate this global/local indexing rule. No general parallel/continuous-field gather convention is assumed here. """ width = (degree + 1) ** 3 require(len(full_density) % width == 0 and max(support) < len(full_density) // width, "Full DG density dimensions") active = [element * width + j for element in sorted(support) for j in range(width)] selected = set(active) require(all(v == 0 for i, v in enumerate(full_density) if i not in selected), "Density outside support nonzero") offset, size = int(block["offset"]), int(block["size"]) require(size == len(active) and state[offset:offset + size] == [full_density[i] for i in active], "Accepted/observed density does not bind to its full DG storage") def physical(get, filename, meta): table = keyed(rows(get(filename)), lambda r: (int(r["quadrature"]), r["metric"])) result = {} for q in sorted({q for q, _ in table}): m = {name: float(row["value"]) for (order, name), row in table.items() if order == q} require(not any(math.isinf(v) for v in m.values()), "Infinite physical observation") for key in ("mass", "volume", "binding_energy", "kinetic_energy", "minimum_density", "minimum_enthalpy"): require(math.isfinite(m[key]), "Nonfinite physical " + key) require(m["complete_volume_measurement"] == 1, "Incomplete physical volume observation") close(m["polytropic_index"], float(meta["n"]), "Physical index") close(m["polytropic_constant"], float(meta["K"]), "Physical K") passes = {} for prefix, pressure, defined_key in (("", "pressure_integral", "pressure_from_density_defined_everywhere"), ("enthalpy_", "enthalpy_pressure_integral", "pressure_from_enthalpy_defined_everywhere")): require(m[defined_key] in (0, 1), "Invalid pressure definition flag") for force, energy in (("", "binding_energy"), ("force_", "force_binding_energy")): key = prefix + force + "virial_" if m[defined_key]: terms = 2 * m["kinetic_energy"], m[energy], 3 * m[pressure] expected = math.fsum(terms) / abs(m["binding_energy"]) close(m[key + "signed"], expected, "Virial scalar identity", atol=2e-14) close(m[key + "error"], abs(m[key + "signed"]), "Virial absolute identity") passes[key + "within_1e5"] = m[key + "error"] <= 1e-5 else: require(math.isnan(m[key + "signed"]) and math.isnan(m[key + "error"]), "Undefined virial fabricated") passes[key + "within_1e5"] = None if "mass_absolute_error" in m: close(m["mass_absolute_error"], abs(m["mass"] - float(meta["M"])), "Physical mass error", atol=1e-15) close(m["angular_momentum_absolute_error"], abs(m["angular_momentum_z"] - float(meta["J"])), "Physical J error", atol=1e-15) result[str(q)] = {"metrics": m, "virial_milestones": passes, "undefined_metrics": [k for k, v in m.items() if math.isnan(v)]} return result def audit_closure_rules(table, support): rules = keyed(table, lambda row: int(row["element"])) require(set(rules) == support, "Closure quadrature support differs from the stellar elements") counts = defaultdict(int) for row in rules.values(): require(set(row) == {"element", "quadrature_order", "quadrature_points"}, "Closure rule schema changed") achieved, points = int(row["quadrature_order"]), int(row["quadrature_points"]) require(achieved > 0 and points > 0 and achieved % 2 == 1, "Invalid achieved tensor Gauss rule") per_axis = (achieved + 1) // 2 require(points == per_axis ** 3, "Achieved order and tensor Gauss point count disagree") counts[(achieved, points)] += 1 return {"scope": "Actual IntegrationRule achieved order and point count; the requested policy order is not exported.", "stellar_elements": len(rules), "distinct_achieved_orders": sorted({order for order, _ in counts}), "distinct_point_counts": sorted({points for _, points in counts}), "rules": [{"achieved_order": order, "quadrature_points": points, "points_per_axis": (order + 1) // 2, "elements": count} for (order, points), count in sorted(counts.items())]} def audit_line_search(table, history, initial_residual, strategy): require(strategy == "baseline", "Trial audit currently covers the recorded baseline policy only") # Reviewed against LineSearchTrialDisposition in solver/newton.cppm and # PreparationStageName/line-search dispatch in solver/stellar_context.cppm. dispositions = {0: "accepted", 1: "inadmissible_state", 2: "non_finite_state", 3: "non_finite_residual", 4: "insufficient_decrease", 5: "auxiliary_solve_failure"} stages = {"generated geometry", "gravity", "barotropic closure", "hydrostatic equilibrium", "displacement residual", "pressure force", "gravity displacement force", "rotational displacement force", "displacement composition", "mass normalization", "model specification", "unspecified preparation", "unknown preparation"} grouped, source_counts, disposition_counts, joint_counts = defaultdict(list), defaultdict(int), defaultdict(int), defaultdict(int) saved = [] for raw in table: r = {**numeric(raw, ("source",)), "source": raw["source"]} iteration, trial, code = int(raw["iteration"]), int(raw["trial"]), int(raw["disposition"]) require(iteration >= 0 and trial >= 0 and code in dispositions, "Invalid trial index/disposition") require(math.isfinite(r["alpha"]) and 0 < r["alpha"] <= 1, "Invalid trial step length") require(not saved or (iteration, trial) > (int(saved[-1]["iteration"]), int(saved[-1]["trial"])), "Trial log is not sequential") require(trial == len(grouped[iteration]), "Missing/repeated trial within an iteration") if grouped[iteration]: require(r["alpha"] < grouped[iteration][-1]["alpha"], "Backtracking trial did not shorten the step") require(grouped[iteration][-1]["disposition"] != 0, "Trial recorded after acceptance") source = r["source"] if code == 0: require(source == "", "Accepted trial has a rejection source") elif code == 4: require(source == "residual metric", "Insufficient decrease source mismatch") elif code == 1: require(source in stages or source == "enthalpy polynomial", "Inadmissible-state source mismatch") elif code == 2: require(source in stages or source == "candidate state", "Nonfinite-state source mismatch") elif code == 3: require(source in stages or source in ("normalized residual", "residual metric"), "Nonfinite-residual source mismatch") else: require(False, "Baseline monolithic trial unexpectedly reports an auxiliary solve failure") if code in (0, 4): require(math.isfinite(r["F"]) and r["F"] >= 0 and math.isfinite(r["min_jacobian"]) and r["min_jacobian"] > 0, "Metric-evaluated trial lacks finite residual/valid geometry") elif source != "residual metric": require(math.isnan(r["F"]), "Rejected preparation has a fabricated residual metric") require(math.isfinite(r["min_jacobian"]) or math.isnan(r["min_jacobian"]), "Invalid optional trial geometry") grouped[iteration].append(r) saved.append(r) source_counts[source] += 1 disposition_counts[dispositions[code]] += 1 joint_counts[(source, code)] += 1 history_by_iteration = keyed(history, lambda row: int(row["iteration"])) require(set(grouped) <= set(history_by_iteration), "Trials lack an after-iteration record") previous = initial_residual accepted_alphas = [] for iteration, row in sorted(history_by_iteration.items()): trials = grouped[iteration] require(len(trials) == row["trials"], "Trial/iteration count mismatch") accepted = [r for r in trials if r["disposition"] == 0] require(len(accepted) == int(row["accepted"]), "Trial/iteration acceptance mismatch") for trial in trials: if trial["disposition"] in (0, 4): limit = previous * math.sqrt(1 - 1e-4 * trial["alpha"]) allowance = 128 * EPS * max(previous, trial["F"]) if trial["disposition"] == 0: require(trial["F"] <= limit + allowance, "Accepted baseline trial fails recorded merit inequality") else: require(trial["F"] >= limit - allowance, "Rejected baseline metric satisfies decrease away from rounding") if accepted: require(accepted[0] is trials[-1], "Accepted trial was not the last attempt") close(accepted[0]["alpha"], row["alpha"], "Accepted alpha binding") close(accepted[0]["F"], row["F"], "Accepted residual binding") accepted_alphas.append(row["alpha"]) else: require(row["alpha"] == 0, "Failed iteration has a nonzero accepted step") close(row["F"], previous, "Failed iteration residual restoration") previous = row["F"] return {"scope": "Recorded baseline trial stages and metric inequalities; rejected preparation has no evaluated residual metric. Geometry may shorten the backtracking ladder.", "disposition_enum": {str(code): name for code, name in dispositions.items()}, "total_trials": len(saved), "counts_by_source": dict(sorted(source_counts.items())), "counts_by_disposition": dict(sorted(disposition_counts.items())), "source_disposition_counts": [{"source": source, "disposition": dispositions[code], "disposition_code": code, "count": count} for (source, code), count in sorted(joint_counts.items())], "pressure_force_rejections": sum(r["source"] == "pressure force" and r["disposition"] != 0 for r in saved), "residual_metric_rejections": sum(r["source"] == "residual metric" and r["disposition"] != 0 for r in saved), "residual_metric_insufficient_decrease_rejections": sum(r["source"] == "residual metric" and r["disposition"] == 4 for r in saved), "minimum_accepted_alpha": min(accepted_alphas) if accepted_alphas else None, "trials": saved} def audit_case(campaign, name): get = lambda filename: campaign.artifact(name, filename) record = campaign.run(name) meta = metadata(get("metadata.txt")) scope = metadata(get("density_order_scope.txt")) require(scope["source_fields_bitwise_identical"] == "1" and scope["trial_preparations"] == "0" and scope["mpi_ranks"] == "1", "Observer changed source fields/state") n, K = float(meta["n"]), float(meta["K"]) for key in ("n", "K"): close(float(scope[key]), float(meta[key]), "Scope model") degree = int(meta["density_order"]) require(degree in (2, 3, 4) and degree == int(scope["source_density_degree"]), "Source degree binding") for key, expected in (("enthalpy_order", 3), ("gravity_potential_order", 2), ("gravity_flux_order", 2), ("displacement_order", 3)): require(int(meta[key]) == expected, "Non-density field order changed: " + key) settings = campaign.provenance["build_settings"] require(int(settings["MEAN_FIELD_DENSITY_POLYNOMIAL_ORDER:STRING"]) == degree and int(settings["MEAN_FIELD_UNIFORM_POLYNOMIAL_ORDER_INCREMENT:STRING"]) == 0, "Build/order metadata contradiction") frozen = scope["projection_controls"] == "1" require(meta["density_order_projection"] == str(int(frozen)) and meta["density_order_final"] == str(int(not frozen)), "Control mode mismatch") require(meta["surface_prescription"] == "harmonics" and int(meta["surface_degree"]) == 4 and int(meta["surface_maximum_order"]) == -1, "Surface prescription changed in density-only study") audit_selfchecks(rows(get("density_order_root_selfchecks.csv"))) controls = {(c, d) for c in CONTROLS for d in (2, 3, 4)} if frozen else {("stored_density", degree)} measures, support, totals, intervals = audit_measures(get, controls, n, K) closure_rule_file = "production_closure_rules.csv" if frozen else "final_closure_rules.csv" closure_rules = audit_closure_rules(rows(get(closure_rule_file)), support) closure_rules["source_artifact"] = closure_rule_file error_quadrature = audit_error_quadrature(get, controls, support, totals, frozen) timing_rows = keyed(rows(get("density_order_timings.csv")), lambda r: r["phase"]) require(set(timing_rows) == {"setup_and_projection", "independent_error_quadrature", "region_observation_and_output", "observer_total"}, "Observer timing scope changed") timing = {phase: float(row["seconds"]) for phase, row in timing_rows.items()} require(all(math.isfinite(v) and v >= 0 for v in timing.values()), "Invalid observer timing") close(math.fsum(value for phase, value in timing.items() if phase != "observer_total"), timing["observer_total"], "Observer phase-time partition", atol=1e-8) layout = rows(get("layout.csv")) blocks = keyed(layout, lambda r: (r["kind"], r["block"])) state = vector(get("observed_state.txt" if frozen else "accepted_state.txt")) for kind in ("value", "residual"): ordered = sorted([r for r in layout if r["kind"] == kind], key=lambda r: int(r["offset"])) end = 0 for r in ordered: require(int(r["offset"]) == end and int(r["size"]) > 0 and math.isfinite(float(r["scale"])), "Invalid block layout") end += int(r["size"]) require(end == len(state), "Layout/state dimensions") require(int(blocks[("value", "density")]["size"]) == len(support) * (degree + 1) ** 3 and int(blocks[("residual", "barotropic_closure")]["size"]) == len(support) * (degree + 1) ** 3, "Active density/closure dof count") result = {"artifact_passed": True, "kind": "frozen" if frozen else "coupled", "n": n, "K": K, "J": float(meta["J"]), "density_degree": degree, "state_size": len(state), "stellar_elements": len(support), "block_sizes": {f"{kind}/{block}": int(row["size"]) for (kind, block), row in blocks.items()}, "measures": measures, "interval_audit": intervals, "unsliced_error_quadrature": error_quadrature, "closure_quadrature": closure_rules, "observer_timings": timing, "run": {k: record[k] for k in ("exit_code", "seconds", "utc_elapsed_seconds", "utc_minus_monotonic_seconds", "child_cpu_seconds")}} if frozen: require(record["exit_code"] == 0, "Frozen observer did not complete") invariant = {r["metric"]: float(r["value"]) for r in rows(get("density_invariants.csv"))} common.finite(invariant) require(all(invariant[k] == 1 for k in ("accepted_physical_unchanged", "accepted_normalized_unchanged", "accepted_residual_unchanged")), "Frozen accepted-state invariant") result["invariants"] = invariant result["projection"] = audit_projection(get, support, n, K, totals) result["physical"] = physical(get, "density_physical.csv", meta) # The runner binds all replay metadata. Add an exact state-content test, # which also distinguishes a replay from the observer's untouched seed. inputs = record["input_state_hashes"] previous = [Path(path) for path in inputs if Path(path).name == "accepted_state.txt"] if previous: require(len(previous) == 1 and state == vector(previous[0]), "Frozen snapshot changed") density = vector(get("observed_density_coefficients.txt")) bind_density(state, density, blocks[("value", "density")], degree, support) else: result["initial_physical"] = physical(get, "initial_physical.csv", meta) result["physical"] = physical(get, "final_physical.csv", meta) density = vector(get("density_order_stored_density_coefficients.txt")) bind_density(state, density, blocks[("value", "density")], degree, support) result["density_storage_size"] = len(density) history = [numeric(row, ()) for row in rows(get("newton.csv"))] accepted = [row for row in history if row["accepted"] == 1] require(len(accepted) == int(meta["accepted_steps"]) and bool(int(meta["converged"])) == (record["exit_code"] == 0), "Convergence/exit/accepted-count contradiction") for i, row in enumerate(history): require(row["iteration"] == i and row["accepted"] in (0, 1), "Iteration ordering") for key in ("F", "seconds", "linear_seconds", "prepare_seconds", "refresh_seconds"): require(math.isfinite(row[key]) and row[key] >= 0, "Invalid iteration metric") if row["accepted"]: require(0 < row["alpha"] <= 1 and row["linear_status"] == 0 and row["linear_true"] <= float(meta["linear_tol"]) * (1 + 1e-10), "Accepted step without reported linear gate") if row["linear_rhs_norm"] > 0: close(row["linear_initial_over_rhs"], row["linear_initial_residual"] / row["linear_rhs_norm"], "Initial linear defect ratio") final_f = result["physical"]["18"]["metrics"]["normalized_bordered_residual"] initial_f = result["initial_physical"]["18"]["metrics"]["normalized_bordered_residual"] result["line_search"] = audit_line_search(rows(get("trials.csv")), history, initial_f, meta["strategy"]) if history: close(history[-1]["F"], final_f, "Final residual/log binding") threshold = max(float(meta["atol"]), float(meta["rtol"]) * initial_f) if int(meta["converged"]): require(final_f <= threshold * (1 + 1e-10), "Reported nonlinear convergence threshold fails") result["newton"] = {"converged": bool(int(meta["converged"])), "failure": meta.get("failure"), "accepted_steps": len(accepted), "initial_residual": initial_f, "final_residual": final_f, "threshold": threshold, "iterations": history, "iteration_seconds": math.fsum(r["seconds"] for r in history), "linear_iterations": sum(int(r["linear_iterations"]) for r in history)} fine = totals[("stored_density", degree, 30)] physical18 = result["physical"]["18"]["metrics"] for region_key, physical_key in (("mass", "mass"), ("volume", "volume")): close(fine[region_key], physical18[physical_key], "Sliced/physical " + region_key, rtol=2e-9, atol=1e-13) # Position factors raise the normal polynomial degree beyond the # eight-point sliced rule's exactness. Tangential refinement does not # test that error, so moment differences are observations, not identities. radius = float(meta["R_reference"]) require(math.isfinite(radius) and radius > 0, "Invalid reference radius") moment_comparisons = {"scope": "Order14 normal slicing versus unsliced q18; first moments and inertia may have normal quadrature error."} values = [("inertia_z", physical18["moment_of_inertia_z"], 2)] for j, axis in enumerate("xyz"): values.append((f"first_moment_{axis}", physical18[f"center_of_mass_{j}"] * physical18["mass"], 1)) for key, reference, power in values: require(math.isfinite(fine[key]) and math.isfinite(reference), "Nonfinite moment observation") difference = fine[key] - reference mass_radius_scale = abs(physical18["mass"]) * radius ** power moment_comparisons[key] = {"sliced": fine[key], "physical_q18": reference, "difference": difference, "difference_over_mass_radius_power": difference / mass_radius_scale if mass_radius_scale else None} result["sliced_vs_physical_moments"] = moment_comparisons return result class SelfTests(unittest.TestCase): def test_trial_sources_and_sequential_binding(self): trial_rows = [dict(iteration="0", trial="0", alpha="1", disposition="1", source="pressure force", F="nan", min_jacobian="nan"), dict(iteration="0", trial="1", alpha="0.5", disposition="4", source="residual metric", F="1.1", min_jacobian=".8"), dict(iteration="0", trial="2", alpha="0.25", disposition="0", source="", F=".8", min_jacobian=".9")] history = [dict(iteration=0., trials=3., accepted=1., alpha=.25, F=.8)] result = audit_line_search(trial_rows, history, 1., "baseline") self.assertEqual(result["pressure_force_rejections"], 1) self.assertEqual(result["residual_metric_rejections"], 1) self.assertEqual(result["minimum_accepted_alpha"], .25) for index, key, value in ((0, "F", "1"), (1, "F", ".7"), (1, "trial", "2"), (2, "alpha", ".5"), (2, "source", "pressure force"), (2, "F", ".7")): altered = [dict(row) for row in trial_rows] altered[index][key] = value with self.subTest(key=key), self.assertRaises(AssertionError): audit_line_search(altered, history, 1., "baseline") failed = [dict(iteration=0., trials=0., accepted=0., alpha=0., F=1.)] self.assertIsNone(audit_line_search([], failed, 1., "baseline")["minimum_accepted_alpha"]) def test_actual_closure_rule_support_and_tensor_counts(self): rows = [{"element": "2", "quadrature_order": "19", "quadrature_points": "1000"}, {"element": "5", "quadrature_order": "21", "quadrature_points": "1331"}] result = audit_closure_rules(rows, {2, 5}) self.assertEqual(result["distinct_achieved_orders"], [19, 21]) self.assertEqual(result["distinct_point_counts"], [1000, 1331]) for altered in (rows[:1], rows + [rows[0]], [{**rows[0], "element": "3"}, rows[1]], [{**rows[0], "quadrature_order": "20"}, rows[1]], [{**rows[0], "quadrature_order": "0"}, rows[1]], [{**rows[0], "quadrature_points": "1001"}, rows[1]]): with self.assertRaises(AssertionError): audit_closure_rules(altered, {2, 5}) def test_separate_enthalpy_observations(self): self.assertEqual(enthalpy_observations([], "h"), {"samples": 0, "minimum_enthalpy": None, "negative_enthalpy_samples": 0}) self.assertEqual(enthalpy_observations([{"h": -1e-18}, {"h": 0.}, {"h": 1.}], "h"), {"samples": 3, "minimum_enthalpy": -1e-18, "negative_enthalpy_samples": 1}) with self.assertRaises(AssertionError): enthalpy_observations([{"witness_h": math.nan}], "witness_h") def test_independent_roots(self): cases = ((["0", ".1875", "-1", "1", "0"], [(.25, .75)]), ([".025634765625", "-.34375", "1.34375", "-2", "1"], [(.125, .375), (.625, .875)]), ([".03515625", "-.375", "1.375", "-2", "1"], [])) for c, expected in cases: self.assertEqual(independent_partition(c)[0], expected) def test_resolved_narrow_quartic(self): delta = 2. ** -20 c = [.25 - delta ** 2, -1, 1.25 - delta ** 2, -1, 1] actual = independent_partition([str(v) for v in c])[0] self.assertEqual(len(actual), 1) for got, expected in zip(actual[0], (.5 - delta, .5 + delta)): close(got, expected, "Narrow quartic root", atol=1e-11) def test_density_state_binding(self): full = [0.] * 81 full[27:54] = [i / 10 for i in range(27)] state = [99., *full[27:54], 88.] block = {"offset": "1", "size": "27"} bind_density(state, full, block, 2, {1}) for target, position in ((full, 0), (state, 3)): original = target[position] target[position] = .12345 with self.assertRaises(AssertionError): bind_density(state, full, block, 2, {1}) target[position] = original def test_unsliced_error_binding(self): import csv with tempfile.TemporaryDirectory() as directory: path = Path(directory) / "density_order_error_quadrature.csv" records = [dict(control="stored_density", degree=3, element=0, attribute=1, quadrature_request=q, quadrature_achieved=q+1, quadrature_points=(q//2+1)**3, volume=1, mass=.25, target_mass=.25, eos_squared_integral=.001, signed_eos_integral=0, minimum_density=0, minimum_enthalpy=.1, negative_enthalpy_samples=0, eos_defined_everywhere_sampled=1) for q in (26, 30)] def write(): with path.open("w", newline="") as stream: writer = csv.DictWriter(stream, fieldnames=list(records[0])) writer.writeheader() writer.writerows(records) totals = {("stored_density", 3, 30): {"volume": 1., "mass": .25, "eos_squared_integral": .002}} write() result = audit_error_quadrature(lambda _: path, {("stored_density", 3)}, {0}, totals, False) self.assertEqual(result["Q3/stored_density"]["sliced_eos_squared_minus_q30"], .001) for key, value in (("quadrature_achieved", 29), ("signed_eos_integral", .1), ("eos_defined_everywhere_sampled", 0), ("mass", .3)): original = records[1][key] records[1][key] = value write() with self.subTest(key=key), self.assertRaises(AssertionError): audit_error_quadrature(lambda _: path, {("stored_density", 3)}, {0}, totals, False) records[1][key] = original def fixture(self): node = common.gauss_rule(16)[3][0] return {k: str(v) for k, v in dict(control="stored_density", degree=3, element=0, axis=0, orientation=1, u=node, v=node, t_left=.25, t_right=.75, c0=0, c1=.1875, c2=-1, c3=1, c4=0, negative_volume=.5, negative_mass=1/96, minimum_density=-.033009119061284215, witness_h=1, reconstruction_tolerance=2048 * EPS * 2.1875, resolved_negative_witness=1).items()} def test_valid_and_corrupted_interval(self): r = self.fixture() critical = independent_partition([r[f"c{i}"] for i in range(5)])[1] c = [float(r[f"c{i}"]) for i in range(5)] r["minimum_density"] = str(min(poly_value(c, x) for x in critical if .25 < x < .75)) nodes = [x for x, _ in common.gauss_rule(16)] interval_line([r], nodes) for key, value in (("t_right", ".8"), ("negative_mass", "2"), ("minimum_density", "-.5"), ("resolved_negative_witness", "0"), ("witness_h", "nan"), ("u", ".5")): with self.subTest(key=key), self.assertRaises(AssertionError): interval_line([{**r, key: value}], nodes) def test_hash_corruption(self): with tempfile.TemporaryDirectory() as directory: folder = Path(directory) (folder / "case").mkdir() item = folder / "case/value.txt" item.write_text("original") campaign = object.__new__(Campaign) campaign.path, campaign.ledger = folder, {} campaign.runs = {"case": {"output_hashes": {"value.txt": digest(item)}}} campaign.artifact("case", "value.txt") item.write_text("changed") with self.assertRaises(AssertionError): campaign.artifact("case", "value.txt") def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("campaigns", nargs="*", type=Path) parser.add_argument("--output", type=Path) parser.add_argument("--self-test", action="store_true") args = parser.parse_args() if args.self_test: test = unittest.TextTestRunner(verbosity=2).run(unittest.defaultTestLoader.loadTestsFromTestCase(SelfTests)) if not test.wasSuccessful(): raise SystemExit(1) if not args.campaigns: return if not args.campaigns or args.output is None: parser.error("Supply one or more completed campaigns and --output") report = {"artifact_passed": True, "auditor_sha256": digest(Path(__file__)), "scope": __doc__, "limitations": ["Local projection M,b and mapped weights are not exported.", "Complete negative partitions are checked only on lines with exported negative intervals; no global positivity claim.", "Shallow unresolved witnesses remain recorded; root position tests have coefficient-scaled allowances.", "Tangential convergence does not establish accuracy of the fixed eight-point normal rule; sliced first moments and inertia are compared observationally with q18.", "Density order alone changes among field spaces; production closure quadrature may also change under the unchanged order-aware policy. Only actual exported rules are reported.", "Coupled timing and convergence gates are producer observations, not independent equation evaluations."], "campaigns": {}, "consumed_sha256": {str(Path(common.__file__).resolve()): digest(Path(common.__file__))}} comparable = defaultdict(list) for path in args.campaigns: campaign = Campaign(path) records = {} for name in campaign.cases: try: records[name] = audit_case(campaign, name) key = records[name]["kind"], records[name]["n"], records[name]["J"] comparable[key].append((str(path), name, records[name])) except (AssertionError, KeyError, ValueError, OSError, StopIteration) as error: records[name] = {"artifact_passed": False, "error": str(error)} report["artifact_passed"] = False report["campaigns"][str(campaign.path)] = {"producer_commit": campaign.provenance["commit"], "mesh_sha256": campaign.provenance["mesh_sha256"], "cases": records} report["consumed_sha256"].update(campaign.ledger) # A cross-order comparison must leave every non-density block dimension # fixed. Different state-vector lengths are expected for density alone. for key, records in comparable.items(): reference = {k: v for k, v in records[0][2]["block_sizes"].items() if k not in ("value/density", "residual/barotropic_closure")} for _, name, record in records[1:]: actual = {k: v for k, v in record["block_sizes"].items() if k not in ("value/density", "residual/barotropic_closure")} if actual != reference: report["artifact_passed"] = False record["cross_order_error"] = "Non-density layout changed in " + str(key) + "/" + name args.output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text(json.dumps(common.clean(report), indent=2, allow_nan=False) + "\n") print("PASS" if report["artifact_passed"] else "FAIL", args.output) raise SystemExit(0 if report["artifact_passed"] else 1) if __name__ == "__main__": main()