From 97204a3253c36edfdbd51cbcb344c1239d772da6 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 09:31:31 -0400 Subject: [PATCH] Add frozen polynomial positivity probes and tighten certificate diagnostics --- ...lateau_positivity_strategy_diagnostics.hpp | 187 ++++++++++++++++++ .../solver/enthalpy_admissibility.cppm | 17 +- tests/solver/enthalpy_admissibility.cpp | 11 ++ 3 files changed, 213 insertions(+), 2 deletions(-) create mode 100644 experiments/plateau_positivity_strategy_diagnostics.hpp diff --git a/experiments/plateau_positivity_strategy_diagnostics.hpp b/experiments/plateau_positivity_strategy_diagnostics.hpp new file mode 100644 index 0000000..b7cdc3d --- /dev/null +++ b/experiments/plateau_positivity_strategy_diagnostics.hpp @@ -0,0 +1,187 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Include after import mean_field. Frozen probes never substitute their base +// into accepted buffers: the context's original seed supplies the safe anchor. +namespace experiment::plateau_diagnostics { +namespace polynomial_probe_detail { +inline const char *Status(mean_field::solver::PolynomialSignStatus status) { + using S = mean_field::solver::PolynomialSignStatus; + switch (status) { + case S::certified_nonnegative: return "certified_nonnegative"; + case S::unresolved: return "unresolved"; + case S::negative: return "negative"; + case S::non_finite: return "non_finite"; + } + return "unknown"; +} +inline const char *RepairStatus(mean_field::solver::EnthalpyDirectionRepairStatus status) { + using S = mean_field::solver::EnthalpyDirectionRepairStatus; + switch (status) { + case S::feasible: return "feasible"; + case S::invalid_base: return "invalid_base"; + case S::invalid_anchor: return "invalid_anchor"; + case S::invalid_direction: return "invalid_direction"; + case S::unresolved: return "unresolved"; + case S::no_descent: return "no_descent"; + } + return "unknown"; +} +inline std::ofstream Open(const std::filesystem::path &path) { + if (std::filesystem::exists(path)) throw std::runtime_error("Refusing to overwrite positivity probe output."); + std::ofstream stream(path); + stream.exceptions(std::ios::badbit | std::ios::failbit); + stream << std::setprecision(17); + return stream; +} +inline void Save(const std::filesystem::path &path,const mfem::Vector &values) { + auto stream = Open(path); stream << values.Size() << '\n'; + for (int i = 0; i < values.Size(); ++i) stream << values(i) << '\n'; +} +inline double Difference(const mfem::Vector &a,const mfem::Vector &b) { + if (a.Size() != b.Size()) return std::numeric_limits::infinity(); + mfem::Vector delta(a); delta -= b; return delta.Norml2(); +} +} + +template +void ProbePolynomialEnthalpy(State &s,const mean_field::fem::FEM &f, + const mfem::Vector &normalizedBase,const mfem::Vector &normalizedDirection, + const std::filesystem::path &output) { + using namespace mean_field; + namespace detail = polynomial_probe_detail; + using Schema = utils::domain::CoreEnvelopeVacuumDomainSchema; + using Domain = field::FieldDomainT; + int ranks = 0; + if (MPI_Comm_size(s.Problem().GetCommunicator(),&ranks) != MPI_SUCCESS || ranks != 1) + throw std::invalid_argument("Frozen polynomial positivity diagnostics require one rank."); + const auto start = std::chrono::steady_clock::now(); + const mfem::Vector acceptedPhysical(s.AcceptedPhysicalState()), acceptedNormalized(s.acceptedNormalizedState), + acceptedResidual(s.acceptedNormalizedResidual), savedCorrection(s.normalizedCorrection), + base(normalizedBase), originalDirection(normalizedDirection); + int preparations = 0; + auto restore = [&]() { + s.trialNormalizedState = base; + ++preparations; + if (!s.PrepareTrial()) throw std::runtime_error("Frozen positivity base failed production preparation."); + }; + auto finish = [&]() { + restore(); s.physicalInverse->Refresh(); s.normalizedCorrection = savedCorrection; + }; + try { + restore(); + const mfem::Vector baseResidual(s.trialNormalizedResidual); + using Problem = std::remove_cvref_t; + using Form = typename Problem::FormType; + constexpr int hIndex = utils::blocks::type_index_v; + constexpr utils::blocks::value_block block{}; + const int offset = s.Problem().GetManifest().layout().offset(block); + const int size = s.Problem().GetManifest().layout().size(block); + std::vector elements; + for (int e = 0; e < f.mesh->GetNE(); ++e) + if (Schema::template attribute_belongs_to(f.mesh->GetAttribute(e))) elements.push_back(e); + const auto &surface = s.Problem().GetPressureSurfaceRows().reduced_dofs(); + solver::PreparedEnthalpyAdmissibility prepared(*f.enthalpyFes, + field::make_field_dof_map(*f.enthalpyFes),elements, + std::span(surface.GetData(),surface.Size())); + mfem::Vector physicalBase,physicalDirection,h(size),p(size),anchor(size); + s.normalizedOperator->DenormalizeState(base,physicalBase); + s.normalizedOperator->DenormalizeState(originalDirection,physicalDirection); + for (int i = 0; i < size; ++i) { + h(i) = physicalBase(offset+i); p(i) = physicalDirection(offset+i); + anchor(i) = acceptedPhysical(offset+i); + } + auto certificates = detail::Open(output / "polynomial_certificates.csv"); + certificates << "state,status,lower_bound,upper_bound,visited_subcells,maximum_depth,witness_element,x,y,z,seconds\n"; + auto check = [&](const char *name,const mfem::Vector &values) { + const auto before = std::chrono::steady_clock::now(); + const auto report = prepared.Check(values); + certificates << name << ',' << detail::Status(report.status) << ',' << report.lowerBound << ',' + << report.upperBound << ',' << report.visitedSubcells << ',' << report.maximumDepth << ',' + << report.witnessElement << ',' << report.witnessReferencePoint[0] << ',' + << report.witnessReferencePoint[1] << ',' << report.witnessReferencePoint[2] << ',' + << std::chrono::duration(std::chrono::steady_clock::now()-before).count() << '\n'; + certificates.flush(); return report; + }; + check("frozen_base",h); check("context_seed_anchor",anchor); + mfem::Vector fullCandidate(h); fullCandidate += p; + check("original_full_candidate",fullCandidate); + const auto repair = prepared.MinimumRepairFraction(h,anchor,p); + auto summary = detail::Open(output / "polynomial_repair.csv"); + summary << "status,minimum_safe_fraction,base_status,anchor_status,candidate_status\n" + << detail::RepairStatus(repair.status) << ',' << repair.minimumSafeFraction << ',' + << detail::Status(repair.base.status) << ',' << detail::Status(repair.anchor.status) << ',' + << detail::Status(repair.candidate.status) << '\n'; + summary.flush(); + auto trials = detail::Open(output / "polynomial_repair_trials.csv"); + trials << "fraction,alpha,certificate,prepared,reason,stage,eos_code,residual_norm,linear_model_norm,merit_slope,direction_change\n"; + auto evaluate = [&](double fraction,double alpha,const mfem::Vector &direction) { + restore(); + mfem::Vector action(s.Problem().EquationSize()); s.normalizedOperator->Mult(direction,action); + const double slope = baseResidual*action; + mfem::Vector model(action); model *= alpha; model += baseResidual; + mfem::Vector candidate(base); candidate.Add(alpha,direction); + mfem::Vector physicalCandidate,reduced(size); + s.normalizedOperator->DenormalizeState(candidate,physicalCandidate); + for (int i = 0; i < size; ++i) reduced(i) = physicalCandidate(offset+i); + const auto certificate = prepared.Check(reduced); + s.trialNormalizedState = candidate; ++preparations; + const auto preparation = s.PrepareTrial(); + trials << fraction << ',' << alpha << ',' << detail::Status(certificate.status) << ',' + << static_cast(preparation) << ','; + if (preparation) trials << "-1,-1,-1," << s.trialNormalizedResidual.Norml2(); + else trials << int(preparation.error().reason) << ',' << int(preparation.error().stage) << ',' + << int(preparation.error().thermodynamicErrorCode) << ",nan"; + trials << ',' << model.Norml2() << ',' << slope << ',' + << detail::Difference(direction,originalDirection) << '\n'; trials.flush(); + }; + evaluate(0.0,1.0,originalDirection); + if (repair.status == solver::EnthalpyDirectionRepairStatus::feasible) { + mfem::Vector repairedH; + const std::array fractions{repair.minimumSafeFraction, + 0.5*(repair.minimumSafeFraction+1.0),1.0}; + for (std::size_t i = 0; i < fractions.size(); ++i) { + const double fraction = fractions[i]; + const auto report = prepared.ApplyRepairFraction(h,anchor,p,fraction,repairedH); + if (report.status != solver::PolynomialSignStatus::certified_nonnegative) continue; + mfem::Vector repairedPhysical(physicalDirection),repairedNormalized; + for (int j = 0; j < size; ++j) repairedPhysical(offset+j) = repairedH(j); + s.normalizedOperator->NormalizeState(repairedPhysical,repairedNormalized); + for (int j = 0; j < repairedNormalized.Size(); ++j) + if (j < offset || j >= offset+size) repairedNormalized(j) = originalDirection(j); + detail::Save(output / ("polynomial_repaired_direction_"+std::to_string(i)+".txt"),repairedNormalized); + for (double alpha : {1.0,0.5,0.125}) evaluate(fraction,alpha,repairedNormalized); + } + } + finish(); + auto checks = detail::Open(output / "polynomial_restoration.csv"); + checks << "metric,value\naccepted_physical_difference," << detail::Difference(acceptedPhysical,s.AcceptedPhysicalState()) + << "\naccepted_normalized_difference," << detail::Difference(acceptedNormalized,s.acceptedNormalizedState) + << "\naccepted_residual_difference," << detail::Difference(acceptedResidual,s.acceptedNormalizedResidual) + << "\nrestored_residual_difference," << detail::Difference(baseResidual,s.trialNormalizedResidual) + << "\nrestored_correction_difference," << detail::Difference(savedCorrection,s.normalizedCorrection) + << "\npreparations," << preparations << "\nseconds," + << std::chrono::duration(std::chrono::steady_clock::now()-start).count() << '\n'; + if (detail::Difference(acceptedPhysical,s.AcceptedPhysicalState()) != 0.0 || + detail::Difference(acceptedNormalized,s.acceptedNormalizedState) != 0.0 || + detail::Difference(acceptedResidual,s.acceptedNormalizedResidual) != 0.0 || + detail::Difference(baseResidual,s.trialNormalizedResidual) != 0.0) + throw std::runtime_error("Polynomial positivity diagnostics changed an accepted state or failed restoration."); + } catch (...) { + const auto error = std::current_exception(); + finish(); std::rethrow_exception(error); + } +} +} diff --git a/libmeanfield/interface/solver/enthalpy_admissibility.cppm b/libmeanfield/interface/solver/enthalpy_admissibility.cppm index 1bdddef..6c58d9b 100644 --- a/libmeanfield/interface/solver/enthalpy_admissibility.cppm +++ b/libmeanfield/interface/solver/enthalpy_admissibility.cppm @@ -135,6 +135,11 @@ export namespace mean_field::solver { EnthalpyAdmissibilityOptions options = {}) const { enthalpy_detail::ValidateOptions(options); auto coefficients = Transform(nodalValues); + double rootLower = std::numeric_limits::infinity(), rootUpper = -rootLower; + for (const auto &coefficient : coefficients) { + rootLower = std::min(rootLower,coefficient.lo); + rootUpper = std::max(rootUpper,coefficient.hi); + } struct Cell { std::vector coefficients; int depth; std::array lo, hi; }; std::vector pending; pending.push_back({std::move(coefficients), 0, {0,0,0}, {1,1,1}}); @@ -164,8 +169,8 @@ export namespace mean_field::solver { } if (cell.coefficients[index].hi < 0.0) { report.status = PolynomialSignStatus::negative; - report.lowerBound = lower; - report.upperBound = upper; + report.lowerBound = rootLower; + report.upperBound = rootUpper; for (int axis = 0; axis < Dimension; ++axis) report.witnessReferencePoint[axis] = (corner & (1<GetComm())) { + report.status = EnthalpyDirectionRepairStatus::invalid_direction; + report.candidate.status = PolynomialSignStatus::non_finite; + return report; + } mfem::Vector candidate(base); candidate += direction; for (int row : m_zeroTrace) candidate(row) = 0.0; report.candidate = Check(candidate); diff --git a/tests/solver/enthalpy_admissibility.cpp b/tests/solver/enthalpy_admissibility.cpp index 40736b4..0c6cadc 100644 --- a/tests/solver/enthalpy_admissibility.cpp +++ b/tests/solver/enthalpy_admissibility.cpp @@ -138,6 +138,7 @@ TEST_CASE("Prepared admissibility preserves shared true DOFs and rejects an inva const auto base = fixture.Project([](const mfem::Vector &) { return 1.0; }); const auto anchor = fixture.Project([](const mfem::Vector &p) { return 1.0+0.2*p(1); }); const auto direction = fixture.Project([](const mfem::Vector &p) { return -2.0-p(0)+0.5*p(2); }); + const mfem::Vector savedBase(base), savedAnchor(anchor), savedDirection(direction); const auto repair = prepared.MinimumRepairFraction(base,anchor,direction); REQUIRE(repair.status == EnthalpyDirectionRepairStatus::feasible); REQUIRE(repair.minimumSafeFraction > 0.0); @@ -149,6 +150,12 @@ TEST_CASE("Prepared admissibility preserves shared true DOFs and rejects an inva CHECK(prepared.Check(candidate).status == PolynomialSignStatus::certified_nonnegative); for (int i = 0; i < repaired.Size(); ++i) CHECK(repaired(i) == Catch::Approx((1.0-beta)*direction(i)+beta*(anchor(i)-base(i)))); + for (int i = 0; i < base.Size(); ++i) { + CHECK(base(i) == savedBase(i)); CHECK(anchor(i) == savedAnchor(i)); + CHECK(direction(i) == savedDirection(i)); + } + const auto repeated = prepared.MinimumRepairFraction(base,anchor,direction); + CHECK(repeated.minimumSafeFraction == repair.minimumSafeFraction); mfem::Vector invalid(base); invalid *= -1.0; CHECK(prepared.MinimumRepairFraction(invalid,anchor,direction).status == EnthalpyDirectionRepairStatus::invalid_base); @@ -181,6 +188,10 @@ TEST_CASE("Feasible direction repair enforces exact homogeneous trace and reject REQUIRE(prepared.ApplyRepairFraction(base,anchor,direction,1.0,repaired).status == PolynomialSignStatus::certified_nonnegative); for (int row : boundary) CHECK(repaired(row) == 0.0); + auto nonfiniteDirection = direction; + for (int row : boundary) nonfiniteDirection(row) = std::numeric_limits::quiet_NaN(); + CHECK(prepared.MinimumRepairFraction(base,anchor,nonfiniteDirection).status == + EnthalpyDirectionRepairStatus::invalid_direction); auto invalidTrace = base; for (int row : boundary) invalidTrace(row) = 1.0e-20; CHECK(prepared.MinimumRepairFraction(invalidTrace,anchor,direction).status ==