"""Read the frozen AC50 fields, without importing or advancing an MHD model. Python standard library only. From an extracted replay directory: python3 -I -B ../measure_saved_fields.py --bundle . An optional --out must name a new file. Source/expected/field bytes are all verified against the bundle manifest before any numerical record is decoded. """ import argparse import gzip import hashlib import json import math from pathlib import Path def require(condition, message): if not condition: raise ValueError(message) def verify_members(bundle): manifest = json.loads((bundle / 'manifest.json').read_text()) paths = set() for member in manifest['members']: relative = member['path'] require(relative not in paths, f'Duplicate manifest member: {relative}') paths.add(relative) path = (bundle / relative).resolve() require(path.is_relative_to(bundle), f'Member escapes bundle: {relative}') body = path.read_bytes() require(len(body) == member['bytes'], f'Byte count differs: {relative}') require(hashlib.sha256(body).hexdigest() == member['sha256'], f'SHA-256 differs: {relative}') require(len(paths) == 299, 'This guide requires the frozen 299-member manifest') return paths def decode_record(path): return json.loads(gzip.decompress(path.read_bytes())) def field_measurements(field, count): require(len(field) == 4 and all(len(component) == count for component in field), 'A saved field must have four complete vertical arrays') require(all(math.isfinite(value) for component in field for value in component), 'A saved field contains a nonfinite value') radial_velocity, azimuthal_velocity, radial_magnetic, azimuthal_magnetic = field energy = math.fsum(value * value / 2 for component in field for value in component) / count reynolds = math.fsum(first * second for first, second in zip(radial_velocity, azimuthal_velocity)) / count maxwell = -math.fsum(first * second for first, second in zip(radial_magnetic, azimuthal_magnetic)) / count projections = [] for component_index, component in enumerate(field): basis = math.sin if component_index < 2 else math.cos projections.append(2 * math.fsum(value * basis(2 * math.pi * index / count) for index, value in enumerate(component)) / count) return dict(energy=energy, reynolds=reynolds, maxwell=maxwell, stress=reynolds + maxwell, projection=projections) def inspect(bundle): members = verify_members(bundle) # Admission precedes all numerical reads. receipt = decode_record(bundle / 'records/receipt.json.gz') require(receipt['status'] == 'PASS_LINEAR_SPATIAL_CHANNEL_ONLY', 'The original commissioning receipt did not pass') require(len(receipt['cases']) == 11 and len(receipt['gates']) == 45, 'Unexpected original commissioning scope') require(all(gate['passed'] for gate in receipt['gates']), 'An original gate failed') rows, frame_count = [], 0 maximum_energy_difference = maximum_stress_difference = 0. for case in receipt['cases']: measured_frames = [] previous_time = None for reference in case['frames']: relative = Path('records') / case['name'] / Path(reference['path']).name require(relative.as_posix() in members, f'Frame missing from manifest: {relative}') frame = decode_record(bundle / relative) measured = field_measurements(frame['field'], case['count']) require(math.isfinite(frame['time']) and (previous_time is None or frame['time'] > previous_time), f'Nonincreasing saved time: {relative}') energy_difference = abs(measured['energy'] - frame['diagnostics']['energy']) stress_difference = abs(measured['stress'] - frame['diagnostics']['stress']) require(energy_difference <= 1e-18, f'Energy mismatch: {relative}') require(stress_difference <= 1e-18, f'Stress mismatch: {relative}') require(math.isfinite(frame['integratedShearWork']), f'Nonfinite work: {relative}') maximum_energy_difference = max(maximum_energy_difference, energy_difference) maximum_stress_difference = max(maximum_stress_difference, stress_difference) measured_frames.append(dict(time=frame['time'], **measured, stageWork=frame['integratedShearWork'])) previous_time = frame['time'] frame_count += 1 first, last = measured_frames[0], measured_frames[-1] expected = case['expectedReferenceProjection'] error = math.sqrt(math.fsum((actual - target) ** 2 for actual, target in zip(last['projection'], expected)) / math.fsum(target ** 2 for target in expected)) saved_trapezoid_work = case['shear'] * case['rotation'] * math.fsum( (later['time'] - earlier['time']) * (earlier['stress'] + later['stress']) / 2 for earlier, later in zip(measured_frames, measured_frames[1:])) rows.append(dict(name=case['name'], referenceKind=case['referenceKind'], count=case['count'], interval=case['interval'], frames=len(measured_frames), initialEnergy=first['energy'], finalEnergy=last['energy'], finalReynolds=last['reynolds'], finalMaxwell=last['maxwell'], stageIntegratedWork=last['stageWork'], savedFrameTrapezoidWork=saved_trapezoid_work, workResidualOverInitialEnergy=abs(last['energy'] - first['energy'] - last['stageWork']) / first['energy'], independentlyProjectedRelativeError=error)) by_name = {row['name']: row for row in rows} spatial_ratios = [by_name[f'spatial-{coarse}']['independentlyProjectedRelativeError'] / by_name[f'spatial-{fine}']['independentlyProjectedRelativeError'] for coarse, fine in [(64, 128), (128, 256)]] temporal_ratios = [by_name[f'temporal-{coarse}']['independentlyProjectedRelativeError'] / by_name[f'temporal-{fine}']['independentlyProjectedRelativeError'] for coarse, fine in [('0.02', '0.01'), ('0.01', '0.005')]] return dict(status='PASS_READ_ONLY_FIELD_MEASUREMENTS', verifiedMembers=len(members), inspectedFrames=frame_count, newTimeSteps=0, cases=rows, maximumEnergyDifference=maximum_energy_difference, maximumStressDifference=maximum_stress_difference, spatialErrorRatios=spatial_ratios, temporalErrorRatios=temporal_ratios, comparisonNote='Projection targets and RK-stage work are saved original values; ' 'energy, stresses and projections above are remeasured from fields. ' 'The saved-frame trapezoid is a coarser, separate work estimate.') def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--bundle', required=True, type=Path) parser.add_argument('--out', type=Path) arguments = parser.parse_args() result = inspect(arguments.bundle.resolve()) encoded = json.dumps(result, indent=2, allow_nan=False) + '\n' if arguments.out: with arguments.out.open('x') as stream: stream.write(encoded) print(encoded, end='') if __name__ == '__main__': main()