Add matched frozen-state preconditioner benchmarks

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 <class Context> 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<double>(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");

experiments/plateau_preconditioner_benchmark.hpp

@@ -0,0 +1,407 @@
+#pragma once
+
+#include <algorithm>
+#include <chrono>
+#include <cmath>
+#include <concepts>
+#include <cstring>
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <stdexcept>
+#include <string>
+#include <type_traits>
+#include <utility>
+
+// 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<double>(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<std::size_t>(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<long double>(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<double>::infinity());
+}
+
+template <pc::StellarStructureFactorizationPolicy StructurePolicy, bool MaterialLDU = false, class Problem>
+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 <class State>
+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 = [&]<pc::StellarStructureFactorizationPolicy Policy, bool MaterialLDU = false>(
+                         const std::string &role, const std::string &variant) {
+        const auto candidateStart = detail::Clock::now();
+        auto physical = detail::Prepare<Policy, MaterialLDU>(state.Problem());
+        static_assert(pc::PreparedStellarInverseFor<decltype(physical), decltype(state.Problem())>);
+        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<Policy, pc::IndependentStellarSubsystems> && !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<double>::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()<pc::IndependentStellarSubsystems>("baseline", "independent");
+    switch (selected) {
+    case PreconditionerBenchmarkVariant::independent:
+        run.template operator()<pc::IndependentStellarSubsystems>("candidate", options.variant); break;
+    case PreconditionerBenchmarkVariant::material_then_gravity:
+        run.template operator()<pc::MaterialThenGravityTriangular>("candidate", options.variant); break;
+    case PreconditionerBenchmarkVariant::gravity_then_material:
+        run.template operator()<pc::GravityThenMaterialTriangular>("candidate", options.variant); break;
+    case PreconditionerBenchmarkVariant::stellar_ldu:
+        run.template operator()<pc::ApproximateStellarBlockLDU>("candidate", options.variant); break;
+    case PreconditionerBenchmarkVariant::material_ldu:
+        run.template operator()<pc::IndependentStellarSubsystems, 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"
+          << "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