From e557e6b481a49609d4f4da816f3bb4c6db5a7ce1 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 16:16:46 -0400 Subject: [PATCH] Add matched frozen-state preconditioner benchmarks --- experiments/plateau_diagnostic.cpp | 47 +- .../plateau_preconditioner_benchmark.hpp | 407 ++++++++++++++++++ 2 files changed, 452 insertions(+), 2 deletions(-) create mode 100644 experiments/plateau_preconditioner_benchmark.hpp diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 3ea38a4..ca82933 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -18,6 +18,7 @@ import mean_field; #include "plateau_enthalpy_step_diagnostics.hpp" #include "plateau_physical_diagnostics.hpp" #include "plateau_preconditioner_diagnostics.hpp" +#include "plateau_preconditioner_benchmark.hpp" #include "plateau_pressure_obstruction_diagnostics.hpp" #include "plateau_shape_curvature_diagnostics.hpp" #include "plateau_numerical_strategy_diagnostics.hpp" @@ -33,6 +34,8 @@ struct Options { int maxNewton = 12, maxLinear = 200, restart = 40, probeRestart = 40; bool probes = false, noCentral = false; bool strategyProbe = false, positivityProbe = false, positivityCertificate = false; + bool preconditionerProbe = false; + std::string preconditionerVariant = "independent"; std::string strategy = "baseline"; double gravityTolerance = 1e-10; int gravityIterations = 300, gravityRestart = 80; @@ -92,6 +95,10 @@ Options Parse(int argc, char **argv) { o.positivityProbe = true; else if (a == "--positivity-certificate") o.positivityCertificate = true; + else if (a == "--preconditioner-probe") + o.preconditionerProbe = true; + else if (a == "--preconditioner-variant") + o.preconditionerVariant = v(); else if (a == "--pressure-probe") o.pressureProbe = true; else if (a == "--shape-curvature") @@ -127,8 +134,17 @@ Options Parse(int argc, char **argv) { if (o.positivityCertificate && o.probeState.empty()) throw std::invalid_argument("--positivity-certificate requires --probe-state"); if (o.positivityCertificate && (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || !o.probeDirection.empty() || !o.replayCandidate.empty())) + o.strategyProbe || o.positivityProbe || o.preconditionerProbe || !o.probeDirection.empty() || !o.replayCandidate.empty())) throw std::invalid_argument("--positivity-certificate is a certificate-only mode; omit other probe options"); + if (o.preconditionerProbe) { + experiment::plateau_diagnostics::ParsePreconditionerBenchmarkVariant(o.preconditionerVariant); + if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || + o.strategyProbe || o.positivityProbe || o.positivityCertificate || !o.probeDirection.empty() || + !o.replayCandidate.empty() || !o.captureIterations.empty() || o.strategy != "baseline") + throw std::invalid_argument("--preconditioner-probe requires --max-newton 0 and no nonlinear/probe strategies"); + if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0) + throw std::invalid_argument("A preconditioner benchmark requires a finite positive --linear-tol"); + } 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" && @@ -397,14 +413,41 @@ template int RunContext(Context &context, const Options &o, cons metadata << "n=" << o.n << "\nK=" << scales.K << "\nrho_c=" << scales.rho << "\nJ=" << o.angularMomentum << "\nM=1\nR_reference=1\nxi1=" << scales.xi << "\nmu1=" << scales.mu << "\ncentral_phase=" << !o.noCentral << "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart - << "\nfrozen_probe_backend=" << (o.positivityCertificate ? "polynomial_certificate_no_linear_solve" : + << "\nfrozen_probe_backend=" << (o.preconditionerProbe ? "direct_mfem_fgmres_preconditioner_benchmark" : + o.positivityCertificate ? "polynomial_certificate_no_linear_solve" : 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 + << "\npreconditioner_probe=" << o.preconditionerProbe + << "\npreconditioner_variant=" << o.preconditionerVariant << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; metadata.flush(); + if (o.preconditionerProbe) { + metadata << "benchmark_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n'; + if (!o.probeState.empty()) metadata << "probe_state=" << o.probeState << '\n'; + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &) { + if (o.probeState.empty()) { + s.trialNormalizedState = s.acceptedNormalizedState; + } else { + const auto snapshot = LoadVector(o.probeState, s.AcceptedPhysicalState().Size()); + s.normalizedOperator->NormalizeState(snapshot, s.trialNormalizedState); + } + const auto preparationStart = Clock::now(); + if (!s.PrepareTrial()) + throw std::runtime_error("Preconditioner benchmark state is inadmissible to the physical operator"); + metadata << "benchmark_frozen_preparation_seconds=" + << std::chrono::duration(Clock::now() - preparationStart).count() << '\n'; + metadata.flush(); + experiment::plateau_diagnostics::BenchmarkPreconditioners(s, o.output, + {.variant = o.preconditionerVariant, .restart = o.restart, + .maximumIterations = o.maxLinear, .relativeTolerance = o.linearTolerance}); + }); + // A missed linear gate is a measured outcome, not failure to execute + // the experiment. The CSV retains both recurrence and true gate results. + return 0; + } if (!o.probeState.empty()) { if (o.maxNewton != 0) throw std::invalid_argument("--probe-state requires --max-newton 0"); diff --git a/experiments/plateau_preconditioner_benchmark.hpp b/experiments/plateau_preconditioner_benchmark.hpp new file mode 100644 index 0000000..cc58f17 --- /dev/null +++ b/experiments/plateau_preconditioner_benchmark.hpp @@ -0,0 +1,407 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Include after import mean_field. Numerical recipes remain typed; the runtime +// selection below belongs only to this single-rank experimental driver. +namespace experiment::plateau_diagnostics { +enum class PreconditionerBenchmarkVariant { + independent, + material_then_gravity, + gravity_then_material, + stellar_ldu, + material_ldu +}; + +inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const std::string &name) { + if (name == "independent") return PreconditionerBenchmarkVariant::independent; + if (name == "material-then-gravity") return PreconditionerBenchmarkVariant::material_then_gravity; + 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; + throw std::invalid_argument("Unknown preconditioner benchmark variant: " + name); +} + +struct PreconditionerBenchmarkOptions final { + std::string variant{"independent"}; + int restart{160}; + int maximumIterations{200}; + double relativeTolerance{0.03}; +}; + +namespace preconditioner_benchmark_detail { +namespace pc = mean_field::preconditioning; +using Clock = std::chrono::steady_clock; + +inline double Seconds(Clock::time_point start) { + return std::chrono::duration(Clock::now() - start).count(); +} + +inline std::ofstream Open(const std::filesystem::path &path) { + std::ofstream stream(path); + stream.exceptions(std::ios::badbit | std::ios::failbit); + stream << std::setprecision(17); + return stream; +} + +inline void Save(const std::filesystem::path &path, const mfem::Vector &values) { + auto stream = Open(path); + stream << values.Size() << '\n'; + values.Print(stream, 1); +} + +inline bool Identical(const mfem::Vector &left, const mfem::Vector &right) { + return left.Size() == right.Size() && (left.Size() == 0 || + std::memcmp(left.GetData(), right.GetData(), + static_cast(left.Size()) * sizeof(mfem::real_t)) == 0); +} + +inline bool Finite(const mfem::Vector &values) { + for (int i = 0; i < values.Size(); ++i) + if (!std::isfinite(values(i))) return false; + return true; +} + +inline double BlockNorm(const mfem::Vector &values, int offset, int size) { + long double sum = 0.0L; + for (int i = offset; i < offset + size; ++i) + sum += static_cast(values(i)) * values(i); + return std::sqrt(sum); +} + +inline double Relative(double numerator, double denominator) { + return denominator > 0.0 ? numerator / denominator : + (numerator == 0.0 ? 0.0 : std::numeric_limits::infinity()); +} + +template +auto Prepare(const Problem &problem) { + const auto defaults = pc::stellarStructureBlock(problem); + const auto material = [&]() { + if constexpr (MaterialLDU) { + const auto &original = defaults.materialSurfaceComponent(); + return pc::materialSurfaceBlock(problem, original.materialBackend(), original.surfaceBackend(), + pc::ApproximateMaterialSurfaceLDU{}, original.surfaceSurrogate(), original.diagonalOptions()); + } else { + return defaults.materialSurfaceComponent(); + } + }(); + auto structure = pc::stellarStructureBlock(problem, material, defaults.gravityComponent(), 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( + problem, std::move(structure), pc::backend::DenseDirect{})); +} + +// Instrumenting J changes its C++ address, whereas the normalization wrapper +// deliberately permits only its original Jacobian identity. Preserve that +// binding contract while allowing MFEM to apply the counted forwarding J. +class InstrumentedBinding final : public mfem::Solver { +public: + InstrumentedBinding(mfem::Solver &inverse, const mfem::Operator &original, + const mfem::Operator &instrumented) + : mfem::Solver(inverse.Height(), inverse.Width()), m_inverse(inverse), + m_original(original), m_instrumented(instrumented) {} + + void SetOperator(const mfem::Operator &operation) override { + if (&operation != &m_instrumented) + throw std::invalid_argument("Benchmark preconditioner received a different counted Jacobian."); + m_inverse.SetOperator(m_original); + } + + void Mult(const mfem::Vector &input, mfem::Vector &output) const override { + m_inverse.Mult(input, output); + } + +private: + mfem::Solver &m_inverse; + const mfem::Operator &m_original; + const mfem::Operator &m_instrumented; +}; +} // namespace preconditioner_benchmark_detail + +template +void BenchmarkPreconditioners(State &state, const std::filesystem::path &output, + const PreconditionerBenchmarkOptions &options) { + namespace detail = preconditioner_benchmark_detail; + namespace pc = mean_field::preconditioning; + namespace solver = mean_field::solver; + const auto selected = ParsePreconditionerBenchmarkVariant(options.variant); + int ranks = 0; + if (MPI_Comm_size(state.Problem().GetCommunicator(), &ranks) != MPI_SUCCESS || ranks != 1) + 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()) + throw std::invalid_argument("Frozen preconditioner benchmark controls or preparation are invalid."); + + const mfem::Vector acceptedPhysical(state.AcceptedPhysicalState()), + acceptedNormalized(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual), + correction(state.normalizedCorrection), frozenState(state.trialNormalizedState), + residual(state.trialNormalizedResidual), candidatePhysical(state.candidatePhysicalState); + const auto generation = state.Problem().GetPreparationGeneration(); + const auto &normalization = state.normalizedOperator->GetNormalization(); + const mfem::Vector stateFactors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors()); + mfem::Vector rhs(residual); + rhs *= -1.0; + const double rhsNorm = rhs.Norml2(); + const double threshold = options.relativeTolerance * rhsNorm; + if (!detail::Finite(rhs) || !std::isfinite(rhsNorm) || rhsNorm <= 0.0 || !std::isfinite(threshold)) + throw std::invalid_argument("A benchmark needs a finite, nonzero frozen Newton right-hand side."); + mfem::Vector physicalBase, physicalRhs; + state.normalizedOperator->DenormalizeState(frozenState, physicalBase); + normalization.DenormalizeResidual(rhs, physicalRhs); + detail::Save(output / "benchmark_physical_base.txt", physicalBase); + detail::Save(output / "benchmark_normalized_base.txt", frozenState); + detail::Save(output / "benchmark_normalized_rhs.txt", rhs); + detail::Save(output / "benchmark_state_factors.txt", stateFactors); + detail::Save(output / "benchmark_residual_factors.txt", residualFactors); + auto layout = detail::Open(output / "layout.csv"); + layout << "kind,block,offset,size,scale\n"; + for (const auto &block : state.Problem().GetManifest().residualBlocks()) + layout << "residual," << block.stableId << ',' << block.offset << ',' << block.size << ',' << block.scale << '\n'; + for (const auto &block : state.Problem().GetManifest().valueBlocks()) + layout << "value," << block.stableId << ',' << block.offset << ',' << block.size << ',' << block.scale << '\n'; + + mfem::Vector direction(rhs.Size()), initialAction(rhs.Size()); + for (int i = 0; i < direction.Size(); ++i) + direction(i) = std::sin(0.019 * (i + 1)) + 0.25 * std::cos(0.071 * (i + 1)); + direction /= direction.Norml2(); + state.normalizedOperator->Mult(direction, initialAction); + if (!detail::Finite(initialAction)) + throw std::runtime_error("The frozen operator invariant probe is nonfinite."); + detail::Save(output / "benchmark_invariant_direction.txt", direction); + detail::Save(output / "benchmark_invariant_action.txt", initialAction); + + auto invariants = detail::Open(output / "preconditioner_invariants.csv"); + invariants << "role,metric,value\n"; + const auto verify = [&](const std::string &role) { + bool valid = true; + const auto check = [&](const char *name, bool result) { + invariants << role << ',' << name << ',' << result << '\n'; + valid = valid && result; + }; + check("accepted_physical_bitwise_identical", detail::Identical(acceptedPhysical, state.AcceptedPhysicalState())); + check("accepted_normalized_bitwise_identical", detail::Identical(acceptedNormalized, state.acceptedNormalizedState)); + check("accepted_residual_bitwise_identical", detail::Identical(acceptedResidual, state.acceptedNormalizedResidual)); + check("correction_bitwise_identical", detail::Identical(correction, state.normalizedCorrection)); + check("trial_state_bitwise_identical", detail::Identical(frozenState, state.trialNormalizedState)); + check("trial_residual_bitwise_identical", detail::Identical(residual, state.trialNormalizedResidual)); + check("candidate_physical_bitwise_identical", detail::Identical(candidatePhysical, state.candidatePhysicalState)); + check("state_factors_bitwise_identical", detail::Identical(stateFactors, normalization.StateFactors())); + check("residual_factors_bitwise_identical", detail::Identical(residualFactors, normalization.ResidualFactors())); + check("preparation_generation_unchanged", generation == state.Problem().GetPreparationGeneration()); + check("normalized_operator_current", state.normalizedOperator->IsPrepared()); + mfem::Vector action(rhs.Size()); + state.normalizedOperator->Mult(direction, action); + check("jacobian_action_bitwise_identical", detail::Identical(initialAction, action)); + action -= initialAction; + invariants << role << ",jacobian_action_difference_norm," << action.Norml2() << '\n'; + invariants.flush(); + if (!valid) + throw std::runtime_error("Frozen preconditioner benchmark changed a context buffer, normalization, or operator."); + }; + + auto reference = detail::Open(output / "preconditioner_reference.csv"); + reference << "metric,value\n"; + auto start = detail::Clock::now(); + const auto referenceRefresh = state.physicalInverse->Refresh(); + reference << "context_inverse_refresh_seconds," << detail::Seconds(start) + << "\ncontext_inverse_refresh_did_work," << referenceRefresh.DidAnyWork() << '\n'; + mfem::Vector referenceCorrection(rhs.Size()); + start = detail::Clock::now(); + state.normalizedInverse->Mult(rhs, referenceCorrection); + reference << "context_inverse_application_seconds," << detail::Seconds(start) + << "\ncontext_inverse_correction_norm," << referenceCorrection.Norml2() << '\n'; + reference.flush(); + if (!detail::Finite(referenceCorrection)) + throw std::runtime_error("The context reference inverse action is nonfinite."); + detail::Save(output / "context_reference_normalized_correction.txt", referenceCorrection); + verify("context_reference"); + + auto summary = detail::Open(output / "preconditioner_summary.csv"); + summary << "role,variant,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," + "krylov_setup_seconds,krylov_seconds,true_check_seconds,setup_plus_solve_seconds," + "setup_solve_and_true_check_seconds,measured_candidate_seconds," + "jacobian_applications,jacobian_seconds,jacobian_maximum_seconds," + "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 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"; + auto blocks = detail::Open(output / "preconditioner_blocks.csv"); + blocks << "role,variant,block,normalized_rhs_norm,normalized_defect_norm,physical_rhs_norm," + "physical_defect_norm,block_relative,global_relative_contribution,squared_defect_fraction\n"; + auto correctionBlocks = detail::Open(output / "preconditioner_correction_blocks.csv"); + correctionBlocks << "role,variant,block,normalized_norm,physical_norm\n"; + auto history = detail::Open(output / "preconditioner_history.csv"); + history << "role,variant,iteration,reported_norm,final\n"; + + const auto run = [&]( + const std::string &role, const std::string &variant) { + const auto candidateStart = detail::Clock::now(); + auto physical = detail::Prepare(state.Problem()); + static_assert(pc::PreparedStellarInverseFor); + const double constructionSeconds = detail::Seconds(candidateStart); + auto sectionStart = detail::Clock::now(); + const auto refresh = physical.Refresh(); + const double refreshSeconds = detail::Seconds(sectionStart); + sectionStart = detail::Clock::now(); + auto normalized = state.normalizedOperator->MakeScaledPreconditioner(physical); + const double bindingSeconds = detail::Seconds(sectionStart); + + mfem::Vector inverseAction(rhs.Size()), inverseJacobianAction(rhs.Size()); + sectionStart = detail::Clock::now(); + normalized.Mult(rhs, inverseAction); + const double applicationSeconds = detail::Seconds(sectionStart); + sectionStart = detail::Clock::now(); + state.normalizedOperator->Mult(inverseAction, inverseJacobianAction); + const double actionSeconds = detail::Seconds(sectionStart); + const double actionNorm = inverseJacobianAction.Norml2(); + inverseJacobianAction -= rhs; + mfem::Vector referenceDifference(inverseAction); + referenceDifference -= referenceCorrection; + const double referenceRelative = detail::Relative(referenceDifference.Norml2(), + std::max(inverseAction.Norml2(), referenceCorrection.Norml2())); + constexpr bool requireEquivalent = std::same_as && !MaterialLDU; + actions << role << ',' << variant << ',' << inverseAction.Norml2() << ',' << actionNorm << ',' + << inverseJacobianAction.Norml2() << ',' << inverseJacobianAction.Norml2() / rhsNorm << ',' + << applicationSeconds << ',' << actionSeconds << ',' << referenceRelative << ',' + << requireEquivalent << '\n'; + actions.flush(); + if (requireEquivalent && (!std::isfinite(referenceRelative) || + referenceRelative > 256.0 * std::numeric_limits::epsilon())) + throw std::runtime_error("Fresh independent inverse differs from the context's default inverse."); + + sectionStart = detail::Clock::now(); + solver::InstrumentedOperator countedJacobian(*state.normalizedOperator); + detail::InstrumentedBinding binding(normalized, *state.normalizedOperator, countedJacobian); + solver::InstrumentedPreconditioner countedInverse(binding); + solver::ResidualHistoryMonitor monitor; + mfem::FGMRESSolver krylov(state.Problem().GetCommunicator()); + krylov.SetPreconditioner(countedInverse); + krylov.SetOperator(countedJacobian); + krylov.SetMonitor(monitor); + krylov.SetKDim(options.restart); + krylov.SetRelTol(0.0); + krylov.SetAbsTol(threshold); + krylov.SetMaxIter(options.maximumIterations); + krylov.SetPrintLevel(-1); + krylov.iterative_mode = false; + mfem::Vector solution(rhs.Size()), defect(rhs.Size()); + solution = 0.0; + const double krylovSetupSeconds = detail::Seconds(sectionStart); + countedJacobian.ResetStatistics(); + countedInverse.ResetStatistics(); + sectionStart = detail::Clock::now(); + krylov.Mult(rhs, solution); + const double krylovSeconds = detail::Seconds(sectionStart); + const auto jacobianStats = countedJacobian.GetStatistics(); + const auto inverseStats = countedInverse.GetStatistics(); + const auto lifecycle = countedInverse.GetLifecycleStatistics(); + + // The recurrence's convergence flag is reported, but only this fresh + // original-operator application determines the linear accuracy gate. + sectionStart = detail::Clock::now(); + state.normalizedOperator->Mult(solution, defect); + defect -= rhs; + const double defectNorm = defect.Norml2(); + const bool finite = detail::Finite(solution) && detail::Finite(defect) && std::isfinite(defectNorm); + const bool gate = finite && defectNorm <= threshold; + const double trueCheckSeconds = detail::Seconds(sectionStart); + const double candidateSeconds = detail::Seconds(candidateStart); + mfem::Vector physicalSolution, physicalDefect; + normalization.DenormalizeState(solution, physicalSolution); + normalization.DenormalizeResidual(defect, physicalDefect); + const double setupAndSolve = constructionSeconds + bindingSeconds + krylovSetupSeconds + krylovSeconds; + summary << role << ',' << variant << ',' << options.restart << ',' << options.maximumIterations << ',' + << options.relativeTolerance << ',' << rhsNorm << ',' << threshold << ',' + << krylov.GetInitialNorm() << ',' << krylov.GetFinalNorm() << ',' << krylov.GetNumIterations() << ',' + << inverseStats.applications << ',' << krylov.GetConverged() << ',' << defectNorm << ',' + << defectNorm / rhsNorm << ',' << gate << ',' << solution.Norml2() << ',' << physicalSolution.Norml2() << ',' + << constructionSeconds << ',' << refreshSeconds << ',' << refresh.DidAnyWork() << ',' << bindingSeconds << ',' + << krylovSetupSeconds << ',' << krylovSeconds << ',' << trueCheckSeconds << ',' << setupAndSolve << ',' + << setupAndSolve + trueCheckSeconds << ',' << candidateSeconds << ',' + << jacobianStats.applications << ',' << jacobianStats.totalSeconds << ',' << jacobianStats.maximumSeconds << ',' + << inverseStats.applications << ',' << inverseStats.totalSeconds << ',' << inverseStats.maximumSeconds << ',' + << lifecycle.setups << ',' << lifecycle.refreshes << ',' << lifecycle.setupSeconds << ',' + << lifecycle.refreshSeconds << ",1," << + std::max(0.0, krylovSeconds - jacobianStats.totalSeconds - inverseStats.totalSeconds) << '\n'; + summary.flush(); + for (const auto &block : state.Problem().GetManifest().residualBlocks()) { + const double b = detail::BlockNorm(rhs, block.offset, block.size); + const double r = detail::BlockNorm(defect, block.offset, block.size); + blocks << role << ',' << variant << ',' << block.stableId << ',' << b << ',' << r << ',' + << detail::BlockNorm(physicalRhs, block.offset, block.size) << ',' + << detail::BlockNorm(physicalDefect, block.offset, block.size) << ',' + << detail::Relative(r, b) << ',' << r / rhsNorm << ',' + << (defectNorm > 0.0 ? std::pow(r / defectNorm, 2) : 0.0) << '\n'; + } + for (const auto &block : state.Problem().GetManifest().valueBlocks()) + correctionBlocks << role << ',' << variant << ',' << block.stableId << ',' + << detail::BlockNorm(solution, block.offset, block.size) << ',' + << detail::BlockNorm(physicalSolution, block.offset, block.size) << '\n'; + for (const auto &sample : monitor.GetHistory()) + history << role << ',' << variant << ',' << sample.iteration << ',' << sample.reportedNorm << ',' + << sample.final << '\n'; + blocks.flush(); + correctionBlocks.flush(); + history.flush(); + detail::Save(output / (role + "_normalized_correction.txt"), solution); + detail::Save(output / (role + "_physical_correction.txt"), physicalSolution); + detail::Save(output / (role + "_normalized_linear_defect.txt"), defect); + detail::Save(output / (role + "_physical_linear_defect.txt"), physicalDefect); + verify(role); + std::cout << "Preconditioner " << role << '/' << variant << " iterations=" << inverseStats.applications + << " true_relative=" << defectNorm / rhsNorm << " true_gate=" << gate + << " krylov_seconds=" << krylovSeconds + << " setup_solve_and_true_check_seconds=" << setupAndSolve + trueCheckSeconds << std::endl; + }; + + run.template operator()("baseline", "independent"); + switch (selected) { + case PreconditionerBenchmarkVariant::independent: + run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::material_then_gravity: + run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::gravity_then_material: + run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::stellar_ldu: + run.template operator()("candidate", options.variant); break; + case PreconditionerBenchmarkVariant::material_ldu: + 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" + << "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" + << "setup_plus_solve_seconds includes construction, normalization binding, Krylov setup, and Krylov solve.\n" + << "setup_solve_and_true_check_seconds adds the separately measured original-J true-residual gate.\n" + << "measured_candidate_seconds also includes diagnostic refresh and the inverse/Jv comparison.\n" + << "Common FEM/context construction and frozen-state preparation occur before this helper.\n" + << "The helper performs no trial preparation or nonlinear update; WithState restores its accepted seed afterward.\n" + << "Block norms are coefficient norms in the stated normalized/physical coordinates, not spatial error norms.\n" + << "A linear gate pass does not establish nonlinear acceptance, EOS admissibility, or virial qualification.\n"; +} +} // namespace experiment::plateau_diagnostics