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

Add bounded fully coupled second-order trial correction and frozen screen

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-18 09:41:00 -04:00
parent 2944dd5019
commit 0d87b8b919
14 changed files with 2442 additions and 35 deletions
Download Patch File Download Diff File Expand all files Collapse all files

127
experiments/plateau_diagnostic.cpp
View File

@@ -24,6 +24,8 @@ import mean_field;
#include "plateau_shape_curvature_diagnostics.hpp"
#include "plateau_numerical_strategy_diagnostics.hpp"
#include "plateau_nonlinear_step_diagnostics.hpp"
#include "plateau_coupled_trial_diagnostics.hpp"
#include "plateau_coupled_trial_physical.hpp"
#include "plateau_gravity_application_diagnostics.hpp"
#include "plateau_positivity_strategy_diagnostics.hpp"
#include "plateau_polynomial_certificate.hpp"
@@ -37,6 +39,9 @@ struct Options {
bool probes = false, noCentral = false;
bool strategyProbe = false, positivityProbe = false, positivityCertificate = false;
bool preconditionerProbe = false, nonlinearStepProbe = false;
bool coupledTrialProbe = false, coupledSeedControl = false;
double coupledTolerance = .03, coupledMaximumRatio = 1.;
int coupledIterations = 160, coupledRestart = 160, coupledAttempts = 1;
std::string preconditionerVariant = "independent";
int schurCorrectionSweeps = 1;
double schurRelaxation = 0.5;
@@ -103,6 +108,20 @@ Options Parse(int argc, char **argv) {
o.preconditionerProbe = true;
else if (a == "--nonlinear-step-probe")
o.nonlinearStepProbe = true;
else if (a == "--coupled-trial-probe")
o.coupledTrialProbe = true;
else if (a == "--coupled-seed-control")
o.coupledSeedControl = true;
else if (a == "--coupled-tol")
o.coupledTolerance = std::stod(v());
else if (a == "--coupled-iterations")
o.coupledIterations = std::stoi(v());
else if (a == "--coupled-restart")
o.coupledRestart = std::stoi(v());
else if (a == "--coupled-max-ratio")
o.coupledMaximumRatio = std::stod(v());
else if (a == "--coupled-attempts")
o.coupledAttempts = std::stoi(v());
else if (a == "--preconditioner-variant")
o.preconditionerVariant = v();
else if (a == "--schur-sweeps")
@@ -139,17 +158,17 @@ Options Parse(int argc, char **argv) {
throw std::invalid_argument("Invalid options");
if (o.probeRestart != 40 && o.probeRestart != 160)
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()))
if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.coupledTrialProbe) && (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.preconditionerProbe || o.nonlinearStepProbe || !o.probeDirection.empty() || !o.replayCandidate.empty()))
o.strategyProbe || o.positivityProbe || o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || !o.probeDirection.empty() || !o.replayCandidate.empty()))
throw std::invalid_argument("--positivity-certificate is a certificate-only mode; omit other probe options");
if (o.preconditionerProbe) {
experiment::plateau_diagnostics::ParsePreconditionerBenchmarkVariant(o.preconditionerVariant);
if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace ||
o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.nonlinearStepProbe || !o.probeDirection.empty() ||
o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.nonlinearStepProbe || o.coupledTrialProbe || !o.probeDirection.empty() ||
!o.replayCandidate.empty() || !o.captureIterations.empty() || o.strategy != "baseline")
throw std::invalid_argument("--preconditioner-probe requires --max-newton 0 and no nonlinear/probe strategies");
if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0)
@@ -157,25 +176,37 @@ Options Parse(int argc, char **argv) {
}
if (o.nonlinearStepProbe) {
if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace ||
o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe ||
o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || o.coupledTrialProbe ||
!o.probeDirection.empty() || !o.replayCandidate.empty() || !o.captureIterations.empty() ||
o.strategy != "baseline")
throw std::invalid_argument("--nonlinear-step-probe requires --max-newton 0 and no other probes or strategy; it compares raw and curvature paths itself");
if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0)
throw std::invalid_argument("A nonlinear-step probe requires a finite positive --linear-tol");
}
if (o.coupledTrialProbe) {
if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace ||
o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe ||
o.nonlinearStepProbe || !o.replayCandidate.empty() || !o.captureIterations.empty() ||
o.coupledAttempts != 1 || (o.strategy != "coupled" && o.strategy != "curvature-coupled"))
throw std::invalid_argument("--coupled-trial-probe requires --max-newton 0 and strategy coupled or curvature-coupled, without other probes");
if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.)
throw std::invalid_argument("A coupled trial probe requires a finite positive outer --linear-tol");
} else if (o.coupledSeedControl) {
throw std::invalid_argument("--coupled-seed-control requires --coupled-trial-probe");
}
if (!o.preconditionerProbe) {
if (o.preconditionerVariant != "independent" && o.preconditionerVariant != "gravity-schur" &&
o.preconditionerVariant != "stellar-ldu" && o.preconditionerVariant != "stellar-ldu-gravity-schur")
throw std::invalid_argument("Newton preconditioner variants are independent, gravity-schur, stellar-ldu, or stellar-ldu-gravity-schur");
if (o.preconditionerVariant != "independent") {
if (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe ||
o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe) || !o.probeDirection.empty() ||
o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe) ||
(!o.probeDirection.empty() && !o.coupledTrialProbe) ||
!o.replayCandidate.empty())
throw std::invalid_argument("A custom Newton preconditioner cannot be combined with frozen or certificate probes");
if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" &&
o.strategy != "positivity-curvature")
throw std::invalid_argument("Custom Newton preconditioners support baseline, positivity, curvature, or positivity-curvature");
o.strategy != "positivity-curvature" && o.strategy != "coupled" && o.strategy != "curvature-coupled")
throw std::invalid_argument("Unsupported custom Newton preconditioner strategy");
if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0 || o.schurCorrectionSweeps < 0 ||
!std::isfinite(o.schurRelaxation) || o.schurRelaxation <= 0.0)
throw std::invalid_argument("Invalid custom Newton preconditioner controls");
@@ -184,10 +215,14 @@ Options Parse(int argc, char **argv) {
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" &&
o.strategy != "reduced" && o.strategy != "positivity-curvature" && o.strategy != "positivity-reduced")
o.strategy != "reduced" && o.strategy != "positivity-curvature" && o.strategy != "positivity-reduced" &&
o.strategy != "coupled" && o.strategy != "curvature-coupled")
throw std::invalid_argument("Unknown numerical strategy");
if (!std::isfinite(o.gravityTolerance) || o.gravityTolerance <= 0.0 || o.gravityIterations < 1 || o.gravityRestart < 1)
throw std::invalid_argument("Invalid gravity solve controls");
if (!std::isfinite(o.coupledTolerance) || o.coupledTolerance <= 0. || o.coupledIterations < 1 ||
o.coupledRestart < 1 || !std::isfinite(o.coupledMaximumRatio) || o.coupledMaximumRatio <= 0. || o.coupledAttempts < 1)
throw std::invalid_argument("Invalid coupled correction controls");
if (std::filesystem::exists(o.output))
throw std::runtime_error("Output directory already exists");
std::filesystem::create_directories(o.output);
@@ -453,6 +488,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
<< "\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.preconditionerProbe ? "direct_mfem_fgmres_preconditioner_benchmark" :
o.coupledTrialProbe ? "supplied_outer_direction_verified_with_full_J_and_auxiliary_fgmres" :
o.nonlinearStepProbe ? "prepared_fgmres_cold_with_independent_true_gate" :
o.positivityCertificate ? "polynomial_certificate_no_linear_solve" :
o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve")
@@ -462,6 +498,13 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
<< "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart
<< "\npreconditioner_probe=" << o.preconditionerProbe
<< "\nnonlinear_step_probe=" << o.nonlinearStepProbe
<< "\ncoupled_trial_probe=" << o.coupledTrialProbe
<< "\ncoupled_seed_control=" << o.coupledSeedControl
<< "\ncoupled_tolerance=" << o.coupledTolerance
<< "\ncoupled_iterations=" << o.coupledIterations
<< "\ncoupled_restart=" << o.coupledRestart
<< "\ncoupled_maximum_ratio=" << o.coupledMaximumRatio
<< "\ncoupled_attempts_per_iteration=" << o.coupledAttempts
<< "\npreconditioner_variant=" << o.preconditionerVariant
<< "\ncontext_preconditioner_variant=" << (o.preconditionerProbe ? "independent" : o.preconditionerVariant)
<< "\ncontext_gravity_schur_sweeps="
@@ -470,6 +513,32 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
<< "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation
<< "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n';
metadata.flush();
if (o.coupledTrialProbe) {
metadata << "coupled_trial_input=physical_snapshot_and_supplied_normalized_direction\n"
<< "coupled_trial_admissibility=sampled_enthalpy_without_repair\n"
<< "coupled_trial_alphas_are_independent=1\n"
<< "probe_state=" << o.probeState << "\nprobe_direction=" << o.probeDirection << '\n';
metadata.flush();
solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
mfem::Vector base;
s.normalizedOperator->NormalizeState(physical, base);
const auto direction = LoadVector(o.probeDirection, base.Size());
experiment::plateau_diagnostics::ProbeCoupledTrial(s, base, direction, o.output, {
.outerLinearTolerance = o.linearTolerance,
.gravity = {.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14,
.maximumIterations = o.gravityIterations},
.backend = {.restartLength = o.gravityRestart, .printLevel = -1}},
.coupled = {.control = {.relativeTolerance = o.coupledTolerance, .absoluteTolerance = 1e-14,
.maximumIterations = o.coupledIterations},
.backend = {.restartLength = o.coupledRestart, .printLevel = -1},
.maximumCorrectionRatio = o.coupledMaximumRatio},
.gravityPreconditionPath = o.strategy == "curvature-coupled",
.seedControl = o.coupledSeedControl});
experiment::plateau_diagnostics::MeasureCoupledTrialSnapshots(s, f, o.output, o.n, scales.K);
});
return 0;
}
if (o.nonlinearStepProbe) {
metadata << "nonlinear_step_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot")
<< "\nnonlinear_step_paths=raw,curvature\nnonlinear_step_positivity_repair=0\n";
@@ -600,6 +669,12 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
blocklog << "iteration,stage,block,normalized_norm,physical_norm,max_abs_normalized\n";
auto triallog = File(o.output / "trials.csv");
triallog << "iteration,trial,alpha,disposition,source,F,min_jacobian\n";
auto coupledLog = File(o.output / "coupled_trials.csv");
coupledLog << "iteration,trial,alpha,diagnostics_available,status,attempts_total,kept_total,fallbacks_total,"
"iterations_total,preparations_total,restorations_total,solve_seconds_total,verification_seconds_total,"
"preparation_seconds_total,refresh_seconds_total,original_merit,corrected_merit,has_solve,"
"linear_status,linear_iterations,rhs_norm,reported_true_norm,verified_true_norm,threshold,verified_gate,"
"correction_norm,maximum_correction_norm\n";
auto stateConsistency = File(o.output / "trial_state_consistency.csv");
stateConsistency << "iteration,trial,normalized_physical_roundtrip_error,relative_error\n";
mfem::Vector captureBase;
@@ -643,6 +718,23 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
std::cout << "Newton " << e.iteration << " F=" << e.residualNorm << std::endl;
},
[&](const solver::nonlinear::AfterLineSearchTrial &e) {
coupledLog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << (e.numericalStrategies != nullptr);
if (e.numericalStrategies) {
const auto &c = e.numericalStrategies->coupledTrial;
coupledLog << ',' << static_cast<int>(c.status) << ',' << c.attempts << ',' << c.kept << ','
<< c.fallbacks << ',' << c.iterations << ',' << c.preparations << ',' << c.restorations << ','
<< c.solveSeconds << ',' << c.verificationSeconds << ',' << c.preparationSeconds << ','
<< c.refreshSeconds << ',' << c.originalMerit << ',' << c.correctedMerit << ',' << c.currentSolve.has_value();
if (c.currentSolve) {
const auto &r = *c.currentSolve;
coupledLog << ',' << static_cast<int>(r.linearSolve.status) << ',' << r.linearSolve.iterations << ','
<< r.linearSolve.rightHandSideNorm << ',' << r.linearSolve.trueResidualNorm << ','
<< r.verifiedTrueResidualNorm << ',' << r.convergenceThreshold << ',' << r.verifiedGate << ','
<< r.correctionNorm << ',' << r.maximumCorrectionNorm;
} else coupledLog << ",-1,0,nan,nan,nan,nan,0,nan,nan";
} else coupledLog << ",-1,0,0,0,0,0,0,0,0,0,0,nan,nan,0,-1,0,nan,nan,nan,nan,0,nan,nan";
coupledLog << '\n';
coupledLog.flush();
if (e.candidatePhysicalState.size() != e.candidateNormalizedState.size() ||
e.candidateNormalizedState.size() != static_cast<size_t>(stateFactors.Size()))
throw std::runtime_error("Trial observer exposed incompatible state representations");
@@ -731,6 +823,12 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14,
.maximumIterations = o.gravityIterations},
.backend = {.restartLength = o.gravityRestart, .printLevel = -1}};
const nl::CoupledCurvatureCorrection coupled{
.correction = {.control = {.relativeTolerance = o.coupledTolerance, .absoluteTolerance = 1e-14,
.maximumIterations = o.coupledIterations},
.backend = {.restartLength = o.coupledRestart, .printLevel = -1},
.maximumCorrectionRatio = o.coupledMaximumRatio},
.maximumAttemptsPerIteration = o.coupledAttempts};
auto run = [&](auto numerics) {
auto solve = solver::make(context, nl::Newton(options, nl::NormalizedResidualMetric{}, numerics), observer);
return solve.evaluate();
@@ -738,6 +836,8 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
auto report = [&] {
if (o.strategy == "positivity") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::MonolithicGravity{}});
if (o.strategy == "curvature") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}});
if (o.strategy == "coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::MonolithicGravity{}, coupled});
if (o.strategy == "curvature-coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}, coupled});
if constexpr (FullDefaultStack)
if (o.strategy == "reduced") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::EliminatedGravity{gravity}});
if (o.strategy == "positivity-curvature") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::GravityCurvatureCorrection{gravity}});
@@ -770,6 +870,17 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
metric("initial_gravity_residual", strategyDiagnostics.initialGravityResidual);
metric("completed_initial_gravity_residual", strategyDiagnostics.completedInitialGravityResidual);
metric("last_lifted_relative_true_residual", strategyDiagnostics.lastLiftedRelativeTrueResidual);
const auto &coupledDiagnostics = strategyDiagnostics.coupledTrial;
metric("coupled_attempts", coupledDiagnostics.attempts);
metric("coupled_kept", coupledDiagnostics.kept);
metric("coupled_fallbacks", coupledDiagnostics.fallbacks);
metric("coupled_iterations", coupledDiagnostics.iterations);
metric("coupled_preparations", coupledDiagnostics.preparations);
metric("coupled_restorations", coupledDiagnostics.restorations);
metric("coupled_solve_seconds", coupledDiagnostics.solveSeconds);
metric("coupled_verification_seconds", coupledDiagnostics.verificationSeconds);
metric("coupled_preparation_seconds", coupledDiagnostics.preparationSeconds);
metric("coupled_refresh_seconds", coupledDiagnostics.refreshSeconds);
converged = report.converged();
metadata << "converged=" << converged << "\naccepted_steps=" << report.completedNonlinearIterations() << '\n';
if (!converged)