245 lines
14 KiB
Python
245 lines
14 KiB
Python
#!/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()
|