Wire typed experimental preconditioners into Newton diagnostics

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 <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,8 +700,10 @@ 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 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 constexpr (FullDefaultStack)
                 if (o.strategy == "positivity-reduced") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::EliminatedGravity{gravity}});
             return run(nl::NewtonNumerics{});
         }();
@@ -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(),
+    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}}));
-    return RunContext(context, o, scales);
+        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);

experiments/plateau_newton_preconditioner.hpp

@@ -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