From bc3c39dcfdc80a94960ddf0861296fe3ed0d1526 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 16:23:13 -0400 Subject: [PATCH] Benchmark isolated and coupled gravity Schur refinements --- experiments/plateau_diagnostic.cpp | 10 ++- .../plateau_preconditioner_benchmark.hpp | 70 ++++++++++++++++--- ...preconditioner_schur_cases_2026-09-18.json | 56 +++++++++++++++ 3 files changed, 125 insertions(+), 11 deletions(-) create mode 100644 experiments/preconditioner_schur_cases_2026-09-18.json diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index ca82933..37415d3 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/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 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 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. diff --git a/experiments/plateau_preconditioner_benchmark.hpp b/experiments/plateau_preconditioner_benchmark.hpp index cc58f17..7f0bbc8 100644 --- a/experiments/plateau_preconditioner_benchmark.hpp +++ b/experiments/plateau_preconditioner_benchmark.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -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::infinity()); } -template -auto Prepare(const Problem &problem) { +template +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 = [&]( + const auto run = [&]( const std::string &role, const std::string &variant) { const auto candidateStart = detail::Clock::now(); - auto physical = detail::Prepare(state.Problem()); + auto physical = detail::Prepare(state.Problem(), options); static_assert(pc::PreparedStellarInverseFor); 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 && !MaterialLDU; + const bool requireEquivalent = std::same_as && !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()("candidate", options.variant); break; case PreconditionerBenchmarkVariant::material_ldu: run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::gravity_schur: + run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::stellar_ldu_gravity_schur: + run.template operator()("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" diff --git a/experiments/preconditioner_schur_cases_2026-09-18.json b/experiments/preconditioner_schur_cases_2026-09-18.json new file mode 100644 index 0000000..264f95f --- /dev/null +++ b/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" + ] + } +]