#!/usr/bin/env python3 """Offline Q3 interpolation and bounded interior sampling; no solver or input mutation.""" import argparse import csv import itertools import json import math from pathlib import Path def read_rows(path): with path.open() as stream: return list(csv.DictReader(stream)) def weights(nodes, x): return [math.prod((x-other)/(node-other) for other in nodes if other != node) for node in nodes] def evaluate(axes, coefficients, point): w = [weights(axis, x) for axis, x in zip(axes, point)] return math.fsum(coefficients[i][j][k]*w[0][i]*w[1][j]*w[2][k] for i, j, k in itertools.product(range(4), repeat=3)) def extrema_record(): return {'minimum': math.inf, 'minimum_reference_point': None, 'negative_samples': 0, 'negative_beyond_roundoff_samples': 0} def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('campaign_root', nargs='?', type=Path, default=Path(__file__).resolve().parent) args = parser.parse_args() root = args.campaign_root.expanduser().resolve() trajectory = {int(r['iteration']): r for r in read_rows(root/'captures/n1_J0p10_capture/newton.csv')} case_ids = [('n1_J0p10_onset_between_bounds', 2), ('n1_J0p10_late_surface_trace', 11)] results = [] for name, iteration in case_ids: candidates = list(root.glob(f'**/{name}/surface_trace_nodes.csv')) if len(candidates) != 1: raise ValueError(f'Expected one nodal export for {name}') folder = candidates[0].parent raw = read_rows(candidates[0]) nodal = [r for r in raw if r['family'] == 'pressure'] closure = [r for r in raw if r['family'] == 'closure'] limits = read_rows(folder/'surface_trace_limits.csv') if len(nodal) != 64 or any(int(r['fe_order']) != 3 or int(r['node_count']) != 64 for r in nodal): raise ValueError(f'{name} is not the expected 64-node cubic element') if len({r['element'] for r in nodal + closure}) != 1: raise ValueError('This bounded follow-up expects one shared limiting element per captured state') for field in ('local_node', 'ip_x', 'ip_y', 'ip_z', 'h', 'dh'): if [r[field] for r in nodal] != [r[field] for r in closure]: raise ValueError(f'Family exports disagree: {name} {field}') axes = [sorted({float(r['ip_'+axis]) for r in nodal}) for axis in 'xyz'] if any(len(axis) != 4 or axis[0] != 0 or axis[-1] != 1 for axis in axes): raise ValueError('Expected four distinct interpolation coordinates per unit-cube axis') coordinates = [tuple(float(r['ip_'+axis]) for axis in 'xyz') for r in nodal] if set(coordinates) != set(itertools.product(*axes)) or len(set(coordinates)) != 64: raise ValueError('Export is not a complete unique tensor-product grid') fields = [] for field in ('h', 'dh'): values = [[[0.0 for _ in range(4)] for _ in range(4)] for _ in range(4)] for r, point in zip(nodal, coordinates): i, j, k = [axis.index(x) for axis, x in zip(axes, point)] values[i][j][k] = float(r[field]) fields.append(values) scale = max(abs(float(r[field])) for r in nodal for field in ('h', 'dh')) tolerance = 128*math.ulp(1.0)*scale checks = [] node_error = max(abs(evaluate(axes, coefficients, point)-float(r[field])) for field, coefficients in zip(('h', 'dh'), fields) for r, point in zip(nodal, coordinates)) if node_error > tolerance: raise ValueError(f'Node interpolation validation failed: {node_error}') for limit in limits: point = tuple(float(limit['ip_'+axis]) for axis in 'xyz') check = {'family': limit['family'], 'reference_point': point} for field, coefficients in zip(('h', 'dh'), fields): value = evaluate(axes, coefficients, point) expected = float(limit['full_'+field]) check[field] = value check[field+'_production'] = expected check[field+'_absolute_disagreement'] = abs(value-expected) if abs(value-expected) > tolerance: raise ValueError(f'Production-point interpolation validation failed: {check}') checks.append(check) pressure_limit = next(r for r in limits if r['family'] == 'pressure') pressure_x, pressure_y = [float(pressure_limit['ip_'+axis]) for axis in 'xy'] surface_line = [] for z in (1e-6, .001, .005, float(pressure_limit['ip_z'])): point = (pressure_x, pressure_y, z) h, dh = [evaluate(axes, values, point) for values in fields] surface_line.append({'reference_point': point, 'h': h, 'dh': dh}) zero_faces = [] for dim, axis in enumerate('xyz'): for side in (0.0, 1.0): face = [r for r in nodal if float(r['ip_'+axis]) == side] if len(face) == 16 and all(float(r['h']) == 0 and float(r['dh']) == 0 for r in face): zero_faces.append((dim, side)) if zero_faces != [(2, 0.0)]: raise ValueError(f'Expected the exported exact-zero face z=0; found {zero_faces}') grids = [] for dim, axis in enumerate('xyz'): values = {i/64 for i in range(1, 64)} | set(axes[dim][1:-1]) values.update(float(r['ip_'+axis]) for r in limits) if dim == 2: values.update(10**(-12+i/4) for i in range(45)) # 1e-12 through 1e-1. grids.append(sorted(values)) union = float(nodal[0]['alpha_union']) alphas = [('base', 0.0), ('accepted_alpha', float(trajectory[iteration]['alpha'])), ('near_below_union', .99*union), ('near_above_union', 1.01*union)] if iteration == 11: alphas.append(('unit_step', 1.0)) extrema = {label: extrema_record() for label, _ in alphas} directional_bound = {'minimum_nonnegative_base_ratio': math.inf, 'reference_point': None, 'h': None, 'dh': None, 'base_negative_and_direction_decreasing_samples': 0} wz = [weights(axes[2], z) for z in grids[2]] # Tensor contractions avoid a 64-term sum at every point. This is an # independent offline interpolation, validated above before use. for x in grids[0]: wx = weights(axes[0], x) yz = [[[math.fsum(values[i][j][k]*wx[i] for i in range(4)) for k in range(4)] for j in range(4)] for values in fields] for y in grids[1]: wy = weights(axes[1], y) z_values = [[math.fsum(values[j][k]*wy[j] for j in range(4)) for k in range(4)] for values in yz] for z, weights_z in zip(grids[2], wz): h, dh = [math.fsum(values[k]*weights_z[k] for k in range(4)) for values in z_values] point = (x, y, z) for label, alpha in alphas: value = h+alpha*dh record = extrema[label] if value < record['minimum']: record['minimum'] = value record['minimum_reference_point'] = point record['negative_samples'] += value < 0 record['negative_beyond_roundoff_samples'] += value < -tolerance if h < 0 and dh < 0: directional_bound['base_negative_and_direction_decreasing_samples'] += 1 if h >= 0 and dh < 0 and -h/dh < directional_bound['minimum_nonnegative_base_ratio']: directional_bound.update(minimum_nonnegative_base_ratio=-h/dh, reference_point=point, h=h, dh=dh) base_admissible = extrema['base']['negative_samples'] == 0 for label, alpha in alphas: record = extrema[label] h, dh = [evaluate(axes, values, record['minimum_reference_point']) for values in fields] error = abs(h+alpha*dh-record['minimum']) record['minimum_direct_sum_verification_error'] = error if error > tolerance: raise ValueError(f'Tensor contraction/direct-sum disagreement at minimum: {error}') directional_bound['sampled_base_admissible'] = base_admissible directional_bound['sampled_admissible_step_bound'] = directional_bound['minimum_nonnegative_base_ratio'] if base_admissible else None result = {'case': name, 'iteration': iteration, 'element': int(nodal[0]['element']), 'fe_order': 3, 'geometry_type': int(nodal[0]['geometry_type']), 'node_count': 64, 'axes': axes, 'zero_faces': [{'axis': 'z', 'side': 0}], 'node_interpolation_max_error': node_error, 'production_point_checks': checks, 'roundoff_scale_estimate': tolerance, 'pressure_xy_surface_line': surface_line, 'grid': {'axes': grids, 'shape': [len(g) for g in grids], 'sample_count': math.prod(map(len, grids)), 'strictly_interior': True, 'log_surface_distances': [1e-12, 1e-1]}, 'union_production_bound': union, 'alphas': dict(alphas), 'sampled_extrema': extrema, 'directional_bound': directional_bound} results.append(result) report = {'schema_version': 1, 'method': 'Independent tensor Q3 nodal interpolation and finite interior sampling', 'limitations': ['No global positivity certificate or rigorous polynomial minimum.', 'Reported minima are over the specified finite reference-coordinate grids.', 'Direction coefficients include saved-direction reconstruction roundoff; base h does not.', 'The roundoff scale is descriptive, not a certified interpolation error bound.'], 'cases': results} (root/'pressure_polynomial_checks.json').write_text(json.dumps(report, indent=2, allow_nan=False)+'\n') lines = ['# Offline enthalpy polynomial sampling', '', 'Both captured limiting elements have 64 distinct tensor-product nodes and cubic enthalpy order. ' 'The interpolation uses only exported physical nodal h and dh; it launches no solver and changes no inputs. ' 'The complete z=0 nodal face has exactly zero h and dh in both cases.', '', '| Iteration | Node reconstruction error | Max production-point discrepancy | Grid points |', '|---:|---:|---:|---:|'] for c in results: error = max(r[key] for r in c['production_point_checks'] for key in ('h_absolute_disagreement', 'dh_absolute_disagreement')) lines.append(f"| {c['iteration']} | {c['node_interpolation_max_error']:.9g} | {error:.9g} | {c['grid']['sample_count']} |") lines += ['', 'Each grid is strictly inside the reference cube: uniform spacing 1/64, exported interior node ' 'coordinates and production-limit coordinates, plus z distances 1e-12 through 1e-1 in quarter decades ' 'from the exactly zero face. Minima below are sampled minima, not global polynomial minima.', '', '| Iteration | State | Alpha | Sampled minimum h | Reference (x,y,z) | Negative samples |', '|---:|---|---:|---:|---|---:|'] for c in results: for label, record in c['sampled_extrema'].items(): point = ', '.join(f'{x:.9g}' for x in record['minimum_reference_point']) lines.append(f"| {c['iteration']} | {label} | {c['alphas'][label]:.9g} | {record['minimum']:.9g} | ({point}) | {record['negative_samples']} |") for c in results: bound = c['directional_bound'] if not bound['sampled_base_admissible']: lines += ['', f"Iteration {c['iteration']} already has negative enthalpy in the saved **base** on this grid. " 'A positivity-preserving step bound from that base is therefore undefined; the JSON separately ' 'reports the smallest nonnegative-base ratio only as a diagnostic. This base conclusion does ' 'not depend on reconstructing the Newton direction.'] else: lines += ['', f"Iteration {c['iteration']} has no negative base samples on this grid. The sampled directional " f"bound is {bound['sampled_admissible_step_bound']:.12g}, versus production-rule bound {c['union_production_bound']:.12g}. " 'This remains a sampled condition, not a guarantee of positivity everywhere.'] lines += ['', 'An independent direct 64-term sum also reproduces every contracted-grid minimum to roundoff. ' 'Along the pressure-limit x,y line, the offline base-field values are:', '', '| Iteration | z reference coordinate | Base h |', '|---:|---:|---:|'] for c in results: for row in c['pressure_xy_surface_line']: lines.append(f"| {c['iteration']} | {row['reference_point'][2]:.12g} | {row['h']:.12g} |") lines += ['', 'The production-point validation establishes consistency with the exported full-field values at ' 'both pressure and closure limit points to roundoff. Dense sampling can prove the existence of negative ' 'values when their size clearly exceeds numerical disagreement; absence of negative samples cannot ' 'prove positivity throughout the element. These data diagnose pointwise enthalpy, not virial balance ' 'or the origin of the selected Newton direction.', '', 'Reproduce with `python3 sample_pressure_polynomial.py /path/to/campaign_artifacts`. ' 'Machine-readable grids, locations, validation errors and sample counts are in pressure_polynomial_checks.json.'] (root/'pressure_polynomial_analysis.md').write_text('\n'.join(lines)+'\n') print(json.dumps([{'iteration': c['iteration'], 'base_minimum': c['sampled_extrema']['base']['minimum'], 'sample_count': c['grid']['sample_count'], 'base_admissible': c['directional_bound']['sampled_base_admissible']} for c in results])) if __name__ == '__main__': main()