6.9 KiB
Polynomial enthalpy admissibility and feasible direction repair
This branch adds a sufficient polynomial nonnegativity certificate and an interchangeable enthalpy direction repair. It does not alter the EOS or declare a constrained residual minimum an equilibrium.
TensorBernsteinPolynomial<Dimension> encodes spatial dimension in its type and accepts polynomial order through the mesh's runtime nodal points. A cached one-dimensional Lagrange-to-Bernstein transformation is applied by tensor contractions. Every nontrivial arithmetic operation in the transformation and de Casteljau subdivision encloses its result using nextafter; a polynomial is certified only when all coefficient lower bounds in every leaf are nonnegative. The implementation assumes ordinary IEEE binary64 arithmetic, gradual underflow, and no unsafe fast-math transformations. The underlying mathematical FE is the Lagrange polynomial at the supplied floating-point nodal coordinates.
A negative control coefficient is not a negative polynomial witness. Such cells are refined. A negative upper interval bound at a tensor corner proves a negative polynomial value. Cells that remain ambiguous at the work budget return unresolved, including some truly nonnegative polynomials with tangential interior zeros. No negative tolerance is accepted as nonnegative. Exact endpoint interpolation preserves identically zero tensor faces without manufacturing a small negative error interval.
The production adapter supports conforming scalar nodal H1 hexahedra with closed Gauss-Lobatto, closed uniform, or closed Gauss-Legendre points. Unsupported element types are rejected explicitly. Its reduced field map and MFEM prolongation preserve shared true DOFs; no elementwise coefficient clipping occurs. The prepared object assumes the finite-element space and topology remain unchanged. Mesh adaptation requires a new object.
Given a certified base h, a certified anchor a, and proposed direction p, the repaired full candidate is (1-beta)*(h+p)+beta*a. Thus the altered enthalpy direction is (1-beta)*p+beta*(a-h). Zero-trace reduced rows are fixed to zero. The adapter constructs a sufficient convex cone from a subdivision on which the anchor has nonnegative Bernstein bounds, and computes the largest required fraction across all coefficients and MPI ranks. This is the minimum sufficient fraction for those conservative coefficient inequalities, not the least fraction for the complete nonnegative polynomial set. A caller can choose any larger fraction through one, then verify the actual rounded full candidate. The original direction is left unchanged when verification fails.
This changes the enthalpy direction while allowing the stellar strategy to retain the other field directions. It is a feasible direction repair, not a constrained Newton optimizer. Full residual acceptance, actual descent checks, and the original physical convergence requirements remain necessary. Invalid or unresolved bases are rejected explicitly; this mechanism is not a feasibility-restoration solver. A certified checkpoint or initial seed is required as anchor. A frozen repaired full candidate and certified base define a nonnegative exact-arithmetic segment for backtracking; every floating-point candidate must nevertheless pass the same check.
Collective adapter calls must occur in the same order on all ranks. Sign status, polynomial bounds, and repair fraction are globally combined. Witness element indices, coordinates, and work counters remain rank-local. Unequal repair fractions and work budgets are rejected collectively.
Tests cover nodal-positive/interior-negative cubics, a positive polynomial with negative coarse Bernstein coefficients, exact zero faces, very small negative pockets, unresolved tangential zeros, arbitrary runtime degree, nonfinite input, finite-element ordering, true-DOF sharing, invalid bases, immutable failure output, and prescribed trace handling. The compiled production module passed serial tests and supported low-level two-rank tests. Full stellar MPI operation remains unvalidated because the existing logical radial surface interpolation is restricted to one rank.
PreparedStellarEnthalpyStrategy<State> resolves the enthalpy block from the compiled form type and captures the accepted physical enthalpy at the start of each evaluate() as its anchor. The anchor remains fixed during that evaluation. Every full-state check denormalizes the state; pointwise positivity is not invariant under a general diagonal Riesz normalization. The strategy minimizes the full normalized linear-residual norm over the sufficient one-dimensional repair interval, verifies the actual normalized root candidate after rounding, and rejects a non-descent result. Unrelated blocks of the normalized correction are retained bit for bit. It does not prepare an operator or change an accepted state. A valid unaltered full candidate costs one Jacobian action for its true merit slope; an altered direction costs three actions (family model plus the actual rounded chosen direction). The anchor's certificate is cached only when coefficients agree exactly on every rank.
A standalone extraction of the scalar/tensor certification engine was compiled with the system C++20 compiler and exercised before module integration. It certified a strictly positive quadratic with negative coarse coefficients and a three-dimensional zero-face bubble; it proved a negative pocket of depth 1e-10; it returned unresolved for an exact tangential zero. The subsequent integrated frozen probes certified the earlier rotating base, detected a negative full trial, and rejected the already-negative late rotating base with invalid_base while preserving accepted buffers. The n=1, J=0 polynomial-policy control converged and passed all eight virial checks; feasible repaired rotating candidates in the frozen probe did not pass full residual acceptance.
Arithmetic requirements are enforced at compilation: __FAST_MATH__ and a positive __FINITE_MATH_ONLY__ are rejected, and static assertions require IEC 559 binary64 with subnormal support. Runtime floating-point modes must retain gradual underflow; the application must not enable flush-to-zero after construction. The existing Release build was inspected: -O3 -g -DNDEBUG -std=c++23 -flto=thin -arch arm64; no unsafe floating-point flags were present in the generated build commands.
The stellar strategy explicitly supports only state-independent pressure-to-specific-enthalpy surface relations. That dependency is checked from the compiled relation type. Its constructor collectively rejects nonzero prescribed surface pressure with an unsupported-prescription error before checking seed coefficients, and verifies that the EOS's enthalpy reference satisfies h(P=0)=0. This prevents a valid positive-pressure model from being mislabeled an invalid numerical base. The focused regression validates zero pressure and collective rejection when only one MPI rank supplies positive pressure.