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.csvandshape_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.