Files
MeanField/experiments/plateau_shape_curvature_README.md

8.4 KiB

Shape and field Taylor decomposition

plateau_shape_curvature_diagnostics.hpp is an experiment-only observer, included after import mean_field;. It also includes the existing sampled enthalpy-bound helper. The root driver supplies an actual frozen normalized base and selected normalized correction:

experiment::plateau_diagnostics::ProbeShapeCurvature(
    state, finiteElements, polytropicIndex, normalizedBase, normalizedDirection,
    outputDirectory);

Use a distinct output directory for each state/direction. The inputs are copied before state workspaces change, so aliases of those workspaces are safe. The caller must use this context's normalization and matching EOS index; the helper does not reconstruct a Newton correction from a damped trial. If the driver does that reconstruction, its alpha and roundoff uncertainty belong in the enclosing provenance, particularly when the accepted step is tiny.

Experiment and preparation budget

The actual value manifest must contain exactly one nonempty block whose stable ID starts with surface. Only that block enters p_shape; all remaining value coordinates, including scalar multipliers and angular velocity, enter p_fields. No value/residual offsets or field sizes are hard-coded.

The helper prepares the exact supplied base immediately before each of its three Jacobian actions. It records repeated base residuals and checks p=p_shape+p_fields and Jp=Jp_shape+Jp_fields. These measurements are made before any step-size gate, even when the direction's enthalpy boundary makes curvature differencing unusable.

At that same prepared base, it measures geometry limits for all three directions using the production preflight rules, and scans the exact production closure and pressure enthalpy quadrature for full and fields directions separately. Their enthalpy bounds should agree because the shape direction changes no enthalpy coefficient. The shape-only bound is therefore infinite for an admissible base. The common maximum alpha is the minimum of one, all geometry-safe bounds, and 0.8 times the finite full/fields enthalpy crossings. The geometric estimator uses its existing safety fraction 0.9.

Six common scales are used: alpha_max times 1, 0.1, 0.01, 0.001, 0.0001, 0.00001. Each scale prepares full, shape-only, and fields-only trials. Every trial vector is constructed directly from the copied base; trials never accumulate. The ordinary budget is 22 preparations (three for J actions, 18 trials, final base restoration), three Jv, and zero preconditioner applications or linear solves. An unusable common cap skips all trial preparations, retaining the action/bound diagnostics. Exceptions attempt one additional base restoration before rethrow.

For each residual block and direction it records the vector Taylor remainder

R_d(alpha) = F(x + alpha p_d) - F(x) - alpha Jp_d.

The mixed term is formed by vector subtraction before taking its norm:

R_mixed = R_full - R_shape - R_fields.

For each block the mixed file also records the linear prediction L=F_base+alpha Jp_full, its norm, the measured full-trial residual norm, and long-double-accumulated signed products <L,R_full>, <L,R_shape>, <L,R_fields>, and <L,R_mixed>. These distinguish curvature amplitude from its effect on the merit. The file checks both scalar identities

||F_full_trial||² = ||L||² + 2<L,R_full> + ||R_full||²
<L,R_full> = <L,R_shape> + <L,R_fields> + <L,R_mixed>.

The errors of these identities are exported, not forced to zero. Individual signed projections remain available when their own trial is valid, even if another direction at that alpha is inadmissible; mixed projections require all three. Summing signed products across residual blocks gives the full normalized merit decomposition. A large norm alone does not establish an adverse merit contribution.

It is not the difference of three norms. Six remainder cosines and the sum of component norms divided by the full norm show cancellation between terms. All three trials must be admissible before a mixed result exists; otherwise mixed values are NaN. Individual valid trials remain available even if another trial at that alpha fails.

At fixed shape, the two gravity residual blocks should be affine in the remaining independent fields. The fields rows therefore test an arithmetic-level gravity remainder. The helper records that measurement without forcing it to zero. Closure, pressure/mechanics, rotation and scalar constraints need not be affine in these remaining coordinates; the same expectation does not apply to every residual block.

Arithmetic limits

The state addition error is measured explicitly as (fl(x+alpha p)-x)-alpha p. A nonzero direction's increment is flagged resolved only if its expected norm exceeds 64 epsilon (1+||x||), its actual increment is nonzero, and its addition error is below 1% of the expected increment. These are diagnostic numerical heuristics, not physical tolerances. Tiny components may be unresolved even when the whole direction is resolved.

Per-block residual roundoff scales use the larger of eight times the observed pre-action base-restoration drift and 128 epsilon (1+||F_base,block||+alpha||Jp_block||). Mixed scales sum the three component scales. The above_roundoff_scale flag also requires admissibility and a resolved state increment. This heuristic is not a certified bound on matrix-free assembly error: terms may cancel before the residual is formed. Use the alpha trend, the fields-only gravity measurement, and the final restoration CSV together. If final restoration drift exceeds the pre-action drift, increase the interpretation floor accordingly; previously written rows are not retroactively rescaled.

If the common alpha cap is zero/nonfinite, its largest full increment is below the update threshold, or the Jacobian split exceeds its recorded roundoff heuristic, common_step_regime_usable=0 and no curvature is inferred. Below-floor remainders are unresolved, not evidence of zero physical curvature. Inadmissible trials and zero-denominator cosines/relative errors stay NaN, never fabricated zeros. Root norms remain in the supplied normalization; manifest scale metadata is not used as a replacement normalization.

alpha_squared_valid requires a positive finite floating-point alpha². If it underflows or is nonfinite, scaled curvature is NaN and the resolved-curvature flag is false. Invalid denominators or overflow in scaled ratios also produce NaN; a nonpositive/nonfinite alpha itself is skipped and recorded in the summary. Absolute measured remainders may still be available when squared scaling is not.

Outputs and lifecycle

  • shape_curvature_summary.csv: base/direction sizes and norms, split checks, arithmetic thresholds, common cap, preparation count, and accepted-vector invariance checks.
  • shape_curvature_preparations.csv: every attempted preparation with exact available rejection reason/stage, names, meaningful EOS code, and the reported minimum mapping determinant. The production payload does not identify an offending element or quadrature point; the helper does not invent one.
  • shape_curvature_restorations.csv: per-block residual differences on each repeated base preparation, including final restoration.
  • shape_curvature_directions.csv and shape_curvature_actions.csv: norms split with the actual value/residual manifests and Jacobian reconstruction.
  • shape_curvature_bounds.csv: geometry and sampled pressure/closure limits, including the enthalpy scanner's limiting sample coordinates.
  • shape_curvature_remainders.csv: absolute R, R/alpha², derivative error R/alpha, relative derivative error, arithmetic flags, and actual state-addition errors for every block/direction/alpha.
  • shape_curvature_mixed.csv: the mixed vector's norm/scaling, cosines, linear prediction and trial norms, signed merit projections, and identity errors.

No trial is committed and no accepted vector is changed. The helper leaves the supplied frozen base prepared and restores the entry normalized correction, physical correction, and volume-displacement-direction workspaces. Preparation revisions advance, so a caller that subsequently applies an inverse must refresh it; the enclosing WithState hook normally restores its accepted state and inverse. The helper performs no preconditioner refresh or policy change itself. Its rejection/event files are flushed incrementally. No physics, normalization, production API, driver, or build configuration is changed by this helper.