From 370dbbf73bd976f11cb80893ec8d6a4842e1faa3 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 16:37:32 -0400 Subject: [PATCH] Compare cached coupling and checked coarse preconditioner variants --- CMakeLists.txt | 1 + .../plateau_preconditioner_benchmark.hpp | 89 +++++++++++++-- .../preconditioner_candidates_2026-09-18.json | 102 ++++++++++++++++++ 3 files changed, 185 insertions(+), 7 deletions(-) create mode 100644 experiments/preconditioner_candidates_2026-09-18.json diff --git a/CMakeLists.txt b/CMakeLists.txt index 776bc21..42f85da 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -372,6 +372,7 @@ add_executable(tests tests/solver/stellar_equilibrium_architecture_internal.cpp tests/solver/stellar_equilibrium_runtime.cpp tests/solver/enthalpy_admissibility.cpp + tests/solver/plateau_coupled_coarse.cpp tests/solver/gravity_subsystem.cpp tests/solver/gravity_compact_applications.cpp tests/solver/gravity_reduction.cpp diff --git a/experiments/plateau_preconditioner_benchmark.hpp b/experiments/plateau_preconditioner_benchmark.hpp index 7f0bbc8..b8125fa 100644 --- a/experiments/plateau_preconditioner_benchmark.hpp +++ b/experiments/plateau_preconditioner_benchmark.hpp @@ -11,6 +11,9 @@ #include #include #include +#include + +#include "plateau_coupled_coarse.hpp" #include #include #include @@ -26,7 +29,9 @@ enum class PreconditionerBenchmarkVariant { stellar_ldu, material_ldu, gravity_schur, - stellar_ldu_gravity_schur + stellar_ldu_gravity_schur, + coupled_shape_coarse, + gravity_schur_coarse }; inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const std::string &name) { @@ -37,6 +42,8 @@ inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const 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; + if (name == "coupled-shape-coarse") return PreconditionerBenchmarkVariant::coupled_shape_coarse; + if (name == "gravity-schur-coarse") return PreconditionerBenchmarkVariant::gravity_schur_coarse; throw std::invalid_argument("Unknown preconditioner benchmark variant: " + name); } @@ -262,6 +269,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 coarseReport = detail::Open(output / "preconditioner_coarse.csv"); + coarseReport << "role,variant,stage,metric,value\n"; auto work = detail::Open(output / "preconditioner_work.csv"); work << "role,variant,stage,metric,count\n"; auto actions = detail::Open(output / "preconditioner_actions.csv"); @@ -276,12 +285,12 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, history << "role,variant,iteration,reported_norm,final\n"; const auto run = [&]( + bool RefinedGravity = false, bool UseCoarse = false>( const std::string &role, const std::string &variant) { const auto candidateStart = detail::Clock::now(); auto physical = detail::Prepare(state.Problem(), options); static_assert(pc::PreparedStellarInverseFor); - const double constructionSeconds = detail::Seconds(candidateStart); + double constructionSeconds = detail::Seconds(candidateStart); const auto recordWork = [&](const char *stage) { const auto &structure = physical.GetGroupedPreconditioner().GetStructurePreconditioner(); const auto &stellar = structure.GetFactorization().GetStatistics(); @@ -308,11 +317,70 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, sectionStart = detail::Clock::now(); auto normalized = state.normalizedOperator->MakeScaledPreconditioner(physical); const double bindingSeconds = detail::Seconds(sectionStart); + using NormalizedOperation = std::remove_reference_t; + using Coarse = PreparedCoupledShapeCoarse; + std::unique_ptr coarse; + mfem::Solver *applicationInverse = &normalized; + const auto recordCoarse = [&](const char *stage) { + if (!coarse) return; + const auto &stats = coarse->GetStatistics(); + const auto metric = [&](const char *name, auto value) { + coarseReport << role << ',' << variant << ',' << stage << ',' << name << ',' << value << '\n'; + }; + metric("input_modes", stats.inputModes); + metric("independent_modes", stats.independentModes); + metric("active_modes", stats.activeModes); + metric("fallback", static_cast(stats.fallback)); + metric("minimum_overlap_lower_bound", stats.minimumOverlapLowerBound); + metric("reciprocal_overlap_condition", stats.reciprocalOverlapCondition); + metric("setup_jacobian_applications", stats.setupOperatorApplications); + metric("setup_inverse_applications", stats.setupInverseApplications); + metric("setup_seconds", stats.lastSetupSeconds); + metric("persistent_vector_bytes", stats.persistentVectorBytes); + metric("applications", stats.applications); + metric("inverse_applications", stats.inverseApplications); + metric("failed_applications", stats.failedApplications); + coarseReport.flush(); + }; + if constexpr (UseCoarse) { + sectionStart = detail::Clock::now(); + const auto modes = MakeCoupledShapeModes(state.Problem(), physical); + const double basisSeconds = detail::Seconds(sectionStart); + coarse = std::make_unique(*state.normalizedOperator, normalized, modes, + state.Problem().GetCommunicator()); + const double completeCoarseSeconds = detail::Seconds(sectionStart); + coarseReport << role << ',' << variant << ",constructed,physical_inverse_construction_seconds," + << constructionSeconds << '\n'; + constructionSeconds += completeCoarseSeconds; + applicationInverse = coarse.get(); + coarseReport << role << ',' << variant << ",constructed,basis_construction_seconds," << basisSeconds << '\n' + << role << ',' << variant << ",constructed,complete_coarse_construction_seconds," + << completeCoarseSeconds << '\n'; + recordCoarse("constructed"); + // These extra checks validate the represented FEM coarse identity; + // they are diagnostics, separate from construction and solve timing. + double maximumRelativeError = 0.0; + for (const auto &source : coarse->GetSourceBasis()) { + mfem::Vector value(rhs.Size()), image(rhs.Size()), recovered(rhs.Size()); + normalized.Mult(source, value); + state.normalizedOperator->Mult(value, image); + coarse->Mult(image, recovered); + recovered -= value; + maximumRelativeError = std::max(maximumRelativeError, + detail::Relative(recovered.Norml2(), value.Norml2())); + } + coarseReport << role << ',' << variant << ",validation,maximum_exact_coarse_relative_error," + << maximumRelativeError << '\n'; + coarseReport.flush(); + if (!std::isfinite(maximumRelativeError) || maximumRelativeError > 1.e-8) + throw std::runtime_error("The actual FEM coarse inverse failed its exact-subspace check."); + } mfem::Vector inverseAction(rhs.Size()), inverseJacobianAction(rhs.Size()); sectionStart = detail::Clock::now(); - normalized.Mult(rhs, inverseAction); + applicationInverse->Mult(rhs, inverseAction); const double applicationSeconds = detail::Seconds(sectionStart); + detail::Save(output / (role + "_inverse_rhs_action.txt"), inverseAction); sectionStart = detail::Clock::now(); state.normalizedOperator->Mult(inverseAction, inverseJacobianAction); const double actionSeconds = detail::Seconds(sectionStart); @@ -322,7 +390,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, referenceDifference -= referenceCorrection; const double referenceRelative = detail::Relative(referenceDifference.Norml2(), std::max(inverseAction.Norml2(), referenceCorrection.Norml2())); - const bool requireEquivalent = std::same_as && !MaterialLDU && + const bool requireEquivalent = std::same_as && !MaterialLDU && !UseCoarse && (!RefinedGravity || options.schurCorrectionSweeps == 0); actions << role << ',' << variant << ',' << inverseAction.Norml2() << ',' << actionNorm << ',' << inverseJacobianAction.Norml2() << ',' << inverseJacobianAction.Norml2() / rhsNorm << ',' @@ -335,7 +403,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, sectionStart = detail::Clock::now(); solver::InstrumentedOperator countedJacobian(*state.normalizedOperator); - detail::InstrumentedBinding binding(normalized, *state.normalizedOperator, countedJacobian); + detail::InstrumentedBinding binding(*applicationInverse, *state.normalizedOperator, countedJacobian); solver::InstrumentedPreconditioner countedInverse(binding); solver::ResidualHistoryMonitor monitor; mfem::FGMRESSolver krylov(state.Problem().GetCommunicator()); @@ -352,12 +420,14 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, solution = 0.0; const double krylovSetupSeconds = detail::Seconds(sectionStart); recordWork("before_krylov"); + recordCoarse("before_krylov"); countedJacobian.ResetStatistics(); countedInverse.ResetStatistics(); sectionStart = detail::Clock::now(); krylov.Mult(rhs, solution); const double krylovSeconds = detail::Seconds(sectionStart); recordWork("after_krylov"); + recordCoarse("after_krylov"); const auto jacobianStats = countedJacobian.GetStatistics(); const auto inverseStats = countedInverse.GetStatistics(); const auto lifecycle = countedInverse.GetLifecycleStatistics(); @@ -437,12 +507,17 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, run.template operator()("candidate", options.variant); break; case PreconditionerBenchmarkVariant::stellar_ldu_gravity_schur: run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::coupled_shape_coarse: + run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::gravity_schur_coarse: + 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" << "the same normalization, the same saved RHS=-F, and zero initial corrections.\n" << "The fresh independent inverse must match the existing context inverse on that RHS.\n" - << "construction_seconds includes the complete physical preconditioner construction.\n" + << "construction_seconds includes physical inverse and any selected coarse/basis construction.\n" + << "Coarse exact-subspace validation is diagnostic work outside setup/solve totals.\n" << "unchanged_state_refresh_seconds measures an immediate Refresh; it is not a changed-state refresh benchmark.\n" << "Each candidate receives one diagnostic inverse/Jv application before the timed cold Krylov solve.\n" << "Application timers/counts in the summary include only Krylov work; the independent true check is separate.\n" diff --git a/experiments/preconditioner_candidates_2026-09-18.json b/experiments/preconditioner_candidates_2026-09-18.json new file mode 100644 index 0000000..9dd66f4 --- /dev/null +++ b/experiments/preconditioner_candidates_2026-09-18.json @@ -0,0 +1,102 @@ +[ + { + "id": "n1p5_cached_stellar_ldu", + "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" + ] + }, + { + "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" + ] + }, + { + "id": "n1p5_coupled_shape_coarse", + "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", + "coupled-shape-coarse" + ] + } +]