#include <algorithm>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <memory>
#include <stdexcept>
#include <string>
#include <cstdint>
#include <utility>
#include <vector>
#include <mfem.hpp>

import mean_field;

namespace {
    using namespace mean_field;
    using Clock = std::chrono::steady_clock;
    using Apply = std::function<void(int, mfem::Vector &)>;

    struct Options {
        std::string mesh{"sandbox.smesh"};
        std::filesystem::path output{"system_isolation_results"};
        std::string system{"hdiv-displacement"};
        std::string implementation{"both"};
        std::string check{"full"};
        int refinements{0}, samples{5}, iterations{3}, warmup{1}, directions{3};
        int prepareSamples{3}, prepareIterations{1};
        double deformation{0.7}, finiteDifferenceStep{1.0e-5}, polytropicIndex{1.0};
        bool help{false};
    };

    int Integer(const std::string &value) {
        std::size_t consumed = 0;
        const int result = std::stoi(value, &consumed);
        if (consumed != value.size()) throw std::invalid_argument("Invalid integer: " + value);
        return result;
    }

    double Real(const std::string &value) {
        std::size_t consumed = 0;
        const double result = std::stod(value, &consumed);
        if (consumed != value.size() || !std::isfinite(result))
            throw std::invalid_argument("Invalid finite number: " + value);
        return result;
    }

    bool IsLinearDerivativeSystem(const std::string &name);

    Options Parse(int argc, char **argv) {
        Options options;
        for (int index = 1; index < argc; ++index) {
            const std::string argument = argv[index];
            auto value = [&]() -> std::string {
                if (++index == argc) throw std::invalid_argument("Missing value after " + argument);
                return argv[index];
            };
            if (argument == "--mesh") options.mesh = value();
            else if (argument == "--output") options.output = value();
            else if (argument == "--system") options.system = value();
            else if (argument == "--implementation") options.implementation = value();
            else if (argument == "--check") options.check = value();
            else if (argument == "--refinements") options.refinements = Integer(value());
            else if (argument == "--samples") options.samples = Integer(value());
            else if (argument == "--iterations") options.iterations = Integer(value());
            else if (argument == "--warmup") options.warmup = Integer(value());
            else if (argument == "--directions") options.directions = Integer(value());
            else if (argument == "--prepare-samples") options.prepareSamples = Integer(value());
            else if (argument == "--prepare-iterations") options.prepareIterations = Integer(value());
            else if (argument == "--deformation") options.deformation = Real(value());
            else if (argument == "--fd-step") options.finiteDifferenceStep = Real(value());
            else if (argument == "--polytropic-index") options.polytropicIndex = Real(value());
            else if (argument == "--help") options.help = true;
            else throw std::invalid_argument("Unknown argument: " + argument);
        }
        if (options.system != "hdiv-displacement" && options.system != "hdiv-mass" &&
            options.system != "mass-normalization" && options.system != "angular-momentum" &&
            options.system != "hydrostatic-displacement" && options.system != "stellar-derivatives" &&
            options.system != "all" && options.system != "runtime-optimizations" &&
            options.system != "geometry-preparation" && options.system != "normalization-diagonals" &&
            options.system != "geometry-inspection" && options.system != "gravity-source-preparation" &&
            options.system != "final-optimizations" && options.system != "gravity-source-displacement" &&
            options.system != "hydrostatic-preparation" && options.system != "gravity-schur-preparation" &&
            !IsLinearDerivativeSystem(options.system))
            throw std::invalid_argument("Unknown system; use --help for the available systems.");
        if (options.implementation != "both" && options.implementation != "optimized" && options.implementation != "reference")
            throw std::invalid_argument("Implementation must be both, optimized, or reference.");
        if (options.check != "full" && options.check != "reference" && options.check != "none")
            throw std::invalid_argument("Check must be full, reference, or none.");
        if (options.refinements < 0 || options.samples < 1 || options.iterations < 1 || options.warmup < 0 ||
            options.directions < 1 || options.prepareSamples < 0 || options.prepareIterations < 1 ||
            options.finiteDifferenceStep <= 0.0 || options.polytropicIndex <= 0.0)
            throw std::invalid_argument("Invalid sample count, refinement, or finite-difference step.");
        return options;
    }

    std::ofstream File(const std::filesystem::path &path) {
        std::ofstream stream(path);
        stream.exceptions(std::ios::badbit | std::ios::failbit);
        stream << std::setprecision(17);
        return stream;
    }

    mfem::Vector Displacement(const fem::FEM &f, double scale, double phase = 0.0) {
        mfem::ParGridFunction field(f.displacementFes.get());
        mfem::VectorFunctionCoefficient coefficient(3, [=](const mfem::Vector &x, mfem::Vector &u) {
            u.SetSize(3);
            u(0) = scale * (0.04 * x(0) + 0.01 * x(1) * x(2) + phase * (0.012 + 0.007 * x(2)));
            u(1) = scale * (-0.03 * x(1) + 0.008 * x(0) * x(2) + phase * (0.009 * x(0) - 0.003));
            u(2) = scale * (0.02 * x(2) - 0.006 * x(0) * x(1) - phase * 0.011 * x(1));
        });
        field.ProjectCoefficient(coefficient);
        mfem::Vector result;
        field.GetTrueDofs(result);
        return result;
    }

    mfem::Vector Flux(int size, double phase) {
        mfem::Vector result(size);
        for (int i = 0; i < size; ++i)
            result(i) = std::sin(0.37 * (i + 1) + phase) + 0.31 * std::cos(0.19 * (i + 1) - 0.5 * phase);
        return result;
    }

    double Norm(const mfem::Vector &vector) {
        for (int i = 0; i < vector.Size(); ++i)
            if (!std::isfinite(vector(i))) throw std::runtime_error("An operator returned non-finite output.");
        const double result = vector.Norml2();
        if (!std::isfinite(result)) throw std::runtime_error("An output norm is non-finite.");
        return result;
    }

    double RelativeError(const mfem::Vector &actual, const mfem::Vector &expected) {
        if (actual.Size() != expected.Size()) throw std::runtime_error("Comparison vector sizes differ.");
        mfem::Vector difference(actual);
        difference -= expected;
        return Norm(difference) / std::max({Norm(actual), Norm(expected), 1.0e-30});
    }

    double Median(std::vector<double> values) {
        std::sort(values.begin(), values.end());
        const std::size_t middle = values.size() / 2;
        return values.size() % 2 ? values[middle] : 0.5 * (values[middle - 1] + values[middle]);
    }

    struct CheckReport {
        std::ofstream stream;
        bool passed{true};
        void Add(const std::string &name, double error, double tolerance) {
            const bool okay = std::isfinite(error) && error <= tolerance;
            stream << name << ',' << error << ',' << tolerance << ',' << okay << '\n';
            stream.flush();
            passed = passed && okay;
            std::cout << "Check " << name << ": " << error << (okay ? " PASS" : " FAIL") << std::endl;
        }
    };

    // All systems use this timing driver. Construction, checks, warmup and output
    // reductions are excluded. Alternate order to reduce systematic A/B bias.
    void Benchmark(const Options &options, const std::string &system,
                   const std::string &referenceName, const Apply &optimized, const Apply &reference,
                   std::ofstream &timings, std::ofstream &summary, double extraPrepareSeconds = 0.0) {
        const bool runOptimized = options.implementation != "reference";
        const bool runReference = options.implementation != "optimized";
        mfem::Vector output;
        for (int iteration = 0; iteration < options.warmup; ++iteration) {
            if (runOptimized) { optimized(iteration % options.directions, output); Norm(output); }
            if (runReference) { reference(iteration % options.directions, output); Norm(output); }
        }
        std::vector<double> optimizedTimes, referenceTimes;
        auto sample = [&](const std::string &name, const Apply &apply, std::vector<double> &times, int index) {
            const auto start = Clock::now();
            for (int iteration = 0; iteration < options.iterations; ++iteration)
                apply(iteration % options.directions, output);
            const double seconds = std::chrono::duration<double>(Clock::now() - start).count();
            const double perCall = seconds / options.iterations;
            times.push_back(perCall);
            timings << system << ',' << name << ',' << index << ',' << options.iterations << ','
                    << seconds << ',' << perCall << ',' << Norm(output) << '\n';
            timings.flush();
            std::cout << system << '/' << name << " sample " << index + 1 << ": " << perCall << " s/apply" << std::endl;
        };
        for (int index = 0; index < options.samples; ++index) {
            if (index % 2 == 0 && runReference) sample(referenceName, reference, referenceTimes, index);
            if (runOptimized) sample("optimized", optimized, optimizedTimes, index);
            if (index % 2 != 0 && runReference) sample(referenceName, reference, referenceTimes, index);
        }
        if (runOptimized) summary << system << ".optimized_median_seconds=" << Median(optimizedTimes) << '\n';
        if (runReference) summary << system << '.' << referenceName << "_median_seconds=" << Median(referenceTimes) << '\n';
        if (runOptimized && runReference) {
            const double speedup = Median(referenceTimes) / Median(optimizedTimes);
            summary << system << ".reference_over_optimized=" << speedup << '\n';
            if (extraPrepareSeconds > 0.0) {
                const double saving = Median(referenceTimes) - Median(optimizedTimes);
                summary << system << ".row_prepare_seconds=" << extraPrepareSeconds << '\n';
                summary << system << ".row_prepare_break_even_applications="
                        << (saving > 0.0 ? std::ceil(extraPrepareSeconds / saving) : -1.0) << '\n';
                for (const int applications : {1, 10, 25})
                    summary << system << ".reference_over_optimized_including_row_prepare_" << applications << '='
                            << applications * Median(referenceTimes) /
                               (extraPrepareSeconds + applications * Median(optimizedTimes)) << '\n';
            }
            std::cout << system << " reference/optimized median ratio: " << speedup << std::endl;
        }
        summary.flush();
    }

    // State/context updates happen before either timer. Each timed operation
    // must rebuild once; changing revisions prevents accidental no-op timings.
    void BenchmarkPreparation(
        const Options &options, const std::string &system, const std::function<void()> &advance,
        const std::function<void()> &optimized, const std::function<void()> &reference,
        const std::function<void(mfem::Vector &)> &optimizedResidual,
        const std::function<void(mfem::Vector &)> &referenceResidual,
        std::ofstream &timings, std::ofstream &summary
    ) {
        if (options.prepareSamples == 0) return;
        const bool runOptimized = options.implementation != "reference";
        const bool runReference = options.implementation != "optimized";
        std::vector<double> optimizedTimes, referenceTimes;
        mfem::Vector residual;
        for (int sample = 0; sample < options.prepareSamples; ++sample) {
            double fastSeconds = 0.0, referenceSeconds = 0.0;
            auto measure = [&](const std::function<void()> &operation, double &seconds) {
                const auto start = Clock::now();
                operation();
                seconds += std::chrono::duration<double>(Clock::now() - start).count();
            };
            for (int iteration = 0; iteration < options.prepareIterations; ++iteration) {
                advance();
                if (sample % 2 == 0 && runReference) measure(reference, referenceSeconds);
                if (runOptimized) measure(optimized, fastSeconds);
                if (sample % 2 != 0 && runReference) measure(reference, referenceSeconds);
            }
            auto record = [&](const std::string &name, double seconds, auto &values, const auto &getResidual) {
                getResidual(residual); // Residual copy/norm is outside preparation timing.
                const double perCall = seconds / options.prepareIterations;
                values.push_back(perCall);
                timings << system << ',' << name << ',' << sample << ',' << options.prepareIterations << ','
                        << seconds << ',' << perCall << ',' << Norm(residual) << '\n';
                std::cout << system << '/' << name << " sample " << sample + 1 << ": "
                          << perCall << " s/prepare" << std::endl;
            };
            if (runOptimized) record("optimized", fastSeconds, optimizedTimes, optimizedResidual);
            if (runReference) record("reference", referenceSeconds, referenceTimes, referenceResidual);
            timings.flush();
        }
        if (runOptimized) summary << system << ".optimized_median_seconds=" << Median(optimizedTimes) << '\n';
        if (runReference) summary << system << ".reference_median_seconds=" << Median(referenceTimes) << '\n';
        if (runOptimized && runReference)
            summary << system << ".reference_over_optimized=" << Median(referenceTimes) / Median(optimizedTimes) << '\n';
        summary.flush();
    }

    mfem::Vector ScalarField(mfem::ParFiniteElementSpace &space, double offset, double phase) {
        mfem::ParGridFunction field(&space);
        mfem::FunctionCoefficient coefficient([=](const mfem::Vector &x) {
            return offset + 0.09 * std::sin(0.71 * x(0) + phase) +
                   0.05 * std::cos(0.63 * x(1) - phase) + 0.02 * x(2);
        });
        field.ProjectCoefficient(coefficient);
        mfem::Vector result;
        field.GetTrueDofs(result);
        return result;
    }

    // Mixed relative/absolute error also handles zero scalar columns without
    // amplifying harmless cancellation into an artificial relative failure.
    double ScaledError(const mfem::Vector &actual, const mfem::Vector &expected) {
        if (actual.Size() != expected.Size()) throw std::runtime_error("Comparison vector sizes differ.");
        mfem::Vector difference(actual);
        difference -= expected;
        return Norm(difference) / std::max({1.0, Norm(actual), Norm(expected)});
    }

    int RunStellarDerivatives(const Options &options, fem::FEM &f, std::ofstream &metadata,
                             std::ofstream &summary, std::ofstream &timings, CheckReport &checks) {
        const bool group = options.system == "all" || options.system == "stellar-derivatives";
        const bool runMass = group || options.system == "mass-normalization";
        const bool runAngular = group || options.system == "angular-momentum";
        const bool runHydro = group || options.system == "hydrostatic-displacement";
        using GravityContext = operators::context::gravity_field::GravityFieldLinearizationContext;
        std::unique_ptr<GravityContext> gravity;
        std::unique_ptr<operators::PreparedMassNormalizationOperator> mass;
        std::unique_ptr<operators::PreparedAngularMomentumOperator> angular;
        std::unique_ptr<operators::PreparedHydrostaticEquilibriumOperator> hydro;
        const auto constructStart = Clock::now();
        if (runMass || runAngular) gravity = std::make_unique<GravityContext>(f, *f.domainMapperStateless);
        if (runMass) mass = std::make_unique<operators::PreparedMassNormalizationOperator>(f, *f.domainMapperStateless, *gravity);
        if (runAngular) angular = std::make_unique<operators::PreparedAngularMomentumOperator>(
            f, *f.domainMapperStateless, *gravity,
            models::compileConstraint(integral::FixedAngularMomentum(
                {.Jtotal = dimensions::AngularMomentumValue{0.41}, .axis = {1.0, -2.0, 3.0}, .center = {0.12, -0.07, 0.03}})));
        if (runHydro) hydro = std::make_unique<operators::PreparedHydrostaticEquilibriumOperator>(f, *f.domainMapperStateless);
        summary << "stellar.operator_construction_seconds=" << std::chrono::duration<double>(Clock::now() - constructStart).count() << '\n';

        const auto baseDensity = ScalarField(*f.densityFes, 0.94, 0.31);
        const auto baseDisplacement = Displacement(f, options.deformation);
        const auto enthalpy = ScalarField(*f.enthalpyFes, 0.73, 0.27);
        const auto potential = ScalarField(*f.gravityPotentialFes, -0.38, 0.61);
        mfem::Vector flux(f.gravityFluxFes->GetTrueVSize());
        flux = 0.0;
        operators::MassNormalizationDependencies md{
            .discretization = {1, 1}, .density = {2, 1}, .displacement = {3, 1}, .targetMass = {4, 1}};
        operators::AngularMomentumDependencies ad{
            .discretization = {1, 1}, .density = {2, 1}, .displacement = {3, 1}, .rotation = {5, 1}};
        operators::context::hydrostatic::HydrostaticEquilibriumDependencies hd{
            .discretization = {1, 1}, .enthalpy = {6, 1}, .gravityPotential = {7, 1},
            .displacement = {3, 1}, .rotation = {5, 1}, .bernoulliConstant = {8, 1}};
        double targetMass = 1.23;
        auto prepare = [&](const mfem::Vector &density, const mfem::Vector &displacement, double omega,
                           bool geometryChanged, bool densityChanged, bool rotationChanged,
                           const std::string &timingLabel) {
            if (geometryChanged) ++md.displacement.revision;
            if (densityChanged) ++md.density.revision;
            if (rotationChanged) ++ad.rotation.revision;
            ad.density.revision = md.density.revision;
            ad.displacement.revision = md.displacement.revision;
            hd.displacement.revision = md.displacement.revision;
            hd.rotation.revision = ad.rotation.revision;
            auto time = [&](const std::string &name, auto &&operation) {
                const auto start = Clock::now();
                operation();
                if (!timingLabel.empty()) summary << name << '.' << timingLabel << "_seconds="
                    << std::chrono::duration<double>(Clock::now() - start).count() << '\n';
            };
            if (gravity) time("stellar.shared_gravity_prepare", [&] {
                gravity->Prepare(
                    {.density = gravity->GetDensityMap().gather(density),
                     .displacement = gravity->GetDisplacementMap().gather(displacement),
                     .gravity_gradient = gravity->GetGravityGradientMap().gather(flux),
                     .gravity_potential = gravity->GetGravityPotentialMap().gather(potential)},
                    {.discretization = {1}, .displacement = {md.displacement.revision},
                     .density = {md.density.revision}, .gravity_gradient = {1}, .gravity_potential = {1}});
            });
            if (mass) time("mass-normalization.prepare", [&] { mass->Prepare({.targetMass = targetMass}, md); });
            if (angular) time("angular-momentum.prepare", [&] { angular->Prepare(omega, ad); });
            if (hydro) time("hydrostatic-displacement.prepare", [&] {
                mfem::Vector velocity(3), center(3);
                velocity(0) = omega / std::sqrt(14.0); velocity(1) = -2.0 * velocity(0); velocity(2) = 3.0 * velocity(0);
                center(0) = 0.12; center(1) = -0.07; center(2) = 0.03;
                hydro->Prepare({.enthalpy = hydro->GetEnthalpyMap().gather(enthalpy),
                                .gravityPotential = hydro->GetGravityPotentialMap().gather(potential),
                                .displacement = hydro->GetDisplacementMap().gather(displacement),
                                .bernoulliConstant = 0.39}, hd, physics::RigidRotation(velocity, center));
            });
        };
        prepare(baseDensity, baseDisplacement, 0.63, true, true, true, "initial");
        std::vector<mfem::Vector> densityDirections, fullDensityDirections, displacementDirections,
                                  hydroDirections, fullDisplacementDirections;
        std::vector<double> omegaDirections;
        for (int i = 0; i < options.directions; ++i) {
            fullDensityDirections.push_back(ScalarField(*f.densityFes, 0.17, 0.4 + i));
            fullDisplacementDirections.push_back(Displacement(f, 0.43, 0.2 + 0.7 * i));
            // Include a direction with variation at the element scale as well
            // as smooth shear/translation directions. It is not a base state.
            if (i == 1) fullDisplacementDirections.back().Add(0.002, Flux(f.displacementFes->GetTrueVSize(), 0.37));
            if (gravity) {
                densityDirections.push_back(gravity->GetDensityMap().gather(fullDensityDirections.back()));
                displacementDirections.push_back(gravity->GetDisplacementMap().gather(fullDisplacementDirections.back()));
            }
            if (hydro) hydroDirections.push_back(hydro->GetDisplacementMap().gather(fullDisplacementDirections.back()));
            omegaDirections.push_back(-0.13 + 0.07 * i);
        }
        const Apply massFast = [&](int i, mfem::Vector &out) { mass->ApplyCompleteJacobianAction(densityDirections[i], displacementDirections[i], out); };
        const Apply massReference = [&](int i, mfem::Vector &out) { mass->ApplyCompleteJacobianActionReference(densityDirections[i], displacementDirections[i], out); };
        const Apply angularFast = [&](int i, mfem::Vector &out) { angular->ApplyCompleteJacobianAction(densityDirections[i], displacementDirections[i], omegaDirections[i], out); };
        const Apply angularReference = [&](int i, mfem::Vector &out) { angular->ApplyCompleteJacobianActionReference(densityDirections[i], displacementDirections[i], omegaDirections[i], out); };
        const Apply hydroFast = [&](int i, mfem::Vector &out) { hydro->ApplyDisplacementJacobianAction(hydroDirections[i], out); };
        const Apply hydroReference = [&](int i, mfem::Vector &out) { hydro->ApplyDisplacementJacobianActionCachedReference(hydroDirections[i], out); };
        const Apply hydroGeneralReference = [&](int i, mfem::Vector &out) { hydro->ApplyDisplacementJacobianActionReference(hydroDirections[i], out); };
        auto compare = [&](const std::string &name, const Apply &fast, const Apply &reference) {
            for (int i = 0; i < options.directions; ++i) {
                mfem::Vector a, b;
                reference(i, b);
                fast(i, a);
                checks.Add(name + "_direction" + std::to_string(i), ScaledError(a, b), 2.0e-11);
            }
        };
        if (options.check != "none") {
            const auto validationStart = Clock::now();
            if (mass) checks.Add("mass_rows_initially_lazy", mass->GetJacobianRowPreparationCount(), 0.0);
            if (angular) checks.Add("angular_rows_initially_lazy", angular->GetJacobianRowPreparationCount(), 0.0);
            const int states = options.check == "full" ? 4 : 1;
            for (int state = 0; state < states; ++state) {
                const auto density = state >= 2 ? ScalarField(*f.densityFes, 1.13, 0.82) : baseDensity;
                const auto displacement = state == 0 && states > 1 ? Displacement(f, 0.0) : baseDisplacement;
                prepare(density, displacement, state == 3 ? -0.27 : 0.63, state < 2, state != 3, true, "");
                const std::string label = "state" + std::to_string(state);
                const auto massCount = mass ? mass->GetJacobianRowPreparationCount() : 0;
                const auto angularCount = angular ? angular->GetJacobianRowPreparationCount() : 0;
                if (mass) {
                    compare(label + "_mass", massFast, massReference);
                    checks.Add(label + "_mass_row_refresh", mass->GetJacobianRowPreparationCount() == massCount + (state != 3) ? 0.0 : 1.0, 0.0);
                }
                if (angular) {
                    compare(label + "_angular", angularFast, angularReference);
                    checks.Add(label + "_angular_row_refresh", angular->GetJacobianRowPreparationCount() == angularCount + (state != 3) ? 0.0 : 1.0, 0.0);
                }
                if (hydro) {
                    compare(label + "_hydro_cached", hydroFast, hydroReference);
                    compare(label + "_hydro_general", hydroFast, hydroGeneralReference);
                    mfem::Vector zero(hydroDirections[0].Size()), action;
                    zero = 0.0;
                    hydro->ApplyDisplacementJacobianAction(zero, action);
                    checks.Add(label + "_hydro_zero", Norm(action), 0.0);
                }
            }
            prepare(baseDensity, baseDisplacement, 0.63, true, true, true, "");
            mfem::Vector massAnalytic, angularAnalytic, hydroAnalytic;
            if (mass) massFast(0, massAnalytic);
            if (angular) angularFast(0, angularAnalytic);
            if (hydro) hydroFast(0, hydroAnalytic);
            if (options.check == "full") {
                // Check individual scalar columns, including zero input, so
                // cancellation in a combined direction cannot hide an error.
                if (gravity) {
                    mfem::Vector zeroDensity(densityDirections[0].Size()), zeroDisplacement(displacementDirections[0].Size());
                    zeroDensity = 0.0; zeroDisplacement = 0.0;
                    for (int column = 0; column < 4; ++column) {
                        const auto &rho = column == 0 ? densityDirections[0] : zeroDensity;
                        const auto &u = column == 1 ? displacementDirections[0] : zeroDisplacement;
                        const double omega = column == 2 ? omegaDirections[0] : 0.0;
                        mfem::Vector a, b;
                        if (mass) {
                            mass->ApplyCompleteJacobianAction(rho, u, a);
                            mass->ApplyCompleteJacobianActionReference(rho, u, b);
                            checks.Add("mass_column" + std::to_string(column), ScaledError(a, b), 2.0e-11);
                        }
                        if (angular) {
                            angular->ApplyCompleteJacobianAction(rho, u, omega, a);
                            angular->ApplyCompleteJacobianActionReference(rho, u, omega, b);
                            checks.Add("angular_column" + std::to_string(column), ScaledError(a, b), 2.0e-11);
                        }
                    }
                }
                if (mass) {
                    mfem::Vector rhoDual, displacementDual;
                    mass->ApplyCompleteJacobianTransposeAction(0.73, rhoDual, displacementDual);
                    mfem::Vector expected(1), actual(1);
                    expected(0) = 0.73 * massAnalytic(0);
                    actual(0) = rhoDual * densityDirections[0] + displacementDual * displacementDirections[0];
                    checks.Add("mass_transpose_identity", ScaledError(actual, expected), 2.0e-11);
                    const auto count = mass->GetJacobianRowPreparationCount();
                    targetMass = 1.37; ++md.targetMass.revision;
                    mass->Prepare({.targetMass = targetMass}, md);
                    mfem::Vector afterTarget;
                    massFast(0, afterTarget);
                    checks.Add("mass_target_only_action", ScaledError(afterTarget, massAnalytic), 2.0e-11);
                    checks.Add("mass_target_only_reuses_rows", mass->GetJacobianRowPreparationCount() == count ? 0.0 : 1.0, 0.0);
                }
                for (const double step : {options.finiteDifferenceStep, 0.5 * options.finiteDifferenceStep}) {
                    mfem::Vector massPlus, angularPlus, hydroPlus, massMinus, angularMinus, hydroMinus;
                    for (const int sign : {1, -1}) {
                        mfem::Vector rho(baseDensity), displacement(baseDisplacement);
                        rho.Add(sign * step, fullDensityDirections[0]);
                        displacement.Add(sign * step, fullDisplacementDirections[0]);
                        // Scalar FD varies rho/u/Omega together. Hydro FD keeps
                        // all fields and rotation fixed except displacement.
                        prepare(rho, displacement, 0.63, true, true, true, "");
                        if (mass) mass->BuildResidual(sign == 1 ? massPlus : massMinus);
                        if (angular) {
                            ++ad.rotation.revision;
                            angular->Prepare(0.63 + sign * step * omegaDirections[0], ad);
                            angular->BuildResidual(sign == 1 ? angularPlus : angularMinus);
                        }
                        if (hydro) hydro->BuildResidual(sign == 1 ? hydroPlus : hydroMinus);
                    }
                    auto difference = [&](const std::string &name, const mfem::Vector &analytic, mfem::Vector &plus, const mfem::Vector &minus) {
                        plus -= minus; plus /= 2.0 * step;
                        checks.Add(name + "_finite_difference_" + std::to_string(step), ScaledError(analytic, plus), 2.0e-6);
                    };
                    if (mass) difference("mass", massAnalytic, massPlus, massMinus);
                    if (angular) difference("angular", angularAnalytic, angularPlus, angularMinus);
                    if (hydro) difference("hydro", hydroAnalytic, hydroPlus, hydroMinus);
                }
            }
            summary << "stellar.validation_seconds=" << std::chrono::duration<double>(Clock::now() - validationStart).count() << '\n';
            if (!checks.passed) { summary << "stellar.checks_passed=0\nstellar.benchmark_skipped=1\n"; return 3; }
        }
        prepare(baseDensity, baseDisplacement, 0.63, true, true, true, "repeat");
        auto rows = [&](const std::string &name, auto &operation) {
            if (!operation) return 0.0;
            const auto start = Clock::now();
            const bool supported = operation->PrepareJacobianRows();
            const double seconds = std::chrono::duration<double>(Clock::now() - start).count();
            metadata << name << ".prepared_rows_supported=" << supported << '\n';
            summary << name << ".explicit_row_prepare_seconds=" << seconds << '\n';
            return seconds;
        };
        const double massRowSeconds = rows("mass-normalization", mass);
        const double angularRowSeconds = rows("angular-momentum", angular);
        if (hydro) {
            metadata << "hydrostatic-displacement.fast_elements=" << hydro->GetFastDisplacementElementCount() << '\n'
                     << "hydrostatic-displacement.tensor_elements=" << hydro->GetTensorDisplacementElementCount() << '\n'
                     << "hydrostatic-displacement.reference=previous_cached_pointwise\n";
            std::size_t nq = 0;
            hydro->VisitMappedGeometryRules([&](int, const mfem::IntegrationRule &rule) { nq += rule.GetNPoints(); });
            metadata << "hydrostatic-displacement.quadrature_points=" << nq << '\n';
        }
        metadata << "stellar.check_error_metric=norm_difference/max(1,norm_actual,norm_reference)\n";
        metadata.flush();
        summary << "stellar.checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
        if (mass) Benchmark(options, "mass-normalization", "prepared_reference", massFast, massReference, timings, summary, massRowSeconds);
        if (angular) Benchmark(options, "angular-momentum", "prepared_reference", angularFast, angularReference, timings, summary, angularRowSeconds);
        if (hydro) Benchmark(options, "hydrostatic-displacement", "cached_pointwise_reference", hydroFast, hydroReference, timings, summary);
        return 0;
    }


    bool IsLinearDerivativeSystem(const std::string &name) {
        return name == "application-optimizations" || name == "linear-derivatives" || name == "barotropic-closure" ||
            name == "barotropic-enthalpy" || name == "barotropic-displacement" ||
            name == "pressure-displacement" || name == "pressure-enthalpy" ||
            name == "gravity-force" || name == "rotation-force" || name == "preparation-optimizations";
    }

    int RunLinearDerivatives(const Options &options, fem::FEM &f, std::ofstream &metadata,
                             std::ofstream &summary, std::ofstream &timings, CheckReport &checks) {
        const bool preparationGroup = options.system == "preparation-optimizations";
        const bool group = preparationGroup || options.system == "all" || options.system == "linear-derivatives";
        const bool applicationGroup = options.system == "application-optimizations";
        const bool runClosure = applicationGroup || (group && !preparationGroup) || options.system.starts_with("barotropic-");
        const bool runPressure = group || options.system.starts_with("pressure-");
        const bool runGravity = group || options.system == "gravity-force";
        const bool runRotation = applicationGroup || group || options.system == "rotation-force";
        using Schema = utils::domain::CoreEnvelopeVacuumDomainSchema;
        const auto densityMap = field::make_field_dof_map<field::Density, Schema>(*f.densityFes);
        const auto enthalpyMap = field::make_field_dof_map<field::Enthalpy, Schema>(*f.enthalpyFes);
        const auto displacementMap = field::make_field_dof_map<field::Displacement, Schema>(*f.displacementFes);
        const eos::Polytrope equationOfState(options.polytropicIndex, 1.0);
        using GravityContext = operators::context::gravity_field::GravityFieldLinearizationContext;
        const auto constructStart = Clock::now();
        std::unique_ptr<GravityContext> gravityContext;
        std::unique_ptr<operators::PreparedBarotropicClosureOperator> closure;
        std::unique_ptr<operators::PreparedPressureForceOperator> pressure, pressureReference;
        std::unique_ptr<operators::PreparedGravityDisplacementForceOperator> gravity, gravityReference;
        std::unique_ptr<operators::PreparedRotationalDisplacementForceOperator> rotation, rotationReference;
        if (runClosure) closure = std::make_unique<operators::PreparedBarotropicClosureOperator>(
            f, *f.domainMapperStateless, equationOfState);
        if (runPressure) {
            pressure = std::make_unique<operators::PreparedPressureForceOperator>(f, *f.domainMapperStateless, equationOfState);
            pressureReference = std::make_unique<operators::PreparedPressureForceOperator>(f, *f.domainMapperStateless, equationOfState);
        }
        if (runGravity) {
            gravityContext = std::make_unique<GravityContext>(f, *f.domainMapperStateless);
            gravity = std::make_unique<operators::PreparedGravityDisplacementForceOperator>(
                f, *f.domainMapperStateless, *gravityContext);
            gravityReference = std::make_unique<operators::PreparedGravityDisplacementForceOperator>(
                f, *f.domainMapperStateless, *gravityContext);
        }
        if (runRotation) {
            rotation = std::make_unique<operators::PreparedRotationalDisplacementForceOperator>(f, *f.domainMapperStateless);
            rotationReference = std::make_unique<operators::PreparedRotationalDisplacementForceOperator>(f, *f.domainMapperStateless);
        }
        summary << "linear.operator_construction_seconds=" <<
            std::chrono::duration<double>(Clock::now() - constructStart).count() << '\n';

        const auto baseDensity = densityMap.gather(ScalarField(*f.densityFes, 0.94, 0.31));
        const auto baseEnthalpy = enthalpyMap.gather(ScalarField(*f.enthalpyFes, 0.83, 0.27));
        const auto baseDisplacement = displacementMap.gather(Displacement(f, options.deformation));
        mfem::Vector baseFlux;
        mfem::Vector potential;
        if (gravityContext) {
            baseFlux = gravityContext->GetGravityGradientMap().gather(Flux(f.gravityFluxFes->GetTrueVSize(), 0.29));
            baseFlux *= 0.04;
            potential = gravityContext->GetGravityPotentialMap().gather(
                ScalarField(*f.gravityPotentialFes, -0.38, 0.61));
        }
        std::vector<mfem::Vector> rhoDirections, hDirections, uDirections, gDirections;
        for (int i = 0; i < options.directions; ++i) {
            rhoDirections.push_back(densityMap.gather(ScalarField(*f.densityFes, 0.17, 0.4 + i)));
            hDirections.push_back(enthalpyMap.gather(ScalarField(*f.enthalpyFes, 0.11, 0.6 + i)));
            auto u = Displacement(f, 0.43, 0.2 + 0.7 * i);
            if (i == 1) u.Add(0.002, Flux(f.displacementFes->GetTrueVSize(), 0.37));
            uDirections.push_back(displacementMap.gather(u));
            if (gravityContext) {
                gDirections.push_back(gravityContext->GetGravityGradientMap().gather(
                    Flux(f.gravityFluxFes->GetTrueVSize(), 0.81 + i)));
                gDirections.back() *= 0.03;
            }
        }
        mfem::Vector zeroRho(baseDensity.Size()), zeroH(baseEnthalpy.Size()),
                     zeroU(baseDisplacement.Size()), zeroG(baseFlux.Size());
        zeroRho = 0.0; zeroH = 0.0; zeroU = 0.0; zeroG = 0.0;
        // The last direction is reserved for checks and is never timed.
        rhoDirections.push_back(zeroRho); hDirections.push_back(zeroH); uDirections.push_back(zeroU);
        if (gravityContext) gDirections.push_back(zeroG);

        std::uint64_t densityRevision = 1, enthalpyRevision = 1, displacementRevision = 1,
                      fluxRevision = 1, rotationRevision = 1;
        const bool needPreparationReferences = options.check != "none" || options.implementation != "optimized";
        bool prepareReferences = needPreparationReferences;
        auto prepare = [&](const mfem::Vector &rho, const mfem::Vector &h, const mfem::Vector &u,
                           const mfem::Vector &g, double omega, unsigned changed, const std::string &label) {
            if (changed & 1) ++densityRevision;
            if (changed & 2) ++enthalpyRevision;
            if (changed & 4) ++displacementRevision;
            if (changed & 8) ++fluxRevision;
            if (changed & 16) ++rotationRevision;
            auto time = [&](const std::string &name, auto &&operation) {
                const auto start = Clock::now(); operation();
                if (!label.empty()) summary << name << ".prepare." << label << "_seconds="
                    << std::chrono::duration<double>(Clock::now() - start).count() << '\n';
            };
            if (gravityContext) time("linear.shared-gravity", [&] {
                gravityContext->Prepare(
                    {.density = rho, .displacement = u, .gravity_gradient = g, .gravity_potential = potential},
                    {.discretization = {1}, .displacement = {displacementRevision}, .density = {densityRevision},
                     .gravity_gradient = {fluxRevision}, .gravity_potential = {1}});
            });
            if (closure) time("barotropic-closure", [&] {
                closure->Prepare({.density = rho, .enthalpy = h, .displacement = u},
                    {.discretization = {1, 1}, .density = {2, densityRevision}, .enthalpy = {3, enthalpyRevision},
                     .displacement = {4, displacementRevision}});
            });
            if (pressure) time("pressure-displacement", [&] {
                pressure->Prepare({.enthalpy = h, .displacement = u},
                    {.discretization = {1, 1}, .enthalpy = {3, enthalpyRevision}, .displacement = {4, displacementRevision}});
            });
            if (pressureReference && prepareReferences) time("pressure-displacement.reference", [&] {
                pressureReference->PrepareCachedReference({.enthalpy = h, .displacement = u},
                    {.discretization = {1, 1}, .enthalpy = {3, enthalpyRevision}, .displacement = {4, displacementRevision}});
            });
            if (gravity) time("gravity-force", [&] { gravity->Prepare(); });
            if (gravityReference && prepareReferences) time("gravity-force.reference", [&] { gravityReference->PrepareReference(); });
            if (rotation) time("rotation-force", [&] {
                mfem::Vector velocity(3), center(3);
                velocity(0) = omega / std::sqrt(14.0); velocity(1) = -2.0 * velocity(0); velocity(2) = 3.0 * velocity(0);
                center(0) = 0.12; center(1) = -0.07; center(2) = 0.03;
                rotation->Prepare({.density = rho, .displacement = u},
                    {.discretization = {1, 1}, .density = {2, densityRevision},
                     .displacement = {4, displacementRevision}, .rotation = {5, rotationRevision}},
                    physics::RigidRotation(velocity, center));
            });
            if (rotationReference && prepareReferences) time("rotation-force.reference", [&] {
                mfem::Vector velocity(3), center(3);
                velocity(0) = omega / std::sqrt(14.0); velocity(1) = -2.0 * velocity(0); velocity(2) = 3.0 * velocity(0);
                center(0) = 0.12; center(1) = -0.07; center(2) = 0.03;
                rotationReference->PrepareReference({.density = rho, .displacement = u},
                    {.discretization = {1, 1}, .density = {2, densityRevision},
                     .displacement = {4, displacementRevision}, .rotation = {5, rotationRevision}},
                    physics::RigidRotation(velocity, center));
            });
        };
        prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 31, "initial");
        if (gravity) metadata << "gravity_tensor_elements=" << gravity->GetTensorElementCount() << '\n';

        struct System {
            std::string name;
            Apply fast, reference;
            std::function<void(mfem::Vector &)> residual;
            unsigned columns; // density=1, enthalpy=2, displacement=4, gravity=8
        };
        std::vector<System> systems;
        auto add = [&](System system) {
            if (applicationGroup || group || options.system == system.name) systems.push_back(std::move(system));
        };
        if (closure) {
            auto residual = [&](mfem::Vector &out) { closure->BuildResidual(out); };
            add({"barotropic-closure",
                [&](int i, mfem::Vector &out) { closure->Mult(rhoDirections[i], hDirections[i], uDirections[i], out); },
                [&](int i, mfem::Vector &out) { closure->MultCachedReference(rhoDirections[i], hDirections[i], uDirections[i], out); },
                residual, 7});
            add({"barotropic-enthalpy",
                [&](int i, mfem::Vector &out) { closure->ApplyEnthalpyJacobianAction(hDirections[i], out); },
                [&](int i, mfem::Vector &out) { closure->ApplyEnthalpyJacobianAction(hDirections[i], out, false); },
                residual, 2});
            add({"barotropic-displacement",
                [&](int i, mfem::Vector &out) { closure->ApplyDisplacementJacobianAction(uDirections[i], out); },
                [&](int i, mfem::Vector &out) { closure->ApplyDisplacementJacobianAction(uDirections[i], out, false); },
                residual, 4});
        }
        if (pressure) add({"pressure-displacement",
            [&](int i, mfem::Vector &out) { pressure->ApplyDisplacementJacobianAction(uDirections[i], out); },
            [&](int i, mfem::Vector &out) { pressure->ApplyDisplacementJacobianActionCachedReference(uDirections[i], out); },
            [&](mfem::Vector &out) { pressure->BuildResidual(out); }, 4});
        if (pressure) add({"pressure-enthalpy",
            [&](int i, mfem::Vector &out) { pressure->ApplyEnthalpyJacobianAction(hDirections[i], out); },
            [&](int i, mfem::Vector &out) { pressureReference->ApplyEnthalpyJacobianAction(hDirections[i], out); },
            [&](mfem::Vector &out) { pressure->BuildResidual(out); }, 2});
        if (gravity) add({"gravity-force",
            [&](int i, mfem::Vector &out) { gravity->ApplyCompleteJacobianAction(rhoDirections[i], uDirections[i], gDirections[i], out); },
            [&](int i, mfem::Vector &out) { gravity->ApplyCompleteJacobianActionCachedReference(rhoDirections[i], uDirections[i], gDirections[i], out); },
            [&](mfem::Vector &out) { gravity->BuildResidual(out); }, 13});
        if (rotation) add({"rotation-force",
            [&](int i, mfem::Vector &out) { rotation->ApplyCompleteJacobianAction(rhoDirections[i], uDirections[i], out); },
            [&](int i, mfem::Vector &out) { rotation->ApplyCompleteJacobianActionCachedReference(rhoDirections[i], uDirections[i], out); },
            [&](mfem::Vector &out) { rotation->BuildResidual(out); }, 5});

        auto compare = [&](const std::string &label, const mfem::Vector &a, const mfem::Vector &b) {
            checks.Add(label, ScaledError(a, b), 2.0e-11);
            if (std::max(Norm(a), Norm(b)) > 1.0e-14)
                checks.Add(label + "_relative", RelativeError(a, b), 2.0e-10);
        };
        if (options.check != "none") {
            const auto start = Clock::now();
            const int states = options.check == "full" ? 5 : 1;
            for (int state = 0; state < states; ++state) {
                const auto rho = state >= 2 ? densityMap.gather(ScalarField(*f.densityFes, 1.13, 0.82)) : baseDensity;
                const auto h = state >= 2 ? enthalpyMap.gather(ScalarField(*f.enthalpyFes, 0.91, 0.42)) : baseEnthalpy;
                const auto &u = state == 0 && states > 1 ? zeroU : baseDisplacement;
                auto g = baseFlux;
                if (state >= 3) g *= -0.71;
                const unsigned revisions = state == 0 ? 31u : state == 1 ? 4u : state == 2 ? 3u : state == 3 ? 24u : 16u;
                const double omega = state == 4 ? 0.0 : state == 3 ? -0.27 : 0.63;
                prepare(rho, h, u, g, omega, revisions, "");
                auto compareResidual = [&](const std::string &name, auto &fast, auto &reference) {
                    if (!fast) return;
                    mfem::Vector a, b;
                    fast->BuildResidual(a); reference->BuildResidual(b);
                    compare("state" + std::to_string(state) + "_" + name + "_prepared_residual", a, b);
                };
                compareResidual("pressure", pressure, pressureReference);
                compareResidual("gravity", gravity, gravityReference);
                compareResidual("rotation", rotation, rotationReference);
                if (pressure) {
                    for (int i = 0; i <= options.directions; ++i) {
                        mfem::Vector a, b;
                        pressure->ApplyEnthalpyJacobianAction(hDirections[i], a);
                        pressureReference->ApplyEnthalpyJacobianAction(hDirections[i], b);
                        compare("state" + std::to_string(state) + "_pressure_assembled_enthalpy_" + std::to_string(i), a, b);
                    }
                }
                for (const auto &system : systems) {
                    for (int i = 0; i <= options.directions; ++i) {
                        mfem::Vector a, b;
                        system.fast(i, a); system.reference(i, b);
                        compare("state" + std::to_string(state) + "_" + system.name + "_direction" + std::to_string(i), a, b);
                        if (system.name == "pressure-displacement") {
                            pressure->ApplyDisplacementJacobianActionReference(uDirections[i], b);
                            compare("state" + std::to_string(state) + "_pressure_general_" + std::to_string(i), a, b);
                        }
                        if (system.name == "gravity-force") {
                            gravity->ApplyCompleteJacobianActionReference(rhoDirections[i], uDirections[i], gDirections[i], b);
                            compare("state" + std::to_string(state) + "_gravity_general_" + std::to_string(i), a, b);
                        }
                        if (system.name == "rotation-force") {
                            rotation->ApplyCompleteJacobianActionReference(rhoDirections[i], uDirections[i], b);
                            compare("state" + std::to_string(state) + "_rotation_general_" + std::to_string(i), a, b);
                        }
                        if (system.name.starts_with("barotropic-")) {
                            closure->MultReference(system.columns & 1 ? rhoDirections[i] : zeroRho,
                                system.columns & 2 ? hDirections[i] : zeroH,
                                system.columns & 4 ? uDirections[i] : zeroU, b);
                            compare("state" + std::to_string(state) + "_" + system.name + "_general_" + std::to_string(i), a, b);
                        }
                        if (i == options.directions) checks.Add(system.name + "_zero", Norm(a), 0.0);
                    }
                }
            }
            prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 31, "");
            if (options.check == "full") {
                // Reference preparation was checked at every state above.
                // Finite differences below need only the production residual.
                prepareReferences = false;
                // Compare separate columns against the old complete action.
                // Individual gravity/rotation column methods remain independent
                // stateless kernels and serve as an additional cross-check.
                for (unsigned column : {1u, 2u, 4u, 8u}) {
                    const auto &rho = column == 1 ? rhoDirections[0] : zeroRho;
                    const auto &h = column == 2 ? hDirections[0] : zeroH;
                    const auto &u = column == 4 ? uDirections[0] : zeroU;
                    const auto &g = column == 8 && gravity ? gDirections[0] : zeroG;
                    mfem::Vector a, b, independent;
                    if (closure) {
                        closure->Mult(rho, h, u, a); closure->MultReference(rho, h, u, b);
                        compare("closure_column" + std::to_string(column), a, b);
                    }
                    if (gravity && column != 2) {
                        gravity->ApplyCompleteJacobianAction(rho, u, g, a);
                        gravity->ApplyCompleteJacobianActionReference(rho, u, g, b);
                        compare("gravity_column" + std::to_string(column), a, b);
                        if (column == 1) gravity->ApplyDensityJacobianAction(rho, independent);
                        if (column == 4) gravity->ApplyDisplacementJacobianAction(u, independent);
                        if (column == 8) gravity->ApplyGravityGradientJacobianAction(g, independent);
                        compare("gravity_stateless_column" + std::to_string(column), a, independent);
                    }
                    if (rotation && (column == 1 || column == 4)) {
                        rotation->ApplyCompleteJacobianAction(rho, u, a);
                        rotation->ApplyCompleteJacobianActionReference(rho, u, b);
                        compare("rotation_column" + std::to_string(column), a, b);
                        if (column == 1) rotation->ApplyDensityJacobianAction(rho, independent);
                        else rotation->ApplyDisplacementJacobianAction(u, independent);
                        compare("rotation_stateless_column" + std::to_string(column), a, independent);
                    }
                }
                for (const auto &system : systems) {
                    // The complete closure FD plus independent column checks
                    // covers its two couplings in a grouped run.
                    if (group && (system.name == "barotropic-enthalpy" || system.name == "barotropic-displacement")) continue;
                    prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 31, "");
                    mfem::Vector analytic;
                    system.fast(0, analytic);
                    for (double step : {options.finiteDifferenceStep, 0.5 * options.finiteDifferenceStep}) {
                        mfem::Vector plus, minus;
                        for (int sign : {1, -1}) {
                            auto rho = baseDensity, h = baseEnthalpy, u = baseDisplacement, g = baseFlux;
                            if (system.columns & 1) rho.Add(sign * step, rhoDirections[0]);
                            if (system.columns & 2) h.Add(sign * step, hDirections[0]);
                            if (system.columns & 4) u.Add(sign * step, uDirections[0]);
                            if (system.columns & 8) g.Add(sign * step, gDirections[0]);
                            prepare(rho, h, u, g, 0.63, 31, "");
                            system.residual(sign == 1 ? plus : minus);
                        }
                        plus -= minus; plus /= 2.0 * step;
                        checks.Add(system.name + "_finite_difference_" + std::to_string(step),
                                   ScaledError(analytic, plus), 2.0e-6);
                        if (Norm(analytic) > 1.0e-12)
                            checks.Add(system.name + "_finite_difference_relative_" + std::to_string(step),
                                       RelativeError(analytic, plus), 2.0e-6);
                    }
                }
            }
            prepareReferences = needPreparationReferences;
            summary << "linear.validation_seconds=" << std::chrono::duration<double>(Clock::now() - start).count() << '\n';
            if (!checks.passed) { summary << "linear.checks=failed\n"; return 3; }
        }
        prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 31, "repeat");
        // Unchanged dependencies must preserve the prepared state and results.
        const auto count = [&] {
            return (closure ? closure->GetPreparationCount() : 0) +
                (pressure ? pressure->GetResidualPreparationCount() : 0) +
                (gravity ? gravity->GetResidualPreparationCount() : 0) +
                (rotation ? rotation->GetResidualPreparationCount() : 0);
        };
        const auto beforeNoop = count();
        prepare(baseDensity, baseEnthalpy, baseDisplacement, baseFlux, 0.63, 0, "unchanged");
        if (options.check != "none") {
            checks.Add("linear_unchanged_dependencies_reuse", count() == beforeNoop ? 0.0 : 1.0, 0.0);
            if (!checks.passed) { summary << "linear.checks=failed\n"; return 3; }
        }
        if (gravity) metadata << "gravity-force.fast_elements=" << gravity->GetFastElementCount() << '\n';
        if (rotation) metadata << "rotation-force.fast_elements=" << rotation->GetFastElementCount()
            << "\nrotation_tensor_elements=" << rotation->GetTensorElementCount() << '\n';
        if (pressure) metadata << "pressure_tensor_elements=" << pressure->GetTensorElementCount() << '\n';
        if (closure) metadata << "closure_tensor_elements=" << closure->GetTensorElementCount() << '\n';
        metadata << "linear.polytropic_index=" << options.polytropicIndex << "\nlinear.polytropic_constant=1\n"
                 << "linear.check_error_metric=norm_difference/max(1,norm_actual,norm_reference)\n"
                 << "linear.relative_check_metric=norm_difference/max(norm_actual,norm_reference,1e-30)\n";
        metadata.flush();
        summary << "linear.checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
        for (const auto &system : systems)
            Benchmark(options, system.name, (system.name == "gravity-force" || system.name == "rotation-force" || system.name == "pressure-displacement" || system.name.starts_with("barotropic-")) ?
                          "cached_pointwise_reference" : "prepared_reference",
                      system.fast, system.reference, timings, summary);

        // Use the same physical state for both paths, with fresh dependency
        // revisions on each pair. Shared gravity context preparation is excluded.
        auto requireRebuild = [](auto &operation, auto &&prepareOperation) {
            const auto before = operation.GetResidualPreparationCount();
            prepareOperation();
            if (operation.GetResidualPreparationCount() != before + 1)
                throw std::runtime_error("Preparation benchmark timed an unchanged-state reuse.");
        };
        if (pressure) {
            auto dependencies = operators::context::pressure_force::PressureForceDependencies{
                .discretization = {1, 1}, .enthalpy = {3, enthalpyRevision}, .displacement = {4, displacementRevision}};
            const auto state = operators::context::pressure_force::PressureForceStateView{
                .enthalpy = baseEnthalpy, .displacement = baseDisplacement};
            BenchmarkPreparation(options, "pressure-preparation",
                [&] { dependencies.enthalpy.revision = ++enthalpyRevision;
                      dependencies.displacement.revision = ++displacementRevision; },
                [&] { requireRebuild(*pressure, [&] { pressure->Prepare(state, dependencies); }); },
                [&] { requireRebuild(*pressureReference, [&] { pressureReference->PrepareCachedReference(state, dependencies); }); },
                [&](mfem::Vector &out) { pressure->BuildResidual(out); },
                [&](mfem::Vector &out) { pressureReference->BuildResidual(out); }, timings, summary);
        }
        if (gravity) {
            BenchmarkPreparation(options, "gravity-force-preparation",
                [&] {
                    gravityContext->Prepare(
                        {.density = baseDensity, .displacement = baseDisplacement,
                         .gravity_gradient = baseFlux, .gravity_potential = potential},
                        {.discretization = {1}, .displacement = {++displacementRevision},
                         .density = {++densityRevision}, .gravity_gradient = {++fluxRevision}, .gravity_potential = {1}});
                },
                [&] { requireRebuild(*gravity, [&] { gravity->Prepare(); }); },
                [&] { requireRebuild(*gravityReference, [&] { gravityReference->PrepareReference(); }); },
                [&](mfem::Vector &out) { gravity->BuildResidual(out); },
                [&](mfem::Vector &out) { gravityReference->BuildResidual(out); }, timings, summary);
        }
        if (rotation) {
            auto dependencies = operators::context::rotational_displacement_force::RotationalDisplacementForceDependencies{
                .discretization = {1, 1}, .density = {2, densityRevision},
                .displacement = {4, displacementRevision}, .rotation = {5, rotationRevision}};
            const auto state = operators::context::rotational_displacement_force::RotationalDisplacementForceStateView{
                .density = baseDensity, .displacement = baseDisplacement};
            mfem::Vector velocity(3), center(3);
            velocity(0) = 0.63 / std::sqrt(14.0); velocity(1) = -2.0 * velocity(0); velocity(2) = 3.0 * velocity(0);
            center(0) = 0.12; center(1) = -0.07; center(2) = 0.03;
            const physics::RigidRotation rigidRotation(velocity, center);
            BenchmarkPreparation(options, "rotation-force-preparation",
                [&] { dependencies.density.revision = ++densityRevision;
                      dependencies.displacement.revision = ++displacementRevision;
                      dependencies.rotation.revision = ++rotationRevision; },
                [&] { requireRebuild(*rotation, [&] { rotation->Prepare(state, dependencies, rigidRotation); }); },
                [&] { requireRebuild(*rotationReference, [&] { rotationReference->PrepareReference(state, dependencies, rigidRotation); }); },
                [&](mfem::Vector &out) { rotation->BuildResidual(out); },
                [&](mfem::Vector &out) { rotationReference->BuildResidual(out); }, timings, summary);
        }
        return 0;
    }

    int RunNormalizationDiagonals(
        const Options &options, fem::FEM &f, std::ofstream &summary, std::ofstream &timings, CheckReport &checks
    ) {
        mfem::Array<int> stellarMarker(f.mesh->attributes.Max());
        stellarMarker = 0;
        for (int e = 0; e < f.mesh->GetNE(); ++e)
            if (!f.domainMapperStateless->IsCompactifiedElement(*f.mesh->GetElementTransformation(e)))
                stellarMarker[f.mesh->GetAttribute(e)-1] = 1;
        const auto surfaceMap = field::make_stellar_surface_scalar_dof_map(*f.surfaceDeformationFes);
        struct Diagonal { std::string name; std::function<mfem::Vector(bool)> assemble; };
        const std::vector<Diagonal> systems{
            {"normalization-hdiv", [&](bool ref) { return normalization::AssembleReferenceHDivMassDiagonal(*f.gravityFluxFes, ref); }},
            {"normalization-potential", [&](bool ref) { return normalization::AssembleReferenceScalarMassDiagonal(*f.gravityPotentialFes, nullptr, ref); }},
            {"normalization-density", [&](bool ref) { return normalization::AssembleReferenceScalarMassDiagonal(*f.densityFes, &stellarMarker, ref); }},
            {"normalization-enthalpy", [&](bool ref) { return normalization::AssembleReferenceScalarMassDiagonal(*f.enthalpyFes, &stellarMarker, ref); }},
            {"normalization-surface", [&](bool ref) { return normalization::AssembleReferenceSurfaceMassDiagonal(f, surfaceMap, ref); }}
        };
        for (const auto &system : systems) {
            mfem::Vector fast, reference;
            if (options.check != "none") {
                const auto start = Clock::now();
                fast = system.assemble(false);
                summary << system.name << ".first_optimized_seconds=" <<
                    std::chrono::duration<double>(Clock::now()-start).count() << '\n';
                reference = system.assemble(true);
                summary << system.name << ".optimized_diagonal_norm=" << Norm(fast) << '\n'
                        << system.name << ".reference_diagonal_norm=" << Norm(reference) << '\n';
                checks.Add(system.name + "_exact_true_diagonal", RelativeError(fast, reference), 2.0e-11);
                if (system.name == "normalization-surface") {
                    int nonpositiveEntries = 0;
                    for (int i = 0; i < fast.Size(); ++i)
                        if (fast(i) <= 0.0) ++nonpositiveEntries;
                    checks.Add(system.name + "_positive_surface_diagonal", nonpositiveEntries, 0.0);
                }
            }
            if (!checks.passed) { summary << "checks_passed=0\nbenchmark_skipped=1\n"; return 3; }
            BenchmarkPreparation(options, system.name, [] {},
                [&] { fast = system.assemble(false); }, [&] { reference = system.assemble(true); },
                [&](mfem::Vector &out) { out = fast; }, [&](mfem::Vector &out) { out = reference; }, timings, summary);
        }
        return 0;
    }

    int RunGeometryPreparation(
        const Options &options, fem::FEM &f, std::ofstream &metadata, std::ofstream &summary,
        std::ofstream &timings, CheckReport &checks
    ) {
        const auto &mapper = *f.domainMapperStateless;
        struct RestoreCacheMode {
            const mapping::DomainMapper &mapper;
            bool enabled;
            ~RestoreCacheMode() { mapper.SetGeometryCacheEnabled(enabled); }
        } restore{mapper, mapper.IsGeometryCacheEnabled()};
        operators::PreparedMappedHDivMassOperator fast(f, mapper), reference(f, mapper);
        const auto flux = Flux(f.gravityFluxFes->GetTrueVSize(), 0.37);
        const auto reducedFlux = fast.GetFluxMap().gather(flux);
        const auto direction = Displacement(f, 0.21, 0.32);
        const std::vector<mfem::Vector> states{
            fast.GetDisplacementMap().gather(Displacement(f, 0.0)),
            fast.GetDisplacementMap().gather(Displacement(f, options.deformation)),
            fast.GetDisplacementMap().gather(Displacement(f, 0.5*options.deformation, 0.4))
        };
        int state = 1;
        auto prepare = [&](auto &op, bool cache) {
            mapper.SetGeometryCacheEnabled(cache);
            const auto count = op.GetPreparationCount();
            op.Prepare(states[state]);
            if (op.GetPreparationCount() != count + 1)
                throw std::runtime_error("Geometry preparation benchmark did not rebuild its operator.");
        };
        mapper.ClearGeometryCache();
        const auto start = Clock::now();
        prepare(fast, true);
        summary << "geometry-preparation.cold_cache_prepare_seconds=" <<
            std::chrono::duration<double>(Clock::now()-start).count() << '\n';
        if (options.check != "none") {
            for (state = 0; state < static_cast<int>(states.size()); ++state) {
                prepare(fast, true); prepare(reference, false);
                mfem::Vector a, b;
                fast.Mult(reducedFlux, a); reference.Mult(reducedFlux, b);
                checks.Add("geometry_state" + std::to_string(state) + "_mass", RelativeError(a,b), 2.0e-11);
                fast.MultDisplacementVariationTrue(flux, direction, a);
                reference.MultDisplacementVariationTrue(flux, direction, b);
                checks.Add("geometry_state" + std::to_string(state) + "_derivative", RelativeError(a,b), 2.0e-11);
            }
        }
        if (!checks.passed) { summary << "checks_passed=0\nbenchmark_skipped=1\n"; return 3; }
        state = 0;
        const auto before = mapper.GetGeometryCacheStatistics();
        BenchmarkPreparation(options, "geometry-preparation",
            [&] { state = (state + 1) % static_cast<int>(states.size()); },
            [&] { prepare(fast, true); }, [&] { prepare(reference, false); },
            [&](mfem::Vector &out) { fast.Mult(reducedFlux, out); },
            [&](mfem::Vector &out) { reference.Mult(reducedFlux, out); }, timings, summary);
        const auto after = mapper.GetGeometryCacheStatistics();
        metadata << "geometry_cache_points=" << after.points << "\ngeometry_cache_hits=" << after.hits
                 << "\ngeometry_cache_misses=" << after.misses << '\n';
        if (options.check != "none" && options.prepareSamples > 0 && options.implementation != "reference") {
            checks.Add("geometry_warm_preparation_new_fixed_points", double(after.misses-before.misses), 0.0);
            checks.Add("geometry_warm_preparation_reuses_fixed_points", after.hits > before.hits ? 0.0 : 1.0, 0.0);
        }
        return checks.passed ? 0 : 3;
    }

    int RunGeometryInspection(
        const Options &options, fem::FEM &f, std::ofstream &metadata,
        std::ofstream &summary, std::ofstream &timings, CheckReport &checks
    ) {
        fem::PreparedGeometryInspector inspector;
        std::vector<mfem::Vector> states{
            Displacement(f, 0.0), Displacement(f, options.deformation),
            Displacement(f, 0.5 * options.deformation, 0.4)
        };
        const auto start = Clock::now();
        double fast = inspector.Inspect(*f.displacementFes, states[1]), reference = 0.0;
        summary << "geometry-inspection.cold_cache_seconds="
                << std::chrono::duration<double>(Clock::now()-start).count() << '\n';
        if (options.check != "none") {
            for (int i = 0; i < static_cast<int>(states.size()); ++i) {
                fast = inspector.Inspect(*f.displacementFes, states[i]);
                reference = inspector.Inspect(*f.displacementFes, states[i], true);
                checks.Add("geometry_inspection_state" + std::to_string(i),
                           std::abs(fast-reference) / std::max(1.0, std::abs(reference)), 2.0e-11);
                checks.Add("geometry_inspection_state" + std::to_string(i) + "_positive",
                           fast > 0.0 ? 0.0 : 1.0, 0.0);
            }
            if (options.check == "full") {
                mfem::ParGridFunction folded(f.displacementFes.get());
                mfem::VectorFunctionCoefficient fold(3, [](const mfem::Vector &x, mfem::Vector &u) {
                    u.SetSize(3); u = 0.0; u(0) = -2.0*x(0);
                });
                folded.ProjectCoefficient(fold);
                mfem::Vector foldedTrue; folded.GetTrueDofs(foldedTrue);
                fast = inspector.Inspect(*f.displacementFes, foldedTrue);
                reference = inspector.Inspect(*f.displacementFes, foldedTrue, true);
                checks.Add("geometry_inspection_fold_matches_reference", std::abs(fast-reference), 2.0e-11);
                checks.Add("geometry_inspection_fold_rejected", fast < 0.0 ? 0.0 : 1.0, 0.0);
                auto nonfinite = states[1];
                nonfinite(0) = std::numeric_limits<double>::quiet_NaN();
                fast = inspector.Inspect(*f.displacementFes, nonfinite);
                reference = inspector.Inspect(*f.displacementFes, nonfinite, true);
                checks.Add("geometry_inspection_nonfinite_rejected",
                           !std::isfinite(fast) && !std::isfinite(reference) ? 0.0 : 1.0, 0.0);
            }
        }
        if (!checks.passed) { summary << "checks_passed=0\nbenchmark_skipped=1\n"; return 3; }
        int state = 0;
        const auto builds = inspector.GetGeometryBuildCount();
        BenchmarkPreparation(options, "geometry-inspection",
            [&] { state = (state + 1) % static_cast<int>(states.size()); },
            [&] { fast = inspector.Inspect(*f.displacementFes, states[state]); },
            [&] { reference = inspector.Inspect(*f.displacementFes, states[state], true); },
            [&](mfem::Vector &out) { out.SetSize(1); out(0) = fast; },
            [&](mfem::Vector &out) { out.SetSize(1); out(0) = reference; }, timings, summary);
        metadata << "geometry_inspection_tensor_elements=" << inspector.GetTensorElementCount()
                 << "\ngeometry_inspection_geometry_builds=" << inspector.GetGeometryBuildCount() << '\n';
        if (options.check != "none")
            checks.Add("geometry_inspection_warm_cache_reuse", double(inspector.GetGeometryBuildCount()-builds), 0.0);
        summary << "geometry-inspection.checks=" << (options.check == "none" ? "not_run" : checks.passed ? "passed" : "failed") << '\n';
        return checks.passed ? 0 : 3;
    }

    int RunGravitySourcePreparation(
        const Options &options, fem::FEM &f, std::ofstream &summary,
        std::ofstream &timings, CheckReport &checks
    ) {
        operators::PreparedMappedGravitySourceOperator fast(f, *f.domainMapperStateless),
                                                      reference(f, *f.domainMapperStateless);
        const auto rhoTrue = ScalarField(*f.densityFes, 0.94, 0.31);
        const auto rho = fast.GetDensityMap().gather(rhoTrue);
        const auto potential = fast.GetPotentialMap().gather(ScalarField(*f.gravityPotentialFes, -0.38, 0.61));
        const auto direction = Displacement(f, 0.21, 0.32);
        const std::vector<mfem::Vector> states{
            fast.GetDisplacementMap().gather(Displacement(f, 0.0)),
            fast.GetDisplacementMap().gather(Displacement(f, options.deformation)),
            fast.GetDisplacementMap().gather(Displacement(f, 0.5 * options.deformation, 0.4))
        };
        int state = 1;
        auto prepare = [&](auto &op, bool oldWeights) {
            const auto count = op.GetPreparationCount();
            if (oldWeights) op.PrepareReference(states[state]); else op.Prepare(states[state]);
            if (op.GetPreparationCount() != count + 1)
                throw std::runtime_error("Gravity-source preparation benchmark did not rebuild.");
        };
        const auto start = Clock::now();
        prepare(fast, false);
        summary << "gravity-source-preparation.first_optimized_seconds="
                << std::chrono::duration<double>(Clock::now()-start).count() << '\n';
        if (options.check != "none") {
            for (state = 0; state < static_cast<int>(states.size()); ++state) {
                prepare(fast, false); prepare(reference, true);
                mfem::Vector a, b;
                fast.Mult(rho, a); reference.Mult(rho, b);
                checks.Add("source_state" + std::to_string(state) + "_mass", RelativeError(a,b), 2.0e-11);
                fast.MultTranspose(potential, a); reference.MultTranspose(potential, b);
                checks.Add("source_state" + std::to_string(state) + "_transpose", RelativeError(a,b), 2.0e-11);
                fast.MultDisplacementVariationTrue(rhoTrue, direction, a);
                reference.MultDisplacementVariationTrue(rhoTrue, direction, b);
                checks.Add("source_state" + std::to_string(state) + "_derivative", RelativeError(a,b), 2.0e-11);
            }
        }
        if (!checks.passed) { summary << "checks_passed=0\nbenchmark_skipped=1\n"; return 3; }
        state = 0;
        BenchmarkPreparation(options, "gravity-source-preparation",
            [&] { state = (state + 1) % static_cast<int>(states.size()); },
            [&] { prepare(fast, false); }, [&] { prepare(reference, true); },
            [&](mfem::Vector &out) { fast.Mult(rho, out); },
            [&](mfem::Vector &out) { reference.Mult(rho, out); }, timings, summary);
        summary << "gravity-source-preparation.checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
        return 0;
    }


    int RunSourceDisplacement(const Options &options, fem::FEM &f, std::ofstream &metadata,
        std::ofstream &summary, std::ofstream &timings, CheckReport &checks) {
        operators::PreparedMappedGravitySourceOperator source(f, *f.domainMapperStateless);
        const auto rho = ScalarField(*f.densityFes, 0.94, 0.31);
        const auto reducedRho = source.GetDensityMap().gather(rho);
        const auto base = Displacement(f, options.deformation);
        std::vector<mfem::Vector> directions;
        for (int i = 0; i < options.directions; ++i) {
            directions.push_back(Displacement(f, 0.43, 0.2+0.7*i));
            if (i == 1) directions.back().Add(0.002, Flux(f.displacementFes->GetTrueVSize(), 0.37));
        }
        auto prepare = [&](const mfem::Vector &u) { source.Prepare(source.GetDisplacementMap().gather(u)); };
        const Apply fast = [&](int i, mfem::Vector &out) { source.MultDisplacementVariationTrue(rho, directions[i], out); };
        const Apply reference = [&](int i, mfem::Vector &out) {
            source.MultDisplacementVariationTrueReference(rho, directions[i], out);
        };
        if (options.check != "none") {
            for (int state = 0; state < (options.check == "full" ? 3 : 1); ++state) {
                const auto u = state == 0 ? base : Displacement(f, state == 1 ? 0.0 : 0.5*options.deformation, 0.3);
                prepare(u);
                for (int i = 0; i < options.directions; ++i) {
                    mfem::Vector a,b;
                    fast(i,a); reference(i,b);
                    const auto name = "source_displacement_state"+std::to_string(state)+"_direction"+std::to_string(i);
                    checks.Add(name, RelativeError(a,b), 2.0e-11);
                    if (options.check == "full" && i == 0) {
                        operators::kernels::apply_mapped_source_variation(f, *f.domainMapperStateless, rho, u, directions[i], b);
                        checks.Add(name+"_stateless", RelativeError(a,b), 2.0e-11);
                    }
                }
                mfem::Vector zero(directions[0]); zero=0.0;
                mfem::Vector a;
                source.MultDisplacementVariationTrue(rho,zero,a);
                checks.Add("source_displacement_zero_state"+std::to_string(state),Norm(a),0.0);
            }
            prepare(base);
            if (options.check == "full") {
                mfem::Vector analytic; fast(0,analytic);
                for (double step : {options.finiteDifferenceStep,0.5*options.finiteDifferenceStep}) {
                    mfem::Vector up(base),um(base),a,b;
                    up.Add(step,directions[0]); um.Add(-step,directions[0]);
                    source.PreparePrimal(source.GetDisplacementMap().gather(up)); source.Mult(reducedRho,a);
                    source.PreparePrimal(source.GetDisplacementMap().gather(um)); source.Mult(reducedRho,b);
                    a-=b; a/=2.0*step;
                    auto expected=source.GetPotentialMap().gather(analytic);
                    checks.Add("source_displacement_fd_"+std::to_string(step),RelativeError(a,expected),2.0e-6);
                }
            }
        }
        if (!checks.passed) return 3;
        prepare(base);
        metadata << "source_tensor_variation_elements=" << source.GetTensorVariationElementCount() << '\n';
        Benchmark(options,"gravity-source-displacement","cached_pointwise_reference",fast,reference,timings,summary);
        return 0;
    }

    int RunHydrostaticPreparation(const Options &options, fem::FEM &f, std::ofstream &metadata,
        std::ofstream &summary, std::ofstream &timings, CheckReport &checks) {
        operators::PreparedHydrostaticEquilibriumOperator fast(f,*f.domainMapperStateless),reference(f,*f.domainMapperStateless);
        auto h = fast.GetEnthalpyMap().gather(ScalarField(*f.enthalpyFes,0.73,0.27));
        auto phi = fast.GetGravityPotentialMap().gather(ScalarField(*f.gravityPotentialFes,-0.38,0.61));
        const auto dh = fast.GetEnthalpyMap().gather(ScalarField(*f.enthalpyFes,0.13,0.41));
        const auto dp = fast.GetGravityPotentialMap().gather(ScalarField(*f.gravityPotentialFes,0.09,0.71));
        const auto du = fast.GetDisplacementMap().gather(Displacement(f,0.31,0.4));
        const std::vector<mfem::Vector> states{
            fast.GetDisplacementMap().gather(Displacement(f,options.deformation)),
            fast.GetDisplacementMap().gather(Displacement(f,0.0)),
            fast.GetDisplacementMap().gather(Displacement(f,0.5*options.deformation,0.3))};
        mfem::Vector velocity(3),center(3); center=0.0;
        velocity(0)=0.17; velocity(1)=-0.21; velocity(2)=0.39;
        const physics::RigidRotation rotation(velocity,center);
        operators::context::hydrostatic::HydrostaticEquilibriumDependencies dependencies{
            .discretization={1,1},.enthalpy={2,1},.gravityPotential={3,1},
            .displacement={4,1},.rotation={5,1},.bernoulliConstant={6,1}};
        int state=0;
        auto prepare = [&](auto &op,bool old) {
            const auto view=operators::context::hydrostatic::HydrostaticEquilibriumStateView{
                .enthalpy=h,.gravityPotential=phi,.displacement=states[state],.bernoulliConstant=0.39};
            const auto before=op.GetResidualPreparationCount();
            if (old) op.PrepareReference(view,dependencies,rotation); else op.Prepare(view,dependencies,rotation);
            if (op.GetResidualPreparationCount()!=before+1)
                throw std::runtime_error("Hydrostatic preparation benchmark did not rebuild.");
        };
        const auto start=Clock::now(); prepare(fast,false);
        summary << "hydrostatic-preparation.first_optimized_seconds=" << std::chrono::duration<double>(Clock::now()-start).count() << '\n';
        if (options.check!="none") {
            for (state=0;state<(options.check=="full"?3:1);++state) {
                ++dependencies.displacement.revision;
                ++dependencies.enthalpy.revision; ++dependencies.gravityPotential.revision;
                h*=1.01; phi*=0.99;
                prepare(fast,false); prepare(reference,true);
                auto compare=[&](const std::string &name,auto &&action) {
                    mfem::Vector a,b; action(fast,a); action(reference,b);
                    checks.Add("hydro_prepare_state"+std::to_string(state)+"_"+name,RelativeError(a,b),2.0e-11);
                };
                compare("residual",[](auto &op,auto &out){op.BuildResidual(out);});
                compare("enthalpy",[&](auto &op,auto &out){op.ApplyEnthalpyJacobianAction(dh,out);});
                compare("potential",[&](auto &op,auto &out){op.ApplyGravityPotentialJacobianAction(dp,out);});
                compare("bernoulli",[&](auto &op,auto &out){op.ApplyBernoulliConstantJacobianAction(0.13,out);});
                compare("displacement",[&](auto &op,auto &out){op.ApplyDisplacementJacobianAction(du,out);});
                compare("complete",[&](auto &op,auto &out){op.ApplyCompleteJacobianAction(dh,dp,0.13,du,out);});
                compare("diagonal",[](auto &op,auto &out){op.AssembleEnthalpyJacobianDiagonal(out);});
                const auto view=operators::context::hydrostatic::HydrostaticEquilibriumStateView{
                    .enthalpy=h,.gravityPotential=phi,.displacement=states[state],.bernoulliConstant=0.39};
                checks.Add("hydro_prepare_unchanged_state"+std::to_string(state),
                    fast.Prepare(view,dependencies,rotation).DidAnyWork()?1.0:0.0,0.0);
            }
        }
        if (options.check == "full") {
            state = 0;
            ++dependencies.displacement.revision;
            prepare(fast, false);
            mfem::Vector analytic;
            fast.ApplyCompleteJacobianAction(dh, dp, 0.13, du, analytic);
            for (double step : {options.finiteDifferenceStep, 0.5*options.finiteDifferenceStep}) {
                mfem::Vector hp(h), hm(h), pp(phi), pm(phi), up(states[0]), um(states[0]), a, b;
                hp.Add(step, dh); hm.Add(-step, dh);
                pp.Add(step, dp); pm.Add(-step, dp);
                up.Add(step, du); um.Add(-step, du);
                ++dependencies.enthalpy.revision;
                ++dependencies.gravityPotential.revision;
                ++dependencies.displacement.revision;
                ++dependencies.bernoulliConstant.revision;
                fast.Prepare({.enthalpy=hp, .gravityPotential=pp, .displacement=up,
                              .bernoulliConstant=0.39+step*0.13}, dependencies, rotation);
                fast.BuildResidual(a);
                ++dependencies.enthalpy.revision;
                ++dependencies.gravityPotential.revision;
                ++dependencies.displacement.revision;
                ++dependencies.bernoulliConstant.revision;
                fast.Prepare({.enthalpy=hm, .gravityPotential=pm, .displacement=um,
                              .bernoulliConstant=0.39-step*0.13}, dependencies, rotation);
                fast.BuildResidual(b);
                a -= b; a /= 2.0*step;
                checks.Add("hydro_prepare_complete_fd_"+std::to_string(step), RelativeError(a, analytic), 2.0e-6);
            }
            mfem::ParGridFunction folded(f.displacementFes.get());
            mfem::VectorFunctionCoefficient fold(3, [](const mfem::Vector &x, mfem::Vector &u) {
                u = x; u *= -2.0;
            });
            folded.ProjectCoefficient(fold);
            mfem::Vector foldedTrue; folded.GetTrueDofs(foldedTrue);
            const auto rejectedDisplacement = fast.GetDisplacementMap().gather(foldedTrue);
            ++dependencies.displacement.revision;
            const auto rejectedView = operators::context::hydrostatic::HydrostaticEquilibriumStateView{
                .enthalpy=h, .gravityPotential=phi, .displacement=rejectedDisplacement, .bernoulliConstant=0.39};
            const auto rejected = fast.TryPrepare(rejectedView, dependencies, rotation);
            const auto referenceRejected = reference.TryPrepareReference(rejectedView, dependencies, rotation);
            checks.Add("hydro_prepare_fold_rejected",
                !rejected && !referenceRejected &&
                rejected.error().mappingStatus == mapping::MappingStatus::non_positive_determinant &&
                rejected.error().mappingStatus == referenceRejected.error().mappingStatus ? 0.0 : 1.0, 0.0);
            checks.Add("hydro_prepare_rejected_state_unavailable",
                !fast.IsPrepared() && !reference.IsPrepared() ? 0.0 : 1.0, 0.0);
        }
        if(!checks.passed)return 3;
        state=0;
        ++dependencies.displacement.revision;
        ++dependencies.enthalpy.revision; ++dependencies.gravityPotential.revision;
        ++dependencies.bernoulliConstant.revision;
        prepare(fast,false); prepare(reference,true);
        metadata << "hydro_preparation_tensor_elements=" << fast.GetTensorDisplacementElementCount() << '\n';
        BenchmarkPreparation(options,"hydrostatic-preparation",
            [&]{state=(state+1)%static_cast<int>(states.size());++dependencies.displacement.revision;},
            [&]{prepare(fast,false);},[&]{prepare(reference,true);},
            [&](auto &out){fast.BuildResidual(out);},[&](auto &out){reference.BuildResidual(out);},timings,summary);
        return 0;
    }

    int RunGravitySchurPreparation(const Options &options,fem::FEM &f,std::ofstream &metadata,
        std::ofstream &summary,std::ofstream &timings,CheckReport &checks) {
        preconditioning::PreparedGravityDivergenceCache cache;
        operators::PreparedMappedHDivMassOperator mass(f,*f.domainMapperStateless);
        std::vector<mfem::Vector> diagonals(3);
        for(int i=0;i<3;++i) {
            mass.PreparePrimal(mass.GetDisplacementMap().gather(Displacement(f,0.5*i*options.deformation,0.2)));
            mass.AssembleTrueDiagonal(diagonals[i]);
        }
        const auto x=Flux(f.gravityPotentialFes->GetTrueVSize(),0.37);
        std::unique_ptr<mfem::HypreParMatrix> fast,reference;
        int state=0;
        auto prepareFast=[&]{fast=preconditioning::assembleGravityPotentialSchurSurrogate(f,diagonals[state],&cache);};
        auto prepareReference=[&]{reference=preconditioning::assembleGravityPotentialSchurSurrogate(f,diagonals[state]);};
        const auto start=Clock::now(); prepareFast();
        summary << "gravity-schur-preparation.cold_cache_seconds=" << std::chrono::duration<double>(Clock::now()-start).count() << '\n';
        auto compare=[&](const std::string &label){
            mfem::Vector a(fast->Height()),b(reference->Height());
            fast->Mult(x,a);reference->Mult(x,b);checks.Add(label,RelativeError(a,b),2.0e-11);
        };
        if(options.check!="none") {
            for(state=0;state<3;++state){prepareFast();prepareReference();compare("schur_state"+std::to_string(state));}
            // An exact scaling check also works with --deformation 0, and catches
            // both a stale Schur matrix and accidental cumulative scaling of B.
            const auto originalDiagonal = diagonals[0];
            state = 0; prepareFast();
            mfem::Vector originalAction(fast->Height()), changedAction(fast->Height());
            fast->Mult(x, originalAction);
            diagonals[0] *= 2.0; prepareFast(); prepareReference(); compare("schur_changed_mass_diagonal");
            fast->Mult(x, changedAction);
            originalAction *= 0.5;
            checks.Add("schur_inverse_mass_scaling", RelativeError(changedAction, originalAction), 2.0e-11);
            diagonals[0] = originalDiagonal; prepareFast(); prepareReference(); compare("schur_restored_mass_diagonal");
            checks.Add("divergence_shared_across_geometry",cache.GetBuildCount()==1?0.0:1.0,0.0);
            state=0;
            if(options.check=="full") {
                cache.Clear();prepareFast();prepareReference();compare("schur_after_clear");
                checks.Add("divergence_clear_rebuild",cache.GetBuildCount()==2?0.0:1.0,0.0);
                f.mesh->NodesUpdated();prepareFast();prepareReference();compare("schur_after_nodes_update");
                checks.Add("divergence_nodes_update_rebuild",cache.GetBuildCount()==3?0.0:1.0,0.0);
            }
        }
        if(!checks.passed)return 3;
        state=0;const auto builds=cache.GetBuildCount();
        BenchmarkPreparation(options,"gravity-schur-preparation",
            [&]{state=(state+1)%3;},prepareFast,prepareReference,
            [&](auto &out){out.SetSize(fast->Height());fast->Mult(x,out);},
            [&](auto &out){out.SetSize(reference->Height());reference->Mult(x,out);},timings,summary);
        if(options.check!="none")checks.Add("divergence_warm_rebuilds",double(cache.GetBuildCount()-builds),0.0);
        metadata << "fixed_divergence_assemblies=" << cache.GetBuildCount() << '\n';
        return checks.passed?0:3;
    }

    int Run(const Options &options) {
        if (options.help) {
            std::cout << "system_isolation_profile_harness\n"
                "  --system hdiv-displacement|hdiv-mass|mass-normalization|angular-momentum|\n"
                "           hydrostatic-displacement|stellar-derivatives|barotropic-closure|barotropic-enthalpy|\n"
                "           barotropic-displacement|pressure-displacement|pressure-enthalpy|gravity-force|rotation-force|\n"
                "           linear-derivatives|preparation-optimizations|runtime-optimizations|\n"
                "           geometry-preparation|normalization-diagonals|application-optimizations|\n"
                "           final-optimizations|gravity-source-displacement|hydrostatic-preparation|gravity-schur-preparation|\n"
                "           geometry-inspection|gravity-source-preparation|all (default hdiv-displacement)\n"
                "  --implementation both|optimized|reference (default both)\n"
                "  --check full|reference|none (default full)\n"
                "  --mesh sandbox.smesh --output NEW_DIRECTORY --refinements 0\n"
                "  --samples 5 --iterations 3 --warmup 1 --directions 3\n"
                "  --prepare-samples 3 --prepare-iterations 1 (prepare-samples 0 skips preparation timing)\n"
                "  --deformation 0.7 --fd-step 1e-5 --polytropic-index 1\n"
                "Single MPI rank, CPU. No Newton solve. Exit 0: requested checks pass;\n"
                "2: execution error; 3: check failure. Existing output directories are refused.\n"
                "hdiv-displacement compares the optimized and retained prepared derivatives.\n"
                "hdiv-mass compares the new CPU application with retained MFEM partial assembly.\n"
                "Full checks also compare the independent stateless mass kernel.\n";
            return 0;
        }
        if (!std::filesystem::create_directory(options.output))
            throw std::invalid_argument("Output directory already exists: " + options.output.string());
        auto metadata = File(options.output / "metadata.txt");
        auto summary = File(options.output / "summary.txt");
        auto timings = File(options.output / "timings.csv");
        CheckReport checks{File(options.output / "checks.csv")};
        checks.stream << "check,normalized_error,tolerance,passed\n";
        timings << "system,implementation,sample,applications,total_seconds,seconds_per_application,last_output_norm\n";
        metadata << "mesh=" << std::filesystem::absolute(options.mesh).string()
                 << "\nsystem=" << options.system << "\nimplementation=" << options.implementation
                 << "\ncheck=" << options.check << "\nrefinements=" << options.refinements
                 << "\nsamples=" << options.samples << "\niterations=" << options.iterations
                 << "\nwarmup=" << options.warmup << "\ndirections=" << options.directions
                 << "\nprepare_samples=" << options.prepareSamples << "\nprepare_iterations=" << options.prepareIterations
                 << "\npreparation_benchmark=forced_revisions_shared_gravity_excluded"
                 << "\npreparation_reference=original_matrix_assembly_and_separate_force_residual_traversals"
                 << "\nhydrostatic_reference=previous_cached_pointwise"
                 << "\nhdiv_displacement_reference=previous_optimized_pointwise"
                 << "\ngravity_force_reference=previous_optimized_pointwise"
                 << "\ngeometry_preparation_reference=same_operator_with_fixed_mesh_cache_disabled"
                 << "\nhdiv_mass_reference=previous_mfem_partial_assembly"
                 << "\nclosure_reference=previous_optimized_pointwise"
                 << "\nrotation_reference=previous_optimized_pointwise"
                 << "\ngeometry_inspection_reference=original_gridfunction_gradient_same_points"
                 << "\ngravity_source_preparation_reference=previous_mesh_weight_recomputation"
                 << "\npressure_displacement_reference=previous_optimized_pointwise"
                 << "\npressure_preparation_reference=previous_blocked_matrix_assembly_and_pointwise_fields"
                 << "\ngravity_source_displacement_reference=previous_pointwise"
                 << "\nhydrostatic_preparation_reference=previous_pointwise_fields_and_dense_matrix_assembly"
                 << "\ngravity_schur_preparation_reference=fixed_divergence_reassembled_every_call"
                 << "\ndeformation=" << options.deformation << "\nfd_step=" << options.finiteDifferenceStep
                 << "\nmpi_ranks=1\ndevice=cpu\ncompiler=" << __VERSION__
                 << "\nuniform_order_increment=" << MEAN_FIELD_UNIFORM_POLYNOMIAL_ORDER_INCREMENT << '\n';
        metadata.flush();

        utils::Args arguments;
        arguments.mesh_file = options.mesh;
        const auto setupStart = Clock::now();
        auto f = fem::setup_fem(arguments.mesh_file, arguments, options.refinements);
        if (!f.okay() || f.mesh->Dimension() != 3) throw std::runtime_error("A valid 3D discretization is required.");
        summary << "fem_setup_seconds=" << std::chrono::duration<double>(Clock::now() - setupStart).count() << '\n';
        int vacuumElements = 0;
        for (int element = 0; element < f.mesh->GetNE(); ++element)
            vacuumElements += f.domainMapperStateless->IsCompactifiedElement(*f.mesh->GetElementTransformation(element));
        metadata << "elements=" << f.mesh->GetNE() << "\nstellar_elements=" << f.mesh->GetNE() - vacuumElements
                 << "\nvacuum_elements=" << vacuumElements << "\nflux_true_dofs=" << f.gravityFluxFes->GetTrueVSize()
                 << "\ndisplacement_true_dofs=" << f.displacementFes->GetTrueVSize()
                 << "\ndensity_true_dofs=" << f.densityFes->GetTrueVSize()
                 << "\nenthalpy_true_dofs=" << f.enthalpyFes->GetTrueVSize() << '\n';

        const bool finalGroup = options.system == "final-optimizations";
        if (finalGroup || options.system == "gravity-source-displacement" ||
            options.system == "hydrostatic-preparation" || options.system == "gravity-schur-preparation") {
            int result=0;
            if(finalGroup) {
                auto pressureOptions=options;pressureOptions.system="pressure-displacement";
                result=RunLinearDerivatives(pressureOptions,f,metadata,summary,timings,checks);
                if(result)return result;
            }
            if(finalGroup || options.system=="gravity-source-displacement") {
                result=RunSourceDisplacement(options,f,metadata,summary,timings,checks);
                if(result)return result;
            }
            if(finalGroup || options.system=="hydrostatic-preparation") {
                result=RunHydrostaticPreparation(options,f,metadata,summary,timings,checks);
                if(result)return result;
            }
            if(finalGroup || options.system=="gravity-schur-preparation") {
                result=RunGravitySchurPreparation(options,f,metadata,summary,timings,checks);
                if(result)return result;
            }
            summary << "final.checks=" << (options.check=="none"?"not_run":"passed") << '\n';
            std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
            return 0;
        }
        if (options.system == "all") {
            int result = RunSourceDisplacement(options, f, metadata, summary, timings, checks);
            if (result != 0) return result;
            result = RunHydrostaticPreparation(options, f, metadata, summary, timings, checks);
            if (result != 0) return result;
            result = RunGravitySchurPreparation(options, f, metadata, summary, timings, checks);
            if (result != 0) return result;
        }
        const bool runtimeGroup = options.system == "runtime-optimizations";
        const bool applicationGroup = options.system == "application-optimizations";
        if (applicationGroup || options.system == "all" || options.system == "geometry-inspection") {
            const int result = RunGeometryInspection(options, f, metadata, summary, timings, checks);
            if (result != 0) return result;
            if (options.system == "geometry-inspection") {
                std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
                return 0;
            }
        }
        if (applicationGroup || options.system == "all" || options.system == "gravity-source-preparation") {
            const int result = RunGravitySourcePreparation(options, f, summary, timings, checks);
            if (result != 0) return result;
            if (options.system == "gravity-source-preparation") {
                std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
                return 0;
            }
        }
        if (runtimeGroup || options.system == "all" || options.system == "normalization-diagonals") {
            const int result = RunNormalizationDiagonals(options, f, summary, timings, checks);
            if (result != 0) return result;
            if (options.system == "normalization-diagonals") {
                summary << "checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
                std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
                return 0;
            }
        }
        if (runtimeGroup || options.system == "all" || options.system == "geometry-preparation") {
            const int result = RunGeometryPreparation(options, f, metadata, summary, timings, checks);
            if (result != 0) return result;
            if (options.system == "geometry-preparation") {
                summary << "checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
                std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
                return 0;
            }
        }

        if (IsLinearDerivativeSystem(options.system) && !applicationGroup) {
            int result = RunLinearDerivatives(options, f, metadata, summary, timings, checks);
            if (result == 0 && options.system == "preparation-optimizations") {
                auto hydroOptions = options;
                hydroOptions.system = "hydrostatic-displacement";
                result = RunStellarDerivatives(hydroOptions, f, metadata, summary, timings, checks);
            }
            std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
            return result;
        }
        if (!runtimeGroup && !applicationGroup && options.system != "hdiv-displacement" && options.system != "hdiv-mass" && options.system != "all") {
            const int result = RunStellarDerivatives(options, f, metadata, summary, timings, checks);
            std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
            return result;
        }

        const auto constructStart = Clock::now();
        operators::PreparedMappedHDivMassOperator mass(f, *f.domainMapperStateless);
        summary << "operator_construction_seconds=" << std::chrono::duration<double>(Clock::now() - constructStart).count() << '\n';
        const mfem::Vector base = Displacement(f, options.deformation);
        std::vector<mfem::Vector> fluxes, directions, reducedFluxes;
        for (int index = 0; index < options.directions; ++index) {
            fluxes.push_back(Flux(f.gravityFluxFes->GetTrueVSize(), 0.31 + 0.47 * index));
            directions.push_back(Displacement(f, 0.43, 0.2 + 0.7 * index));
            if (index == 1) directions.back().Add(0.002, Flux(f.displacementFes->GetTrueVSize(), 0.37));
            reducedFluxes.push_back(mass.GetFluxMap().gather(fluxes.back()));
        }
        const bool derivative = options.system != "hdiv-mass" && !applicationGroup;
        auto prepare = [&](const mfem::Vector &state) {
            const auto reduced = mass.GetDisplacementMap().gather(state);
            if (derivative) mass.Prepare(reduced); else mass.PreparePrimal(reduced);
        };
        auto start = Clock::now();
        prepare(base);
        summary << "initial_prepare_seconds=" << std::chrono::duration<double>(Clock::now() - start).count() << '\n';
        metadata << "fast_mass_elements=" << mass.GetFastMassElementCount() << '\n';
        metadata << "fast_variation_elements=" << mass.GetFastVariationElementCount()
                 << "\ntensor_variation_elements=" << mass.GetTensorVariationElementCount() << '\n';
        std::size_t quadraturePoints = 0;
        mass.VisitMappedGeometryRules([&](int, const mfem::IntegrationRule &rule) { quadraturePoints += rule.GetNPoints(); });
        metadata << "variation_quadrature_points=" << quadraturePoints << '\n';
        metadata.flush();

        const Apply fastDerivative = [&](int index, mfem::Vector &out) {
            mass.MultDisplacementVariationTrue(fluxes[index], directions[index], out);
        };
        const Apply referenceDerivative = [&](int index, mfem::Vector &out) {
            mass.MultDisplacementVariationTrueCachedReference(fluxes[index], directions[index], out);
        };
        const Apply fastMass = [&](int index, mfem::Vector &out) { mass.Mult(reducedFluxes[index], out); };
        const Apply referenceMass = [&](int index, mfem::Vector &out) { mass.MultReference(reducedFluxes[index], out); };

        if (options.check != "none") {
            const auto validationStart = Clock::now();
            // Reprepare distinct geometries to expose stale geometry caches.
            const std::vector<double> scales = options.check == "full" ? std::vector<double>{0.0, options.deformation}
                                                                       : std::vector<double>{options.deformation};
            for (std::size_t state = 0; state < scales.size(); ++state) {
                const auto displacement = Displacement(f, scales[state]);
                prepare(displacement);
                for (int index = 0; index < options.directions; ++index) {
                    mfem::Vector actual, expected;
                    const std::string label = "state" + std::to_string(state) + "_direction" + std::to_string(index);
                    if (derivative) {
                        fastDerivative(index, actual);
                        referenceDerivative(index, expected);
                        checks.Add(label + "_cached_derivative", RelativeError(actual, expected), 2.0e-11);
                        mass.MultDisplacementVariationTrueReference(fluxes[index], directions[index], expected);
                        checks.Add(label + "_general_prepared_derivative", RelativeError(actual, expected), 2.0e-11);
                        if (options.check == "full" && index == 0) {
                            operators::kernels::apply_mapped_hdiv_mass_variation(
                                f, *f.domainMapperStateless, fluxes[index], displacement, directions[index], expected);
                            checks.Add(label + "_stateless_derivative", RelativeError(actual, expected), 2.0e-11);
                        }
                    }
                    if (options.system != "hdiv-displacement" && !runtimeGroup) {
                        fastMass(index, actual);
                        referenceMass(index, expected);
                        checks.Add(label + "_mfem_pa_mass", RelativeError(actual, expected), 2.0e-11);
                        if (options.check == "full" && index == 0) {
                            mfem::Vector full;
                            operators::kernels::apply_mapped_hdiv_mass(f, *f.domainMapperStateless, fluxes[index], displacement, full);
                            expected = mass.GetFluxMap().gather(full);
                            checks.Add(label + "_stateless_mass", RelativeError(actual, expected), 2.0e-11);
                        }
                    }
                }
                if (options.system != "hdiv-displacement" && !runtimeGroup) {
                    const std::string label = "state" + std::to_string(state) + "_mass";
                    mfem::Vector other = mass.GetFluxMap().gather(
                        Flux(f.gravityFluxFes->GetTrueVSize(), 1.37));
                    mfem::Vector firstAction, otherAction, combinedAction, expected;
                    mass.Mult(reducedFluxes[0], firstAction);
                    mass.Mult(other, otherAction);
                    const double xy = reducedFluxes[0] * otherAction;
                    const double yx = other * firstAction;
                    checks.Add(label + "_symmetry", std::abs(xy-yx) / std::max({1.0, std::abs(xy), std::abs(yx)}), 2.0e-11);
                    checks.Add(label + "_positive_energy", reducedFluxes[0] * firstAction > 0.0 ? 0.0 : 1.0, 0.0);
                    mfem::Vector combined(reducedFluxes[0]);
                    combined *= 0.7;
                    combined.Add(-0.4, other);
                    mass.Mult(combined, combinedAction);
                    expected = firstAction;
                    expected *= 0.7;
                    expected.Add(-0.4, otherAction);
                    checks.Add(label + "_linearity", RelativeError(combinedAction, expected), 2.0e-11);
                    combined = 0.0;
                    mass.Mult(combined, combinedAction);
                    checks.Add(label + "_zero", combinedAction.Norml2(), 0.0);
                    mass.MultReference(reducedFluxes[0], expected);
                    mass.Mult(reducedFluxes[0], combinedAction);
                    checks.Add(label + "_reference_restores_fast_path", RelativeError(combinedAction, firstAction), 0.0);
                }
            }
            if (derivative && options.check == "full") {
                mfem::Vector analyticTrue;
                fastDerivative(0, analyticTrue);
                const auto analytic = mass.GetFluxMap().gather(analyticTrue);
                operators::PreparedMappedHDivMassOperator differenceMass(f, *f.domainMapperStateless);
                for (const double step : {options.finiteDifferenceStep, 0.5 * options.finiteDifferenceStep}) {
                    mfem::Vector plus(base), minus(base), plusAction, minusAction;
                    plus.Add(step, directions[0]);
                    minus.Add(-step, directions[0]);
                    differenceMass.PreparePrimal(differenceMass.GetDisplacementMap().gather(plus));
                    differenceMass.Mult(reducedFluxes[0], plusAction);
                    differenceMass.PreparePrimal(differenceMass.GetDisplacementMap().gather(minus));
                    differenceMass.Mult(reducedFluxes[0], minusAction);
                    plusAction -= minusAction;
                    plusAction /= 2.0 * step;
                    checks.Add("finite_difference_" + std::to_string(step), RelativeError(analytic, plusAction), 2.0e-6);
                }
            }
            summary << "validation_seconds=" << std::chrono::duration<double>(Clock::now() - validationStart).count() << '\n';
            if (!checks.passed) {
                summary << "checks_passed=0\nbenchmark_skipped=1\n";
                return 3;
            }
        }

        start = Clock::now();
        prepare(base);
        summary << "repeat_prepare_seconds=" << std::chrono::duration<double>(Clock::now() - start).count() << '\n';
        summary << "checks=" << (options.check == "none" ? "not_run" : "passed") << '\n';
        summary.flush();
        if (derivative) Benchmark(options, "hdiv-displacement", "cached_pointwise_reference", fastDerivative, referenceDerivative, timings, summary);
        if (options.system != "hdiv-displacement" && !runtimeGroup) Benchmark(options, "hdiv-mass", "mfem_pa_reference", fastMass, referenceMass, timings, summary);
        if (applicationGroup) {
            auto linearOptions = options;
            linearOptions.prepareSamples = 0; // preparation measurements above isolate this batch
            const int result = RunLinearDerivatives(linearOptions, f, metadata, summary, timings, checks);
            if (result != 0) return result;
        }
        if (runtimeGroup) {
            auto gravityOptions = options;
            gravityOptions.system = "gravity-force";
            gravityOptions.prepareSamples = 0; // this group measures new application and geometry changes
            const int result = RunLinearDerivatives(gravityOptions, f, metadata, summary, timings, checks);
            if (result != 0) return result;
        }
        if (options.system == "all") {
            const int result = RunStellarDerivatives(options, f, metadata, summary, timings, checks);
            if (result != 0) return result;
            const int linearResult = RunLinearDerivatives(options, f, metadata, summary, timings, checks);
            if (linearResult != 0) return linearResult;
        }
        std::cout << "Results: " << std::filesystem::absolute(options.output) << std::endl;
        return 0;
    }
}

int main(int argc, char **argv) {
    // Parse before MPI initialization so --help does not require a runtime.
    Options options;
    try {
        options = Parse(argc, argv);
        if (options.help) return Run(options);
    } catch (const std::exception &error) {
        std::cerr << error.what() << '\n';
        return 2;
    }
    mfem::Mpi::Init(argc, argv);
    int result = 0;
    try {
        int ranks = 0;
        MPI_Comm_size(MPI_COMM_WORLD, &ranks);
        if (ranks != 1) throw std::invalid_argument("This harness currently requires exactly one MPI rank.");
        mfem::Device device("cpu");
        std::cout << std::setprecision(12);
        result = Run(options);
    } catch (const std::exception &error) {
        std::cerr << "Isolation harness error: " << error.what() << std::endl;
        result = 2;
    }
    mfem::Mpi::Finalize();
    return result;
}