From 9e0e1bc42af6cf11986db113dfbeed694fdf4132 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 09:28:01 -0400 Subject: [PATCH] Integrate polynomial feasibility with typed stellar enthalpy strategy --- ...OMIAL_ENTHALPY_ADMISSIBILITY_2026-09-17.md | 4 + .../solver/enthalpy_admissibility.cppm | 21 +- .../solver/stellar_enthalpy_strategy.cppm | 192 ++++++++++++++++++ 3 files changed, 215 insertions(+), 2 deletions(-) create mode 100644 libmeanfield/interface/solver/stellar_enthalpy_strategy.cppm diff --git a/experiments/POLYNOMIAL_ENTHALPY_ADMISSIBILITY_2026-09-17.md b/experiments/POLYNOMIAL_ENTHALPY_ADMISSIBILITY_2026-09-17.md index 352cd46..ecda41e 100644 --- a/experiments/POLYNOMIAL_ENTHALPY_ADMISSIBILITY_2026-09-17.md +++ b/experiments/POLYNOMIAL_ENTHALPY_ADMISSIBILITY_2026-09-17.md @@ -15,3 +15,7 @@ This changes the enthalpy direction while allowing the stellar strategy to retai 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. Integration tests are intended for both one and two MPI ranks. + +`PreparedStellarEnthalpyStrategy` resolves the enthalpy block from the compiled form type and captures the initial physical enthalpy as its anchor. 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. Production module compilation and MPI execution remain part of the integration campaign. diff --git a/libmeanfield/interface/solver/enthalpy_admissibility.cppm b/libmeanfield/interface/solver/enthalpy_admissibility.cppm index fc5ba41..1bdddef 100644 --- a/libmeanfield/interface/solver/enthalpy_admissibility.cppm +++ b/libmeanfield/interface/solver/enthalpy_admissibility.cppm @@ -37,12 +37,16 @@ export namespace mean_field::solver { }; enum class EnthalpyDirectionRepairStatus { - feasible, invalid_base, invalid_anchor, invalid_direction, unresolved + feasible, invalid_base, invalid_anchor, invalid_direction, unresolved, no_descent }; struct EnthalpyDirectionRepairReport final { EnthalpyDirectionRepairStatus status{EnthalpyDirectionRepairStatus::unresolved}; double minimumSafeFraction{1.0}; + double chosenFraction{0.0}; + double meritDirectionalDerivative{std::numeric_limits::quiet_NaN()}; + double linearModelResidualNorm{std::numeric_limits::quiet_NaN()}; + int jacobianApplications{0}; PolynomialSignReport base; PolynomialSignReport anchor; PolynomialSignReport candidate; @@ -415,7 +419,17 @@ export namespace mean_field::solver { if (report.base.status != PolynomialSignStatus::certified_nonnegative || !ZeroTrace(base)) { report.status = EnthalpyDirectionRepairStatus::invalid_base; return report; } - report.anchor = Check(anchor); + // The anchor can be reused only for identical coefficients. All + // ranks agree on the cache path before any later scatter/reduction. + bool cached = m_cachedAnchor.Size() == anchor.Size() && m_hasCachedAnchor; + if (cached) + for (int i = 0; i < anchor.Size(); ++i) cached = cached && m_cachedAnchor(i) == anchor(i); + if (enthalpy_detail::CollectiveMaximum(cached ? 0 : 1,m_space->GetComm())) { + m_cachedAnchorReport = Check(anchor); + m_cachedAnchor = anchor; + m_hasCachedAnchor = true; + } + report.anchor = m_cachedAnchorReport; if (report.anchor.status != PolynomialSignStatus::certified_nonnegative || !ZeroTrace(anchor)) { report.status = EnthalpyDirectionRepairStatus::invalid_anchor; return report; } @@ -504,5 +518,8 @@ export namespace mean_field::solver { EnthalpyAdmissibilityOptions m_options; std::map,std::unique_ptr>> m_bases; std::vector m_elements; + mutable mfem::Vector m_cachedAnchor; + mutable PolynomialSignReport m_cachedAnchorReport; + mutable bool m_hasCachedAnchor{false}; }; } diff --git a/libmeanfield/interface/solver/stellar_enthalpy_strategy.cppm b/libmeanfield/interface/solver/stellar_enthalpy_strategy.cppm new file mode 100644 index 0000000..034fc5d --- /dev/null +++ b/libmeanfield/interface/solver/stellar_enthalpy_strategy.cppm @@ -0,0 +1,192 @@ +module; + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +export module mean_field:solver.stellar_enthalpy_strategy; +export import :solver.enthalpy_admissibility; +export import :normalization.stellar_equilibrium; +export import :utils.domain; + +export namespace mean_field::solver { + /* + * The solver owns this workspace only when a positivity policy is selected. + * Certification uses physical h: arbitrary positive diagonal normalization + * need not preserve the polynomial represented by its normalized entries. + */ + template + class PreparedStellarEnthalpyStrategy final { + using Problem = std::remove_cvref_t().Problem())>; + using Form = typename Problem::FormType; + using Schema = utils::domain::CoreEnvelopeVacuumDomainSchema; + using Domain = field::FieldDomainT; + static constexpr int enthalpyIndex = utils::blocks::type_index_v< + utils::blocks::enthalpy::specific::value,typename Form::value_blocks>; + static_assert(enthalpyIndex >= 0, "Polynomial positivity requires the compiled specific-enthalpy block."); + static constexpr utils::blocks::value_block enthalpyBlock{}; + public: + explicit PreparedStellarEnthalpyStrategy(State &state, EnthalpyAdmissibilityOptions options = {}) + : m_state(&state), + m_offset(state.Problem().GetManifest().layout().offset(enthalpyBlock)), + m_size(state.Problem().GetManifest().layout().size(enthalpyBlock)) { + const auto &fem = equilibrium::detail::StellarEquilibriumProblemFactory::FiniteElementModel(state.Problem()); + if (!fem.mesh || !fem.enthalpyFes) + throw std::invalid_argument("Stellar enthalpy certification requires a complete enthalpy FE space."); + std::vector elements; + for (int e = 0; e < fem.mesh->GetNE(); ++e) + if (Schema::template attribute_belongs_to(fem.mesh->GetAttribute(e))) elements.push_back(e); + const auto &surface = state.Problem().GetPressureSurfaceRows().reduced_dofs(); + m_prepared = std::make_unique( + *fem.enthalpyFes,field::make_field_dof_map(*fem.enthalpyFes),elements, + std::span(surface.GetData(),surface.Size()),options); + if (m_prepared->Size() != m_size) + throw std::invalid_argument("Compiled enthalpy block disagrees with the supported field map."); + ExtractEnthalpy(state.AcceptedPhysicalState(),m_anchor); + m_anchorReport = m_prepared->Check(m_anchor); + } + + [[nodiscard]] const PolynomialSignReport &AnchorReport() const noexcept { return m_anchorReport; } + + [[nodiscard]] PolynomialSignReport CheckState(const mfem::Vector &fullNormalized) const { + m_state->normalizedOperator->DenormalizeState(fullNormalized,m_physicalWork); + ExtractEnthalpy(m_physicalWork,m_enthalpyWork); + return m_prepared->Check(m_enthalpyWork); + } + + [[nodiscard]] EnthalpyDirectionRepairReport RepairDirection(State &state) { + if (&state != m_state) throw std::invalid_argument("A prepared enthalpy strategy belongs to one context."); + ExtractEnthalpy(state.AcceptedPhysicalState(),m_base); + state.normalizedOperator->DenormalizeState(state.normalizedCorrection,m_physicalDirection); + ExtractEnthalpy(m_physicalDirection,m_originalDirection); + auto report = m_prepared->MinimumRepairFraction(m_base,m_anchor,m_originalDirection); + if (report.status != EnthalpyDirectionRepairStatus::feasible) return report; + + const auto &surface = state.Problem().GetPressureSurfaceRows().reduced_dofs(); + bool traceUnchanged = true; + for (int row : surface) traceUnchanged = traceUnchanged && m_originalDirection(row) == 0.0; + int localChange = traceUnchanged ? 0 : 1, globalChange = 0; + if (MPI_Allreduce(&localChange,&globalChange,1,MPI_INT,MPI_MAX,state.Problem().GetCommunicator()) != MPI_SUCCESS) + throw std::runtime_error("Could not resolve collective enthalpy trace preservation."); + if (report.minimumSafeFraction == 0.0 && globalChange == 0) { + m_candidateState = state.acceptedNormalizedState; m_candidateState += state.normalizedCorrection; + report.candidate = CheckState(m_candidateState); + if (report.candidate.status != PolynomialSignStatus::certified_nonnegative) { + report.status = EnthalpyDirectionRepairStatus::unresolved; + return report; + } + m_chosenAction.SetSize(state.Problem().EquationSize()); + state.normalizedOperator->Mult(state.normalizedCorrection,m_chosenAction); + report.jacobianApplications = 1; + m_linearResidual = state.acceptedNormalizedResidual; m_linearResidual += m_chosenAction; + report.meritDirectionalDerivative = Dot(state.acceptedNormalizedResidual,m_chosenAction); + report.linearModelResidualNorm = std::sqrt(Dot(m_linearResidual,m_linearResidual)); + if (!std::isfinite(report.meritDirectionalDerivative) || report.meritDirectionalDerivative >= 0.0) + report.status = EnthalpyDirectionRepairStatus::no_descent; + return report; + } + m_zeroFractionDirection = state.normalizedCorrection; + m_physicalWork = m_physicalDirection; + for (int row : surface) m_physicalWork(m_offset+row) = 0.0; + state.normalizedOperator->NormalizeState(m_physicalWork,m_zeroFractionDirection); + // Preserve unrelated blocks bit for bit. A normalize/denormalize + // round trip is not an excuse to perturb their accepted direction. + CopyOtherBlocks(state.normalizedCorrection,m_zeroFractionDirection); + m_unitFractionDirection = state.normalizedCorrection; + m_physicalWork = m_physicalDirection; + for (int i = 0; i < m_size; ++i) m_physicalWork(m_offset+i) = m_anchor(i)-m_base(i); + for (int row : surface) m_physicalWork(m_offset+row) = 0.0; + state.normalizedOperator->NormalizeState(m_physicalWork,m_unitFractionDirection); + CopyOtherBlocks(state.normalizedCorrection,m_unitFractionDirection); + m_familyDirection = m_unitFractionDirection; m_familyDirection -= m_zeroFractionDirection; + m_firstAction.SetSize(state.Problem().EquationSize()); + m_familyAction.SetSize(state.Problem().EquationSize()); + state.normalizedOperator->Mult(m_zeroFractionDirection,m_firstAction); + state.normalizedOperator->Mult(m_familyDirection,m_familyAction); + report.jacobianApplications = 2; + m_linearResidual = state.acceptedNormalizedResidual; m_linearResidual += m_firstAction; + const double denominator = Dot(m_familyAction,m_familyAction); + const double numerator = Dot(m_familyAction,m_linearResidual); + double fraction = report.minimumSafeFraction; + if (std::isfinite(denominator) && denominator > 0.0 && std::isfinite(numerator)) + fraction = std::clamp(-numerator/denominator,report.minimumSafeFraction,1.0); + + auto candidate = m_prepared->ApplyRepairFraction(m_base,m_anchor,m_originalDirection,fraction,m_repaired); + if (candidate.status != PolynomialSignStatus::certified_nonnegative && fraction != 1.0) { + fraction = 1.0; + candidate = m_prepared->ApplyRepairFraction(m_base,m_anchor,m_originalDirection,fraction,m_repaired); + } + if (candidate.status != PolynomialSignStatus::certified_nonnegative) { + report.status = EnthalpyDirectionRepairStatus::unresolved; + report.candidate = candidate; + return report; + } + m_physicalWork = m_physicalDirection; + for (int i = 0; i < m_size; ++i) m_physicalWork(m_offset+i) = m_repaired(i); + state.normalizedOperator->NormalizeState(m_physicalWork,m_chosenDirection); + CopyOtherBlocks(state.normalizedCorrection,m_chosenDirection); + + // State normalization and root addition introduce another rounding + // boundary. Certify exactly what backtracking will construct. + m_candidateState = state.acceptedNormalizedState; m_candidateState += m_chosenDirection; + candidate = CheckState(m_candidateState); + if (candidate.status != PolynomialSignStatus::certified_nonnegative) { + report.status = EnthalpyDirectionRepairStatus::unresolved; + report.candidate = candidate; + return report; + } + m_chosenAction.SetSize(state.Problem().EquationSize()); + state.normalizedOperator->Mult(m_chosenDirection,m_chosenAction); + ++report.jacobianApplications; + m_linearResidual = state.acceptedNormalizedResidual; m_linearResidual += m_chosenAction; + report.meritDirectionalDerivative = Dot(state.acceptedNormalizedResidual,m_chosenAction); + report.linearModelResidualNorm = std::sqrt(Dot(m_linearResidual,m_linearResidual)); + report.chosenFraction = fraction; + report.candidate = candidate; + // A feasible family can lack descent. Leave the original direction + // available for diagnostics rather than silently accepting that case. + if (!std::isfinite(report.meritDirectionalDerivative) || report.meritDirectionalDerivative >= 0.0) { + report.status = EnthalpyDirectionRepairStatus::no_descent; + return report; + } + state.normalizedCorrection = m_chosenDirection; + return report; + } + + private: + void ExtractEnthalpy(const mfem::Vector &physical,mfem::Vector &enthalpy) const { + if (physical.Size() != m_state->Problem().StateSize()) + throw std::invalid_argument("Stellar positivity received an incompatible root state."); + enthalpy.SetSize(m_size); + for (int i = 0; i < m_size; ++i) enthalpy(i) = physical(m_offset+i); + } + void CopyOtherBlocks(const mfem::Vector &source,mfem::Vector &target) const { + for (int i = 0; i < m_offset; ++i) target(i) = source(i); + for (int i = m_offset+m_size; i < target.Size(); ++i) target(i) = source(i); + } + [[nodiscard]] double Dot(const mfem::Vector &left,const mfem::Vector &right) const { + const double local = left*right; + double global = 0.0; + if (MPI_Allreduce(&local,&global,1,MPI_DOUBLE,MPI_SUM,m_state->Problem().GetCommunicator()) != MPI_SUCCESS) + throw std::runtime_error("Could not evaluate the collective feasible-direction model."); + return global; + } + State *m_state; + int m_offset,m_size; + std::unique_ptr m_prepared; + mfem::Vector m_anchor,m_base,m_originalDirection,m_repaired; + PolynomialSignReport m_anchorReport; + mutable mfem::Vector m_physicalWork,m_enthalpyWork; + mfem::Vector m_physicalDirection,m_zeroFractionDirection,m_unitFractionDirection,m_familyDirection; + mfem::Vector m_firstAction,m_familyAction,m_chosenAction,m_linearResidual,m_chosenDirection,m_candidateState; + }; +}