From fcb7833844e1636572f03c001e7fe51243789349 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 17:05:44 -0400 Subject: [PATCH] Wire typed experimental preconditioners into Newton diagnostics --- experiments/plateau_diagnostic.cpp | 64 ++++++++++++++--- experiments/plateau_newton_preconditioner.hpp | 69 +++++++++++++++++++ 2 files changed, 123 insertions(+), 10 deletions(-) create mode 100644 experiments/plateau_newton_preconditioner.hpp diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 37415d3..1663f48 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -19,6 +19,7 @@ import mean_field; #include "plateau_physical_diagnostics.hpp" #include "plateau_preconditioner_diagnostics.hpp" #include "plateau_preconditioner_benchmark.hpp" +#include "plateau_newton_preconditioner.hpp" #include "plateau_pressure_obstruction_diagnostics.hpp" #include "plateau_shape_curvature_diagnostics.hpp" #include "plateau_numerical_strategy_diagnostics.hpp" @@ -151,6 +152,23 @@ Options Parse(int argc, char **argv) { if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0) throw std::invalid_argument("A preconditioner benchmark requires a finite positive --linear-tol"); } + if (!o.preconditionerProbe) { + if (o.preconditionerVariant != "independent" && o.preconditionerVariant != "gravity-schur" && + o.preconditionerVariant != "stellar-ldu") + throw std::invalid_argument("Newton preconditioner variants are independent, gravity-schur, or stellar-ldu"); + if (o.preconditionerVariant != "independent") { + if (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || + o.positivityProbe || o.positivityCertificate || !o.probeState.empty() || !o.probeDirection.empty() || + !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"); + 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"); + } + } 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" && @@ -410,7 +428,8 @@ void FrozenProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vec } } -template int RunContext(Context &context, const Options &o, const Scales &scales) { +template +int RunContext(Context &context, const Options &o, const Scales &scales) { std::vector rows, values; mfem::Vector residualFactors, stateFactors; std::vector centralRows; @@ -428,10 +447,14 @@ template int RunContext(Context &context, const Options &o, cons << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart << "\npreconditioner_probe=" << o.preconditionerProbe << "\npreconditioner_variant=" << o.preconditionerVariant + << "\ncontext_preconditioner_variant=" << (o.preconditionerProbe ? "independent" : o.preconditionerVariant) + << "\ncontext_gravity_schur_sweeps=" + << (!o.preconditionerProbe && o.preconditionerVariant == "gravity-schur" ? o.schurCorrectionSweeps : 0) << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; metadata.flush(); - if (o.preconditionerProbe) { + // Custom Newton recipes do not instantiate the nine-variant frozen suite. + if constexpr (FullDefaultStack) if (o.preconditionerProbe) { metadata << "benchmark_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n'; if (!o.probeState.empty()) metadata << "probe_state=" << o.probeState << '\n'; solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &) { @@ -456,7 +479,7 @@ template int RunContext(Context &context, const Options &o, cons // the experiment. The CSV retains both recurrence and true gate results. return 0; } - if (!o.probeState.empty()) { + if constexpr (FullDefaultStack) if (!o.probeState.empty()) { if (o.maxNewton != 0) throw std::invalid_argument("--probe-state requires --max-newton 0"); metadata << "probe_state=" << o.probeState << '\n'; @@ -677,9 +700,11 @@ template int RunContext(Context &context, const Options &o, cons 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 == "reduced") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::EliminatedGravity{gravity}}); + 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}}); - if (o.strategy == "positivity-reduced") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::EliminatedGravity{gravity}}); + if constexpr (FullDefaultStack) + if (o.strategy == "positivity-reduced") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::EliminatedGravity{gravity}}); return run(nl::NewtonNumerics{}); }(); const auto &strategyDiagnostics = report.diagnostics().numericalStrategies; @@ -723,7 +748,7 @@ template int RunContext(Context &context, const Options &o, cons sh << (j ? "," : "") << schur(i, j); sh << '\n'; } - if (o.probes) { + if constexpr (FullDefaultStack) if (o.probes) { experiment::plateau_diagnostics::InspectPreconditioner(s, o.output); FrozenProbes(s, f, o); } @@ -733,10 +758,29 @@ template int RunContext(Context &context, const Options &o, cons template int Setup(Model model, fem::FEM f, const Options &o, const Scales &scales) { auto d = equilibrium::makeStellarDiscretization( std::move(f), normalization::PhysicalRieszDiagonal{dimensions::LengthValue{1.}, utils::G}); - auto context = solver::makeContext(std::move(model), std::move(d), preconditioning::makePreconditioner(), - solver::linear::FGMRES({.restartLength = o.restart, .printLevel = -1}), - seed::LaneEmden({.centralDensity = dimensions::DensityValue{scales.rho}})); - return RunContext(context, o, scales); + const auto withRecipe = [&](auto recipe) { + using Recipe = decltype(recipe); + using Problem = equilibrium::StellarEquilibriumProblem; + static_assert(preconditioning::StellarPreconditionerRuntimeAvailableFor); + using Inverse = preconditioning::PreparedStellarInverseType; + static_assert(solver::LinearPreconditionerApplicationContractAvailable); + static_assert(solver::linearPreconditionerApplicationContract == + preconditioning::ApplicationContract::stationary_linear); + auto context = solver::makeContext(std::move(model), std::move(d), std::move(recipe), + solver::linear::FGMRES({.restartLength = o.restart, .printLevel = -1}), + seed::LaneEmden({.centralDensity = dimensions::DensityValue{scales.rho}})); + constexpr bool fullDefaultStack = std::same_as; + return RunContext(context, o, scales); + }; + if (o.preconditionerProbe || o.preconditionerVariant == "independent") + return withRecipe(preconditioning::makePreconditioner()); + if (o.preconditionerVariant == "gravity-schur") + return withRecipe(experiment::plateau_diagnostics::GravitySchurNewtonPreconditioner{ + {.correctionSweeps = o.schurCorrectionSweeps, .relaxation = o.schurRelaxation}}); + if (o.preconditionerVariant == "stellar-ldu") + return withRecipe(experiment::plateau_diagnostics::StellarLDUNewtonPreconditioner{ + {.correctionSweeps = 0, .relaxation = 1.0}}); + throw std::logic_error("An unsupported Newton preconditioner bypassed option validation."); } int Run(const Options &o) { const auto scales = Derive(o.n); diff --git a/experiments/plateau_newton_preconditioner.hpp b/experiments/plateau_newton_preconditioner.hpp new file mode 100644 index 0000000..405a6d7 --- /dev/null +++ b/experiments/plateau_newton_preconditioner.hpp @@ -0,0 +1,69 @@ +#pragma once + +#include +#include +#include +#include + +#include "plateau_preconditioner_benchmark.hpp" + +// Include after import mean_field. These owning recipes belong to the +// experiment driver; the public physics-facing solver factory stays unchanged. +namespace experiment::plateau_diagnostics { + template + concept NewtonBenchmarkPolicy = + (std::same_as && RefinedGravity) || + (std::same_as && !RefinedGravity); + + template + requires NewtonBenchmarkPolicy + class NewtonBenchmarkPreconditioner final + : public mean_field::preconditioning::StellarPreconditionerPrescriptionTag { + public: + explicit NewtonBenchmarkPreconditioner(mean_field::preconditioning::GravityConsistentSchurLDU gravity) + : m_gravity(gravity) { + if (gravity.correctionSweeps < 0 || !std::isfinite(gravity.relaxation) || gravity.relaxation <= 0.0) + throw std::invalid_argument("Newton benchmark gravity controls must be finite and valid."); + } + + [[nodiscard]] mean_field::preconditioning::GravityConsistentSchurLDU GravityPolicy() const noexcept { + return m_gravity; + } + + private: + // A recipe can outlive command-line parsing. Keep numeric controls by + // value; prepared objects acquire their problem binding only via ADL. + mean_field::preconditioning::GravityConsistentSchurLDU m_gravity; + }; + + template < + class Policy, + bool RefinedGravity, + mean_field::preconditioning::DefaultStellarPreconditionerAvailableFor Problem> + requires NewtonBenchmarkPolicy< + Policy, + RefinedGravity> + [[nodiscard]] auto prepareStellarPreconditioner( + NewtonBenchmarkPreconditioner< + Policy, + RefinedGravity> prescription, + const Problem &problem + ) { + PreconditionerBenchmarkOptions options; + options.schurCorrectionSweeps = prescription.GravityPolicy().correctionSweeps; + options.schurRelaxation = prescription.GravityPolicy().relaxation; + // Returning the benchmark's concrete prepared type lets the existing + // context own refreshes and rebuilt scalar-border completion unchanged. + return preconditioner_benchmark_detail::Prepare(problem, options); + } + + using GravitySchurNewtonPreconditioner = + NewtonBenchmarkPreconditioner; + using StellarLDUNewtonPreconditioner = + NewtonBenchmarkPreconditioner; + + static_assert(mean_field::preconditioning::StellarPreconditionerPrescription); + static_assert(mean_field::preconditioning::StellarPreconditionerPrescription); + static_assert(std::is_trivially_copyable_v); + static_assert(std::is_trivially_copyable_v); +} // namespace experiment::plateau_diagnostics