#!/usr/bin/env python3 """Audit frozen density observations and summarize their physical location. This checks saved arithmetic and provenance, not a second FEM implementation. The projection and production closure are independently assembled in C++. """ import argparse import csv import hashlib import json import math from pathlib import Path import numpy as np def digest(path): return hashlib.sha256(path.read_bytes()).hexdigest() def rows(path): with path.open() as stream: return list(csv.DictReader(stream)) def vector(path): value = np.loadtxt(path) if value.ndim != 1 or int(value[0]) != len(value) - 1 or not np.all(np.isfinite(value)): raise ValueError(f"Invalid saved vector: {path}") return value[1:] def require(condition, message): if not condition: raise AssertionError(message) def close(left, right, label, rtol=2e-10, atol=1e-15): require(math.isfinite(left) and math.isfinite(right) and math.isclose(left, right, rel_tol=rtol, abs_tol=atol), f"{label}: {left} != {right}") def require_complete_hashes(record): required = {'metadata.txt', 'density_summary.csv', 'density_elements.csv', 'density_radial_bins.csv', 'density_physical.csv', 'negative_density_samples.csv', 'negative_density_projection_elements.csv', 'production_closure_rules.csv', 'production_closure_residual_full.txt', 'projection_production_closure_residual_full.txt', 'observed_density_coefficients.txt', 'density_invariants.csv', 'density_element_topology.csv', 'density_surface_traces.csv'} required.update(f'negative_density_projection_{control}_coefficients.txt' for control in ('production', 'q18', 'q26')) require(required <= record['output_hashes'].keys(), 'Run record omits a consumed artifact hash.') def clean(value): if isinstance(value, dict): return {str(k): clean(v) for k, v in value.items()} if isinstance(value, list): return [clean(v) for v in value] if isinstance(value, float) and not math.isfinite(value): return None return value def numeric(row, omit=()): return {key: float(value) for key, value in row.items() if key not in omit} def audit_case(folder): meta = dict(line.split('=', 1) for line in (folder / 'metadata.txt').read_text().splitlines() if '=' in line) n, K = float(meta['n']), float(meta['K']) summary = {(int(r['quadrature']), r['field']): numeric(r, ('quadrature', 'field')) for r in rows(folder / 'density_summary.csv')} element_rows = rows(folder / 'density_elements.csv') for (order, field), totals in summary.items(): matching = [numeric(r, ('quadrature', 'field', 'element', 'attribute', 'surface_layer')) for r in element_rows if int(r['quadrature']) == order and r['field'] == field] for key in ('samples', 'negative_samples', 'negative_h_samples', 'joint_negative_samples', 'volume', 'mass', 'negative_mass', 'negative_volume', 'negative_h_volume', 'closure_squared_integral'): close(math.fsum(r[key] for r in matching), totals[key], f"element sum {order}/{field}/{key}") close(min(r['minimum_density'] for r in matching), totals['minimum_density'], 'element minimum') close(max(r['maximum_density'] for r in matching), totals['maximum_density'], 'element maximum') radial = [numeric(r, ('quadrature', 'field', 'radial_bin_left', 'radial_bin_right')) for r in rows(folder / 'density_radial_bins.csv') if int(r['quadrature']) == order and r['field'] == field] for key in ('volume', 'mass', 'negative_mass', 'negative_volume'): close(math.fsum(r[key] for r in radial), totals[key], f"radial sum {order}/{field}/{key}") physical = {r['metric']: float(r['value']) for r in rows(folder / 'density_physical.csv')} observed18, observed26 = summary[(18, 'state')], summary[(26, 'state')] for key, physical_key in [('minimum_density', 'minimum_density'), ('minimum_h', 'minimum_enthalpy'), ('volume', 'volume'), ('mass', 'mass'), ('negative_mass', 'negative_density_mass_magnitude'), ('negative_samples', 'negative_density_samples'), ('negative_h_samples', 'negative_enthalpy_samples')]: close(observed18[key], physical[physical_key], f"physical observer {key}") close(observed18['negative_volume'] / observed18['volume'], physical['negative_density_volume_fraction'], 'physical negative volume fraction') negative = rows(folder / 'negative_density_samples.csv') require(len(negative) == observed26['negative_samples'], 'Sparse negative sample count does not bind to complete observation.') for r in negative: rho, h = float(r['rho']), float(r['h']) require(rho < 0 and int(r['quadrature']) == 26, 'Unexpected negative sample record.') close(float(r['radius']), math.sqrt(sum(float(r[d]) ** 2 for d in ('x', 'y', 'z'))), 'negative sample radius') close(float(r['rho_from_h']), (h / ((n + 1) * K)) ** n if h > 0 else 0., 'negative sample EOS') close(math.fsum(float(r['weight']) for r in negative), observed26['negative_volume'], 'negative export volume') close(math.fsum(-float(r['weight']) * float(r['rho']) for r in negative), observed26['negative_mass'], 'negative export mass') projection_rows = rows(folder / 'negative_density_projection_elements.csv') production = {int(r['element']): r for r in projection_rows if r['control'] == 'production'} actual_rules = {int(r['element']): r for r in rows(folder / 'production_closure_rules.csv')} require(production.keys() == actual_rules.keys(), 'Projection support differs from the actual closure.') for element, r in production.items(): rule = actual_rules[element] # MFEM may report the achieved odd Gauss degree for an even requested # degree. Equal points and this exact one-degree alias are intentional. requested, achieved = int(r['quadrature_order']), int(rule['quadrature_order']) require(achieved in (requested, requested + 1) and r['quadrature_points'] == rule['quadrature_points'], f"Production quadrature mismatch in element {element}.") actual = vector(folder / 'production_closure_residual_full.txt') reconstructed = vector(folder / 'projection_production_closure_residual_full.txt') require(actual.shape == reconstructed.shape, 'Full closure residual layout mismatch.') difference = actual - reconstructed rhs_norm = math.sqrt(math.fsum(float(r['rhs_norm']) ** 2 for r in production.values())) allowance = 2048 * np.finfo(float).eps * rhs_norm defect_norm = float(np.linalg.norm(difference)) require(defect_norm <= allowance, f"Independent closure residual mismatch {defect_norm}, allowance {allowance}.") saved = vector(folder / 'observed_density_coefficients.txt') coefficient_comparison = {} for control in ('production', 'q18', 'q26'): projected = vector(folder / f'negative_density_projection_{control}_coefficients.txt') require(projected.shape == saved.shape, 'Projected coefficient layout mismatch.') delta = projected - saved matching = [r for r in projection_rows if r['control'] == control] maximum = float(np.max(np.abs(delta))) close(maximum, max(float(r['coefficient_maximum_difference']) for r in matching), 'coefficient error summary') close(float(np.linalg.norm(delta)), math.sqrt(math.fsum(float(r['coefficient_difference_norm']) ** 2 for r in matching)), 'coefficient difference norm assembled from disjoint elements') for row in matching: require(float(row['constant_coefficient_error_maximum']) < 1e-10 and float(row['constant_sample_error_maximum']) < 1e-10, 'Independent constant reproduction failed.') coefficient_comparison[control] = dict(maximum_absolute_difference=maximum, l2_coefficient_difference=float(np.linalg.norm(delta)), relative_l2_coefficient_difference=float(np.linalg.norm(delta) / np.linalg.norm(saved))) invariants = {r['metric']: float(r['value']) for r in rows(folder / 'density_invariants.csv')} require(all(invariants[k] == 1 for k in ('accepted_physical_unchanged', 'accepted_normalized_unchanged', 'accepted_residual_unchanged')), 'Accepted state changed during observation.') topology = {int(r['element']): r for r in rows(folder / 'density_element_topology.csv')} affected = {int(r['element']) for r in negative} layers, attributes = {}, {} for r in negative: element, weight = int(r['element']), float(r['weight']) require(r['surface_layer'] == topology[element]['surface_layer'], 'Negative sample topology mismatch.') for table, key in ((layers, int(r['surface_layer'])), (attributes, int(r['attribute']))): item = table.setdefault(key, {'volume': 0., 'negative_mass': 0.}) item['volume'] += weight item['negative_mass'] -= weight * float(r['rho']) traces = rows(folder / 'density_surface_traces.csv') surface = [r for r in traces if float(r['inward_coordinate']) == 0.] for r in traces: require(all(math.isfinite(float(r[key])) for key in r), 'Nonfinite surface trace observation.') h = float(r['h']) close(float(r['rho_from_h']), (h / ((n + 1) * K)) ** n if h > 0 else 0., 'surface trace EOS') result = dict(n=n, K=K, source=meta['density_source'], normalized_residual=invariants['prepared_residual_norm'], state={str(q): summary[(q, 'state')] for q in (18, 26)}, controls={field: {str(q): summary[(q, field)] for q in (18, 26)} for field in ('eos_pointwise', 'production', 'q18', 'q26')}, coefficient_comparison=coefficient_comparison, closure_residual_check=dict(actual_norm=float(np.linalg.norm(actual)), independent_difference_norm=defect_norm, arithmetic_allowance=allowance, rhs_norm=rhs_norm), affected_elements=len(affected), stellar_elements=len(topology), surface_elements=sum(int(r['surface_layer']) == 0 for r in topology.values()), negative_volume_by_surface_layer=layers, negative_volume_by_attribute=attributes, density_dof_minimum=min(float(r['density_dof_minimum']) for r in topology.values()), surface_trace=dict(samples=len(surface), negative_samples=sum(float(r['rho']) < 0 for r in surface), minimum_density=min(float(r['rho']) for r in surface), maximum_abs_enthalpy=max(abs(float(r['h'])) for r in surface)), full_trace_minimum_density=min(float(r['rho']) for r in traces), physical_q18=physical, artifact_passed=True) return clean(result) def plots(campaign, output): import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt output.mkdir(parents=True, exist_ok=True) fig, axes = plt.subplots(1, 2, figsize=(11, 4.4), constrained_layout=True) for ax, name, title in zip(axes, ('n1p5_equilibrium', 'n3_seed'), ('n = 1.5 converged state', 'n = 3 initial state')): table = rows(campaign / name / 'density_surface_traces.csv') centers = [r for r in table if float(r['tangent_u']) == .5 and float(r['tangent_v']) == .5] witness = min(centers, key=lambda r: float(r['rho'])) line = [r for r in centers if r['face'] == witness['face']] line.sort(key=lambda r: float(r['distance_from_face'])) x = [float(r['distance_from_face']) for r in line] for key, label, style in [('rho_from_h', 'Positive EOS density from h', 'k-'), ('rho', 'Stored quadratic density', 'o-'), ('rho_projected_production', 'Independent closure projection', '--')]: ax.plot(x, [float(r[key]) for r in line], style, label=label, markersize=3) ax.axhline(0., color='0.55', linewidth=.8) ax.set_title(title) ax.set_xlabel('Distance inward from surface / reference radius') ax.set_ylabel('Density (model units)') depth_limit = .08 if name == 'n1p5_equilibrium' else .5 ax.set_xlim(0., x[-1] * depth_limit) shown = [float(r[k]) for r in line if float(r['inward_coordinate']) <= depth_limit for k in ('rho', 'rho_from_h', 'rho_projected_production')] span = max(shown) - min(shown) ax.set_ylim(min(shown) - .08 * span, max(shown) + .08 * span) ax.grid(alpha=.18) axes[0].legend(fontsize=8) fig.savefig(output / 'surface_density_traces.png', dpi=180) plt.close(fig) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('campaign', type=Path) parser.add_argument('--output', type=Path, required=True) parser.add_argument('--plots', type=Path) args = parser.parse_args() campaign = args.campaign.resolve() provenance = json.loads((campaign / 'provenance.json').read_text()) require(digest(campaign / 'plateau_diagnostic') == provenance['binary_sha256'], 'Producer binary changed.') require(digest(Path(provenance['mesh'])) == provenance['mesh_sha256'], 'Input mesh changed.') results = {} for case in provenance['cases']: name = case['id'] run = json.loads((campaign / (name + '_run.json')).read_text()) require(run['exit_code'] == 0 and run['inputs_unchanged'], f'Incomplete frozen diagnostic: {name}') require_complete_hashes(run) for relative, value in run['output_hashes'].items(): require(digest(campaign / name / relative) == value, f'Changed output: {name}/{relative}') for path, value in run['input_state_hashes'].items(): require(digest(Path(path)) == value, f'Changed input: {path}') results[name] = audit_case(campaign / name) payload = dict(source_commit=provenance['commit'], artifact_passed=True, auditor_sha256=digest(Path(__file__)), scope='Saved algebra, projection/production agreement and provenance; sampled spatial localization is not a sign certificate.', cases=results) args.output.write_text(json.dumps(clean(payload), indent=2, allow_nan=False) + '\n') if args.plots: plots(campaign, args.plots) print('PASS', args.output) if __name__ == '__main__': main()