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

Wire typed experimental preconditioners into Newton diagnostics

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 17:05:44 -04:00
parent dc8d4b0501
commit fcb7833844
2 changed files with 123 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

64
experiments/plateau_diagnostic.cpp
View File

@@ -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 <class Context> int RunContext(Context &context, const Options &o, const Scales &scales) {
template <bool FullDefaultStack, class Context>
int RunContext(Context &context, const Options &o, const Scales &scales) {
std::vector<operators::RootBlockDescriptor> rows, values;
mfem::Vector residualFactors, stateFactors;
std::vector<int> centralRows;
@@ -428,10 +447,14 @@ template <class Context> 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 <class Context> 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 <class Context> 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 <class Context> 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 <class Context> int RunContext(Context &context, const Options &o, cons
template <class Model> 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<Model, decltype(d)>;
static_assert(preconditioning::StellarPreconditionerRuntimeAvailableFor<Recipe, Problem>);
using Inverse = preconditioning::PreparedStellarInverseType<Recipe, Problem>;
static_assert(solver::LinearPreconditionerApplicationContractAvailable<Inverse>);
static_assert(solver::linearPreconditionerApplicationContract<Inverse> ==
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<Recipe, preconditioning::DefaultStellarPreconditioner>;
return RunContext<fullDefaultStack>(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);

69
experiments/plateau_newton_preconditioner.hpp Normal file
View File

@@ -0,0 +1,69 @@
#pragma once
#include <cmath>
#include <concepts>
#include <stdexcept>
#include <type_traits>
#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 <class Policy, bool RefinedGravity>
concept NewtonBenchmarkPolicy =
(std::same_as<Policy, mean_field::preconditioning::IndependentStellarSubsystems> && RefinedGravity) ||
(std::same_as<Policy, mean_field::preconditioning::ApproximateStellarBlockLDU> && !RefinedGravity);
template <class Policy, bool RefinedGravity>
requires NewtonBenchmarkPolicy<Policy, RefinedGravity>
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<Policy, false, RefinedGravity>(problem, options);
}
using GravitySchurNewtonPreconditioner =
NewtonBenchmarkPreconditioner<mean_field::preconditioning::IndependentStellarSubsystems, true>;
using StellarLDUNewtonPreconditioner =
NewtonBenchmarkPreconditioner<mean_field::preconditioning::ApproximateStellarBlockLDU, false>;
static_assert(mean_field::preconditioning::StellarPreconditionerPrescription<GravitySchurNewtonPreconditioner>);
static_assert(mean_field::preconditioning::StellarPreconditionerPrescription<StellarLDUNewtonPreconditioner>);
static_assert(std::is_trivially_copyable_v<GravitySchurNewtonPreconditioner>);
static_assert(std::is_trivially_copyable_v<StellarLDUNewtonPreconditioner>);
} // namespace experiment::plateau_diagnostics