tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add read-only polynomial endpoint certification mode

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 10:51:07 -04:00
parent b21fa40249
commit 7e1ab69fcd
2 changed files with 112 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
experiments/plateau_diagnostic.cpp
View File

@@ -23,6 +23,7 @@ import mean_field;
#include "plateau_numerical_strategy_diagnostics.hpp"
#include "plateau_gravity_application_diagnostics.hpp"
#include "plateau_positivity_strategy_diagnostics.hpp"
#include "plateau_polynomial_certificate.hpp"
using namespace mean_field;
namespace {
@@ -31,7 +32,7 @@ struct Options {
double n = 1., angularMomentum = 0., atol = 0., rtol = 1e-8, linearTolerance = .03;
int maxNewton = 12, maxLinear = 200, restart = 40, probeRestart = 40;
bool probes = false, noCentral = false;
bool strategyProbe = false, positivityProbe = false;
bool strategyProbe = false, positivityProbe = false, positivityCertificate = false;
std::string strategy = "baseline";
double gravityTolerance = 1e-10;
int gravityIterations = 300, gravityRestart = 80;
@@ -89,6 +90,8 @@ Options Parse(int argc, char **argv) {
o.strategyProbe = true;
else if (a == "--positivity-probe")
o.positivityProbe = true;
else if (a == "--positivity-certificate")
o.positivityCertificate = true;
else if (a == "--pressure-probe")
o.pressureProbe = true;
else if (a == "--shape-curvature")
@@ -121,6 +124,11 @@ Options Parse(int argc, char **argv) {
throw std::invalid_argument("--probe-restart must be 40 or 160");
if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe) && (o.probeState.empty() || o.probeDirection.empty()))
throw std::invalid_argument("Focused probes require --probe-state and --probe-direction");
if (o.positivityCertificate && o.probeState.empty())
throw std::invalid_argument("--positivity-certificate requires --probe-state");
if (o.positivityCertificate && (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace ||
o.strategyProbe || o.positivityProbe || !o.probeDirection.empty() || !o.replayCandidate.empty()))
throw std::invalid_argument("--positivity-certificate is a certificate-only mode; omit other probe options");
if ((!o.probeDirection.empty() || !o.replayCandidate.empty()) && o.probeState.empty())
throw std::invalid_argument("A supplied direction or candidate requires --probe-state");
if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" &&
@@ -389,7 +397,8 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
metadata << "n=" << o.n << "\nK=" << scales.K << "\nrho_c=" << scales.rho << "\nJ=" << o.angularMomentum
<< "\nM=1\nR_reference=1\nxi1=" << scales.xi << "\nmu1=" << scales.mu << "\ncentral_phase=" << !o.noCentral
<< "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart
<< "\nfrozen_probe_backend=" << (o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve")
<< "\nfrozen_probe_backend=" << (o.positivityCertificate ? "polynomial_certificate_no_linear_solve" :
o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve")
<< "\nlinear_tol=" << o.linearTolerance << "\nmax_linear=" << o.maxLinear
<< "\nmax_newton=" << o.maxNewton << "\natol=" << o.atol << "\nrtol=" << o.rtol
<< "\nstrategy=" << o.strategy << "\ngravity_tolerance=" << o.gravityTolerance
@@ -405,10 +414,18 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
if (!o.replayCandidate.empty())
metadata << "replay_candidate=" << o.replayCandidate << '\n';
metadata << "pressure_probe=" << o.pressureProbe << "\nshape_curvature=" << o.shapeCurvature
<< "\nsurface_trace=" << o.surfaceTrace << '\n';
<< "\nsurface_trace=" << o.surfaceTrace << "\npositivity_certificate=" << o.positivityCertificate << '\n';
metadata.flush();
solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
const auto snapshot = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
if (o.positivityCertificate) {
// A legacy endpoint may already violate admissibility. Its
// polynomial must be inspectable without preparing that state.
mfem::Vector normalizedSnapshot;
s.normalizedOperator->NormalizeState(snapshot,normalizedSnapshot);
experiment::plateau_diagnostics::CertifyPolynomialEnthalpy(s,normalizedSnapshot,o.output);
return;
}
s.normalizedOperator->NormalizeState(snapshot, s.trialNormalizedState);
auto prep = s.PrepareTrial();
if (!prep)

92
experiments/plateau_polynomial_certificate.hpp Normal file
View File

@@ -0,0 +1,92 @@
#pragma once
#include <chrono>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <stdexcept>
// Include after import mean_field. This callback only reads the supplied state;
// its WithState owner performs the usual accepted-state restoration afterward.
namespace experiment::plateau_diagnostics {
namespace polynomial_certificate_detail {
inline const char *Status(mean_field::solver::PolynomialSignStatus status) {
using Status = mean_field::solver::PolynomialSignStatus;
switch (status) {
case Status::certified_nonnegative: return "certified_nonnegative";
case Status::unresolved: return "unresolved";
case Status::negative: return "negative";
case Status::non_finite: return "non_finite";
}
return "unknown";
}
inline std::ofstream Open(const std::filesystem::path &path) {
if (std::filesystem::exists(path))
throw std::runtime_error("Refusing to overwrite a polynomial endpoint certificate.");
std::ofstream stream(path);
stream.exceptions(std::ios::badbit | std::ios::failbit);
stream << std::setprecision(17);
return stream;
}
inline bool Identical(const mfem::Vector &left,const mfem::Vector &right) {
return left.Size() == right.Size() && (left.Size() == 0 ||
std::memcmp(left.GetData(),right.GetData(),static_cast<std::size_t>(left.Size())*sizeof(double)) == 0);
}
}
template <class State>
void CertifyPolynomialEnthalpy(State &state,const mfem::Vector &fullNormalized,
const std::filesystem::path &output) {
namespace detail = polynomial_certificate_detail;
using Clock = std::chrono::steady_clock;
int ranks = 0;
if (MPI_Comm_size(state.Problem().GetCommunicator(),&ranks) != MPI_SUCCESS || ranks != 1)
throw std::invalid_argument("Polynomial endpoint output requires exactly one MPI rank.");
const mfem::Vector acceptedPhysical(state.AcceptedPhysicalState()),
acceptedNormalized(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual),
correction(state.normalizedCorrection), trialState(state.trialNormalizedState),
trialResidual(state.trialNormalizedResidual), candidatePhysical(state.candidatePhysicalState);
const auto generation = state.Problem().GetPreparationGeneration();
const auto start = Clock::now();
mean_field::solver::PreparedStellarEnthalpyStrategy<State> prepared(state);
const auto certificateStart = Clock::now();
const auto report = prepared.CheckState(fullNormalized);
const double checkSeconds = std::chrono::duration<double>(Clock::now()-certificateStart).count();
const double workspaceSeconds = std::chrono::duration<double>(certificateStart-start).count();
auto certificate = detail::Open(output / "polynomial_endpoint_certificate.csv");
certificate << "status,lower_bound,upper_bound,visited_subcells,maximum_depth,witness_element,x,y,z,certificate_seconds,workspace_seconds\n"
<< detail::Status(report.status) << ',' << report.lowerBound << ',' << report.upperBound << ','
<< report.visitedSubcells << ',' << report.maximumDepth << ',' << report.witnessElement << ','
<< report.witnessReferencePoint[0] << ',' << report.witnessReferencePoint[1] << ','
<< report.witnessReferencePoint[2] << ',' << checkSeconds << ',' << workspaceSeconds << '\n';
const bool sameAcceptedPhysical = detail::Identical(acceptedPhysical,state.AcceptedPhysicalState());
const bool sameAcceptedNormalized = detail::Identical(acceptedNormalized,state.acceptedNormalizedState);
const bool sameAcceptedResidual = detail::Identical(acceptedResidual,state.acceptedNormalizedResidual);
const bool sameCorrection = detail::Identical(correction,state.normalizedCorrection);
const bool sameTrialState = detail::Identical(trialState,state.trialNormalizedState);
const bool sameTrialResidual = detail::Identical(trialResidual,state.trialNormalizedResidual);
const bool sameCandidate = detail::Identical(candidatePhysical,state.candidatePhysicalState);
const bool samePreparation = generation == state.Problem().GetPreparationGeneration();
auto invariants = detail::Open(output / "polynomial_endpoint_invariants.csv");
invariants << "metric,value\naccepted_physical_bitwise_identical," << sameAcceptedPhysical
<< "\naccepted_normalized_bitwise_identical," << sameAcceptedNormalized
<< "\naccepted_residual_bitwise_identical," << sameAcceptedResidual
<< "\ncorrection_bitwise_identical," << sameCorrection
<< "\ntrial_state_bitwise_identical," << sameTrialState
<< "\ntrial_residual_bitwise_identical," << sameTrialResidual
<< "\ncandidate_physical_bitwise_identical," << sameCandidate
<< "\npreparation_generation_unchanged," << samePreparation
<< "\nhelper_trial_preparations,0\nhelper_jacobian_applications,0\nhelper_linear_solves,0\n";
auto scope = detail::Open(output / "polynomial_endpoint_scope.txt");
scope << "The helper constructs a positivity workspace and certifies only the supplied normalized state.\n"
<< "It does not prepare that state, apply a Jacobian, repair a direction, or solve any equations.\n"
<< "The existing WithState wrapper restores the accepted context after this callback returns;\n"
<< "that standard restoration is outside the helper counters and timings reported here.\n"
<< "Bounds are polynomial enclosures, not measured minimum/maximum values. Witness coordinates\n"
<< "identify a proved negative point when status is negative; unresolved is never a certificate.\n";
if (!(sameAcceptedPhysical && sameAcceptedNormalized && sameAcceptedResidual && sameCorrection &&
sameTrialState && sameTrialResidual && sameCandidate && samePreparation))
throw std::runtime_error("Polynomial endpoint certification changed a context buffer or preparation.");
}
}