Benchmark isolated and coupled gravity Schur refinements

experiments/plateau_diagnostic.cpp

@@ -36,6 +36,8 @@ struct Options {
     bool strategyProbe = false, positivityProbe = false, positivityCertificate = false;
     bool preconditionerProbe = false;
     std::string preconditionerVariant = "independent";
+    int schurCorrectionSweeps = 1;
+    double schurRelaxation = 0.5;
     std::string strategy = "baseline";
     double gravityTolerance = 1e-10;
     int gravityIterations = 300, gravityRestart = 80;
@@ -99,6 +101,10 @@ Options Parse(int argc, char **argv) {
             o.preconditionerProbe = true;
         else if (a == "--preconditioner-variant")
             o.preconditionerVariant = v();
+        else if (a == "--schur-sweeps")
+            o.schurCorrectionSweeps = std::stoi(v());
+        else if (a == "--schur-relaxation")
+            o.schurRelaxation = std::stod(v());
         else if (a == "--pressure-probe")
             o.pressureProbe = true;
         else if (a == "--shape-curvature")
@@ -422,6 +428,7 @@ 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
+             << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation
              << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n';
     metadata.flush();
     if (o.preconditionerProbe) {
@@ -442,7 +449,8 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
             metadata.flush();
             experiment::plateau_diagnostics::BenchmarkPreconditioners(s, o.output,
                 {.variant = o.preconditionerVariant, .restart = o.restart,
-                 .maximumIterations = o.maxLinear, .relativeTolerance = o.linearTolerance});
+                 .maximumIterations = o.maxLinear, .relativeTolerance = o.linearTolerance,
+                 .schurCorrectionSweeps = o.schurCorrectionSweeps, .schurRelaxation = o.schurRelaxation});
         });
         // A missed linear gate is a measured outcome, not failure to execute
         // the experiment. The CSV retains both recurrence and true gate results.

experiments/plateau_preconditioner_benchmark.hpp

@@ -5,6 +5,7 @@
 #include <cmath>
 #include <concepts>
 #include <cstring>
+#include <cstdint>
 #include <filesystem>
 #include <fstream>
 #include <iomanip>
@@ -23,7 +24,9 @@ enum class PreconditionerBenchmarkVariant {
     material_then_gravity,
     gravity_then_material,
     stellar_ldu,
-    material_ldu
+    material_ldu,
+    gravity_schur,
+    stellar_ldu_gravity_schur
 };
 
 inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const std::string &name) {
@@ -32,6 +35,8 @@ inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const
     if (name == "gravity-then-material") return PreconditionerBenchmarkVariant::gravity_then_material;
     if (name == "stellar-ldu") return PreconditionerBenchmarkVariant::stellar_ldu;
     if (name == "material-ldu") return PreconditionerBenchmarkVariant::material_ldu;
+    if (name == "gravity-schur") return PreconditionerBenchmarkVariant::gravity_schur;
+    if (name == "stellar-ldu-gravity-schur") return PreconditionerBenchmarkVariant::stellar_ldu_gravity_schur;
     throw std::invalid_argument("Unknown preconditioner benchmark variant: " + name);
 }
 
@@ -40,6 +45,8 @@ struct PreconditionerBenchmarkOptions final {
     int restart{160};
     int maximumIterations{200};
     double relativeTolerance{0.03};
+    int schurCorrectionSweeps{1};
+    double schurRelaxation{0.5};
 };
 
 namespace preconditioner_benchmark_detail {
@@ -87,8 +94,9 @@ inline double Relative(double numerator, double denominator) {
         (numerator == 0.0 ? 0.0 : std::numeric_limits<double>::infinity());
 }
 
-template <pc::StellarStructureFactorizationPolicy StructurePolicy, bool MaterialLDU = false, class Problem>
+template <pc::StellarStructureFactorizationPolicy StructurePolicy, bool MaterialLDU = false,
-auto Prepare(const Problem &problem) {
+          bool RefinedGravity = false, class Problem>
+auto Prepare(const Problem &problem, const PreconditionerBenchmarkOptions &options) {
     const auto defaults = pc::stellarStructureBlock(problem);
     const auto material = [&]() {
         if constexpr (MaterialLDU) {
@@ -99,7 +107,17 @@ auto Prepare(const Problem &problem) {
             return defaults.materialSurfaceComponent();
         }
     }();
-    auto structure = pc::stellarStructureBlock(problem, material, defaults.gravityComponent(), StructurePolicy{});
+    auto gravity = [&]() {
+        if constexpr (RefinedGravity) {
+            const auto &original = defaults.gravityComponent();
+            return pc::GravityFieldBlock(original.massInverseBackend(), original.potentialSchurBackend(),
+                pc::GravityConsistentSchurLDU{.correctionSweeps = options.schurCorrectionSweeps,
+                                              .relaxation = options.schurRelaxation});
+        } else {
+            return defaults.gravityComponent();
+        }
+    }();
+    auto structure = pc::stellarStructureBlock(problem, material, std::move(gravity), StructurePolicy{});
     // Each candidate needs its own completed border columns and scalar Schur;
     // reusing the default border would confound the factorization comparison.
     return pc::prepare(problem, pc::specificationBorderBlock(
@@ -145,7 +163,8 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         throw std::invalid_argument("Frozen preconditioner benchmarks require one MPI rank.");
     if (options.restart <= 0 || options.maximumIterations <= 0 ||
         !std::isfinite(options.relativeTolerance) || options.relativeTolerance <= 0.0 ||
-        !state.normalizedOperator->IsPrepared())
+        options.schurCorrectionSweeps < 0 || !std::isfinite(options.schurRelaxation) ||
+        options.schurRelaxation <= 0.0 || !state.normalizedOperator->IsPrepared())
         throw std::invalid_argument("Frozen preconditioner benchmark controls or preparation are invalid.");
 
     const mfem::Vector acceptedPhysical(state.AcceptedPhysicalState()),
@@ -233,7 +252,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
     verify("context_reference");
 
     auto summary = detail::Open(output / "preconditioner_summary.csv");
-    summary << "role,variant,restart,maximum_iterations,relative_tolerance,rhs_norm,absolute_target,"
+    summary << "role,variant,schur_correction_sweeps,schur_relaxation,restart,maximum_iterations,relative_tolerance,rhs_norm,absolute_target,"
                "reported_initial_norm,reported_final_norm,reported_iterations,iterations,solver_converged,"
                "true_residual_norm,true_relative,true_gate_pass,normalized_correction_norm,physical_correction_norm,"
                "construction_seconds,unchanged_state_refresh_seconds,refresh_did_work,normalization_binding_seconds,"
@@ -243,6 +262,8 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
                "preconditioner_applications,preconditioner_seconds,preconditioner_maximum_seconds,"
                "binding_setups,binding_refreshes,binding_setup_seconds,binding_refresh_seconds,"
                "true_check_jacobian_applications,krylov_other_seconds\n";
+    auto work = detail::Open(output / "preconditioner_work.csv");
+    work << "role,variant,stage,metric,count\n";
     auto actions = detail::Open(output / "preconditioner_actions.csv");
     actions << "role,variant,correction_norm,jacobian_action_norm,defect_norm,relative_defect,"
                "application_seconds,jacobian_seconds,context_action_relative_difference,context_equivalence_required\n";
@@ -254,12 +275,33 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
     auto history = detail::Open(output / "preconditioner_history.csv");
     history << "role,variant,iteration,reported_norm,final\n";
 
-    const auto run = [&]<pc::StellarStructureFactorizationPolicy Policy, bool MaterialLDU = false>(
+    const auto run = [&]<pc::StellarStructureFactorizationPolicy Policy, bool MaterialLDU = false,
+                         bool RefinedGravity = false>(
                          const std::string &role, const std::string &variant) {
         const auto candidateStart = detail::Clock::now();
-        auto physical = detail::Prepare<Policy, MaterialLDU>(state.Problem());
+        auto physical = detail::Prepare<Policy, MaterialLDU, RefinedGravity>(state.Problem(), options);
         static_assert(pc::PreparedStellarInverseFor<decltype(physical), decltype(state.Problem())>);
         const double constructionSeconds = detail::Seconds(candidateStart);
+        const auto recordWork = [&](const char *stage) {
+            const auto &structure = physical.GetGroupedPreconditioner().GetStructurePreconditioner();
+            const auto &stellar = structure.GetFactorization().GetStatistics();
+            const auto &gravity = structure.GetGravityPreconditioner().GetFactorization().GetStatistics();
+            const auto record = [&](const char *metric, std::uint64_t count) {
+                work << role << ',' << variant << ',' << stage << ',' << metric << ',' << count << '\n';
+            };
+            record("stellar_applications", stellar.applications);
+            record("material_inverse_applications", stellar.materialSurfaceInverseApplications);
+            record("gravity_inverse_applications", stellar.gravityInverseApplications);
+            record("material_to_gravity_applications", stellar.materialToGravityApplications);
+            record("gravity_to_material_applications", stellar.gravityToMaterialApplications);
+            record("gravity_mass_inverse_applications", gravity.massInverseApplications);
+            record("gravity_potential_schur_applications", gravity.potentialSchurApplications);
+            record("gravity_divergence_applications", gravity.divergenceApplications);
+            record("gravity_transpose_divergence_applications", gravity.transposeDivergenceApplications);
+            record("gravity_schur_refinement_sweeps", gravity.schurRefinementSweeps);
+            work.flush();
+        };
+        recordWork("constructed");
         auto sectionStart = detail::Clock::now();
         const auto refresh = physical.Refresh();
         const double refreshSeconds = detail::Seconds(sectionStart);
@@ -280,7 +322,8 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         referenceDifference -= referenceCorrection;
         const double referenceRelative = detail::Relative(referenceDifference.Norml2(),
             std::max(inverseAction.Norml2(), referenceCorrection.Norml2()));
-        constexpr bool requireEquivalent = std::same_as<Policy, pc::IndependentStellarSubsystems> && !MaterialLDU;
+        const bool requireEquivalent = std::same_as<Policy, pc::IndependentStellarSubsystems> && !MaterialLDU &&
+            (!RefinedGravity || options.schurCorrectionSweeps == 0);
         actions << role << ',' << variant << ',' << inverseAction.Norml2() << ',' << actionNorm << ','
                 << inverseJacobianAction.Norml2() << ',' << inverseJacobianAction.Norml2() / rhsNorm << ','
                 << applicationSeconds << ',' << actionSeconds << ',' << referenceRelative << ','
@@ -308,11 +351,13 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         mfem::Vector solution(rhs.Size()), defect(rhs.Size());
         solution = 0.0;
         const double krylovSetupSeconds = detail::Seconds(sectionStart);
+        recordWork("before_krylov");
         countedJacobian.ResetStatistics();
         countedInverse.ResetStatistics();
         sectionStart = detail::Clock::now();
         krylov.Mult(rhs, solution);
         const double krylovSeconds = detail::Seconds(sectionStart);
+        recordWork("after_krylov");
         const auto jacobianStats = countedJacobian.GetStatistics();
         const auto inverseStats = countedInverse.GetStatistics();
         const auto lifecycle = countedInverse.GetLifecycleStatistics();
@@ -331,7 +376,8 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         normalization.DenormalizeState(solution, physicalSolution);
         normalization.DenormalizeResidual(defect, physicalDefect);
         const double setupAndSolve = constructionSeconds + bindingSeconds + krylovSetupSeconds + krylovSeconds;
-        summary << role << ',' << variant << ',' << options.restart << ',' << options.maximumIterations << ','
+        summary << role << ',' << variant << ',' << (RefinedGravity ? options.schurCorrectionSweeps : 0) << ','
+                << (RefinedGravity ? options.schurRelaxation : 0.0) << ',' << options.restart << ',' << options.maximumIterations << ','
                 << options.relativeTolerance << ',' << rhsNorm << ',' << threshold << ','
                 << krylov.GetInitialNorm() << ',' << krylov.GetFinalNorm() << ',' << krylov.GetNumIterations() << ','
                 << inverseStats.applications << ',' << krylov.GetConverged() << ',' << defectNorm << ','
@@ -387,6 +433,10 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         run.template operator()<pc::ApproximateStellarBlockLDU>("candidate", options.variant); break;
     case PreconditionerBenchmarkVariant::material_ldu:
         run.template operator()<pc::IndependentStellarSubsystems, true>("candidate", options.variant); break;
+    case PreconditionerBenchmarkVariant::gravity_schur:
+        run.template operator()<pc::IndependentStellarSubsystems, false, true>("candidate", options.variant); break;
+    case PreconditionerBenchmarkVariant::stellar_ldu_gravity_schur:
+        run.template operator()<pc::ApproximateStellarBlockLDU, false, true>("candidate", options.variant); break;
     }
     auto scope = detail::Open(output / "preconditioner_scope.txt");
     scope << "Both solves use fresh physical inverses, separately completed scalar borders, the same physical J,\n"

experiments/preconditioner_schur_cases_2026-09-18.json

@@ -0,0 +1,56 @@
+[
+  {
+    "id": "n1p5_gravity_schur_k1_w0p5",
+    "arguments": [
+      "--index",
+      "1.5",
+      "--angular-momentum",
+      "0",
+      "--no-central",
+      "--max-newton",
+      "0",
+      "--restart",
+      "160",
+      "--max-linear",
+      "200",
+      "--linear-tol",
+      "0.03",
+      "--probe-state",
+      "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/numerical_strategies_2026-09-17/models_v2/n1p5_J0_curvature/accepted_state.txt",
+      "--preconditioner-probe",
+      "--preconditioner-variant",
+      "gravity-schur",
+      "--schur-sweeps",
+      "1",
+      "--schur-relaxation",
+      "0.5"
+    ]
+  },
+  {
+    "id": "n1p5_stellar_ldu_gravity_schur_k1_w0p5",
+    "arguments": [
+      "--index",
+      "1.5",
+      "--angular-momentum",
+      "0",
+      "--no-central",
+      "--max-newton",
+      "0",
+      "--restart",
+      "160",
+      "--max-linear",
+      "200",
+      "--linear-tol",
+      "0.03",
+      "--probe-state",
+      "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/numerical_strategies_2026-09-17/models_v2/n1p5_J0_curvature/accepted_state.txt",
+      "--preconditioner-probe",
+      "--preconditioner-variant",
+      "stellar-ldu-gravity-schur",
+      "--schur-sweeps",
+      "1",
+      "--schur-relaxation",
+      "0.5"
+    ]
+  }
+]