From 7a07caeee33e2addb897b89682678b4cc41b1f3e Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 08:55:55 -0400 Subject: [PATCH] experiment: record pressure positivity and shape curvature evidence --- .../PRESSURE_AND_SHAPE_CAMPAIGN_2026-09-17.md | 378 ++++++++++++++ .../plateau_analysis/analyze_pressure.py | 475 ++++++++++++++++++ experiments/plateau_analysis/analyze_shape.py | 234 +++++++++ .../plateau_analysis/audit_focused_outputs.py | 143 ++++++ .../plateau_analysis/audit_surface_outputs.py | 167 ++++++ .../sample_pressure_polynomial.py | 218 ++++++++ experiments/plateau_diagnostic_README.md | 13 + 7 files changed, 1628 insertions(+) create mode 100644 experiments/PRESSURE_AND_SHAPE_CAMPAIGN_2026-09-17.md create mode 100644 experiments/plateau_analysis/analyze_pressure.py create mode 100644 experiments/plateau_analysis/analyze_shape.py create mode 100644 experiments/plateau_analysis/audit_focused_outputs.py create mode 100644 experiments/plateau_analysis/audit_surface_outputs.py create mode 100644 experiments/plateau_analysis/sample_pressure_polynomial.py diff --git a/experiments/PRESSURE_AND_SHAPE_CAMPAIGN_2026-09-17.md b/experiments/PRESSURE_AND_SHAPE_CAMPAIGN_2026-09-17.md new file mode 100644 index 0000000..1d71d10 --- /dev/null +++ b/experiments/PRESSURE_AND_SHAPE_CAMPAIGN_2026-09-17.md @@ -0,0 +1,378 @@ +# Pressure-preparation obstruction and shape curvature — September 17, 2026 + +The rotating plateau has a demonstrated local enthalpy-positivity obstruction: +the same interior polynomial sample approaches zero, while the surface trace +remains exactly zero. The late accepted polynomial is already negative between +production quadrature points. Separately, the n=1.5 full step is admissible but +nonlinear shape/field curvature raises its residual to 12.24 times the base. +The measured small-step behavior supports second-order curvature rather than a +first-order Jacobian defect along these directions. No production solver fix +or additional virial-qualified model is claimed. + +## Scope and controlled source + +This campaign adds experiment-only diagnostics to the preserved combined solver +at `95d9fba`. It changes no production operator, admissibility rule, normalization, +preconditioner, solver policy, or physics-facing API. The central-density +constraint remains optional. The working campaign branch is +`research/diagnostic-campaign-20260917`; pressure and shape helpers are retained +on separate hypothesis branches: + +| Purpose | Branch | Retained revision | +|---|---|---| +| Preserved combined solver | `research/combined-20260916` | `95d9fba` | +| Pressure and surface diagnostic helper | `research/pressure-obstruction-20260917` | `c6572d3` | +| Shape/field diagnostic helper | `research/shape-curvature-20260917` | `b24e6e0` | +| Integrated captures and focused probes | `research/diagnostic-campaign-20260917` | tested `fd08f7b`, `9e1b65d`, `4280c53` | + +The primary checkout remains on the diagnostic campaign branch. `main` remains +at `75cc638`. Existing untracked files were preserved. Source additions in this +campaign are confined to `experiments/`. + +The capture executable records source `fd08f7b`. It copies observer spans without +calling the prepared operator during Newton. A first-trial direction is +reconstructed from candidate minus base divided by alpha; the exact candidate +is also saved. A zero candidate-reconstruction error does not prove recovery of +every bit of the original internal correction. Roundoff estimates and direct +candidate replay distinguish that limitation from the original rejection. + +Raw artifacts, copied executables, manifests, input/output hashes and analysis +are retained under: +`/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/pressure_shape_2026-09-17`. + +## Trajectory fidelity + +The full twelve-attempt n=1, J=0.1 replay exactly reproduces every non-timing +Newton column, every trial and block record, the initial/final physical CSVs, +layout, and accepted-state bytes from the September 16 combined run. Its +accepted-state SHA-256 is +`ba55e7295c05f0c9ac5f32490f74d87556502d1744aa42deb7964729b2be2018`. +The n=1.5 replay reproduces the first three guarded attempts exactly, including +all ten trials and all 42 block records. There is no archived three-step endpoint +for an independent final-state comparison with that shortened run. + +Seven selected iterations save finite 178073-entry vectors. Independently +reconstructing each candidate gives zero componentwise round-trip error. The +late n=1 capture nevertheless has a direction-roundoff scale estimate of +1.40534e-8, about 4.88e-7 of its direction norm, because its first trial alpha is +1.90735e-6. The onset capture uses alpha=1 and has an estimate of 2.68e-14. +Exact candidate replay remains authoritative for the historical rejection. + +## Pressure-preparation diagnosis + +The exact saved first candidates from iterations 2, 3, 4, 7 and 11 all reproduce +`thermodynamic_domain / pressure_force / outside_domain` (enum values 2 / 6 / 5). +The production pressure source reaches this branch for finite negative enthalpy. +These are enthalpy-domain failures rather than inverted geometry. No new Krylov +solve is involved in the frozen replays. + +All iteration labels below are zero-based. Physical quantities use the model +units M=R_reference=G=1; h values below are not normalized residual entries. The pressure limiter stays at local +stellar element 459, pressure quadrature point 50; the closure limiter is point +36 of its different production rule on the same element. + +| Iteration | Pressure-point h | Direction dh there | Pressure alpha bound | Closure alpha bound | Live accepted alpha | +|---:|---:|---:|---:|---:|---:| +| 2 | 9.63674e-4 | -1.31441e-2 | 7.33159e-2 | 7.27489e-2 | 6.25e-2 | +| 3 | 1.42165e-4 | -2.11119e-3 | 6.73390e-2 | 7.15585e-2 | 6.25e-2 | +| 4 | 1.02159e-5 | -1.89704e-3 | 5.38520e-3 | 1.01581e-2 | 3.90625e-3 | +| 7 | 7.04574e-8 | -1.86569e-3 | 3.77648e-5 | 4.90600e-3 | 3.05176e-5 | +| 11 | 1.07227e-9 | -1.86493e-3 | 5.74963e-7 | 4.87083e-3 | 4.76837e-7 | + +The pressure step budget falls by **127,514 times**. These are physical enthalpy +values and corrections evaluated with unchanged normalization at the same +production points. For consecutive captured states, the next h +agrees with `h + accepted_alpha * dh` to 2.01e-18. At the late state, the closure +bound is about 8,472 times looser than the pressure bound. The shrinking pressure +h explains the vanishing admissible steps along these measured directions; it +does not yet explain why the computed directions keep decreasing h there. + +The element touches the stellar surface. Its limiting pressure point has +reference Cartesian radius 0.997553 and lies near the negative rotation pole +(about 8.66 degrees from the negative z axis at onset). This is an interior +quadrature point in a surface-adjacent element, not a surface boundary row. +Topology establishes adjacency, not distance to the deformed surface. + +Both onset and late probes bracket the production enthalpy bound with 0.99 and +1.01 times that bound. Enthalpy-only steps reproduce the crossing rejection; +shape-only steps prepare successfully. Across the two 19-trial grids, no +geometry rejection occurs. Affine limiting-point h predictions agree with +independent point evaluations within 1.05e-17. Preparation success alone does +not imply merit acceptance. + +At onset the closure bound is slightly tighter. The original large candidate +crosses both bounds, and pressure is prepared first, so its rejection stage +cannot identify the tighter family. The supplemental candidate at +alpha=0.07303236911639208 lies strictly between the two bounds. Its predicted +closure h is -5.01979e-6 while the limiting pressure h remains +3.72619e-6. +The production replay confirms that sequence: pressure preparation passes and +closure preparation rejects with `thermodynamic_domain / barotropic_closure / +outside_domain` (2 / 3 / 5). This distinguishes the first reported failure stage +from the tightest sampled bound. + +The supplemental trace inspection uses the actual compiled surface carrier map. +All **866 surface enthalpy coefficients and directions are exactly zero** at +onset and late plateau. The exported indices also establish exact zero surface +enthalpy in all five captured physical bases and their five exact first-trial +candidates. Their contributions to h and dh at both limiting interior points +are exactly zero. The surface condition exists and is satisfied here; boundary +leakage does not explain this wall. + +The limiting element is cubic, with 64 enthalpy FE nodes. At onset, a step of +1.01 times the sampled union bound leaves all 64 nodal h values nonnegative, +while closure and pressure interior samples have h=-1.28819e-5 and -2.11e-6, +respectively. At the late state, **even the full unit step leaves all nodes +nonnegative**, but h at the pressure limiting point is -1.86493e-3. This +establishes interior polynomial undershoot along these corrections. Nodal +positivity or exact zero surface trace is insufficient to prevent it. It does +not establish an interpolation implementation bug. The onset full unit step +also makes five nodes negative, so that larger step must not be described as +an exclusively between-node failure. + +A final offline check reconstructs the two local tensor-Q3 polynomials from +those 64 physical FE nodal values. It reproduces all nodal values exactly and +both production limiting-point h/dh values within 6.94e-18. The same interior +sampling grid has 495,264 points per inspected state, including uniform spacing +1/64 and logarithmic distances from 1e-12 to 1e-1 from the zero-trace z=0 face. + +| Local polynomial state | Minimum sampled h | Reference-element location | Interpretation | +|---|---:|---|---| +| Onset base | positive on this grid | — | No negativity detected; not a positivity certificate | +| Onset accepted alpha=0.0625 | positive on this grid | — | Same limited conclusion | +| Onset alpha=0.99 times production bound | -2.36406e-5 | (0.0130467, 0.515625, 0.01) | Negative between production samples | +| Late captured base | **-2.87114e-5** | (0.0130467, 0.515625, 0.01) | Existing pointwise violation before the next step | + +The late base has 2,401 negative samples on that grid. This conclusion uses only +the captured base polynomial and is independent of direction-reconstruction +roundoff. A separate line evaluation at the pressure point's x,y coordinates +and reference z=0.005 gives h=-1.54933e-5. An independent 60-digit Decimal +reconstruction gives -1.5493345386780819e-5 at that point; its input hashes, +formula and production-point checks are retained in +`tensor_q3_independent_check.json`. The negative magnitudes are many orders +above interpolation-check discrepancies. They prove the existence of negative +values in this element; the sampled minimum is not a certified global minimum. + +For the onset direction, the finer sampled step bound is 0.0665364 versus the +production union bound 0.0727489. For the late base, a positivity-preserving step +bound from an everywhere-admissible base cannot be asserted because the base +already has negative h. **Passing production quadrature preparation currently +does not ensure internal pointwise enthalpy positivity.** The inspected late +state is nonconverged and already fails the virial milestone; this finding does +not assert that the previously converged nonrotating/low-J controls share the +same defect. No new quadrature or EOS rule was installed. + +## Shape and field curvature + +Let p be the captured correction, p_s its actual surface-deformation block, +and p_f all remaining fields and scalars. The diagnostic prepares the same +frozen base before every Jacobian action and measures + +`R_d(alpha) = F(x + alpha*p_d) - F(x) - alpha*J*p_d`. + +The mixed component is the vector `R_full - R_shape - R_fields`, formed before +any norm. Each case uses six common alphas separated by decades; the largest +is capped by geometry and both enthalpy rules so all three directions can be +compared at the same step. Directions are not renormalized separately. + +| Frozen case | Common alpha | ‖Ftrial‖/‖F0‖ | ‖L‖/‖F0‖ | ‖R‖ / alpha² | Gravity fraction of ‖R‖² | +|---|---:|---:|---:|---:|---:| +| n=1, J=0.1, iteration 2 | 0.058199106 | 0.94307844 | 0.94180757 | 0.062255721 | 99.8390% | +| n=1.5, J=0, iteration 1 | 1 | 12.243161 | 0.02576812 | 0.0021890295 | 99.7747% | +| n=3 constrained seed | 0.12194828 | 1.0237652 | 0.87826014 | 0.11177311 | 90.7812% | + +For n=1.5 at iteration 1, the exact first candidate is preparation-admissible, +with positive minimum Jacobian 0.896116. Its residual norm reproduces the original +merit-rejected trial within 1.43e-17. The full-step linear prediction would leave +only 2.5768% of the starting residual, whereas the actual trial has **12.243 times** +the starting residual. At alpha=0.1, the actual ratio is still 1.02052 versus the +linear prediction 0.900070; at 0.01 it falls to 0.991209. This directly measures +curvature-driven merit rejection separately from the rotating enthalpy wall. + +At this n=1.5 base, gravity-gradient shape and field action norms are 0.115959 +and 0.116055, with cosine -0.999999831; their combined action norm is only +1.16990e-4. Poisson actions likewise nearly cancel (cosine -0.999996055). +Truncating just the shape block would destroy this measured first-order +cancellation. A full-step second-order remainder of 2.18903e-3 survives, with +99.7747% of its squared norm in the two gravity residual blocks. + +The remainder divided by alpha squared stabilizes as alpha decreases, while +resolved gravity remainder slopes approach two. With shape held fixed, gravity +is affine in its other fields in the source, and the measured fields-only gravity +remainders remain at arithmetic scale. The n=1.5 maximum is 2.96e-14, 1.042 times +the diagnostic's heuristic floor; the strict below-floor flag is therefore false, +but that floor is not a certified assembly-error bound. This slight excess is +not evidence of nonlinear fixed-shape gravity. The smallest steps that reach the +arithmetic floor are excluded from slope claims. + +| Frozen case | 2⟨L,R_shape⟩ | 2⟨L,R_fields⟩ | 2⟨L,R_mixed⟩ | ‖R_full‖² | Total curvature cost | +|---|---:|---:|---:|---:|---:| +| n=1, J=0.1, iteration 2 | -1.952454e-08 | -3.042692e-10 | 9.736527e-08 | 4.446556e-08 | 1.220020e-07 | +| n=1.5, J=0, iteration 1 | -7.601828e-10 | 1.587389e-12 | 1.185527e-09 | 4.791850e-06 | 4.792277e-06 | +| n=3 constrained seed | 1.458147e-09 | -9.452314e-08 | -9.285962e-09 | 2.762978e-06 | 2.660627e-06 | + +Here `L = F0 + alpha*Jp` and the exact scalar identity is +`||Ftrial||² = ||L||² + 2 + ||Rfull||²`. +The three signed projection terms add to `2`; the last squared term +contains cross terms between remainder components. Positive values worsen the +linear squared-merit prediction. At the n=1.5 full step, the dominant cost is +the squared remainder itself, not the small signed projections onto the nearly +zero linear prediction. At alpha=0.1, both shape and mixed projections are adverse. +At the admissible rotating onset step, curvature modestly worsens the prediction +but still allows residual decrease; its larger historical steps instead fail +enthalpy preparation. + +These observations support genuine nonlinear shape/field curvature along the +measured directions, with no detected first-order Jacobian inconsistency there. +They do not certify every Jacobian column. The gravity-block aggregates include +both stellar and exterior domains and do not establish where the nonlinear +remainder originates. No physical importance is assigned to exterior accuracy +beyond its influence on the interior solution. + +The n=3 seed retains its explicit central-density constraint. Its saved +restart-160 correction has true relative linear defect 0.0402590212946322, +exactly matching the earlier failed linear solve and exceeding the 0.03 gate. +No new solve was substituted. At the common cap the trial residual is 2.3765% +above the base despite a predicted 12.174% decrease. The mixed component is large, +but its signed projection is favorable; the positive squared-remainder term +outweighs all favorable projections. At one tenth the cap, the residual decreases +and net curvature slightly improves the linear prediction. This prevents a +blanket attribution of every failed step to an adverse mixed projection. + +## Consequences and next experiments + +The two mechanisms should be tested and treated separately. + +1. **Pressure/closure feasibility.** A production-rule enthalpy step cap would + avoid known inadmissible trials. It would not, by itself, remove a wall where + the same h tends to zero while dh stays negative: the existing accepted steps + are already 72.5–92.8% of the sampled union bound. The next discriminating + numerical experiment is a tighter frozen linear solve, with local h-direction + and component defects measured explicitly, to test how strongly the local + direction depends on linear accuracy. A + persistent negative direction would still permit nonlinear-path or basin + effects; it would not prove that no admissible discrete root exists. + Preserve EOS admissibility and use quadrature/interior polynomial checks when + testing a positivity representation or constrained step; clipping nodal values + cannot prevent the demonstrated between-node undershoot. Test interior + polynomial positivity at accepted states as a separate validation check. + Increasing quadrature coverage can detect the violation earlier, but detection + alone is not a convergence remedy. A trial starting from the late base also + needs feasibility restoration if positivity is to be maintained everywhere. + Clipping h would + change the equations and would not establish internal correctness. +2. **Shape/field nonlinear balance.** Test a trial-only gravity completion at + fixed candidate shape and density, solving the linear mixed-gravity subsystem + there before reevaluating every coupled residual. Gravity completion changes the trial path and reduces the gravity defect by + construction. The test is whether that produces a better coupled trial, + interior solution and virial after every residual is reevaluated. + It is a hypothesis, not an implemented production change; other residuals, + pointwise interior checks and virial must all be measured. Shape-only damping + is poorly motivated by the strong first-order cancellation in these cases. +3. **Where the merit is spent.** Before changing normalization or exterior + weighting, partition the dominant gravity remainder by stellar/exterior + contribution and measure its effect on interior hydrostatic balance. These + whole-block probes do not justify demanding exterior accuracy for its own + sake, nor do they justify dropping exterior residuals that influence the star. + +The present data do not single out preconditioner failure for the rotating or +n=1.5 captured steps: their current true linear defects still satisfy the 0.03 +forcing gate. That gate is a global norm criterion, not proof that the local +enthalpy direction is adequate near a vanishing h. The n=3 frozen correction +fails that gate and is reported separately. No production preconditioner, +normalization, EOS admissibility or line-search rule has been altered here. + +## Verification and reproducibility limits + +The driver was built in Release with Homebrew Clang 22.1.8, IPO enabled, +profiling disabled and uniform polynomial order increment zero. Runs were +serial, on one MPI rank, using the same 1216-hex mesh (832 stellar, 384 vacuum) +and the unchanged production quadrature rules. Mesh SHA-256: +`ce11ae99e21e6c3bbfbee402acd2e191c1da0d8261d2227b4203f73f2f337a74`. +Case manifests and copied immutable executables retain the tested sources even +while the experiment branch acquires later diagnostic additions. Performance +was not compared in this campaign. + +Two instrumented trajectory runs reached their requested iteration caps +(expected model exit 3). All eight focused cases and both supplemental cases +completed successfully (exit 0), with matching binary and input/output hashes. +The n=1.5 rejected-candidate residual norms agree with their historical values +to 1.43e-17 and 1.84e-18; the two norm implementations have different reduction +orders, so those numbers are roundoff comparisons rather than byte comparisons. +All captured base residual vectors reproduce exactly. + +All ten frozen cases restore their final base residual exactly. Each shape case +has 22 admissible preparation events, for 66 total; all accepted normalized +states, physical states and residuals are unchanged. Direction partitions are +exact. The largest Jacobian partition error is 8.33e-17. Aggregate signed-merit +and projection identities agree to better than 6.8e-15 relative. All 38 pressure +crossing trials have complete, unique records and no geometry rejection. +Independent audits check vector copies, quadrature-bound identity, CSV coverage, +all surface coefficient statistics and nodal identities. + +Implicit neighboring capture state/residual inputs are covered by the capture +output hashes as well as independent checks; the runner's explicit-input hashes +alone do not cover that auxiliary-read path. Late reconstructed directions retain +the subtraction-roundoff limitation documented above. The late recomputed true +linear defect differs from its historical value by 2.25e-9, while both remain +below 0.03; exact-candidate replay is independent of that reconstruction. + +No production/library or public-header code changed relative to `95d9fba`; +the existing production C++ suite was not rerun for these experiment-only +changes. Validation instead includes compilation of the diagnostic executables, +exact trajectory comparisons, production candidate replays, independently +recomputed CSV identities, input/output hashes, and frozen-state restoration. +A zero diagnostic process exit means the observations completed. It does not +mean a stellar equilibrium converged. + +This campaign establishes no additional virial-qualified model. The captured +rotating endpoint is exactly the previous nonconverged endpoint: its P(h)-based +virial errors remain approximately 4.75e-4 and 8.18e-4, above the required 1e-5. +The n=1.5 capture is intentionally only a three-step prefix; its eventual +convergence and pointwise accuracy are not inferred from these probes. J=0.1 +is an angular-momentum input, not 10% of critical rotation. Critical-rotation +fractions were not computed. + +## Reproduction and retained analyses + +The three case manifests are `plateau_capture_cases_2026-09-17.json`, +`plateau_focused_cases_2026-09-17.json`, and +`plateau_supplemental_cases_2026-09-17.json` in `experiments/`. Their tested source +commits are respectively `fd08f7b`, `9e1b65d`, and `4280c53`. Each artifact folder +contains `provenance.json`, exact per-case commands, case input/output hashes, +logs and the copied executable. The midpoint candidate also retains a derivation +record and SHA-256. No raw state files or binary build products are committed. + +Run a new numerical campaign with a compatible configured build and a fresh +output directory: + +```sh +python3 experiments/run_plateau_campaign.py \ + --source . --build cmake-build-release-homebrew --mesh sandbox.smesh \ + --cases experiments/plateau_focused_cases_2026-09-17.json \ + --cmake /Applications/CLion.app/Contents/bin/cmake/mac/aarch64/bin/cmake \ + --jobs 8 --output /absolute/path/to/new-results +``` + +Use the corresponding tested source checkout for exact source reproduction. +The supplied case manifests reference the retained absolute capture paths; +relocating artifacts requires changing those paths consistently. The runner +checks committed source, rebuilds, records provenance and refuses reused output +directories. Replaying the copied executables with the recorded arguments and a +new output path is another way to retain the exact tested binary. + +The portable analysis scripts under `experiments/plateau_analysis/` take the +artifact root containing `captures/`, `focused/`, and `supplemental/`: + +```sh +python3 experiments/plateau_analysis/analyze_pressure.py /absolute/path/to/artifact-root +python3 experiments/plateau_analysis/analyze_shape.py --campaign /absolute/path/to/artifact-root +python3 experiments/plateau_analysis/audit_focused_outputs.py /absolute/path/to/artifact-root +python3 experiments/plateau_analysis/audit_surface_outputs.py /absolute/path/to/artifact-root +python3 experiments/plateau_analysis/sample_pressure_polynomial.py /absolute/path/to/artifact-root +``` + +They do not run the solver or alter raw case data. They regenerate compact +Markdown/JSON analyses alongside those data. See `pressure_analysis.md`, +`shape_analysis.md`, `pressure_polynomial_analysis.md`, `capture_validation.md`, +and `focused_method_validation.md` +in the retained artifact root for the detailed per-case checks. diff --git a/experiments/plateau_analysis/analyze_pressure.py b/experiments/plateau_analysis/analyze_pressure.py new file mode 100644 index 0000000..457b077 --- /dev/null +++ b/experiments/plateau_analysis/analyze_pressure.py @@ -0,0 +1,475 @@ +#!/usr/bin/env python3 +"""Read completed/partial pressure diagnostics; never launch a solver or modify inputs.""" +import argparse +import csv +import json +import math +from pathlib import Path + +ROOT = Path(__file__).resolve().parent +ITERATIONS = (2, 3, 4, 7, 11) + + +def rows(path): + if not path.exists(): + return [] + with path.open() as stream: + return [r for r in csv.DictReader(stream) if None not in r and all(v is not None for v in r.values())] + + +def metrics(path): + return {r['metric']: float(r['value']) for r in rows(path)} + + +def metadata(path): + if not path.exists(): + return {} + return dict(line.split('=', 1) for line in path.read_text().splitlines() if '=' in line) + + +def numeric(row): + result = {} + for key, value in row.items(): + try: + result[key] = float(value) + except ValueError: + result[key] = value + return result + + +def negative_h_rejection(row): + # Enums at source9e1b65d: thermodynamic_domain=2, pressure_force=6, + # outside_domain=5. The prepared pressure operator emits this for h<0. + return all(float(row.get(k, -999)) == v for k, v in + [('admissible', 0), ('reason', 2), ('stage', 6), ('eos_code', 5)]) + + +def sanitize(value): + if isinstance(value, float) and not math.isfinite(value): + return None + if isinstance(value, dict): + return {k: sanitize(v) for k, v in value.items()} + if isinstance(value, list): + return [sanitize(v) for v in value] + return value + + +def case(iteration): + folder = ROOT / 'focused' / f'n1_J0p10_iteration_{iteration}' + capture = metadata(ROOT / 'captures/n1_J0p10_capture' / f'capture_{iteration}/capture.txt') + bound = metrics(folder / 'frozen_enthalpy_bound.csv') + inputs = metrics(folder / 'frozen_input_checks.csv') + replay = rows(folder / 'replayed_candidate.csv') + trial_rows = rows(folder / 'pressure_obstruction_trials.csv') + limits = [numeric(r) for r in rows(folder / 'pressure_obstruction_limits.csv')] + summary = metrics(folder / 'pressure_obstruction_summary.csv') + live = next((r for r in rows(ROOT / 'captures/n1_J0p10_capture/newton.csv') + if int(r['iteration']) == iteration), None) + result = {'iteration': iteration, 'folder': str(folder), 'capture': capture, + 'bounds': bound, 'input_checks': inputs, 'limits': limits, + 'exact_candidate': numeric(replay[0]) if replay else None, + 'exact_candidate_confirms_pressure_negative_h': negative_h_rejection(replay[0]) if replay else None, + 'restoration_completed': 'final_restoration_residual_difference' in inputs, + 'pressure_trials_completed': len(trial_rows), + 'pressure_trials_expected': summary.get('candidate_trial_count'), + 'trial_checks': [], 'live_newton_step': numeric(live) if live else None} + result['pressure_grid_complete'] = (summary.get('candidate_trial_count') is not None and + len(trial_rows) == summary['candidate_trial_count']) + if capture: + result['direction_roundoff_relative_estimate'] = ( + float(capture['direction_roundoff_scale_estimate']) / float(capture['direction_norm'])) + if bound: + result['closure_to_pressure_bound_ratio'] = bound['closure_alpha_h'] / bound['pressure_alpha_h'] + result['exact_first_trial_over_pressure_bound'] = float(capture['first_trial_alpha']) / bound['pressure_alpha_h'] + if live: + result['accepted_alpha_over_union_bound'] = float(live['alpha']) / bound['alpha_thermodynamic'] + result['tightest_sampled_family'] = ('closure' if bound['closure_alpha_h'] < bound['pressure_alpha_h'] else + 'pressure' if bound['pressure_alpha_h'] < bound['closure_alpha_h'] else 'tie') + midpoint = .5 * (bound['closure_alpha_h'] + bound['pressure_alpha_h']) + result['between_bounds_midpoint'] = { + 'alpha': midpoint, + 'closure_limiting_h_prediction': bound['limiting_h'] + midpoint * bound['limiting_dh'], + 'pressure_limiting_h_prediction': bound['pressure_limiting_h'] + midpoint * bound['pressure_limiting_dh'], + } + for row in trial_rows: + r = numeric(row) + errors = [abs(r[f'{family}_limiting_h_evaluated'] - r[f'{family}_limiting_h_prediction']) + for family in ('pressure', 'closure')] + sign_checks = [((r[f'{family}_limiting_h_evaluated'] < 0) == + (r[f'{family}_limiting_h_prediction'] < 0)) for family in ('pressure', 'closure')] + ratio = r['alpha_over_union_bound'] + r['maximum_limiting_h_evaluation_prediction_absolute_difference'] = max(errors) + r['both_limiting_point_signs_match_prediction'] = all(sign_checks) + r['confirms_pressure_negative_h'] = negative_h_rejection(row) + # Descriptive expectation check for these fixed cases, not a general + # assertion that h>=0 guarantees preparation or merit acceptance. + expected = (r['direction'] == 'shape_only' or ratio < 1) + r['admissibility_matches_enthalpy_only_prediction'] = bool(r['admissible']) == expected + result['trial_checks'].append(r) + if trial_rows: + result['all_limiting_sign_predictions_match'] = all(r['both_limiting_point_signs_match_prediction'] + for r in result['trial_checks']) + result['all_admissibility_matches_enthalpy_only_prediction'] = all( + r['admissibility_matches_enthalpy_only_prediction'] for r in result['trial_checks']) + result['mapping_rejections'] = sum(r['admissible'] == 0 and r['reason'] in (0, 1) + for r in result['trial_checks']) + result['maximum_limiting_h_absolute_disagreement'] = max( + r['maximum_limiting_h_evaluation_prediction_absolute_difference'] for r in result['trial_checks']) + shape_rows = [r for r in result['trial_checks'] if r['direction'] == 'shape_only'] + h_crossing_rows = [r for r in result['trial_checks'] + if r['direction'] == 'enthalpy_only' and r['alpha_over_union_bound'] > 1] + result['shape_only_all_prepared'] = all(r['admissible'] == 1 for r in shape_rows) if len(shape_rows) == 4 else None + result['enthalpy_only_crossing_failures_all_pressure_negative_h'] = ( + all(r['confirms_pressure_negative_h'] for r in h_crossing_rows) if len(h_crossing_rows) == 2 else None) + return result + + +def fmt(value): + return 'pending' if value is None else f'{value:.9g}' + + +def vector_entries(path, indices): + """Read only selected zero-based entries; reject an incompatible saved vector.""" + if not path.exists(): + return None + selected = set(indices) + found = {} + with path.open() as stream: + size = int(next(stream).strip()) + if not selected or min(selected) < 0 or max(selected) >= size: + raise ValueError(f'Carrier indices incompatible with {path}') + count = 0 + for index, line in enumerate(stream): + value = float(line) + count += 1 + if index in selected: + found[index] = value + if count != size or len(found) != len(selected): + raise ValueError(f'Incomplete vector {path}') + return [found[index] for index in indices] + + +def extrema(values): + return {'count': len(values), 'minimum': min(values), 'maximum': max(values), + 'maximum_absolute': max(abs(x) for x in values), + 'l2_norm': math.sqrt(math.fsum(x*x for x in values)), + 'negative_count': sum(x < 0 for x in values), + 'nonzero_count': sum(x != 0 for x in values)} + + +def supplemental(cases): + """Find the two supplemental case IDs without imposing an output-parent name.""" + result = {'cases': [], 'offline_surface_checks': []} + case_ids = [('n1_J0p10_onset_between_bounds', 2), ('n1_J0p10_late_surface_trace', 11)] + reference_rows = None + for name, iteration in case_ids: + candidates = [path.parent for path in ROOT.glob(f'**/{name}/surface_trace_summary.csv')] + if len(candidates) > 1: + raise ValueError(f'Ambiguous duplicate supplemental case: {name}') + if not candidates: + result['cases'].append({'id': name, 'iteration': iteration, 'available': False}) + continue + folder = candidates[0] + carrier_rows = [numeric(r) for r in rows(folder / 'surface_trace_rows.csv')] + summary = metrics(folder / 'surface_trace_summary.csv') + if 'pressure_alpha_h' not in summary: + result['cases'].append({'id': name, 'iteration': iteration, 'available': False}) + continue + contribution_rows = [numeric(r) for r in rows(folder / 'surface_trace_limits.csv')] + node_rows = [numeric(r) for r in rows(folder / 'surface_trace_nodes.csv')] + replay_rows = rows(folder / 'replayed_candidate.csv') + source_case = next(c for c in cases if c['iteration'] == iteration) + details = {'id': name, 'iteration': iteration, 'available': True, 'folder': str(folder), + 'summary': summary, 'contributions': contribution_rows, + 'replayed_candidate': numeric(replay_rows[0]) if replay_rows else None, + 'carrier_rows_complete': len(carrier_rows) == summary.get('surface_row_count'), + 'nodal_checks': [], + 'restoration': metrics(folder / 'frozen_input_checks.csv')} + provenance = folder.parent / 'provenance.json' + if provenance.exists(): + p = json.loads(provenance.read_text()) + details['provenance'] = {key: p.get(key) for key in ('commit', 'branch', 'mesh_sha256')} + if carrier_rows and details['carrier_rows_complete']: + indices = [int(r['root_index']) for r in carrier_rows] + if len(set(indices)) != len(indices): + raise ValueError(f'Duplicate surface indices in {folder}') + if reference_rows is not None and indices != reference_rows: + raise ValueError('Supplemental surface carrier maps differ') + reference_rows = indices + for family in ('closure', 'pressure'): + nodal = [r for r in node_rows if r['family'] == family] + if not nodal: + continue + if len(nodal) != nodal[0]['node_count']: + continue # A live CSV may not yet have all nodal rows. + union = nodal[0]['alpha_union'] + factors = [('base', 0), ('near_below_union', .99*union), + ('near_above_union', 1.01*union), ('unit_step', 1)] + if source_case['capture']: + factors.append(('captured_first_alpha', float(source_case['capture']['first_trial_alpha']))) + if source_case['live_newton_step']: + factors.append(('accepted_alpha', source_case['live_newton_step']['alpha'])) + scale = max(abs(r[k]) for r in nodal for k in ('h', 'dh')) + tolerance = 64 * math.ulp(1.0) * scale + limit = next((r for r in contribution_rows if r['family'] == family), None) + if limit is None: + continue + for label, alpha in factors: + values = [r['h'] + alpha*r['dh'] for r in nodal] + stats = extrema(values) + stats.update({'family': family, 'element': nodal[0]['element'], 'fe_order': nodal[0]['fe_order'], + 'label': label, 'alpha': alpha, 'roundoff_scale_estimate': tolerance, + 'negative_beyond_roundoff_count': sum(v < -tolerance for v in values), + 'limiting_point_h': limit['full_h'] + alpha*limit['full_dh']}) + details['nodal_checks'].append(stats) + if name.endswith('between_bounds') and details['replayed_candidate']: + r = details['replayed_candidate'] + details['replay_confirms_closure_outside_domain'] = all(r.get(k) == v for k, v in + [('admissible', 0), ('reason', 2), ('stage', 3), ('eos_code', 5)]) + details['candidate_derivation'] = json.loads((ROOT / 'supplemental_inputs/between_bounds_derivation.json').read_text()) + result['cases'].append(details) + if reference_rows is not None: + result['surface_root_indices'] = reference_rows + for c in cases: + folder = ROOT / 'captures/n1_J0p10_capture' / f"capture_{c['iteration']}" + base = vector_entries(folder / 'physical_base.txt', reference_rows) + candidate = vector_entries(folder / 'physical_candidate.txt', reference_rows) + if base is None or candidate is None: + continue + alpha = float(c['capture']['first_trial_alpha']) + dh = [(v-b)/alpha for b, v in zip(base, candidate)] + result['offline_surface_checks'].append({'iteration': c['iteration'], 'first_trial_alpha': alpha, + 'base_h': extrema(base), 'exact_candidate_h': extrema(candidate), + 'reconstructed_direction_dh': extrema(dh), + 'base_and_exact_candidate_entire_surface_identically_zero': all(v == 0 for v in base + candidate)}) + return result + + +def main(): + global ROOT + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('campaign_root', nargs='?', type=Path, default=ROOT, + help='Campaign artifact directory (default: directory containing this script).') + ROOT = parser.parse_args().campaign_root.expanduser().resolve() + if not ROOT.is_dir(): + parser.error(f'Campaign directory does not exist: {ROOT}') + cases = [case(i) for i in ITERATIONS] + transitions = [] + for first, second in zip(cases, cases[1:]): + if second['iteration'] != first['iteration'] + 1 or not first['bounds'] or not second['bounds']: + continue + alpha = first['live_newton_step']['alpha'] + for family, prefix in [('closure', ''), ('pressure', 'pressure_')]: + b, c = first['bounds'], second['bounds'] + same = all(b[prefix + key] == c[prefix + key] for key in ['limiting_element', 'limiting_point']) + predicted = b[prefix + 'limiting_h'] + alpha * b[prefix + 'limiting_dh'] + observed = c[prefix + 'limiting_h'] + transitions.append({'from_iteration': first['iteration'], 'to_iteration': second['iteration'], + 'family': family, 'same_limiting_point': same, 'accepted_alpha': alpha, + 'predicted_next_h': predicted, 'observed_next_h': observed, + 'absolute_disagreement': abs(predicted - observed) if same else None}) + supplements = supplemental(cases) + provenance_file = ROOT / 'focused/provenance.json' + source_revision = json.loads(provenance_file.read_text()).get('commit', '9e1b65d') if provenance_file.exists() else '9e1b65d' + report = {'schema_version': 2, 'source_revision': source_revision, + 'interpretation': 'Exact candidate replay is authoritative; direction probes include reconstruction roundoff.', + 'cases': cases, 'consecutive_captured_step_checks': transitions, 'supplemental': supplements} + polynomial_path = ROOT / 'pressure_polynomial_checks.json' + polynomial = json.loads(polynomial_path.read_text()) if polynomial_path.exists() else None + if polynomial: + report['offline_polynomial_artifacts'] = ['pressure_polynomial_checks.json', 'pressure_polynomial_analysis.md'] + (ROOT / 'pressure_checks.json').write_text(json.dumps(sanitize(report), indent=2, allow_nan=False) + '\n') + lines = ['# Pressure obstruction analysis', '', + 'Read-only analysis of the frozen n=1, J=.1 direction captures. This script launches no solver and does not modify case data. ' + f'The original focused source is {source_revision}; iteration numbers are zero-based. ' + 'Supplemental source revisions are recorded per case in pressure_checks.json. Partial cases are explicitly pending.', '', + '| Iteration | Exact first alpha | Pressure bound | Closure bound | Limiting pressure h | dh | Element / point | Exact replay |', + '|---:|---:|---:|---:|---:|---:|---|---|'] + for c in cases: + b = c['bounds'] + loc = f"{int(b['pressure_limiting_element'])} / {int(b['pressure_limiting_point'])}" if b else 'pending' + replay = ('pressure negative h' if c['exact_candidate_confirms_pressure_negative_h'] else + ('other; inspect JSON' if c['exact_candidate'] else 'pending')) + alpha = float(c['capture']['first_trial_alpha']) if c['capture'] else None + lines.append(f"| {c['iteration']} | {fmt(alpha)} | {fmt(b.get('pressure_alpha_h'))} | " + f"{fmt(b.get('closure_alpha_h'))} | {fmt(b.get('pressure_limiting_h'))} | " + f"{fmt(b.get('pressure_limiting_dh'))} | {loc} | {replay} |") + if all(c['exact_candidate_confirms_pressure_negative_h'] for c in cases): + lines += ['', '**All five exact saved first candidates reproduce finite negative enthalpy in pressure preparation.** ' + 'These failures are thermodynamic-domain rejections, not geometry rejections. The saved candidates are ' + 'replayed directly, so this conclusion does not depend on reconstructing the Newton direction.'] + populated = [c for c in cases if c['bounds']] + if populated: + lines += ['', 'The family attaining the smaller sampled bound is ' + ', '.join( + f"iteration {c['iteration']}: {c['tightest_sampled_family']}" for c in populated) + '.', '', + '| Iteration | Live accepted alpha | Accepted alpha / union bound |', '|---:|---:|---:|'] + for c in populated: + lines.append(f"| {c['iteration']} | {fmt(c['live_newton_step']['alpha'])} | {fmt(c['accepted_alpha_over_union_bound'])} |") + if cases[-1]['bounds'] and cases[0]['bounds']: + first, last = cases[0]['bounds'], cases[-1]['bounds'] + lines += ['', f"From iteration 2 to 11 the pressure bound shrinks by " + f"{first['pressure_alpha_h']/last['pressure_alpha_h']:.6g} times. The limiting pressure sample remains " + f"element {int(last['pressure_limiting_element'])}, point {int(last['pressure_limiting_point'])}; " + 'its h approaches zero while dh remains negative. ' + f"At iteration 11 the closure bound is {last['closure_alpha_h']/last['pressure_alpha_h']:.6g} times larger. " + 'The mechanism of the step-length collapse is established along these directions; why the Newton ' + 'directions continue to decrease this h remains unresolved.'] + if transitions: + maximum = max(r['absolute_disagreement'] for r in transitions if r['absolute_disagreement'] is not None) + lines += ['', f"For consecutive captured states at the same limiting point, h_next = h + alpha_accepted*dh " + f"agrees with the independently scanned next base to at most {maximum:.3g} absolute. " + 'The pointwise approach to the enthalpy boundary is therefore directly tied to the accepted steps.'] + lines += ['', '## Crossing controls', '', + '| Iteration | Candidate rows / expected | Shape-only all prepared | h-only crossing rejects negative h | Geometry rejections | Max h prediction/evaluation difference | Restored |', + '|---:|---:|---|---|---:|---:|---|'] + for c in cases: + if c['iteration'] not in (2, 11): + continue + lines.append(f"| {c['iteration']} | {c['pressure_trials_completed']} / {fmt(c['pressure_trials_expected'])} | " + f"{c.get('shape_only_all_prepared') if c.get('shape_only_all_prepared') is not None else 'pending'} | " + f"{c.get('enthalpy_only_crossing_failures_all_pressure_negative_h') if c.get('enthalpy_only_crossing_failures_all_pressure_negative_h') is not None else 'pending'} | " + f"{c.get('mapping_rejections', 'pending')} | {fmt(c.get('maximum_limiting_h_absolute_disagreement'))} | " + f"{c['restoration_completed']} |") + onset = cases[0] + if onset['bounds']: + b = onset['bounds'] + midpoint = onset['between_bounds_midpoint'] + midpoint_case = next(c for c in supplements['cases'] if c['iteration'] == 2) + measured_midpoint = midpoint_case.get('replay_confirms_closure_outside_domain') + midpoint_outcome = ('The supplemental candidate reproduces thermodynamic-domain rejection at the closure stage ' + '(reason=2, stage=3, EOS code=5), confirming pressure preparation passes at this alpha before closure fails.' + if measured_midpoint is True else + 'The supplemental replay has a different outcome; inspect pressure_checks.json.' + if measured_midpoint is False else + 'Its production outcome is pending; the supplemental replay can separate closure from pressure.') + lines += ['', '## Onset interpretation', '', + f"The onset has a tighter **{onset['tightest_sampled_family']}** bound: closure {b['closure_alpha_h']:.12g}, " + f"pressure {b['pressure_alpha_h']:.12g}. The original alpha=1 trial crosses both. " + 'Pressure is prepared before closure, so the recorded pressure-stage rejection identifies the first failing stage at that alpha; ' + 'it does not establish the pressure quadrature as the tightest restriction.', '', + f"A midpoint alpha={midpoint['alpha']:.12g} predicts closure h=" + f"{midpoint['closure_limiting_h_prediction']:.9g} and pressure h=" + f"{midpoint['pressure_limiting_h_prediction']:.9g} at the respective limiting points. " + 'This midpoint is outside the original probe grid. ' + midpoint_outcome, ''] + populated_limits = [c for c in cases if c['limits']] + if populated_limits: + lines += ['', '## Limiting-point locations', '', + '| Iteration | Family | Surface-adjacent element | Reference radius | Physical base (x, y, z) |', + '|---:|---|---|---:|---|'] + for c in populated_limits: + for r in c['limits']: + xyz = ', '.join(fmt(r[f'physical_{d}']) for d in 'xyz') + lines.append(f"| {c['iteration']} | {r['family']} | {int(r['element'])}; " + f"{bool(r['element_touches_stellar_surface'])} | {fmt(r['reference_radius'])} | ({xyz}) |") + lines += ['', '## Surface trace is an enforced equation', '', + 'The runtime isobaric condition is not absent merely because the field registry has no essential constraint. ' + 'The root constructs an enthalpy boundary map on StellarSurface (prepared_stellar_equilibrium.cpp:407), ' + 'replaces surface hydrostatic residual rows (lines 809–810), and replaces their Jacobian rows (904–908). ' + 'The compiled constraint is h−h(P_surface), with identity h Jacobian for fixed pressure ' + '(surface/compiled.cppm:47–59). At P_surface=0 the polytropic EOS gives h=0. ' + 'The radial seed explicitly sets these boundary coefficients to zero ' + '(seed/stellar_equilibrium_projection.cppm:203–208).', '', + 'This is algebraic enforcement, not a per-trial coefficient projection. Source inspection alone does not ' + 'establish exact zero in every captured iterate under an inexact linear solve. Even a zero polynomial ' + 'boundary trace does not guarantee nonnegative enthalpy within the adjacent high-order element. ' + 'The observed near-surface wall therefore does not establish a missing boundary condition or interpolation ' + 'as the cause of the Newton direction.'] + trace_cases = [c for c in supplements['cases'] if c['available']] + if trace_cases: + lines += ['', '## Supplemental surface carrier and nodal checks', '', + 'The actual surface carrier indices come from GetPressureSurfaceRows(), then select entries in each ' + 'saved physical state. Boundary-only fields retain those coefficients and zero every other enthalpy ' + 'coefficient; evaluating this field at a limiting interior point gives the boundary contribution to ' + 'that polynomial value, not a separate trace sample on the stellar surface.', '', + '| Iteration | Carrier rows | Base max abs h | Direction max abs dh | Unit-step max abs h |', + '|---:|---:|---:|---:|---:|'] + for c in trace_cases: + m = c['summary'] + lines.append(f"| {c['iteration']} | {fmt(m.get('surface_row_count'))} | {fmt(m.get('base_h_maximum_absolute'))} | " + f"{fmt(m.get('direction_dh_maximum_absolute'))} | {fmt(m.get('h_plus_dh_maximum_absolute'))} |") + offline = supplements['offline_surface_checks'] + if len(offline) == len(ITERATIONS) and all(c['base_and_exact_candidate_entire_surface_identically_zero'] for c in offline): + lines += ['', '**All actual surface enthalpy coefficients are exactly zero in all five saved bases and all ' + 'five exact first candidates.** Their reconstructed boundary directions are therefore also zero. ' + 'Boundary-row leakage is ruled out for these captured directions, independently of the inexact ' + 'linear-solve possibility from source inspection.'] + else: + lines += ['', 'Offline surface checks (including nonzero and negative counts) are recorded per capture in ' + 'pressure_checks.json; inspect them before inferring exact trace enforcement in captured states.'] + lines += ['', '| Iteration | Family | Full h | Boundary contribution h | Full dh | Boundary contribution dh |', + '|---:|---|---:|---:|---:|---:|'] + for c in trace_cases: + for r in c['contributions']: + lines.append(f"| {c['iteration']} | {r['family']} | {fmt(r['full_h'])} | {fmt(r['boundary_h'])} | " + f"{fmt(r['full_dh'])} | {fmt(r['boundary_dh'])} |") + lines += ['', 'Full physical fields are scattered through the production field adapter, then evaluated at FE ' + 'GetNodes() coordinates. These values do not assume that raw signed local DOF indices are nodal values. ' + 'The same element is intentionally exported separately for each limiting rule family.', '', + '| Iteration | Family | Step | Nodal minimum h | Negative nodes | Limiting interior-point h |', + '|---:|---|---|---:|---:|---:|'] + for c in trace_cases: + for r in c['nodal_checks']: + if r['label'] in ('base', 'near_above_union', 'unit_step'): + lines.append(f"| {c['iteration']} | {r['family']} | {r['label']} | {fmt(r['minimum'])} | " + f"{r['negative_count']} | {fmt(r['limiting_point_h'])} |") + undershoots = [(c['iteration'], r['family']) for c in trace_cases for r in c['nodal_checks'] + if r['label'] == 'near_above_union' and r['minimum'] >= 0 and r['limiting_point_h'] < 0] + if undershoots: + lines += ['', 'At 1.01 times the union bound, nonnegative values at every FE node coexist with negative ' + 'enthalpy at the limiting interior point in ' + ', '.join(f'iteration {i} {family}' for i, family in undershoots) + + '. This demonstrates an interior polynomial undershoot in these trial directions. It does not ' + 'establish an interpolation implementation error or explain why Newton selects this direction.'] + lines += ['', 'Nodal nonnegativity is not a certificate of positivity between nodes. The JSON includes strict ' + 'negative counts and a separate 64-epsilon times nodal-value-scale estimate, which is descriptive ' + 'rather than a certified rounding bound. Equality exactly at the union bound is deliberately not ' + 'used to decide a crossing. The production quadrature union bound also does not certify positivity ' + 'at unsampled points; these results do not establish the minimum enthalpy over an entire element.'] + if polynomial: + lines += ['', '## Offline interior polynomial check', '', + 'A separate independently validated tensor-Q3 reconstruction samples 495,264 strictly interior ' + 'reference points in each exported limiting element, including logarithmic distances to its zero ' + 'face. It reconstructs all 64 nodes exactly and both production limit points to roundoff.'] + for c in polynomial['cases']: + minimum = c['sampled_extrema']['base']['minimum'] + if minimum < -c['roundoff_scale_estimate']: + point = ', '.join(fmt(x) for x in c['sampled_extrema']['base']['minimum_reference_point']) + lines += ['', f"**The saved iteration {c['iteration']} base already has negative enthalpy between " + f"production quadrature points:** sampled h={minimum:.9g} at reference ({point}). " + 'This base-state finding is independent of Newton-direction reconstruction. The surface ' + 'coefficients and all 64 FE node values remain nonnegative; production quadrature positivity ' + 'therefore misses a real polynomial undershoot in this accepted state.'] + else: + lines += ['', f"Iteration {c['iteration']} has no negative base samples on this grid; its sampled " + f"directional bound {c['directional_bound']['sampled_admissible_step_bound']:.12g} is tighter " + f"than the production bound {c['union_production_bound']:.12g}. Absence of negative samples " + 'does not certify positivity everywhere.'] + lines += ['', 'See pressure_polynomial_analysis.md and pressure_polynomial_checks.json for interpolation ' + 'validation, minima, locations and limitations. No solver or preparation rule was changed.'] + lines += ['', '## Scope and interpretation', '', + 'A replay with reason=2, stage=6, EOS code=5 is the production pressure operator’s finite h<0 branch. ' + 'It resolves the old CSV ambiguity with inverted geometry. Missing determinant/F data on rejection are not measured nonfinite residuals.', '', + 'The bound scans the exact production pressure/closure quadrature points using pointwise polynomial contractions. ' + 'The independently re-evaluated h values in the crossing table concern the two frozen limiting points, not a full search for each trial’s most negative point. ' + 'Optimized tensor contractions can differ at roundoff; .99/1.01 factors deliberately avoid equality at the crossing.', '', + 'Element-to-surface adjacency, when available in pressure_obstruction_limits.csv, is an exact topology flag. ' + 'It is not distance to the surface. Physical/reference Cartesian coordinates and element-local integration coordinates are distinct.', '', + 'At iteration 11 the supplied direction is reconstructed from a first-trial difference at alpha=1.9073486328125e-6; ' + 'its relative subtraction-roundoff scale estimate is approximately 4.88e-7, not a certified error bound. ' + 'The exact saved candidate is replayed without that reconstruction and is authoritative for the observed rejection.', '', + 'Preparation success in a frozen probe does not establish line-search merit acceptance, physical equilibrium, or eventual Newton convergence. ' + 'No requested stopping test, enthalpy admissibility rule, or solver was changed. ' + 'The evolving bound is directional: a falling h/dh bound does not by itself identify why the nonlinear correction keeps pointing outward from the admissible set.', '', + 'Machine-readable evidence and all individual trial checks: pressure_checks.json. ' + 'Reproduce this analysis with `python3 analyze_pressure.py /path/to/campaign_artifacts`; ' + 'the optional campaign path defaults to the script directory.'] + (ROOT / 'pressure_analysis.md').write_text('\n'.join(lines) + '\n') + print(json.dumps({'cases_with_exact_replay': sum(c['exact_candidate'] is not None for c in cases), + 'cases_fully_restored': sum(c['restoration_completed'] for c in cases), + 'supplemental_available': sum(c['available'] for c in supplements['cases']), + 'supplemental_fully_restored': sum('final_restoration_residual_difference' in c.get('restoration', {}) + for c in supplements['cases']), + 'candidate_rows': {c['iteration']: c['pressure_trials_completed'] for c in cases}}, sort_keys=True)) + + +if __name__ == '__main__': + main() diff --git a/experiments/plateau_analysis/analyze_shape.py b/experiments/plateau_analysis/analyze_shape.py new file mode 100644 index 0000000..76a4349 --- /dev/null +++ b/experiments/plateau_analysis/analyze_shape.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python3 +"""Read-only shape CSV audit. No solver/library dependency; outputs live elsewhere than raw case files.""" +from __future__ import annotations +import argparse, csv, hashlib, json, math +from pathlib import Path +from collections import defaultdict + +CASES = ('n1_J0p10_iteration_2', 'n1p5_J0_iteration_1', 'n3_J0_seed_R160_curvature') + +def rows(path): + with path.open(newline='') as stream: + reader = csv.DictReader(stream) + result = list(reader) + for r in result: + if None in r or any(v is None for v in r.values()): + raise ValueError(f'Malformed/incomplete CSV row in {path}: {r}') + return result + +def num(r, key): return float(r[key]) +def flag(r, key): return r[key] == '1' +def norm(values): return math.sqrt(math.fsum(v*v for v in values)) +def ratio(a,b): return a/b if math.isfinite(a) and math.isfinite(b) and b > 0 else math.nan + +def clean(obj): + if isinstance(obj, float) and not math.isfinite(obj): return None + if isinstance(obj, dict): return {str(k):clean(v) for k,v in obj.items()} + if isinstance(obj, (tuple,list)): return [clean(v) for v in obj] + return obj + +def signed_aggregate(rs): + out = {} + for kind in ('full','shape','fields','mixed'): + out[kind+'_remainder_norm'] = norm(num(r,kind+'_remainder_norm') for r in rs) + out['linear_dot_'+kind] = math.fsum(num(r,'linear_dot_'+kind) for r in rs) + for name in ('linear_prediction_norm','full_trial_norm','heuristic_roundoff_scale'): + out[name] = norm(num(r,name) for r in rs) + for name in ('merit_identity_error','projection_decomposition_error'): + out[name] = math.fsum(num(r,name) for r in rs) + for a,b in (('full','shape'),('full','fields'),('full','mixed'),('shape','fields'),('shape','mixed'),('fields','mixed')): + dot = 0. + for r in rs: + na,nb = num(r,a+'_remainder_norm'),num(r,b+'_remainder_norm') + if na == 0 or nb == 0: continue + dot += na*nb*num(r,a+'_'+b+'_cosine') + out[a+'_'+b+'_cosine'] = ratio(dot, out[a+'_remainder_norm']*out[b+'_remainder_norm']) + out['linear_curvature_cross_term'] = 2*out['linear_dot_full'] + out['curvature_squared_term'] = out['full_remainder_norm']**2 + out['merit_curvature_total'] = out['linear_curvature_cross_term']+out['curvature_squared_term'] + return out + +def analyze_case(case): + summary_path = case/'shape_curvature_summary.csv' + if not summary_path.exists(): return {'status':'pending'} + summary = {r['metric']:num(r,'value') for r in rows(summary_path)} + required = ('preparations','accepted_normalized_state_change','accepted_residual_change','accepted_physical_state_change') + if not all(k in summary for k in required): return {'status':'incomplete', 'summary':summary} + names = ('preparations','restorations','directions','actions','bounds','remainders','mixed') + data = {name:rows(case/f'shape_curvature_{name}.csv') for name in names} + result = {'status':'complete','summary':summary,'sha256':{p.name:hashlib.sha256(p.read_bytes()).hexdigest() for p in case.glob('shape_curvature_*.csv')}} + input_checks=case/'frozen_input_checks.csv' + if input_checks.exists(): result['input_checks']={r['metric']:num(r,'value') for r in rows(input_checks)} + replay=case/'replayed_candidate.csv' + if replay.exists(): result['replayed_candidate']=rows(replay) + checks = result['checks'] = {} + checks['accepted_vectors_unchanged'] = all(summary[k] == 0 for k in required[1:]) + checks['direction_split_exact'] = summary['direction_split_error_norm'] == 0 + checks['action_split_within_recorded_heuristic'] = summary['action_linearity_within_heuristic'] == 1 + checks['preparation_count_matches_events'] = summary['preparations'] == len(data['preparations']) + checks['all_preparations_admissible'] = all(flag(r,'admissible') for r in data['preparations']) + result['rejections'] = [r for r in data['preparations'] if not flag(r,'admissible')] + checks['all_state_increments_resolved'] = all(flag(r,'state_increment_resolved') for r in data['remainders']) + checks['all_alpha_squared_valid'] = all(flag(r,'alpha_squared_valid') for r in data['remainders']) + result['max_restoration_difference'] = max((num(r,'restoration_difference_norm') for r in data['restorations']),default=0) + result['actions'] = data['actions'] + result['directions'] = data['directions'] + drift = defaultdict(float) + for r in data['restorations']: drift[r['block']] = max(drift[r['block']],num(r,'restoration_difference_norm')) + for r in data['remainders']: + r['effective_floor'] = max(num(r,'heuristic_roundoff_scale'),8*drift[r['block']]) + r['effective_floor_ratio'] = ratio(num(r,'remainder_norm'),r['effective_floor']) + r['robustly_above_floor'] = flag(r,'admissible') and flag(r,'state_increment_resolved') and flag(r,'alpha_squared_valid') and r['effective_floor_ratio'] >= 10 + gravity = lambda r: r['block'] in ('gravity_gradient_relation', 'poisson_balance') + affine = [r for r in data['remainders'] if r['direction']=='fields' and gravity(r)] + result['fixed_shape_gravity'] = { + 'blocks':sorted({r['block'] for r in affine}), + 'max_remainder':max((num(r,'remainder_norm') for r in affine),default=math.nan), + 'max_effective_floor_ratio':max((r['effective_floor_ratio'] for r in affine),default=math.nan), + 'all_below_recorded_or_restoration_floor': bool(affine) and all(r['effective_floor_ratio'] <= 1 for r in affine), + 'rows':affine, + } + by_direction_block = defaultdict(list) + for r in data['remainders']: by_direction_block[r['direction'],r['block']].append(r) + slopes = [] + for (direction,block),rs in by_direction_block.items(): + rs.sort(key=lambda r:num(r,'alpha'),reverse=True) + for high,low in zip(rs,rs[1:]): + if not high['robustly_above_floor'] or not low['robustly_above_floor']: continue + slope = math.log(num(high,'remainder_norm')/num(low,'remainder_norm'))/math.log(num(high,'alpha')/num(low,'alpha')) + slopes.append({'direction':direction,'block':block,'alpha_high':num(high,'alpha'),'alpha_low':num(low,'alpha'),'slope':slope,'high_floor_ratio':high['effective_floor_ratio'],'low_floor_ratio':low['effective_floor_ratio']}) + result['resolved_slopes'] = slopes + by_alpha = defaultdict(list) + mixed_by_block = defaultdict(list) + block_identity_ratios, block_projection_ratios = [], [] + for r in data['mixed']: + by_alpha[num(r,'alpha')].append(r) + mixed_by_block[r['block']].append(r) + if not flag(r,'all_trials_admissible'): continue + den = num(r,'full_trial_norm')**2 + num(r,'linear_prediction_norm')**2 + 2*abs(num(r,'linear_dot_full')) + num(r,'full_remainder_norm')**2 + block_identity_ratios.append(ratio(abs(num(r,'merit_identity_error')),den)) + pden = math.fsum(abs(num(r,'linear_dot_'+k)) for k in ('full','shape','fields','mixed')) + if pden > 0: block_projection_ratios.append(ratio(abs(num(r,'projection_decomposition_error')),pden)) + for block,rs in mixed_by_block.items(): + rs.sort(key=lambda r:num(r,'alpha'),reverse=True) + for high,low in zip(rs,rs[1:]): + floor_high=max(num(high,'heuristic_roundoff_scale'),24*drift[block]) + floor_low=max(num(low,'heuristic_roundoff_scale'),24*drift[block]) + def robust(r,floor): + return all(flag(r,key) for key in ('all_trials_admissible','all_state_increments_resolved','alpha_squared_valid')) and num(r,'mixed_remainder_norm') >= 10*floor + if not robust(high,floor_high) or not robust(low,floor_low): continue + slopes.append({'direction':'mixed','block':block,'alpha_high':num(high,'alpha'),'alpha_low':num(low,'alpha'),'slope':math.log(num(high,'mixed_remainder_norm')/num(low,'mixed_remainder_norm'))/math.log(num(high,'alpha')/num(low,'alpha')),'high_floor_ratio':ratio(num(high,'mixed_remainder_norm'),floor_high),'low_floor_ratio':ratio(num(low,'mixed_remainder_norm'),floor_low)}) + aggregates = [] + identity_ratios,projection_ratios = [],[] + for alpha,rs in sorted(by_alpha.items(),reverse=True): + item = {'alpha':alpha,'all_admissible':all(flag(r,'all_trials_admissible') for r in rs),'all_resolved':all(flag(r,'all_state_increments_resolved') for r in rs)} + if item['all_admissible']: + item.update(signed_aggregate(rs)) + item['base_residual_norm'] = summary['base_residual_norm'] + item['trial_over_base'] = ratio(item['full_trial_norm'],summary['base_residual_norm']) + item['linear_over_base'] = ratio(item['linear_prediction_norm'],summary['base_residual_norm']) + item['full_remainder_over_alpha_squared'] = ratio(item['full_remainder_norm'],alpha**2) + full_rows = [r for r in data['remainders'] if r['direction']=='full' and num(r,'alpha') == alpha] + item['full_effective_floor_norm'] = norm(r['effective_floor'] for r in full_rows) + item['full_effective_floor_ratio'] = ratio(item['full_remainder_norm'],item['full_effective_floor_norm']) + item['gravity'] = signed_aggregate([r for r in rs if gravity(r)]) + item['gravity_full_remainder_squared_fraction'] = ratio(item['gravity']['full_remainder_norm']**2,item['full_remainder_norm']**2) + item['blocks'] = rs + den = item['full_trial_norm']**2+item['linear_prediction_norm']**2+2*abs(item['linear_dot_full'])+item['full_remainder_norm']**2 + item['merit_identity_relative_error'] = ratio(abs(item['merit_identity_error']),den) + pden = math.fsum(abs(item['linear_dot_'+k]) for k in ('full','shape','fields','mixed')) + item['projection_decomposition_relative_error'] = ratio(abs(item['projection_decomposition_error']),pden) + identity_ratios.append(item['merit_identity_relative_error']) + projection_ratios.append(item['projection_decomposition_relative_error']) + aggregates.append(item) + result['alphas'] = aggregates + result['max_merit_identity_relative_error'] = max(identity_ratios,default=math.nan) + result['max_projection_decomposition_relative_error'] = max(projection_ratios,default=math.nan) + result['max_block_merit_identity_relative_error'] = max(block_identity_ratios,default=math.nan) + result['max_block_projection_decomposition_relative_error'] = max(block_projection_ratios,default=math.nan) + result['remainder_rows'] = data['remainders'] + return result + +def make_markdown(results,provenance): + lines = ['# Shape-curvature analysis','',f"Focused executable source `{provenance.get('commit','unknown')}`; binary SHA-256 `{provenance.get('binary_sha256','unknown')}`.",'', + 'This analysis reads the saved full, shape-only, and fields-only residual probes. Norms and projections use the existing normalization. The mixed remainder is a vector subtraction before its norm. No solver was run by this analysis.','', + 'Source expectation: gravity has the form Fg=M(q)g+BᵀPhi and Fp=B g−S(q)rho. With q fixed, it is affine in all remaining fields; shape dependence and products of shape changes with field changes can supply second-order remainders. See [residual assembly](/Users/tboudreaux/Programming/4DSTAR/MeanField/libmeanfield/impl/operators/gravity_field.cpp:415) and [Jacobian shape/source actions](/Users/tboudreaux/Programming/4DSTAR/MeanField/libmeanfield/impl/operators/gravity_field_jacobian.cpp:242).','', + 'Arithmetic checks use the reported preparation history and the final restoration drift. Slopes are reported only when both adjacent residual remainders exceed the larger of the recorded floor and eight times final restoration drift by at least 10×, and both increments are resolved. These floors are heuristic, not certified assembly-error bounds.',''] + lines += ['## Common-cap comparison','', + 'The signed columns below divide squared-merit contributions by each case’s own ‖F0‖². Their sum equals (‖Ftrial‖²−‖L‖²)/‖F0‖². They are evaluated at each direction’s common admissible cap; these are generally different alpha values.','', + '| case | alpha | ‖Ftrial‖/‖F0‖ | 2/F0² | 2/F0² | 2/F0² | ‖Rfull‖²/F0² | total curvature/F0² |', + '|---|---:|---:|---:|---:|---:|---:|---:|'] + for name,result in results.items(): + if result['status'] != 'complete': + lines += [f"| {name} | {result['status']} | | | | | | |"] + continue + if not result['alphas'] or not result['alphas'][0]['all_admissible']: + lines += [f"| {name} | no complete common-cap triple | | | | | | |"] + continue + a=result['alphas'][0]; den=result['summary']['base_residual_norm']**2 + lines += [f"| {name} | {a['alpha']:.7g} | {a['trial_over_base']:.8g} | {2*a['linear_dot_shape']/den:.7g} | {2*a['linear_dot_fields']/den:.7g} | {2*a['linear_dot_mixed']/den:.7g} | {a['curvature_squared_term']/den:.7g} | {a['merit_curvature_total']/den:.7g} |"] + lines += [''] + if all(r['status']=='complete' for r in results.values()): + lines += ['Across these three saved directions, every one of the shape helper’s 66 preparations is admissible and every accepted-vector change and restoration drift is exactly zero. The coupled gravity Taylor remainders have resolved slopes near two, while fields-only gravity remainders stay at approximately 2.5e-14–4.0e-14. This supports geometry–field curvature along the tested directions; no first-derivative inconsistency is visible in the resolved regime. The numerical floor is heuristic: the n1.5 and n3 fields-only maxima exceed it by factors 1.04 and 1.39, which is recorded explicitly.','', + 'The practical mechanisms differ. The rotating n1 direction is blocked at its original full step by pressure admissibility; its sampled safe interval decreases the residual. The n1.5 full step is admissible and is strongly rejected by curvature in the residual metric. The n3 saved direction encounters both an enthalpy restriction and a curvature penalty at the largest tested safe step, but smaller sampled steps improve the residual. Its originating linear solve also missed the 3% gate.',''] + lines += ['## First-order gravity cancellation','', + 'These are cancellations between the two parts of each selected correction. They are not condition-number or symmetry measurements.','', + '| case | block | norm(Jshape) | norm(Jfields) | norm(Jfull) | shape/fields cosine |', + '|---|---|---:|---:|---:|---:|'] + for name,result in results.items(): + if result['status'] != 'complete': continue + for r in result['actions']: + if r['block'] not in ('gravity_gradient_relation','poisson_balance'): continue + lines += [f"| {name} | {r['block']} | {num(r,'shape_action_norm'):.8g} | {num(r,'fields_action_norm'):.8g} | {num(r,'full_action_norm'):.8g} | {num(r,'shape_fields_cosine'):.12g} |"] + lines += [''] + for name,result in results.items(): + lines += ['## '+name,''] + if result['status'] != 'complete': + lines += ['Status: '+result['status']+'.','']; continue + s=result['summary']; checks=result['checks']; aff=result['fixed_shape_gravity'] + if name=='n1_J0p10_iteration_2': + lines += ['**Interpretation:** the sampled admissible interval improves the residual, although mixed curvature makes it worse than the linear prediction. The original full-step pressure rejection is a separate admissibility event. The largest-step mixed signed projection is adverse, while the shape-only signed projection is favorable. The last alpha is limited by residual arithmetic and does not diagnose loss of second-order consistency.',''] + elif name=='n1p5_J0_iteration_1': + lines += ['**Interpretation:** this is an admissible, curvature-driven merit rejection. At alpha=1 the linear prediction almost removes the residual, while the actual residual grows by over 12×. Alpha=.1 remains merit-worse; alpha=.01 improves it. Gravity shape and mixed remainders dominate despite partial cancellation. Their squared remainder term dominates the full-step merit cost; at alpha=.1 both shape and mixed signed projections are adverse. The fields-only gravity maximum exceeds the chosen heuristic floor by only 4.2%, remaining at approximately 3e-14 absolute; that minor heuristic exceedance is not evidence of nonlinear fixed-shape gravity.',''] + elif name=='n3_J0_seed_R160_curvature': + lines += ['**Interpretation:** the largest common safe step is admissible but merit-worse; the next decade improves the residual. The mixed remainder is the largest component and is strongly aligned with the full remainder. Its signed projection onto the linear prediction is nevertheless favorable at the largest alpha, as is the fields-only projection: the adverse total curvature cost comes from the squared remainder norm overwhelming those favorable cross terms. At smaller resolved alphas the net curvature slightly improves the linear prediction. The saved R160 direction retains true linear defect 0.0402590213 and had failed the 0.03 gate; this probe is not a live Newton convergence result. Fields-only gravity stays below 4e-14 absolute, up to 1.393× the heuristic floor.',''] + lines += [f"Frozen ‖F‖={s['base_residual_norm']:.9g}; ‖p‖={s['full_direction_norm']:.9g}; common alpha cap={s['common_alpha_cap']:.9g}. {int(s['preparations'])} preparations. Accepted-vector changes are all zero: **{checks['accepted_vectors_unchanged']}**. Direction split exact: **{checks['direction_split_exact']}**. Jacobian split error {s['action_split_error_norm']:.3e}, recorded scale {s['action_linearity_heuristic_scale']:.3e}; maximum per-block restoration drift {result['max_restoration_difference']:.3e}.",'', + f"All preparations admissible: **{checks['all_preparations_admissible']}**; all increments resolved: **{checks['all_state_increments_resolved']}**. Fields-only gravity maximum remainder {aff['max_remainder']:.3e}, maximum ratio to the effective heuristic floor {aff['max_effective_floor_ratio']:.3e}. Every sampled fields-only gravity remainder lies below that floor: **{aff['all_below_recorded_or_restoration_floor']}**.",'', + f"Maximum aggregate merit-identity relative error {result['max_merit_identity_relative_error']:.3e}; projection-decomposition relative error {result['max_projection_decomposition_relative_error']:.3e}.",'', + '| alpha | ‖Ftrial‖/‖F0‖ | ‖L‖/‖F0‖ | ‖Rfull‖/alpha² | full R / floor | ‖Rs‖ | ‖Rf‖ | ‖Rm‖ | gravity share of ‖Rfull‖² |','|---:|---:|---:|---:|---:|---:|---:|---:|---:|'] + for a in result['alphas']: + if not a['all_admissible']: lines += [f"| {a['alpha']:.6g} | inadmissible | | | | | | | |"] + else: lines += [f"| {a['alpha']:.6g} | {a['trial_over_base']:.8g} | {a['linear_over_base']:.8g} | {a['full_remainder_over_alpha_squared']:.8g} | {a['full_effective_floor_ratio']:.3g} | {a['shape_remainder_norm']:.4e} | {a['fields_remainder_norm']:.4e} | {a['mixed_remainder_norm']:.4e} | {a['gravity_full_remainder_squared_fraction']:.6%} |"] + lines += ['', '| alpha | 2 | 2 | 2 | ‖Rfull‖² | curvature contribution to squared merit |','|---:|---:|---:|---:|---:|---:|'] + for a in result['alphas']: + if a['all_admissible']: lines += [f"| {a['alpha']:.6g} | {2*a['linear_dot_shape']:.7e} | {2*a['linear_dot_fields']:.7e} | {2*a['linear_dot_mixed']:.7e} | {a['curvature_squared_term']:.7e} | {a['merit_curvature_total']:.7e} |"] + lines += ['', 'Positive signed terms increase the squared residual relative to the linear prediction; negative terms decrease it. The quadratic remainder-norm term contains component cross terms, so the signed projections alone are not a unique additive attribution of the entire curvature cost.','', + '| direction | residual block | resolved intervals | slope range | smallest resolved alpha |','|---|---|---:|---:|---:|'] + grouped_slopes=defaultdict(list) + for p in result['resolved_slopes']: grouped_slopes[p['direction'],p['block']].append(p) + for (direction,block),ps in grouped_slopes.items(): + lines += [f"| {direction} | {block} | {len(ps)} | {min(p['slope'] for p in ps):.6f}–{max(p['slope'] for p in ps):.6f} | {min(p['alpha_low'] for p in ps):.3g} |"] + lines += [''] + lines += ['## Discriminating next experiment','', + 'Where the full Jacobian action is small because shape and field actions cancel, suppressing the shape component alone would destroy that first-order balance. A targeted hypothesis is that the current straight Newton trial loses the geometry–gravity balance at second order. Test it by holding each saved trial’s shape, density, and other material/scalar coordinates fixed, solving only its linear gravity subsystem to a documented true defect, and comparing the completed trial against the original one. Record all residual blocks, the new hydrostatic/force response, interior diagnostics and virial; lowering the gravity residual alone is insufficient. This is a proposed diagnostic, not a demonstrated solver improvement or a recommendation to promote a new policy.','', + 'The completed comparison should preserve the original direction and trial as a control, disclose solve cost, and separate the pressure-admissibility obstruction from residual-metric rejection. A linear one-way preconditioner’s previous failure does not test this nonlinear conditional completion. These probes do not justify clipping enthalpy or deleting shape coordinates.','', + 'Before promoting such a change, localize the important remainder and signed-merit contributions by material and exterior element assembly, preserving their shared-DOF cross terms. The existing aggregate gravity norms cannot tell whether an apparently large improvement benefits the interior or mainly removes exterior error. Interior consistency and virial ≤1e-5 remain the promotion criteria.','', + '## Interpretation limits','', + 'A resolved second-order remainder, first-order decay of its difference-quotient error, and arithmetic-level fields-only gravity remainder support geometric/field curvature along these particular directions. They do not prove every Jacobian column correct. Jacobian additivity tests linearity and restoration, not agreement with the nonlinear derivative by itself.','', + 'These residual-block aggregates do not separate interior from exterior elements. They neither locate a physical error nor identify a symmetry or near-null mode. A large correction or curvature may coexist with a weak normalized direction, but this split is not a condition-number measurement. No new model was accepted and these probes do not establish a virial milestone.',''] + return '\n'.join(lines) + +def main(): + ap=argparse.ArgumentParser(description=__doc__) + ap.add_argument('--campaign',type=Path,default=Path(__file__).resolve().parent, + help='Campaign root containing focused/; writes shape_checks.json and shape_analysis.md here. Defaults to the script directory.') + args=ap.parse_args() + campaign=args.campaign; focused=campaign/'focused' + provenance=json.loads((focused/'provenance.json').read_text()) + results={case:analyze_case(focused/case) for case in CASES} + out={'schema_version':1,'source_commit':provenance.get('commit'),'binary_sha256':provenance.get('binary_sha256'),'cases':results} + (campaign/'shape_checks.json').write_text(json.dumps(clean(out),indent=2,allow_nan=False)+'\n') + (campaign/'shape_analysis.md').write_text(make_markdown(results,provenance)) + for name,r in results.items(): print(name,r['status']) + +if __name__=='__main__':main() diff --git a/experiments/plateau_analysis/audit_focused_outputs.py b/experiments/plateau_analysis/audit_focused_outputs.py new file mode 100644 index 0000000..a369dc4 --- /dev/null +++ b/experiments/plateau_analysis/audit_focused_outputs.py @@ -0,0 +1,143 @@ +"""Audit the pressure/shape campaign's completed outputs without running solvers. + +The sole argument is the campaign directory containing captures/ and focused/. +The script reads existing results and writes focused_audit_checks.json there. +Unfinished cases are reported as pending; this is an observation audit, not a +stellar convergence or physical acceptance gate. +""" +import argparse +import csv +import hashlib +import json +import math +from pathlib import Path + +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument("campaign_root", type=Path) +ROOT = parser.parse_args().campaign_root.resolve() +FOCUSED = ROOT / "focused" + + +def rows(path): + with path.open() as source: + result = list(csv.DictReader(source)) + if any(None in r or any(v is None for v in r.values()) for r in result): + raise ValueError(f"Malformed CSV row: {path}") + return result + + +def metrics(path): + data = rows(path) + if len({r["metric"] for r in data}) != len(data): + raise ValueError(f"Duplicate metric: {path}") + return {r["metric"]: float(r["value"]) for r in data} + + +def sha(path): + return hashlib.sha256(path.read_bytes()).hexdigest() + + +def clean(value): + if isinstance(value, dict): + return {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 + + +provenance = json.loads((FOCUSED / "provenance.json").read_text()) +result = { + "source_commit": provenance["commit"], + "binary_hash_matches": sha(FOCUSED / "plateau_diagnostic") == provenance["binary_sha256"], + "completed": [], + "pending": [], +} +for case in json.loads((FOCUSED / "cases.json").read_text()): + case_id = case["id"] + manifest = FOCUSED / f"{case_id}_run.json" + if not manifest.exists(): + result["pending"].append(case_id) + continue + run = json.loads(manifest.read_text()) + if "exit_code" not in run: + result["pending"].append(case_id) + continue + folder = FOCUSED / case_id + data = { + "id": case_id, + "exit_code": run["exit_code"], + "seconds": run["seconds"], + "input_hash_count": len(run["input_state_hashes"]), + "input_hash_mismatches": [str(p) for p, h in run["input_state_hashes"].items() if sha(Path(p)) != h], + "output_hash_count": len(run["output_hashes"]), + "output_hash_mismatches": [p for p, h in run["output_hashes"].items() if sha(folder / p) != h], + "unhashed_data_files": sorted(str(p.relative_to(folder)) for p in folder.rglob("*") if p.is_file() and p.suffix in (".txt", ".csv", ".json") and str(p.relative_to(folder)) not in run["output_hashes"]), + "csv_rows": {p.name: len(rows(p)) for p in sorted(folder.glob("*.csv"))}, + "frozen": metrics(folder / "frozen_input_checks.csv"), + } + args = case["arguments"] + direction = Path(args[args.index("--probe-direction") + 1]) + physical = Path(args[args.index("--probe-state") + 1]) + data["direction_copy_byte_identical"] = direction.read_bytes() == (folder / "probe_normalized_correction.txt").read_bytes() + data["physical_base_copy_byte_identical"] = physical.read_bytes() == (folder / "frozen_physical_base.txt").read_bytes() + if "iteration_" in case_id: + iteration = int(case_id.rsplit("_", 1)[1]) + capture_root = ROOT / "captures" / ("n1p5_J0_capture" if case_id.startswith("n1p5") else "n1_J0p10_capture") + capture_manifest = json.loads((capture_root.parent / f"{capture_root.name}_run.json").read_text()) + data["auxiliary_capture_input_hashes_match"] = all(sha(physical.parent / name) == capture_manifest["output_hashes"][str((physical.parent / name).relative_to(capture_root))] for name in ("normalized_base.txt", "normalized_residual.txt")) + old = next(r for r in rows(capture_root / "newton.csv") if int(r["iteration"]) == iteration) + data["historical_linear_defect_relative"] = float(old["linear_true"]) + trial = next(r for r in rows(capture_root / "trials.csv") if int(r["iteration"]) == iteration and r["trial"] == "0") + data["original_first_trial"] = trial + replay = rows(folder / "replayed_candidate.csv")[0] + data["candidate_replay"] = replay + if trial["disposition"] == "4": + data["replayed_candidate_residual_exact"] = replay["admissible"] == "1" and replay["residual_norm"] == trial["F"] + data["replayed_candidate_residual_difference"] = float(replay["residual_norm"]) - float(trial["F"]) + data["replayed_candidate_minimum_jacobian_exact"] = replay["minimum_jacobian"] == trial["min_jacobian"] + else: + previous = physical.parent + data["historical_linear_defect_relative"] = float(next(r for r in rows(previous / "frozen_linear.csv") if r["restart"] == "160")["true_relative"]) + data["linear_defect_relative_difference"] = data["frozen"]["linear_defect_relative"] - data["historical_linear_defect_relative"] + layout = rows(folder / "layout.csv") + residual_blocks = {r["block"] for r in layout if r["kind"] == "residual"} + value_blocks = {r["block"] for r in layout if r["kind"] == "value"} + nr = sum(r["kind"] == "residual" for r in layout) + nv = sum(r["kind"] == "value" for r in layout) + if "--shape-curvature" in args: + data["shape_summary"] = metrics(folder / "shape_curvature_summary.csv") + expected = {"shape_curvature_directions.csv": 3 * nv, "shape_curvature_actions.csv": nr, "shape_curvature_preparations.csv": 22, "shape_curvature_restorations.csv": 4 * nr, "shape_curvature_remainders.csv": 18 * nr, "shape_curvature_mixed.csv": 6 * nr} + data["shape_row_count_mismatches"] = {f: [data["csv_rows"][f], count] for f, count in expected.items() if data["csv_rows"][f] != count} + restorations = rows(folder / "shape_curvature_restorations.csv") + data["shape_max_restoration_difference"] = max(float(r["restoration_difference_norm"]) for r in restorations) + mixed = rows(folder / "shape_curvature_mixed.csv") + alphas = {r["alpha"] for r in mixed} + remainders = rows(folder / "shape_curvature_remainders.csv") + data["shape_six_alphas"] = len(alphas) == 6 + data["shape_mixed_keys_complete"] = {(r["alpha"], r["block"]) for r in mixed} == {(a, b) for a in alphas for b in residual_blocks} + data["shape_remainder_keys_complete"] = {(r["alpha"], r["direction"], r["block"]) for r in remainders} == {(a, d, b) for a in alphas for d in ("full", "shape", "fields") for b in residual_blocks} + data["shape_direction_keys_complete"] = {(r["direction"], r["value_block"]) for r in rows(folder / "shape_curvature_directions.csv")} == {(d, b) for d in ("full", "shape", "fields") for b in value_blocks} + data["shape_max_absolute_merit_identity_error"] = max(abs(float(r["merit_identity_error"])) for r in mixed) + data["shape_max_absolute_projection_decomposition_error"] = max(abs(float(r["projection_decomposition_error"])) for r in mixed) + data["shape_all_trials_admissible"] = all(r["all_trials_admissible"] == "1" for r in mixed) + data["shape_all_increments_resolved"] = all(r["all_state_increments_resolved"] == "1" for r in mixed) + data["shape_all_alpha_squared_valid"] = all(r["alpha_squared_valid"] == "1" for r in mixed) + if "--pressure-probe" in args: + data["pressure_summary"] = metrics(folder / "pressure_obstruction_summary.csv") + pressure = rows(folder / "pressure_obstruction_trials.csv") + blocks = rows(folder / "pressure_obstruction_blocks.csv") + data["pressure_admissible_trials"] = sum(r["admissible"] == "1" for r in pressure) + data["pressure_expected_block_rows"] = nr * data["pressure_admissible_trials"] + data["pressure_block_rows_complete"] = len(blocks) == data["pressure_expected_block_rows"] + data["pressure_trial_keys_unique"] = len({(r["direction"], r["alpha"]) for r in pressure}) == len(pressure) + data["pressure_block_keys_complete"] = {(r["direction"], r["alpha"], r["block"]) for r in blocks} == {(r["direction"], r["alpha"], b) for r in pressure if r["admissible"] == "1" for b in residual_blocks} + data["pressure_limiting_h_prediction_max_absolute_error"] = { + family: max(abs(float(r[f"{family}_limiting_h_prediction"]) - float(r[f"{family}_limiting_h_evaluated"])) for r in pressure) + for family in ("closure", "pressure") + } + result["completed"].append(data) + +(ROOT / "focused_audit_checks.json").write_text(json.dumps(clean(result), indent=2) + "\n") +print(json.dumps({"completed": [r["id"] for r in result["completed"]], "pending": result["pending"]}, indent=2)) diff --git a/experiments/plateau_analysis/audit_surface_outputs.py b/experiments/plateau_analysis/audit_surface_outputs.py new file mode 100644 index 0000000..73cf3cd --- /dev/null +++ b/experiments/plateau_analysis/audit_surface_outputs.py @@ -0,0 +1,167 @@ +"""Audit completed supplemental surface observations without running solvers. + +Pass the campaign directory containing captures/, focused/, supplemental/ and +supplemental_inputs/. Writes surface_audit_checks.json in that directory. +This records method invariants and sampled values, not physical acceptance. +""" +import argparse +import csv +import hashlib +import json +import math +from pathlib import Path + +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument("campaign_root", type=Path) +ROOT = parser.parse_args().campaign_root.resolve() +SUPPLEMENTAL = ROOT / "supplemental" + + +def sha(path): + return hashlib.sha256(path.read_bytes()).hexdigest() + + +def rows(path): + with path.open() as source: + data = list(csv.DictReader(source)) + if any(None in r or any(v is None for v in r.values()) for r in data): + raise ValueError(f"Malformed CSV row: {path}") + return data + + +def metrics(path): + data = rows(path) + if len({r["metric"] for r in data}) != len(data): + raise ValueError(f"Duplicate metric: {path}") + return {r["metric"]: float(r["value"]) for r in data} + + +def vector(path): + tokens = path.read_text().split() + data = list(map(float, tokens[1:])) + if len(data) != int(tokens[0]) or not all(map(math.isfinite, data)): + raise ValueError(f"Malformed/nonfinite vector: {path}") + return data + + +def clean(value): + if isinstance(value, dict): + return {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 + + +provenance = json.loads((SUPPLEMENTAL / "provenance.json").read_text()) +derivation = json.loads((ROOT / "supplemental_inputs/between_bounds_derivation.json").read_text()) +base = vector(Path(derivation["normalized_base"])) +direction = vector(Path(derivation["direction"])) +candidate = vector(Path(derivation["candidate"])) +if not len(base) == len(direction) == len(candidate): + raise ValueError("Midpoint candidate vector dimensions differ") +alpha = derivation["alpha"] +bounds = metrics(ROOT / "focused/n1_J0p10_iteration_2/frozen_enthalpy_bound.csv") +result = { + "source_commit": provenance["commit"], + "binary_hash_matches": sha(SUPPLEMENTAL / "plateau_diagnostic") == provenance["binary_sha256"], + "midpoint": { + "alpha": alpha, + "vector_size": len(candidate), + "maximum_reconstruction_error": max(abs(b + alpha * d - c) for b, d, c in zip(base, direction, candidate)), + "candidate_hash_matches": sha(Path(derivation["candidate"])) == derivation["candidate_sha256"], + "bounds_match_original_scan": bounds["closure_alpha_h"] == derivation["closure_bound"] and bounds["pressure_alpha_h"] == derivation["pressure_bound"], + "strictly_between_bounds": derivation["closure_bound"] < alpha < derivation["pressure_bound"], + "exact_arithmetic_midpoint": alpha == 0.5 * (derivation["closure_bound"] + derivation["pressure_bound"]), + }, + "completed": [], + "pending": [], +} +for case in json.loads((SUPPLEMENTAL / "cases.json").read_text()): + case_id = case["id"] + manifest = SUPPLEMENTAL / f"{case_id}_run.json" + if not manifest.exists(): + result["pending"].append(case_id) + continue + run = json.loads(manifest.read_text()) + if "exit_code" not in run: + result["pending"].append(case_id) + continue + folder = SUPPLEMENTAL / case_id + args = case["arguments"] + physical_path = Path(args[args.index("--probe-state") + 1]) + direction_path = Path(args[args.index("--probe-direction") + 1]) + capture_index = int(physical_path.parent.name.split("_")[-1]) + original_focused = ROOT / "focused" / f"n1_J0p10_iteration_{capture_index}" + capture_root = physical_path.parent.parent + capture_manifest = json.loads((capture_root.parent / f"{capture_root.name}_run.json").read_text()) + summary = metrics(folder / "surface_trace_summary.csv") + carrier = rows(folder / "surface_trace_rows.csv") + limits = rows(folder / "surface_trace_limits.csv") + nodes = rows(folder / "surface_trace_nodes.csv") + physical_base = vector(physical_path) + normalized_direction = vector(direction_path) + capture_metadata = dict(line.split("=", 1) for line in (physical_path.parent / "capture.txt").read_text().splitlines()) + replay_alpha = alpha if "between_bounds" in case_id else float(capture_metadata["first_trial_alpha"]) + data = { + "id": case_id, + "exit_code": run["exit_code"], + "seconds": run["seconds"], + "input_hash_count": len(run["input_state_hashes"]), + "input_hash_mismatches": [p for p, h in run["input_state_hashes"].items() if sha(Path(p)) != h], + "output_hash_count": len(run["output_hashes"]), + "output_hash_mismatches": [p for p, h in run["output_hashes"].items() if sha(folder / p) != h], + "unhashed_data_files": sorted(str(p.relative_to(folder)) for p in folder.rglob("*") if p.is_file() and p.suffix in (".txt", ".csv", ".json") and str(p.relative_to(folder)) not in run["output_hashes"]), + "auxiliary_capture_input_hashes_match": all(sha(physical_path.parent / name) == capture_manifest["output_hashes"][str((physical_path.parent / name).relative_to(capture_root))] for name in ("normalized_base.txt", "normalized_residual.txt")), + "frozen": metrics(folder / "frozen_input_checks.csv"), + "frozen_residual_byte_identical_to_original_focus": (folder / "frozen_normalized_residual.txt").read_bytes() == (original_focused / "frozen_normalized_residual.txt").read_bytes(), + "enthalpy_bound_byte_identical_to_original_focus": (folder / "frozen_enthalpy_bound.csv").read_bytes() == (original_focused / "frozen_enthalpy_bound.csv").read_bytes(), + "direction_copy_byte_identical": direction_path.read_bytes() == (folder / "probe_normalized_correction.txt").read_bytes(), + "physical_base_copy_byte_identical": physical_path.read_bytes() == (folder / "frozen_physical_base.txt").read_bytes(), + "candidate_replay": rows(folder / "replayed_candidate.csv")[0], + "surface_summary": summary, + "csv_rows": {p.name: len(rows(p)) for p in sorted(folder.glob("*.csv"))}, + "carrier_rows_complete": len(carrier) == int(summary["surface_row_count"]) and {int(r["surface_index"]) for r in carrier} == set(range(len(carrier))), + "carrier_reduced_indices_unique": len({r["reduced_h_index"] for r in carrier}) == len(carrier), + "carrier_root_offsets_correct": all(int(r["root_index"]) == int(summary["h_root_offset"]) + int(r["reduced_h_index"]) for r in carrier), + "carrier_base_matches_physical_input": all(float(r["base_h"]) == physical_base[int(r["root_index"])] for r in carrier), + "all_input_carrier_normalized_direction_coefficients_zero": all(normalized_direction[int(r["root_index"])] == 0 for r in carrier), + "carrier_unit_step_arithmetic_error": max(abs(float(r["base_h"]) + float(r["direction_dh"]) - float(r["h_plus_dh"])) for r in carrier), + "carrier_recomputed": {}, + "carrier_summary_differences": {}, + "limit_families_complete": {r["family"] for r in limits} == {"closure", "pressure"} and len(limits) == 2, + "limit_partition_max_errors": {}, + "nodes_by_family": {}, + } + for name in ("base_h", "direction_dh", "h_plus_dh"): + values = [float(r[name]) for r in carrier] + data["carrier_recomputed"][name] = { + "minimum": min(values), "maximum": max(values), "maximum_absolute": max(map(abs, values)), + "norm": math.sqrt(math.fsum(v * v for v in values)), + "negative_count": sum(v < 0 for v in values), "nonzero_count": sum(v != 0 for v in values), + } + data["carrier_summary_differences"][name] = {stat: value - summary[f"{name}_{stat}"] for stat, value in data["carrier_recomputed"][name].items()} + for suffix in ("h", "dh", "h_plus_dh"): + data["limit_partition_max_errors"][suffix] = max(abs(float(r[f"full_{suffix}"]) - float(r[f"boundary_{suffix}"]) - float(r[f"interior_{suffix}"])) for r in limits) + for limit in limits: + family = limit["family"] + subset = [r for r in nodes if r["family"] == family] + data["nodes_by_family"][family] = { + "rows": len(subset), + "element_matches_limit": all(r["element"] == limit["element"] for r in subset), + "node_indices_complete": len(subset) > 0 and len(subset) == int(subset[0]["node_count"]) and {int(r["local_node"]) for r in subset} == set(range(len(subset))), + "fe_orders": sorted({int(r["fe_order"]) for r in subset}), + "geometry_types": sorted({int(r["geometry_type"]) for r in subset}), + "base_minimum": min(float(r["h"]) for r in subset), + "union_prediction_minimum": min(float(r["h_plus_alpha_union_dh"]) for r in subset), + "unit_step_minimum": min(float(r["h_plus_dh"]) for r in subset), + "replay_alpha": replay_alpha, + "replay_candidate_nodal_prediction_minimum": min(float(r["h"]) + replay_alpha * float(r["dh"]) for r in subset), + "union_prediction_max_arithmetic_error": max(abs(float(r["h"]) + float(r["alpha_union"]) * float(r["dh"]) - float(r["h_plus_alpha_union_dh"])) for r in subset), + "unit_step_max_arithmetic_error": max(abs(float(r["h"]) + float(r["dh"]) - float(r["h_plus_dh"])) for r in subset), + } + result["completed"].append(data) + +(ROOT / "surface_audit_checks.json").write_text(json.dumps(clean(result), indent=2) + "\n") +print(json.dumps({"completed": [r["id"] for r in result["completed"]], "pending": result["pending"]}, indent=2)) diff --git a/experiments/plateau_analysis/sample_pressure_polynomial.py b/experiments/plateau_analysis/sample_pressure_polynomial.py new file mode 100644 index 0000000..3bb21ec --- /dev/null +++ b/experiments/plateau_analysis/sample_pressure_polynomial.py @@ -0,0 +1,218 @@ +#!/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() diff --git a/experiments/plateau_diagnostic_README.md b/experiments/plateau_diagnostic_README.md index cbcf446..c3e4bf3 100644 --- a/experiments/plateau_diagnostic_README.md +++ b/experiments/plateau_diagnostic_README.md @@ -193,3 +193,16 @@ between the closure and pressure enthalpy bounds and compares surface traces at onset and late plateau. The candidate is generated from captured normalized base and direction at the arithmetic midpoint of those two measured bounds; its derivation and hash are retained with the supplemental input. + +The surface trace output also records actual FE nodal h and dh values on each +limiting element, obtained through full physical-field scatters and evaluations +at `GetNodes()`. Shared elements are emitted once per quadrature family. This +allows sign comparisons at arbitrary alpha without confusing coefficient signs +with interior polynomial positivity. The union-bound column is an affine +prediction at the boundary, not a separately prepared candidate. + +The completed campaign and evidence limits are recorded in +`PRESSURE_AND_SHAPE_CAMPAIGN_2026-09-17.md`. Portable offline analyses are in +`plateau_analysis/`; pass the retained artifact root explicitly. The Q3 follow-up +checks a finite interior grid against the exported polynomial and does not +certify positivity everywhere.