Integrate typed Newton strategies with checked trial gravity and feasibility

CMakeLists.txt

@@ -152,6 +152,12 @@ target_sources(mean_field
         libmeanfield/interface/solver/fields.cppm
         libmeanfield/interface/solver/linear_backend.cppm
         libmeanfield/interface/solver/newton.cppm
+    libmeanfield/interface/solver/enthalpy_admissibility.cppm
+    libmeanfield/interface/solver/gravity_subsystem.cppm
+    libmeanfield/interface/solver/gravity_reduction.cppm
+    libmeanfield/interface/solver/stellar_enthalpy_strategy.cppm
+    libmeanfield/interface/solver/newton_numerics.cppm
+    libmeanfield/interface/solver/stellar_newton_strategy.cppm
         libmeanfield/interface/solver/preconditioning_diagnostics.cppm
         libmeanfield/interface/solver/stellar_equilibrium_types.cppm
         libmeanfield/interface/solver/stellar_structure.cppm
@@ -365,6 +371,10 @@ add_executable(tests
         tests/solver/stellar_equilibrium_architecture.cpp
         tests/solver/stellar_equilibrium_architecture_internal.cpp
         tests/solver/stellar_equilibrium_runtime.cpp
+        tests/solver/enthalpy_admissibility.cpp
+        tests/solver/gravity_subsystem.cpp
+        tests/solver/gravity_reduction.cpp
+        tests/solver/stellar_numerical_strategies.cpp
 )
 
 target_link_libraries(tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost)

experiments/plateau_diagnostic.cpp

@@ -20,6 +20,8 @@ import mean_field;
 #include "plateau_preconditioner_diagnostics.hpp"
 #include "plateau_pressure_obstruction_diagnostics.hpp"
 #include "plateau_shape_curvature_diagnostics.hpp"
+#include "plateau_numerical_strategy_diagnostics.hpp"
+#include "plateau_positivity_strategy_diagnostics.hpp"
 
 using namespace mean_field;
 namespace {
@@ -28,6 +30,10 @@ 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;
+    std::string strategy = "baseline";
+    double gravityTolerance = 1e-10;
+    int gravityIterations = 300, gravityRestart = 80;
     bool pressureProbe = false, shapeCurvature = false, surfaceTrace = false;
     std::vector<int> captureIterations;
     std::string mesh = "sandbox.smesh";
@@ -70,6 +76,18 @@ Options Parse(int argc, char **argv) {
             o.probeDirection = v();
         else if (a == "--replay-candidate")
             o.replayCandidate = v();
+        else if (a == "--strategy")
+            o.strategy = v();
+        else if (a == "--gravity-tol")
+            o.gravityTolerance = std::stod(v());
+        else if (a == "--gravity-iterations")
+            o.gravityIterations = std::stoi(v());
+        else if (a == "--gravity-restart")
+            o.gravityRestart = std::stoi(v());
+        else if (a == "--strategy-probe")
+            o.strategyProbe = true;
+        else if (a == "--positivity-probe")
+            o.positivityProbe = true;
         else if (a == "--pressure-probe")
             o.pressureProbe = true;
         else if (a == "--shape-curvature")
@@ -100,10 +118,15 @@ 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.probeState.empty() || o.probeDirection.empty()))
+    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.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")
+        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::filesystem::exists(o.output))
         throw std::runtime_error("Output directory already exists");
     std::filesystem::create_directories(o.output);
@@ -227,6 +250,10 @@ void FocusedProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Ve
         for (const auto &[key, value] : experiment::plateau_diagnostics::MeasureEnthalpyStep(s, f, o.n, &physicalBase))
             bounds << key << ',' << value << '\n';
         bounds.flush();
+        if (o.strategyProbe)
+            experiment::plateau_diagnostics::ProbeNumericalStrategies(s, f, base, direction, o.output);
+        if (o.positivityProbe)
+            experiment::plateau_diagnostics::ProbePolynomialEnthalpy(s, f, base, direction, o.output);
         if (o.surfaceTrace)
             experiment::plateau_diagnostics::ProbeSurfaceTrace(s, f, o.n, base, direction, o.output);
         if (o.pressureProbe)
@@ -357,6 +384,8 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
              << "\nfrozen_probe_backend=" << (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
+             << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart
              << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n';
     metadata.flush();
     if (!o.probeState.empty()) {
@@ -537,18 +566,44 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
         });
     bool converged = false;
     if (o.maxNewton > 0) {
-        auto report = [&]() {
-            auto solve = solver::make(context,
-                                      solver::nonlinear::Newton(solver::nonlinear::NewtonOptions{
-                                          .relativeTolerance = o.rtol,
-                                          .absoluteTolerance = o.atol,
-                                          .maximumIterations = o.maxNewton,
-                                          .linearSolve = {.relativeTolerance = o.linearTolerance,
-                                                          .absoluteTolerance = 0.,
-                                                          .maximumIterations = o.maxLinear}}),
-                                      observer);
+        namespace nl = solver::nonlinear;
+        const nl::NewtonOptions options{
+            .relativeTolerance = o.rtol, .absoluteTolerance = o.atol, .maximumIterations = o.maxNewton,
+            .linearSolve = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0., .maximumIterations = o.maxLinear}};
+        const nl::GravitySolveOptions gravity{
+            .control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14,
+                        .maximumIterations = o.gravityIterations},
+            .backend = {.restartLength = o.gravityRestart}};
+        auto run = [&](auto numerics) {
+            auto solve = solver::make(context, nl::Newton(options, nl::NormalizedResidualMetric{}, numerics), observer);
             return solve.evaluate();
+        };
+        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 (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}});
+            return run(nl::NewtonNumerics{});
         }();
+        const auto &strategyDiagnostics = report.diagnostics().numericalStrategies;
+        auto strategyMetrics = File(o.output / "numerical_strategies.csv");
+        strategyMetrics << "metric,value\n";
+        const auto metric = [&](const char *name, auto value) { strategyMetrics << name << ',' << value << '\n'; };
+        metric("polynomial_checks", strategyDiagnostics.polynomialChecks);
+        metric("direction_repairs", strategyDiagnostics.directionRepairs);
+        metric("last_repair_fraction", strategyDiagnostics.lastRepairFraction);
+        metric("last_merit_directional_derivative", strategyDiagnostics.lastMeritDirectionalDerivative);
+        metric("gravity_solves", strategyDiagnostics.gravitySolves);
+        metric("gravity_iterations", strategyDiagnostics.gravityIterations);
+        metric("gravity_failures", strategyDiagnostics.gravityFailures);
+        metric("gravity_solve_seconds", strategyDiagnostics.gravitySolveSeconds);
+        metric("last_gravity_true_residual", strategyDiagnostics.lastGravityTrueResidual);
+        metric("last_gravity_correction_norm", strategyDiagnostics.lastGravityCorrectionNorm);
+        metric("last_gravity_prepared_defect", strategyDiagnostics.lastGravityPreparedDefect);
+        metric("initial_gravity_residual", strategyDiagnostics.initialGravityResidual);
+        metric("completed_initial_gravity_residual", strategyDiagnostics.completedInitialGravityResidual);
+        metric("last_lifted_relative_true_residual", strategyDiagnostics.lastLiftedRelativeTrueResidual);
         converged = report.converged();
         metadata << "converged=" << converged << "\naccepted_steps=" << report.completedNonlinearIterations() << '\n';
         if (!converged)

libmeanfield/interface/mean_field.cppm

@@ -103,3 +103,10 @@ export import :equilibrium.stellar_discretization;
 export import :operators.stellar_equilibrium_problem;
 export import :seed.stellar_equilibrium_projection;
 export import :operators.stellar_equilibrium_system;
+
+export import :solver.enthalpy_admissibility;
+export import :solver.gravity_subsystem;
+export import :solver.gravity_reduction;
+
+export import :solver.stellar_enthalpy_strategy;
+export import :solver.newton_numerics;

libmeanfield/interface/solver/newton.cppm

@@ -19,6 +19,7 @@ export module mean_field:solver.newton;
 
 export import :deformation.safe_newton_step;
 export import :solver.linear_backend;
+export import :solver.newton_numerics;
 
 export namespace mean_field::solver::nonlinear {
     /*
@@ -143,11 +144,12 @@ export namespace mean_field::solver::nonlinear {
         return {.residualNorm = residualNorm, .merit = 0.5 * residualNorm * residualNorm};
     }
 
-    template <typename Metric = NormalizedResidualMetric>
+    template <typename Metric = NormalizedResidualMetric, NewtonNumericalConfiguration Numerics = NewtonNumerics<>>
         requires std::move_constructible<std::remove_cvref_t<Metric>>
     class Newton final {
     public:
         using MetricType = std::remove_cvref_t<Metric>;
+        using NumericsType = std::remove_cvref_t<Numerics>;
 
         Newton()
             requires std::default_initializable<MetricType>
@@ -167,10 +169,12 @@ export namespace mean_field::solver::nonlinear {
 
         Newton(
             NewtonOptions options,
-            MetricType metric
+            MetricType metric,
+            NumericsType numerics = {}
         )
             : m_options(std::move(options)),
-              m_metric(std::move(metric)) {
+              m_metric(std::move(metric)),
+              m_numerics(std::move(numerics)) {
             m_options.Validate();
         }
 
@@ -182,9 +186,12 @@ export namespace mean_field::solver::nonlinear {
             return m_metric;
         }
 
+        [[nodiscard]] const NumericsType &numerics() const noexcept { return m_numerics; }
+
     private:
         NewtonOptions m_options;
         [[no_unique_address]] MetricType m_metric;
+        [[no_unique_address]] NumericsType m_numerics;
     };
 
     Newton() -> Newton<NormalizedResidualMetric>;
@@ -196,9 +203,12 @@ export namespace mean_field::solver::nonlinear {
         Metric
     ) -> Newton<std::remove_cvref_t<Metric>>;
 
+    template <typename Metric, typename Numerics>
+    Newton(NewtonOptions, Metric, Numerics) -> Newton<std::remove_cvref_t<Metric>, std::remove_cvref_t<Numerics>>;
+
     template <typename Candidate> struct IsNewtonConfiguration : std::false_type { };
 
-    template <typename Metric> struct IsNewtonConfiguration<Newton<Metric>> : std::true_type { };
+    template <typename Metric, typename Numerics> struct IsNewtonConfiguration<Newton<Metric, Numerics>> : std::true_type { };
 
     template <typename Candidate>
     concept NewtonConfiguration = IsNewtonConfiguration<std::remove_cvref_t<Candidate>>::value;
@@ -237,7 +247,8 @@ export namespace mean_field::solver::nonlinear {
         inadmissible_state,
         non_finite_state,
         non_finite_residual,
-        insufficient_decrease
+        insufficient_decrease,
+        auxiliary_solve_failure
     };
 
     struct AfterLineSearchTrial final {

libmeanfield/interface/solver/newton_numerics.cppm

@@ -0,0 +1,83 @@
+module;
+
+#include <concepts>
+#include <type_traits>
+#include <utility>
+
+export module mean_field:solver.newton_numerics;
+
+export import :solver.enthalpy_admissibility;
+export import :solver.linear_backend;
+
+export namespace mean_field::solver::nonlinear {
+    struct SampledEnthalpy final { };
+    struct PolynomialEnthalpy final {
+        EnthalpyAdmissibilityOptions options{};
+    };
+
+    struct MonolithicGravity final { };
+
+    // These tolerances control an auxiliary equation solve, not the physical
+    // model's acceptance threshold. The complete residual remains authoritative.
+    struct GravitySolveOptions final {
+        LinearSolveControl control{.relativeTolerance = 1.0e-10, .absoluteTolerance = 1.0e-14,
+                                   .maximumIterations = 300};
+        linear::FGMRESOptions backend{.restartLength = 80};
+
+        void Validate() const { control.Validate(); backend.Validate(); }
+    };
+
+    struct GravityCurvatureCorrection final {
+        GravitySolveOptions gravity{};
+    };
+
+    struct EliminatedGravity final {
+        GravitySolveOptions gravity{};
+    };
+
+    template <typename Candidate>
+    concept EnthalpyNumerics = std::same_as<Candidate, SampledEnthalpy> ||
+                              std::same_as<Candidate, PolynomialEnthalpy>;
+
+    template <typename Candidate>
+    concept GravityNumerics = std::same_as<Candidate, MonolithicGravity> ||
+                              std::same_as<Candidate, GravityCurvatureCorrection> ||
+                              std::same_as<Candidate, EliminatedGravity>;
+
+    // Types select algorithms; tolerances and the mesh-dependent work sizes
+    // remain runtime data. Adding a strategy cannot silently change a model's
+    // equation set or impose an additional specification such as rho_c.
+    template <EnthalpyNumerics Enthalpy = SampledEnthalpy, GravityNumerics Gravity = MonolithicGravity>
+    struct NewtonNumerics final {
+        using EnthalpyType = Enthalpy;
+        using GravityType = Gravity;
+        [[no_unique_address]] Enthalpy enthalpy{};
+        [[no_unique_address]] Gravity gravity{};
+    };
+
+    template <typename Enthalpy, typename Gravity>
+    NewtonNumerics(Enthalpy, Gravity) -> NewtonNumerics<Enthalpy, Gravity>;
+
+    template <typename> struct IsNewtonNumerics : std::false_type { };
+    template <typename Enthalpy, typename Gravity>
+    struct IsNewtonNumerics<NewtonNumerics<Enthalpy, Gravity>> : std::true_type { };
+    template <typename Candidate>
+    concept NewtonNumericalConfiguration = IsNewtonNumerics<std::remove_cvref_t<Candidate>>::value;
+
+    struct NumericalStrategyDiagnostics final {
+        int polynomialChecks{0};
+        int directionRepairs{0};
+        double lastRepairFraction{0.0};
+        double lastMeritDirectionalDerivative{0.0};
+        int gravitySolves{0};
+        int gravityIterations{0};
+        int gravityFailures{0};
+        double gravitySolveSeconds{0.0};
+        double lastGravityTrueResidual{0.0};
+        double lastGravityCorrectionNorm{0.0};
+        double lastGravityPreparedDefect{0.0};
+        double initialGravityResidual{0.0};
+        double completedInitialGravityResidual{0.0};
+        double lastLiftedRelativeTrueResidual{0.0};
+    };
+}

libmeanfield/interface/solver/stellar_context.cppm

@@ -32,6 +32,7 @@ export import :seed.stellar_equilibrium_projection;
 export import :solver.linear_backend;
 export import :solver.newton;
 export import :solver.stellar_structure;
+export import :solver.stellar_newton_strategy;
 
 export namespace mean_field::solver {
     template <typename Model, typename Discretization, typename PreconditionerPrescription, typename LinearBackend>
@@ -523,6 +524,10 @@ export namespace mean_field::solver {
             std::unique_ptr<PhysicalInverseType> physicalInverse;
             std::unique_ptr<NormalizedInverseType> normalizedInverse;
             std::unique_ptr<LinearBackendType> linearBackend;
+            // Only destruction is erased: strategy construction and all
+            // numerical calls retain their statically checked concrete types.
+            // Declaring this last keeps borrowed dependencies alive on teardown.
+            std::unique_ptr<detail::NewtonStrategyStorage> numericalStrategyStorage;
 
         private:
             [[nodiscard]] static ProblemType *ProblemPointer(const std::shared_ptr<Storage> &candidate) noexcept {
@@ -1209,9 +1214,37 @@ export namespace mean_field::solver {
             const auto &options = m_newton.options();
             options.Validate();
             detail::RequireCollectivelyIdenticalNewtonOptions(options, communicator);
+            detail::RequireCollectivelyIdenticalNumerics(m_newton.numerics(), communicator);
             state.BeginEvaluation();
+            using Strategy = detail::PreparedNewtonStrategy<std::remove_reference_t<decltype(state)>,
+                                                            typename NewtonType::NumericsType>;
+            static_assert(!Strategy::altersDirection ||
+                          std::same_as<typename NewtonType::MetricType, nonlinear::NormalizedResidualMetric>,
+                          "Polynomial direction repair currently requires the normalized residual merit.");
+            state.numericalStrategyStorage.reset();
+            auto preparedStrategy = std::make_unique<Strategy>(state, m_newton.numerics());
+            auto &strategy = *preparedStrategy;
+            state.numericalStrategyStorage = std::move(preparedStrategy);
 
             StellarEquilibriumEvaluationDiagnostics diagnostics{};
+            if (!strategy.CheckState(state.acceptedNormalizedState, diagnostics.numericalStrategies)) {
+                const auto initial = EvaluateMetric(state.acceptedNormalizedResidual);
+                diagnostics.initialResidualNorm = diagnostics.finalResidualNorm = initial.residualNorm;
+                return Failure(state, std::move(diagnostics), StellarEquilibriumFailureReason::inadmissible_state,
+                               "The initial enthalpy polynomial could not be certified nonnegative.");
+            }
+            try {
+                if (!strategy.Initialize(state, diagnostics.numericalStrategies)) {
+                    state.RestoreAccepted();
+                    const auto initial = EvaluateMetric(state.acceptedNormalizedResidual);
+                    diagnostics.initialResidualNorm = diagnostics.finalResidualNorm = initial.residualNorm;
+                    return Failure(state, std::move(diagnostics), StellarEquilibriumFailureReason::linear_solve_failure,
+                                   "Initial gravity completion failed its checked inner solve or prepared residual.");
+                }
+            } catch (...) {
+                state.RestoreAccepted();
+                throw;
+            }
             nonlinear::MetricEvaluation acceptedMetric = EvaluateMetric(state.acceptedNormalizedResidual);
             diagnostics.initialResidualNorm            = acceptedMetric.residualNorm;
             diagnostics.finalResidualNorm              = acceptedMetric.residualNorm;
@@ -1240,9 +1273,7 @@ export namespace mean_field::solver {
 
                 state.linearRightHandSide = state.acceptedNormalizedResidual;
                 state.linearRightHandSide *= -1.0;
-                LinearSolveReport linearReport = state.linearBackend->Solve(
-                    state.linearRightHandSide, state.normalizedCorrection, options.linearSolve
-                );
+                LinearSolveReport linearReport = strategy.Solve(state, options.linearSolve, diagnostics.numericalStrategies);
                 diagnostics.lastLinearSolve = linearReport;
                 diagnostics.totalLinearSolveSeconds += linearReport.solveSeconds;
 
@@ -1275,6 +1306,13 @@ export namespace mean_field::solver {
                     );
                 }
 
+                if (!strategy.CaptureDirection(state, diagnostics.numericalStrategies)) {
+                    NotifyAfter(iteration, nonlinear::IterationDisposition::globalization_failure, false, 0.0, 0,
+                                diagnostics.initialResidualNorm, acceptedMetric, acceptedMetric, linearReport, timings, state);
+                    return Failure(state, std::move(diagnostics), StellarEquilibriumFailureReason::globalization_failure,
+                                   "No certified enthalpy repair with a decreasing linear merit was available.");
+                }
+
                 const nonlinear::MetricEvaluation previousMetric = acceptedMetric;
                 if (options.backtracking.geometryPreflight) {
                     const Clock::time_point geometryPreflightStart   = Clock::now();
@@ -1345,9 +1383,15 @@ export namespace mean_field::solver {
                             trialDisposition = nonlinear::LineSearchTrialDisposition::non_finite_state;
                             rejectionSource  = "candidate state";
                             ++diagnostics.nonFiniteLineSearchTrials;
+                        } else if (!strategy.CheckState(state.trialNormalizedState, diagnostics.numericalStrategies)) {
+                            rejectionReason = StellarEquilibriumFailureReason::inadmissible_state;
+                            rejectionMessage = "The candidate enthalpy polynomial could not be certified nonnegative.";
+                            trialDisposition = nonlinear::LineSearchTrialDisposition::inadmissible_state;
+                            rejectionSource = "enthalpy polynomial";
+                            ++diagnostics.inadmissibleLineSearchTrials;
                         } else {
                             const Clock::time_point preparationStart = Clock::now();
-                            auto preparationResult                   = state.PrepareTrial();
+                            auto preparationResult = strategy.PrepareTrial(state, stepLength, diagnostics.numericalStrategies);
                             preparationSeconds = std::chrono::duration<double>(Clock::now() - preparationStart).count();
                             timings.trialPreparationSeconds += preparationSeconds;
                             diagnostics.totalTrialPreparationSeconds += preparationSeconds;
@@ -1422,6 +1466,12 @@ export namespace mean_field::solver {
                                         minimumJacobianDeterminant = preparationRejection.minimumJacobianDeterminant;
                                     }
                                 }
+                            } else if (!strategy.TrialSolveSucceeded()) {
+                                rejectionReason = StellarEquilibriumFailureReason::linear_solve_failure;
+                                rejectionMessage = "The candidate gravity correction failed its checked solve or prepared residual.";
+                                trialDisposition = nonlinear::LineSearchTrialDisposition::auxiliary_solve_failure;
+                                rejectionSource = "gravity correction";
+                                minimumJacobianDeterminant = preparationResult->minimumJacobianDeterminant;
                             } else if (!detail::AllStellarRanksSatisfy(
                                            AllFinite(state.trialNormalizedResidual), communicator
                                        )) {
@@ -1453,7 +1503,15 @@ export namespace mean_field::solver {
                                     ++diagnostics.nonFiniteLineSearchTrials;
                                 } else if (
                                     detail::AllStellarRanksSatisfy(
-                                        AcceptStep(acceptedMetric, trialMetric, stepLength, options.backtracking),
+                                        [&] {
+                                            if constexpr (Strategy::altersDirection) {
+                                                return trialMetric.merit <= acceptedMetric.merit +
+                                                    options.backtracking.sufficientDecrease * stepLength *
+                                                    diagnostics.numericalStrategies.lastMeritDirectionalDerivative;
+                                            } else {
+                                                return AcceptStep(acceptedMetric, trialMetric, stepLength, options.backtracking);
+                                            }
+                                        }(),
                                         communicator
                                     )
                                 ) {

libmeanfield/interface/solver/stellar_equilibrium_types.cppm

@@ -8,6 +8,7 @@ export module mean_field:solver.stellar_equilibrium_types;
 
 export import :deformation.safe_newton_step;
 export import :solver.linear_backend;
+export import :solver.newton_numerics;
 
 export namespace mean_field::solver {
     enum class StellarEquilibriumFailureReason : std::uint8_t {
@@ -60,5 +61,6 @@ export namespace mean_field::solver {
         double totalRollbackSeconds{0.0};
         std::optional<deformation::LargestSafeNewtonStepSizeEstimate> lastGeometryPreflight;
         std::optional<LinearSolveReport> lastLinearSolve;
+        nonlinear::NumericalStrategyDiagnostics numericalStrategies{};
     };
 } // namespace mean_field::solver

libmeanfield/interface/solver/stellar_newton_strategy.cppm

@@ -0,0 +1,396 @@
+module;
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <concepts>
+#include <limits>
+#include <memory>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+
+#include <mfem.hpp>
+#include <mpi.h>
+
+export module mean_field:solver.stellar_newton_strategy;
+
+export import :solver.newton_numerics;
+export import :solver.gravity_reduction;
+export import :solver.stellar_enthalpy_strategy;
+
+namespace mean_field::solver::detail {
+    struct NewtonStrategyStorage {
+        virtual ~NewtonStrategyStorage() = default;
+    };
+
+    inline double StrategyDot(const mfem::Vector &a, const mfem::Vector &b, const MPI_Comm communicator) {
+        if (a.Size() != b.Size()) { throw std::invalid_argument("Strategy dot product has incompatible vectors."); }
+        double local = a * b;
+        double global = 0.0;
+        if (MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, communicator) != MPI_SUCCESS) {
+            throw std::runtime_error("A numerical strategy could not reduce its vector product.");
+        }
+        return global;
+    }
+
+    inline double StrategyNorm(const mfem::Vector &vector, const MPI_Comm communicator) {
+        return std::sqrt(StrategyDot(vector, vector, communicator));
+    }
+
+    template <nonlinear::NewtonNumericalConfiguration Recipe>
+    void RequireCollectivelyIdenticalNumerics(const Recipe &recipe, const MPI_Comm communicator) {
+        using Enthalpy = typename Recipe::EnthalpyType;
+        using Gravity = typename Recipe::GravityType;
+        std::array<int, 8> integers{};
+        std::array<double, 2> reals{};
+        integers[0] = std::same_as<Enthalpy, nonlinear::PolynomialEnthalpy>;
+        integers[1] = std::same_as<Gravity, nonlinear::GravityCurvatureCorrection> ? 1 :
+                      std::same_as<Gravity, nonlinear::EliminatedGravity> ? 2 : 0;
+        if constexpr (std::same_as<Enthalpy, nonlinear::PolynomialEnthalpy>) {
+            integers[2] = recipe.enthalpy.options.maximumSubdivisionDepth;
+            integers[3] = recipe.enthalpy.options.maximumSubcells;
+        }
+        if constexpr (!std::same_as<Gravity, nonlinear::MonolithicGravity>) {
+            const auto &options = recipe.gravity.gravity;
+            integers[4] = options.control.maximumIterations;
+            integers[5] = options.backend.restartLength;
+            integers[6] = options.backend.printLevel;
+            reals = {options.control.relativeTolerance, options.control.absoluteTolerance};
+        }
+        integers[7] = std::isfinite(reals[0]) && std::isfinite(reals[1]);
+        std::array<int, 8> minimumIntegers{}, maximumIntegers{};
+        std::array<double, 2> minimumReals{}, maximumReals{};
+        if (MPI_Allreduce(integers.data(), minimumIntegers.data(), 8, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS ||
+            MPI_Allreduce(integers.data(), maximumIntegers.data(), 8, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS ||
+            MPI_Allreduce(reals.data(), minimumReals.data(), 2, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS ||
+            MPI_Allreduce(reals.data(), maximumReals.data(), 2, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) {
+            throw std::runtime_error("Newton could not compare numerical strategies across ranks.");
+        }
+        if (minimumIntegers != maximumIntegers || minimumReals != maximumReals ||
+            minimumIntegers[7] == 0) {
+            throw std::invalid_argument("Newton numerical strategies and controls must agree on every rank.");
+        }
+        if constexpr (!std::same_as<Gravity, nonlinear::MonolithicGravity>) { recipe.gravity.gravity.Validate(); }
+    }
+
+    template <typename State, typename Policy> class PreparedGravityStrategy;
+
+    template <typename State>
+    class PreparedGravityStrategy<State, nonlinear::MonolithicGravity> final {
+    public:
+        PreparedGravityStrategy(State &, const nonlinear::MonolithicGravity &) { }
+        bool Initialize(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; }
+        LinearSolveReport Solve(State &state, const LinearSolveControl &control) {
+            return state.linearBackend->Solve(state.linearRightHandSide, state.normalizedCorrection, control);
+        }
+        bool CaptureDirection(State &) { return true; }
+        auto PrepareTrial(State &state, double) { return state.PrepareTrial(); }
+        bool TrialSolveSucceeded() const noexcept { return true; }
+        void CopyDiagnostics(nonlinear::NumericalStrategyDiagnostics &) const { }
+    };
+
+    template <typename State>
+    class StellarGravityWorkspace {
+    public:
+        using Operation = std::remove_cvref_t<decltype(*std::declval<State>().normalizedOperator)>;
+        using Inverse = std::remove_cvref_t<decltype(*std::declval<State>().normalizedInverse)>;
+        using Problem = std::remove_cvref_t<decltype(std::declval<State>().Problem())>;
+        using Subsystem = PreparedGravitySubsystem<Operation, Inverse>;
+
+        StellarGravityWorkspace(State &state, const nonlinear::GravitySolveOptions &options)
+            : m_options(options),
+              subsystem(*state.normalizedOperator, *state.normalizedInverse,
+                        MakeGravityBlockPartition<typename Problem::FormType>(state.Problem().GetManifest().layout()),
+                        state.Problem().GetCommunicator(), linear::FGMRES(options.backend)),
+              gravityResidual(subsystem.GetPartition().GravityResidualSize()),
+              gravityCorrection(subsystem.GetPartition().GravityValueSize()),
+              expectedResidual(subsystem.GetPartition().GravityResidualSize()),
+              gravityAction(subsystem.GetPartition().GravityResidualSize()),
+              desiredResidual(subsystem.GetPartition().GravityResidualSize()),
+              fullCorrection(state.Problem().StateSize()), fullAction(state.Problem().EquationSize()) {
+            m_options.Validate();
+        }
+
+        bool TrialSolveSucceeded() const noexcept { return trialSolveSucceeded; }
+
+        void CopyDiagnostics(nonlinear::NumericalStrategyDiagnostics &report) const {
+            report.gravitySolves = subsystem.GetSolveCount();
+            report.gravityIterations = subsystem.GetIterationCount();
+            report.gravityFailures = subsystem.GetFailureCount();
+            report.gravitySolveSeconds = subsystem.GetSolveSeconds();
+            if (const auto &last = subsystem.GetLastSolveReport(); last.has_value()) {
+                report.lastGravityTrueResidual = last->trueResidualNorm;
+            }
+            report.lastGravityCorrectionNorm = lastCorrectionNorm;
+            report.lastGravityPreparedDefect = lastPreparedDefect;
+            report.lastLiftedRelativeTrueResidual = liftedRelativeDefect;
+        }
+
+    protected:
+        void Refresh(State &state) {
+            // Preparing a trial invalidates both normalized and physical
+            // inverses. Refreshing only the restricted wrapper would hide this.
+            state.physicalInverse->Refresh();
+            subsystem.Refresh();
+        }
+
+        bool SolveCompletion(State &state) {
+            subsystem.GetPartition().GatherGravityResiduals(state.trialNormalizedResidual, gravityResidual);
+            gravityResidual *= -1.0;
+            desiredResidual = 0.0;
+            gravityCorrection = 0.0;
+            lastSolve = subsystem.Solve(gravityResidual, gravityCorrection, m_options.control);
+            return lastSolve.Converged();
+        }
+
+        auto ApplyCorrection(State &state) {
+            const auto &partition = subsystem.GetPartition();
+            partition.InjectGravityValues(gravityCorrection, fullCorrection);
+            state.normalizedOperator->Mult(fullCorrection, fullAction);
+            partition.GatherGravityResiduals(state.trialNormalizedResidual, expectedResidual);
+            partition.GatherGravityResiduals(fullAction, gravityAction);
+            expectedResidual += gravityAction;
+            lastCorrectionNorm = StrategyNorm(gravityCorrection, subsystem.GetCommunicator());
+            partition.AddGravityValues(gravityCorrection, state.trialNormalizedState);
+            state.normalizedOperator->DenormalizeState(state.trialNormalizedState, state.candidatePhysicalState);
+            auto preparation = state.PrepareTrial();
+            if (preparation) {
+                partition.GatherGravityResiduals(state.trialNormalizedResidual, gravityResidual);
+                const double scale = StrategyNorm(gravityResidual, subsystem.GetCommunicator()) +
+                                     StrategyNorm(expectedResidual, subsystem.GetCommunicator()) +
+                                     StrategyNorm(desiredResidual, subsystem.GetCommunicator()) +
+                                     StrategyNorm(gravityAction, subsystem.GetCommunicator());
+                expectedResidual -= gravityResidual;
+                const double disagreement = StrategyNorm(expectedResidual, subsystem.GetCommunicator());
+                gravityResidual -= desiredResidual;
+                const double targetDefect = StrategyNorm(gravityResidual, subsystem.GetCommunicator());
+                const double tolerance = m_options.control.ConvergenceThreshold(lastSolve.rightHandSideNorm) +
+                    256.0 * std::numeric_limits<double>::epsilon() * scale;
+                // G is affine in gravity for fixed retained variables. Verify
+                // that property against a fresh preparation, including its
+                // boundary rows, before trusting a completed candidate.
+                trialSolveSucceeded = std::isfinite(disagreement) && disagreement <= tolerance &&
+                                      std::isfinite(targetDefect) && targetDefect <= tolerance;
+                lastPreparedDefect = targetDefect;
+            }
+            return preparation;
+        }
+
+        nonlinear::GravitySolveOptions m_options;
+        Subsystem subsystem;
+        mfem::Vector gravityResidual, gravityCorrection, expectedResidual, gravityAction, desiredResidual;
+        mfem::Vector fullCorrection, fullAction;
+        LinearSolveReport lastSolve{};
+        bool trialSolveSucceeded{true};
+        double lastCorrectionNorm{0.0};
+        double lastPreparedDefect{0.0};
+        double liftedRelativeDefect{0.0};
+    };
+
+    template <typename State>
+    class PreparedGravityStrategy<State, nonlinear::GravityCurvatureCorrection> final
+        : public StellarGravityWorkspace<State> {
+        using Base = StellarGravityWorkspace<State>;
+    public:
+        PreparedGravityStrategy(State &state, const nonlinear::GravityCurvatureCorrection &policy)
+            : Base(state, policy.gravity), m_curvature(this->subsystem) { }
+        bool Initialize(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; }
+        LinearSolveReport Solve(State &state, const LinearSolveControl &control) {
+            return state.linearBackend->Solve(state.linearRightHandSide, state.normalizedCorrection, control);
+        }
+        bool CaptureDirection(State &state) {
+            this->subsystem.Refresh();
+            m_curvature.CaptureBase(state.acceptedNormalizedResidual, state.normalizedCorrection);
+            return true;
+        }
+        auto PrepareTrial(State &state, const double alpha) {
+            this->trialSolveSucceeded = true;
+            auto preparation = state.PrepareTrial();
+            if (!preparation) { return preparation; }
+            this->Refresh(state);
+            m_curvature.BuildLinearPrediction(alpha, this->desiredResidual);
+            this->lastSolve = m_curvature.SolveCorrection(alpha, state.trialNormalizedResidual,
+                                                         this->gravityCorrection, this->m_options.control);
+            this->trialSolveSucceeded = this->lastSolve.Converged();
+            if (!this->trialSolveSucceeded) { return preparation; }
+            return this->ApplyCorrection(state);
+        }
+    private:
+        PreparedGravityCurvatureCorrection<typename Base::Subsystem> m_curvature;
+    };
+
+    template <typename State>
+    class PreparedGravityStrategy<State, nonlinear::EliminatedGravity> final
+        : public StellarGravityWorkspace<State> {
+        using Base = StellarGravityWorkspace<State>;
+        using ReducedOperator = GravityReducedOperator<typename Base::Subsystem>;
+        using ReducedInverse = GravityReducedPreconditioner<typename Base::Subsystem, typename Base::Inverse>;
+        using Configuration = std::remove_cvref_t<decltype(std::declval<State>().linearBackend->GetConfiguration())>;
+        using Backend = PreparedLinearBackendType<Configuration, ReducedOperator, ReducedInverse>;
+    public:
+        PreparedGravityStrategy(State &state, const nonlinear::EliminatedGravity &policy)
+            : Base(state, policy.gravity), m_operator(this->subsystem, policy.gravity.control),
+              m_inverse(m_operator, *state.normalizedInverse),
+              m_backend(prepareLinearBackend(state.linearBackend->GetConfiguration(), m_operator, m_inverse,
+                                             state.Problem().GetCommunicator())),
+              m_rhs(m_operator.Height()), m_correction(m_operator.Width()) { }
+
+        bool Initialize(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
+            const auto &partition = this->subsystem.GetPartition();
+            partition.GatherGravityResiduals(state.acceptedNormalizedResidual, this->gravityResidual);
+            diagnostics.initialGravityResidual = StrategyNorm(this->gravityResidual, this->subsystem.GetCommunicator());
+            state.trialNormalizedState = state.acceptedNormalizedState;
+            state.trialNormalizedResidual = state.acceptedNormalizedResidual;
+            this->trialSolveSucceeded = this->SolveCompletion(state);
+            if (!this->trialSolveSucceeded) { return false; }
+            const auto preparation = this->ApplyCorrection(state);
+            if (!preparation || !this->trialSolveSucceeded ||
+                !std::isfinite(StrategyNorm(state.trialNormalizedResidual, this->subsystem.GetCommunicator()))) { return false; }
+            partition.GatherGravityResiduals(state.trialNormalizedResidual, this->gravityResidual);
+            diagnostics.completedInitialGravityResidual =
+                StrategyNorm(this->gravityResidual, this->subsystem.GetCommunicator());
+            state.PrepareTrialCommit();
+            state.CommitPreparedTrial(preparation->minimumJacobianDeterminant);
+            this->subsystem.Refresh();
+            return true;
+        }
+
+        LinearSolveReport Solve(State &state, const LinearSolveControl &control) {
+            this->subsystem.Refresh();
+            try {
+                auto inner = m_operator.TryCondenseRightHandSide(state.linearRightHandSide, m_rhs);
+                if (!inner.Converged()) { return inner; }
+                this->subsystem.GetPartition().GatherRetainedValues(state.normalizedCorrection, m_correction);
+                auto report = m_backend.Solve(m_rhs, m_correction, control);
+                if (!report.Converged()) { return report; }
+                inner = m_operator.TryLiftCorrection(m_correction, state.linearRightHandSide, state.normalizedCorrection);
+                if (!inner.Converged()) { return inner; }
+                state.normalizedOperator->Mult(state.normalizedCorrection, this->fullAction);
+                this->fullAction -= state.linearRightHandSide;
+                const auto communicator = this->subsystem.GetCommunicator();
+                const double rhsNorm = StrategyNorm(state.linearRightHandSide, communicator);
+                const double defect = StrategyNorm(this->fullAction, communicator);
+                this->liftedRelativeDefect = rhsNorm > 0.0 ? defect / rhsNorm : (defect == 0.0 ? 0.0 : INFINITY);
+                // Inner Krylov truncation makes Schur actions approximate.
+                // The original equations decide whether the lifted result is
+                // actually an admissible inexact Newton correction.
+                // Keep backend recurrence/initial/true norms in their reduced
+                // coordinates. The separate lifted diagnostic reports the
+                // original-system check; these are different linear systems.
+                if (!std::isfinite(defect) || defect > control.ConvergenceThreshold(rhsNorm)) {
+                    report.status = std::isfinite(defect) ? LinearSolveStatus::backend_failure : LinearSolveStatus::non_finite;
+                }
+                return report;
+            } catch (const GravityReductionFailure &failure) {
+                return failure.GetReport();
+            }
+        }
+
+        bool CaptureDirection(State &state) {
+            // A positivity repair changes retained coordinates. Lift that
+            // repaired tangent again so gravity follows the reduced manifold.
+            this->subsystem.GetPartition().GatherRetainedValues(state.normalizedCorrection, m_correction);
+            this->lastSolve = m_operator.TryLiftCorrection(m_correction, state.linearRightHandSide,
+                                                          state.normalizedCorrection);
+            return this->lastSolve.Converged();
+        }
+
+        auto PrepareTrial(State &state, double) {
+            this->trialSolveSucceeded = true;
+            auto preparation = state.PrepareTrial();
+            if (!preparation) { return preparation; }
+            this->Refresh(state);
+            this->trialSolveSucceeded = this->SolveCompletion(state);
+            if (!this->trialSolveSucceeded) { return preparation; }
+            return this->ApplyCorrection(state);
+        }
+    private:
+        ReducedOperator m_operator;
+        ReducedInverse m_inverse;
+        Backend m_backend;
+        mfem::Vector m_rhs, m_correction;
+    };
+
+    template <typename State, typename Policy> class PreparedEnthalpyStrategy;
+
+    template <typename State>
+    class PreparedEnthalpyStrategy<State, nonlinear::SampledEnthalpy> final {
+    public:
+        PreparedEnthalpyStrategy(State &, const nonlinear::SampledEnthalpy &) { }
+        bool Check(const mfem::Vector &, nonlinear::NumericalStrategyDiagnostics &) { return true; }
+        bool Repair(State &, nonlinear::NumericalStrategyDiagnostics &) { return true; }
+    };
+
+    template <typename State>
+    class PreparedEnthalpyStrategy<State, nonlinear::PolynomialEnthalpy> final {
+    public:
+        PreparedEnthalpyStrategy(State &state, const nonlinear::PolynomialEnthalpy &policy)
+            : m_prepared(state, policy.options) { }
+        bool Check(const mfem::Vector &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
+            ++diagnostics.polynomialChecks;
+            return m_prepared.CheckState(state).status == PolynomialSignStatus::certified_nonnegative;
+        }
+        bool Repair(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
+            const auto report = m_prepared.RepairDirection(state);
+            diagnostics.lastRepairFraction = report.chosenFraction;
+            diagnostics.lastMeritDirectionalDerivative = report.meritDirectionalDerivative;
+            if (report.chosenFraction > 0.0) { ++diagnostics.directionRepairs; }
+            return report.status == EnthalpyDirectionRepairStatus::feasible &&
+                   std::isfinite(report.meritDirectionalDerivative) && report.meritDirectionalDerivative < 0.0;
+        }
+    private:
+        PreparedStellarEnthalpyStrategy<State> m_prepared;
+    };
+
+    template <typename State, nonlinear::NewtonNumericalConfiguration Recipe>
+    class PreparedNewtonStrategy final : public NewtonStrategyStorage {
+    public:
+        using Enthalpy = typename Recipe::EnthalpyType;
+        using Gravity = typename Recipe::GravityType;
+        static constexpr bool altersDirection = std::same_as<Enthalpy, nonlinear::PolynomialEnthalpy>;
+
+        PreparedNewtonStrategy(State &state, const Recipe &recipe)
+            : m_enthalpy(state, recipe.enthalpy), m_gravity(state, recipe.gravity) { }
+
+        bool CheckState(const mfem::Vector &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
+            return m_enthalpy.Check(state, diagnostics);
+        }
+        bool Initialize(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
+            const bool ready = m_gravity.Initialize(state, diagnostics);
+            m_gravity.CopyDiagnostics(diagnostics);
+            return ready;
+        }
+        LinearSolveReport Solve(State &state, const LinearSolveControl &control,
+                                nonlinear::NumericalStrategyDiagnostics &diagnostics) {
+            auto report = m_gravity.Solve(state, control);
+            m_gravity.CopyDiagnostics(diagnostics);
+            return report;
+        }
+        bool CaptureDirection(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
+            if (!m_enthalpy.Repair(state, diagnostics)) { return false; }
+            if constexpr (!std::same_as<Gravity, nonlinear::EliminatedGravity> || altersDirection) {
+                if (!m_gravity.CaptureDirection(state)) { m_gravity.CopyDiagnostics(diagnostics); return false; }
+            }
+            if constexpr (altersDirection && std::same_as<Gravity, nonlinear::EliminatedGravity>) {
+                mfem::Vector action(state.Problem().EquationSize());
+                state.normalizedOperator->Mult(state.normalizedCorrection, action);
+                diagnostics.lastMeritDirectionalDerivative =
+                    StrategyDot(state.acceptedNormalizedResidual, action, state.Problem().GetCommunicator());
+                if (!std::isfinite(diagnostics.lastMeritDirectionalDerivative) ||
+                    diagnostics.lastMeritDirectionalDerivative >= 0.0) { return false; }
+            }
+            m_gravity.CopyDiagnostics(diagnostics);
+            return true;
+        }
+        auto PrepareTrial(State &state, const double alpha, nonlinear::NumericalStrategyDiagnostics &diagnostics) {
+            auto preparation = m_gravity.PrepareTrial(state, alpha);
+            m_gravity.CopyDiagnostics(diagnostics);
+            return preparation;
+        }
+        bool TrialSolveSucceeded() const noexcept { return m_gravity.TrialSolveSucceeded(); }
+    private:
+        PreparedEnthalpyStrategy<State, Enthalpy> m_enthalpy;
+        PreparedGravityStrategy<State, Gravity> m_gravity;
+    };
+}