#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <exception>
#include <iostream>
#include <limits>
#include <map>
#include <numbers>
#include <ranges>
#include <span>
#include <string>
#include <string_view>
#include <utility>
#include <vector>

#include <catch2/catch_test_macros.hpp>
#include <mfem.hpp>
#include <mpi.h>

import experiment;
import mean_field;
import test_helpers;

namespace {
    using Clock               = std::chrono::steady_clock;

    namespace backend         = mean_field::preconditioning::backend;
    namespace preconditioning = mean_field::preconditioning;
    namespace solver          = mean_field::solver;

    struct CandidateDescription final {
        std::string name;
        std::string materialFactorization;
        std::string surfaceCalibration;
        std::string stellarStructureFactorization;
        std::string gravityFactorization;
        std::string specificationBorder;
    };

    struct SetupTimings final {
        double finiteElementSeconds{0.0};
        double laneEmdenCalibrationSeconds{0.0};
        double problemConstructionSeconds{0.0};
        double seedProjectionSeconds{0.0};
        double operatorPreparationSeconds{0.0};
        double manufacturedRightHandSideSeconds{0.0};
    };

    struct NewtonCandidateResult final {
        CandidateDescription candidate;
        mfem::Vector correction;
        mfem::Vector linearAction;
    };

    [[nodiscard]] const char *buildConfiguration() noexcept {
#ifdef NDEBUG
        return "release";
#else
        return "debug";
#endif
    }

    [[nodiscard]] double maximumRankSeconds(
        const Clock::time_point start,
        const MPI_Comm communicator
    ) {
        const double localSeconds = std::chrono::duration<double>(Clock::now() - start).count();
        double maximumSeconds     = 0.0;
        MPI_Allreduce(&localSeconds, &maximumSeconds, 1, MPI_DOUBLE, MPI_MAX, communicator);
        return maximumSeconds;
    }

    [[nodiscard]] double globalNorm(
        const mfem::Vector &vector,
        const MPI_Comm communicator
    ) {
        const double localSquared = vector * vector;
        double globalSquared      = 0.0;
        MPI_Allreduce(&localSquared, &globalSquared, 1, MPI_DOUBLE, MPI_SUM, communicator);
        return std::sqrt(std::max(globalSquared, 0.0));
    }

    void announce(
        const MPI_Comm communicator,
        const std::string &message
    ) {
        int rank = 0;
        MPI_Comm_rank(communicator, &rank);
        if (rank == 0) {
            std::cout << "[P10 full stellar] " << message << std::endl;
        }
    }

    class ArnoldiProgressOperator final : public mfem::Operator {
    public:
        ArnoldiProgressOperator(
            const mfem::Operator &operation,
            const MPI_Comm communicator,
            const int expectedApplications
        )
            : mfem::Operator(
                  operation.Height(),
                  operation.Width()
              ),
              m_operation(&operation),
              m_communicator(communicator),
              m_expectedApplications(expectedApplications) {
        }

        void Mult(
            const mfem::Vector &input,
            mfem::Vector &output
        ) const override {
            m_operation->Mult(input, output);
            ++m_completedApplications;
            if (m_completedApplications == 1 || m_completedApplications == m_expectedApplications ||
                m_completedApplications % 4 == 0) {
                announce(
                    m_communicator, "Arnoldi progress " + std::to_string(m_completedApplications) + "/" +
                                        std::to_string(m_expectedApplications)
                );
            }
        }

    private:
        const mfem::Operator *m_operation;
        MPI_Comm m_communicator;
        int m_expectedApplications;
        mutable int m_completedApplications{0};
    };

    [[nodiscard]] mean_field::operators::StellarEquilibriumDependencies makeDependencies() {
        return {
            .discretization     = {.identity = 111103, .revision = 1},
            .density            = {.identity = 111109, .revision = 1},
            .surfaceDeformation = {.identity = 111119, .revision = 1},
            .gravityGradient    = {.identity = 111121, .revision = 1},
            .gravityPotential   = {.identity = 111127, .revision = 1},
            .enthalpy           = {.identity = 111143, .revision = 1},
            .bernoulliConstant  = {.identity = 111149, .revision = 1},
            .rotation           = {.identity = 111151, .revision = 1},
            .targetMass         = {.identity = 111157, .revision = 1}
        };
    }

    [[nodiscard]] mean_field::physics::RigidRotation zeroRotation() {
        mfem::Vector angularVelocity(3);
        mfem::Vector center(3);
        angularVelocity = 0.0;
        center          = 0.0;
        return {angularVelocity, center};
    }

    [[nodiscard]] mfem::Vector blockBalancedDirection(
        const int size,
        const std::span<const mean_field::operators::RootBlockDescriptor> blocks,
        const double phase,
        const MPI_Comm communicator
    ) {
        mfem::Vector direction(size);
        direction = 0.0;
        for (const auto &block : blocks) {
            REQUIRE(block.offset >= 0);
            REQUIRE(block.size > 0);
            REQUIRE(block.offset + block.size <= size);
            mfem::Vector values(direction.GetData() + block.offset, block.size);
            for (int index = 0; index < values.Size(); ++index) {
                const double ordinal = static_cast<double>(index + 1);
                values(index) =
                    std::sin(0.6180339887498948 * ordinal + phase + static_cast<double>(block.canonicalIndex + 1)) +
                    0.31 * std::cos(0.1732050807568877 * ordinal - phase);
            }
            const double norm = globalNorm(values, communicator);
            REQUIRE(norm > 0.0);
            values /= norm;
        }
        return direction;
    }

    [[nodiscard]] std::map<
        std::string,
        std::string>
    commonParameters(
        const CandidateDescription &candidate,
        const std::string &measurement,
        const int dimension
    ) {
        return {
            {"build_configuration", buildConfiguration()},
            {"candidate", candidate.name},
            {"equation_of_state", "Polytrope(n=3)"},
            {"experiment_schema", "p9_p10_full_stellar_v1"},
            {"gravity_factorization", candidate.gravityFactorization},
            {"linearization_state", "projected_lane_emden"},
            {"manufactured_rhs", "J_times_value_block_balanced_correction"},
            {"material_factorization", candidate.materialFactorization},
            {"measurement", measurement},
            {"mesh_file", test_utils::setup_args().mesh_file},
            {"operator", "canonical_bordered_stellar_jacobian"},
            {"preconditioned_product", "J M^-1"},
            {"residual_arnoldi_seed", "residual_block_balanced"},
            {"root_dimension", std::to_string(dimension)},
            {"rotation", "zero"},
            {"specification_border", candidate.specificationBorder},
            {"stellar_structure_factorization", candidate.stellarStructureFactorization},
            {"surface_calibration", candidate.surfaceCalibration}
        };
    }

    [[nodiscard]] std::map<
        std::string,
        std::string>
    newtonParameters(
        const CandidateDescription &candidate,
        const std::string &measurement,
        const int dimension
    ) {
        auto parameters                      = commonParameters(candidate, measurement, dimension);
        parameters["experiment_schema"]      = "p10_physical_newton_rhs_v1";
        parameters["manufactured_rhs"]       = "none";
        parameters["right_hand_side"]        = "negative_nonlinear_residual";
        parameters["preconditioned_product"] = "not_measured";
        parameters["residual_arnoldi_seed"]  = "not_applicable";
        return parameters;
    }

    void incrementStateRevisions(mean_field::operators::StellarEquilibriumDependencies &dependencies) {
        ++dependencies.density.revision;
        ++dependencies.surfaceDeformation.revision;
        ++dependencies.gravityGradient.revision;
        ++dependencies.gravityPotential.revision;
        ++dependencies.enthalpy.revision;
        ++dependencies.bernoulliConstant.revision;
        ++dependencies.rotation.revision;
        ++dependencies.targetMass.revision;
    }

    [[nodiscard]] const mean_field::operators::RootBlockDescriptor &findBlock(
        const std::span<const mean_field::operators::RootBlockDescriptor> blocks,
        const std::string_view stableId
    ) {
        const auto iterator =
            std::ranges::find(blocks, stableId, &mean_field::operators::RootBlockDescriptor::stableId);
        REQUIRE(iterator != blocks.end());
        return *iterator;
    }

    [[nodiscard]] int firstThresholdIteration(
        const std::vector<solver::IterationResidualMeasurement> &history,
        const double initialNorm,
        const double threshold
    ) {
        if (!std::isfinite(initialNorm) || initialNorm <= 0.0) {
            return -1;
        }
        for (const auto &sample : history) {
            if (std::abs(sample.reportedNorm) / initialNorm <= threshold) {
                return sample.iteration;
            }
        }
        return -1;
    }

    void recordSpectrum(
        const CandidateDescription &candidate,
        const solver::ArnoldiSpectralMeasurement &spectrum,
        const int dimension,
        const double setupSeconds,
        const solver::OperatorApplicationStatistics &jacobianStatistics,
        const solver::OperatorApplicationStatistics &preconditionerStatistics
    ) {
        experiment::record_experiment_result(
            "stellar_preconditioning_p10", candidate.name + "_arnoldi_summary",
            commonParameters(candidate, "arnoldi_summary", dimension),
            {{"preconditioner_setup_seconds_maximum_rank", setupSeconds},
             {"requested_dimension", static_cast<double>(spectrum.requestedDimension)},
             {"achieved_dimension", static_cast<double>(spectrum.achievedDimension)},
             {"invariant_subspace_found", spectrum.invariantSubspaceFound ? 1.0 : 0.0},
             {"operator_applications", static_cast<double>(spectrum.operatorApplications)},
             {"measurement_seconds_maximum_rank", spectrum.measurementSecondsMaximumRank},
             {"operator_application_seconds_maximum_rank", spectrum.operatorApplicationSecondsMaximumRank},
             {"operator_maximum_application_seconds_maximum_rank",
              spectrum.operatorMaximumApplicationSecondsMaximumRank},
             {"nonapplication_seconds_maximum_rank", spectrum.nonApplicationSecondsMaximumRank},
             {"measured_jacobian_applications", static_cast<double>(jacobianStatistics.applications)},
             {"measured_jacobian_application_seconds", jacobianStatistics.totalSeconds},
             {"measured_jacobian_maximum_application_seconds", jacobianStatistics.maximumSeconds},
             {"measured_preconditioner_applications", static_cast<double>(preconditionerStatistics.applications)},
             {"measured_preconditioner_application_seconds", preconditionerStatistics.totalSeconds},
             {"measured_preconditioner_maximum_application_seconds", preconditionerStatistics.maximumSeconds},
             {"projected_condition_proxy", spectrum.projectedConditionProxy},
             {"projected_largest_singular_value", spectrum.projectedLargestSingularValue},
             {"projected_smallest_singular_value", spectrum.projectedSmallestSingularValue},
             {"centroid_real_part", spectrum.centroidRealPart},
             {"centroid_imaginary_part", spectrum.centroidImaginaryPart},
             {"rms_distance_from_one", spectrum.rmsDistanceFromOne},
             {"rms_cluster_radius", spectrum.rmsClusterRadius},
             {"minimum_magnitude", spectrum.minimumMagnitude},
             {"maximum_magnitude", spectrum.maximumMagnitude},
             {"minimum_real_part", spectrum.minimumRealPart},
             {"maximum_real_part", spectrum.maximumRealPart},
             {"maximum_absolute_imaginary_part", spectrum.maximumAbsoluteImaginaryPart},
             {"negative_real_part_count", static_cast<double>(spectrum.negativeRealPartCount)},
             {"converged_ritz_value_count", static_cast<double>(spectrum.convergedRitzValueCount)},
             {"conjugate_pair_defect", spectrum.conjugatePairDefect},
             {"projected_departure_from_normality", spectrum.projectedDepartureFromNormality},
             {"field_of_values_minimum_real_part", spectrum.projectedFieldOfValuesMinimumRealPart},
             {"field_of_values_maximum_real_part", spectrum.projectedFieldOfValuesMaximumRealPart}}
        );

        std::vector<solver::RitzValueMeasurement> ordered = spectrum.ritzValues;
        std::ranges::sort(ordered, [](const auto &left, const auto &right) {
            if (left.realPart != right.realPart) {
                return left.realPart < right.realPart;
            }
            return left.imaginaryPart < right.imaginaryPart;
        });
        for (std::size_t index = 0; index < ordered.size(); ++index) {
            const auto &value = ordered[index];
            experiment::record_experiment_result(
                "stellar_preconditioning_p10", candidate.name + "_ritz_" + std::to_string(index),
                commonParameters(candidate, "ritz_value", dimension),
                {{"ritz_index", static_cast<double>(index)},
                 {"real_part", value.realPart},
                 {"imaginary_part", value.imaginaryPart},
                 {"magnitude", value.magnitude},
                 {"distance_from_one", value.distanceFromOne},
                 {"residual_estimate", value.residualEstimate},
                 {"relative_residual_estimate", value.relativeResidualEstimate},
                 {"converged", value.converged ? 1.0 : 0.0}}
            );
        }
    }

    template <
        typename Preconditioner,
        typename Problem>
    void measureCandidate(
        const CandidateDescription &candidate,
        Preconditioner &inversePreconditioner,
        const double preconditionerSetupSeconds,
        const Problem &problem,
        const mfem::Vector &exactCorrection,
        const mfem::Vector &rightHandSide,
        const mfem::Vector &arnoldiDirection,
        const SetupTimings &setupTimings,
        const MPI_Comm communicator
    ) {
        constexpr int maximumIterations = 36;
        constexpr int restartDimension  = 18;
        constexpr int arnoldiDimension  = 12;

        const mfem::Operator &jacobian  = problem.GetLinearizationOperator();
        solver::InstrumentedOperator instrumentedJacobian(jacobian);
        solver::InstrumentedPreconditioner instrumentedPreconditioner(inversePreconditioner);
        solver::ResidualHistoryMonitor monitor;
        mfem::FGMRESSolver krylov(communicator);
        krylov.SetPreconditioner(instrumentedPreconditioner);
        krylov.SetOperator(instrumentedJacobian);
        krylov.SetMonitor(monitor);
        krylov.SetRelTol(1.0e-8);
        krylov.SetAbsTol(1.0e-12);
        krylov.SetMaxIter(maximumIterations);
        krylov.SetKDim(restartDimension);
        krylov.SetPrintLevel(0);

        mfem::Vector solution(problem.StateSize());
        solution = 0.0;
        announce(communicator, "solving the manufactured system with " + candidate.name);
        const Clock::time_point solveStart = Clock::now();
        krylov.Mult(rightHandSide, solution);
        const double localSolveSeconds             = std::chrono::duration<double>(Clock::now() - solveStart).count();

        const solver::LinearSolveMeasurement solve = solver::measureLinearSolve(
            krylov, jacobian, rightHandSide, solution, problem.GetManifest().residualBlocks(),
            instrumentedJacobian.GetStatistics(), instrumentedPreconditioner.GetStatistics(),
            instrumentedPreconditioner.GetLifecycleStatistics(), monitor, localSolveSeconds, communicator
        );
        mfem::Vector correctionError(solution);
        correctionError -= exactCorrection;
        const double relativeCorrectionError =
            globalNorm(correctionError, communicator) /
            std::max(globalNorm(exactCorrection, communicator), std::numeric_limits<double>::min());

        const double reportedInitial = std::max(std::abs(krylov.GetInitialNorm()), 1.0e-300);
        const int iteration1e2       = firstThresholdIteration(solve.reportedResidualHistory, reportedInitial, 1.0e-2);
        const int iteration1e4       = firstThresholdIteration(solve.reportedResidualHistory, reportedInitial, 1.0e-4);
        const int iteration1e6       = firstThresholdIteration(solve.reportedResidualHistory, reportedInitial, 1.0e-6);
        const int iteration1e8       = firstThresholdIteration(solve.reportedResidualHistory, reportedInitial, 1.0e-8);

        REQUIRE(std::isfinite(solve.directResidual.relativeResidual));
        REQUIRE(std::isfinite(relativeCorrectionError));
        std::map<std::string, double> metrics{
            {"maximum_iterations", static_cast<double>(maximumIterations)},
            {"restart_dimension", static_cast<double>(restartDimension)},
            {"solver_converged", solve.solverConverged ? 1.0 : 0.0},
            {"outer_iterations", static_cast<double>(solve.outerIterations)},
            {"reported_initial_residual_norm", solve.solverReportedInitialNorm},
            {"reported_final_residual_norm", solve.solverReportedFinalNorm},
            {"reported_residual_reduction", solve.solverReportedResidualReduction},
            {"reported_iteration_to_1e-2", static_cast<double>(iteration1e2)},
            {"reported_iteration_to_1e-4", static_cast<double>(iteration1e4)},
            {"reported_iteration_to_1e-6", static_cast<double>(iteration1e6)},
            {"reported_iteration_to_1e-8", static_cast<double>(iteration1e8)},
            {"rhs_norm", solve.directResidual.rightHandSideNorm},
            {"true_residual_norm", solve.directResidual.trueResidualNorm},
            {"true_relative_residual", solve.directResidual.relativeResidual},
            {"relative_correction_error", relativeCorrectionError},
            {"true_residual_digits_per_jacobian_application", solve.trueResidualDigitsReducedPerJacobianApplication},
            {"solve_seconds_maximum_rank", solve.solveSecondsMaximumRank},
            {"jacobian_applications", static_cast<double>(solve.jacobian.applications)},
            {"jacobian_application_seconds", solve.jacobian.totalSeconds},
            {"jacobian_maximum_application_seconds", solve.jacobian.maximumSeconds},
            {"preconditioner_applications", static_cast<double>(solve.inversePreconditioner.applications)},
            {"preconditioner_application_seconds", solve.inversePreconditioner.totalSeconds},
            {"preconditioner_maximum_application_seconds", solve.inversePreconditioner.maximumSeconds},
            {"preconditioner_setups", static_cast<double>(solve.inversePreconditionerLifecycle.setups)},
            {"preconditioner_setup_seconds_in_solver", solve.inversePreconditionerLifecycle.setupSeconds},
            {"preconditioner_setup_seconds_maximum_rank", preconditionerSetupSeconds},
            {"finite_element_setup_seconds", setupTimings.finiteElementSeconds},
            {"lane_emden_calibration_seconds", setupTimings.laneEmdenCalibrationSeconds},
            {"problem_construction_seconds", setupTimings.problemConstructionSeconds},
            {"seed_projection_seconds", setupTimings.seedProjectionSeconds},
            {"operator_preparation_seconds", setupTimings.operatorPreparationSeconds},
            {"manufactured_rhs_seconds", setupTimings.manufacturedRightHandSideSeconds}
        };
        for (const auto &block : solve.directResidual.blocks) {
            const std::string prefix                              = "residual_block." + block.stableId;
            metrics[prefix + ".size"]                             = static_cast<double>(block.size);
            metrics[prefix + ".descriptor_scale"]                 = block.descriptorScale;
            metrics[prefix + ".rhs_norm"]                         = block.rightHandSideNorm;
            metrics[prefix + ".true_norm"]                        = block.trueResidualNorm;
            metrics[prefix + ".block_relative_residual"]          = block.blockRelativeResidual;
            metrics[prefix + ".scaled_rhs_norm"]                  = block.scaledRightHandSideNorm;
            metrics[prefix + ".scaled_true_norm"]                 = block.scaledTrueResidualNorm;
            metrics[prefix + ".global_relative_contribution"]     = block.contributionToGlobalRelativeResidual;
            metrics[prefix + ".fraction_global_squared_residual"] = block.fractionOfGlobalSquaredResidualNorm;
        }
        const double correctionErrorNorm = globalNorm(correctionError, communicator);
        for (const auto &block : problem.GetManifest().valueBlocks()) {
            const mfem::Vector exactBlock(
                const_cast<mfem::real_t *>(exactCorrection.GetData()) + block.offset, block.size
            );
            const mfem::Vector errorBlock(correctionError.GetData() + block.offset, block.size);
            const double exactBlockNorm           = globalNorm(exactBlock, communicator);
            const double errorBlockNorm           = globalNorm(errorBlock, communicator);
            const std::string prefix              = "correction_block." + std::string(block.stableId);
            metrics[prefix + ".size"]             = static_cast<double>(block.size);
            metrics[prefix + ".descriptor_scale"] = block.scale;
            metrics[prefix + ".exact_norm"]       = exactBlockNorm;
            metrics[prefix + ".error_norm"]       = errorBlockNorm;
            metrics[prefix + ".block_relative_error"] =
                errorBlockNorm / std::max(exactBlockNorm, std::numeric_limits<double>::min());
            metrics[prefix + ".scaled_exact_norm"] = exactBlockNorm / block.scale;
            metrics[prefix + ".scaled_error_norm"] = errorBlockNorm / block.scale;
            metrics[prefix + ".fraction_global_squared_error"] =
                correctionErrorNorm > 0.0
                    ? errorBlockNorm * errorBlockNorm / (correctionErrorNorm * correctionErrorNorm)
                    : 0.0;
        }
        experiment::record_experiment_result(
            "stellar_preconditioning_p10", candidate.name + "_linear_solve",
            commonParameters(candidate, "manufactured_linear_solve", problem.StateSize()), std::move(metrics)
        );

        for (std::size_t index = 0; index < solve.reportedResidualHistory.size(); ++index) {
            const auto &sample = solve.reportedResidualHistory[index];
            experiment::record_experiment_result(
                "stellar_preconditioning_p10", candidate.name + "_history_" + std::to_string(index),
                commonParameters(candidate, "fgmres_residual_history", problem.StateSize()),
                {{"history_sample", static_cast<double>(index)},
                 {"iteration", static_cast<double>(sample.iteration)},
                 {"reported_residual_norm", sample.reportedNorm},
                 {"reported_relative_residual", std::abs(sample.reportedNorm) / reportedInitial},
                 {"final_measurement", sample.final ? 1.0 : 0.0}}
            );
        }

        instrumentedJacobian.ResetStatistics();
        instrumentedPreconditioner.ResetStatistics();
        solver::FixedRightPreconditionedOperator product(instrumentedJacobian, instrumentedPreconditioner);
        ArnoldiProgressOperator progress(product, communicator, arnoldiDimension);
        announce(communicator, "measuring " + candidate.name + " with 12-vector Arnoldi");
        const solver::ArnoldiSpectralMeasurement spectrum = solver::measureArnoldiSpectrum(
            progress, arnoldiDirection, communicator,
            {.krylovDimension                  = arnoldiDimension,
             .breakdownRelativeTolerance       = 1.0e-13,
             .ritzConvergenceRelativeTolerance = 1.0e-7,
             .reorthogonalize                  = true}
        );
        REQUIRE(spectrum.achievedDimension > 0);
        recordSpectrum(
            candidate, spectrum, problem.StateSize(), preconditionerSetupSeconds, instrumentedJacobian.GetStatistics(),
            instrumentedPreconditioner.GetStatistics()
        );

        int rank = 0;
        MPI_Comm_rank(communicator, &rank);
        if (rank == 0) {
            std::cout << "[P10 full stellar] " << candidate.name << ": iterations=" << solve.outerIterations
                      << ", converged=" << (solve.solverConverged ? "yes" : "no")
                      << ", true residual=" << solve.directResidual.relativeResidual
                      << ", correction error=" << relativeCorrectionError
                      << ", projected condition=" << spectrum.projectedConditionProxy << '\n';
        }
    }

    template <typename Problem>
    void prepareAndMeasureDefault(
        const CandidateDescription &candidate,
        const Problem &problem,
        const mfem::Vector &exactCorrection,
        const mfem::Vector &rightHandSide,
        const mfem::Vector &arnoldiDirection,
        const SetupTimings &setupTimings,
        const MPI_Comm communicator
    ) {
        const Clock::time_point setupStart = Clock::now();
        auto block                         = preconditioning::makePreconditioner(problem);
        auto prepared                      = preconditioning::prepare(problem, std::move(block));
        const double setupSeconds          = maximumRankSeconds(setupStart, communicator);
        measureCandidate(
            candidate, prepared, setupSeconds, problem, exactCorrection, rightHandSide, arnoldiDirection, setupTimings,
            communicator
        );
    }

    template <
        preconditioning::MaterialSurfaceFactorizationPolicy MaterialPolicy,
        preconditioning::StellarStructureFactorizationPolicy StructurePolicy,
        typename Problem,
        backend::Registered GravityMassBackend = backend::Diagonal>
        requires backend::Compatible<
            GravityMassBackend,
            preconditioning::GravityMassInverseCharacteristics>
    void prepareAndMeasureComposed(
        const CandidateDescription &candidate,
        const Problem &problem,
        const MaterialPolicy materialPolicy,
        const StructurePolicy structurePolicy,
        const preconditioning::MaterialSurfaceDiagonalOptions materialOptions,
        const mfem::Vector &exactCorrection,
        const mfem::Vector &rightHandSide,
        const mfem::Vector &arnoldiDirection,
        const SetupTimings &setupTimings,
        const MPI_Comm communicator,
        GravityMassBackend gravityMassBackend = {},
        const int gravityAmgCycles            = 1
    ) {
        const Clock::time_point setupStart = Clock::now();
        auto material                      = preconditioning::materialSurfaceBlock(
            problem, backend::Diagonal{}, backend::Diagonal{}, materialPolicy, materialOptions
        );
        using FixedAMG = backend::HypreBoomerAMG<backend::FixedCycles>;
        auto gravity   = preconditioning::GravityFieldBlock(
            std::move(gravityMassBackend), FixedAMG{backend::FixedCycles{.cycles = gravityAmgCycles}},
            preconditioning::GravityApproximateLDU{}
        );
        auto structure =
            preconditioning::stellarStructureBlock(problem, std::move(material), std::move(gravity), structurePolicy);
        auto block = preconditioning::specificationBorderBlock(problem, std::move(structure), backend::DenseDirect{});
        auto prepared             = preconditioning::prepare(problem, std::move(block));
        const double setupSeconds = maximumRankSeconds(setupStart, communicator);
        measureCandidate(
            candidate, prepared, setupSeconds, problem, exactCorrection, rightHandSide, arnoldiDirection, setupTimings,
            communicator
        );
    }

    template <
        typename Preconditioner,
        typename Problem>
    [[nodiscard]] NewtonCandidateResult solveNewtonCandidate(
        const CandidateDescription &candidate,
        Preconditioner &inversePreconditioner,
        const double preconditionerSetupSeconds,
        const Problem &problem,
        const mfem::Vector &baseResidual,
        const mfem::Vector &rightHandSide,
        const MPI_Comm communicator
    ) {
        constexpr int maximumIterations = 48;
        constexpr int restartDimension  = 20;

        const mfem::Operator &jacobian  = problem.GetLinearizationOperator();
        solver::InstrumentedOperator instrumentedJacobian(jacobian);
        solver::InstrumentedPreconditioner instrumentedPreconditioner(inversePreconditioner);
        solver::ResidualHistoryMonitor monitor;
        mfem::FGMRESSolver krylov(communicator);
        krylov.SetPreconditioner(instrumentedPreconditioner);
        krylov.SetOperator(instrumentedJacobian);
        krylov.SetMonitor(monitor);
        krylov.SetRelTol(1.0e-8);
        krylov.SetAbsTol(1.0e-12);
        krylov.SetMaxIter(maximumIterations);
        krylov.SetKDim(restartDimension);
        krylov.SetPrintLevel(0);

        mfem::Vector correction(problem.StateSize());
        correction = 0.0;
        announce(communicator, "solving the physical Newton system with " + candidate.name);
        const Clock::time_point solveStart = Clock::now();
        krylov.Mult(rightHandSide, correction);
        const double localSolveSeconds             = std::chrono::duration<double>(Clock::now() - solveStart).count();

        const solver::LinearSolveMeasurement solve = solver::measureLinearSolve(
            krylov, jacobian, rightHandSide, correction, problem.GetManifest().residualBlocks(),
            instrumentedJacobian.GetStatistics(), instrumentedPreconditioner.GetStatistics(),
            instrumentedPreconditioner.GetLifecycleStatistics(), monitor, localSolveSeconds, communicator
        );
        mfem::Vector linearAction(problem.EquationSize());
        jacobian.Mult(correction, linearAction);
        mfem::Vector predictedResidual(baseResidual);
        predictedResidual += linearAction;

        const double baseResidualNorm      = globalNorm(baseResidual, communicator);
        const double correctionNorm        = globalNorm(correction, communicator);
        const double predictedResidualNorm = globalNorm(predictedResidual, communicator);
        REQUIRE(std::isfinite(solve.directResidual.relativeResidual));
        REQUIRE(std::isfinite(correctionNorm));
        REQUIRE(std::isfinite(predictedResidualNorm));

        std::map<std::string, double> metrics{
            {"maximum_iterations", static_cast<double>(maximumIterations)},
            {"restart_dimension", static_cast<double>(restartDimension)},
            {"solver_converged", solve.solverConverged ? 1.0 : 0.0},
            {"outer_iterations", static_cast<double>(solve.outerIterations)},
            {"reported_initial_residual_norm", solve.solverReportedInitialNorm},
            {"reported_final_residual_norm", solve.solverReportedFinalNorm},
            {"reported_residual_reduction", solve.solverReportedResidualReduction},
            {"base_nonlinear_residual_norm", baseResidualNorm},
            {"rhs_norm", solve.directResidual.rightHandSideNorm},
            {"true_linear_residual_norm", solve.directResidual.trueResidualNorm},
            {"true_linear_relative_residual", solve.directResidual.relativeResidual},
            {"predicted_full_step_residual_norm", predictedResidualNorm},
            {"predicted_full_step_relative_residual",
             predictedResidualNorm / std::max(baseResidualNorm, std::numeric_limits<double>::min())},
            {"correction_norm", correctionNorm},
            {"solve_seconds_maximum_rank", solve.solveSecondsMaximumRank},
            {"jacobian_applications", static_cast<double>(solve.jacobian.applications)},
            {"jacobian_application_seconds", solve.jacobian.totalSeconds},
            {"preconditioner_applications", static_cast<double>(solve.inversePreconditioner.applications)},
            {"preconditioner_application_seconds", solve.inversePreconditioner.totalSeconds},
            {"preconditioner_setup_seconds_maximum_rank", preconditionerSetupSeconds}
        };
        for (const auto &block : solve.directResidual.blocks) {
            const std::string prefix                          = "linear_residual_block." + block.stableId;
            metrics[prefix + ".rhs_norm"]                     = block.rightHandSideNorm;
            metrics[prefix + ".true_norm"]                    = block.trueResidualNorm;
            metrics[prefix + ".block_relative_residual"]      = block.blockRelativeResidual;
            metrics[prefix + ".global_relative_contribution"] = block.contributionToGlobalRelativeResidual;
        }
        for (const auto &block : problem.GetManifest().valueBlocks()) {
            const mfem::Vector correctionBlock(correction.GetData() + block.offset, block.size);
            const std::string prefix              = "correction_block." + std::string(block.stableId);
            metrics[prefix + ".norm"]             = globalNorm(correctionBlock, communicator);
            metrics[prefix + ".descriptor_scale"] = block.scale;
            metrics[prefix + ".scaled_norm"]      = metrics[prefix + ".norm"] / block.scale;
        }
        experiment::record_experiment_result(
            "stellar_preconditioning_p10_newton", candidate.name + "_linear_solve",
            newtonParameters(candidate, "physical_newton_linear_solve", problem.StateSize()), std::move(metrics)
        );

        const double reportedInitial = std::max(std::abs(krylov.GetInitialNorm()), 1.0e-300);
        for (std::size_t index = 0; index < solve.reportedResidualHistory.size(); ++index) {
            const auto &sample = solve.reportedResidualHistory[index];
            experiment::record_experiment_result(
                "stellar_preconditioning_p10_newton", candidate.name + "_history_" + std::to_string(index),
                newtonParameters(candidate, "fgmres_residual_history", problem.StateSize()),
                {{"history_sample", static_cast<double>(index)},
                 {"iteration", static_cast<double>(sample.iteration)},
                 {"reported_residual_norm", sample.reportedNorm},
                 {"reported_relative_residual", std::abs(sample.reportedNorm) / reportedInitial},
                 {"final_measurement", sample.final ? 1.0 : 0.0}}
            );
        }

        return {.candidate = candidate, .correction = std::move(correction), .linearAction = std::move(linearAction)};
    }

    template <
        preconditioning::StellarStructureFactorizationPolicy StructurePolicy,
        typename Problem>
    [[nodiscard]] NewtonCandidateResult prepareAndSolveNewtonComposed(
        const CandidateDescription &candidate,
        const Problem &problem,
        const StructurePolicy structurePolicy,
        const preconditioning::MaterialSurfaceDiagonalOptions materialOptions,
        const mfem::Vector &baseResidual,
        const mfem::Vector &rightHandSide,
        const MPI_Comm communicator
    ) {
        const Clock::time_point setupStart = Clock::now();
        auto material                      = preconditioning::materialSurfaceBlock(
            problem, backend::Diagonal{}, backend::Diagonal{}, preconditioning::SurfaceThenMaterialTriangular{},
            materialOptions
        );
        using FixedAMG = backend::HypreBoomerAMG<backend::FixedCycles>;
        auto gravity   = preconditioning::GravityFieldBlock(
            backend::Diagonal{}, FixedAMG{backend::FixedCycles{.cycles = 1}}, preconditioning::GravityApproximateLDU{}
        );
        auto structure =
            preconditioning::stellarStructureBlock(problem, std::move(material), std::move(gravity), structurePolicy);
        auto block = preconditioning::specificationBorderBlock(problem, std::move(structure), backend::DenseDirect{});
        auto prepared             = preconditioning::prepare(problem, std::move(block));
        const double setupSeconds = maximumRankSeconds(setupStart, communicator);
        return solveNewtonCandidate(
            candidate, prepared, setupSeconds, problem, baseResidual, rightHandSide, communicator
        );
    }
} // namespace

TEST_CASE(
    "Full Stellar P9 P10 Canonical Preconditioner Comparison",
    "[preconditioning][diagnostics][experiment][spectrum][p9][p10][full_system]"
) {
    using namespace mean_field;

    constexpr int radialSampleCount            = 4096;
    const Clock::time_point finiteElementStart = Clock::now();
    const utils::Args arguments                = test_utils::setup_args();
    fem::FEM finiteElements                    = fem::setup_fem(arguments.mesh_file, arguments, 0);
    REQUIRE(finiteElements.okay());
    const MPI_Comm communicator = finiteElements.mesh->GetComm();
    int communicatorSize        = 0;
    MPI_Comm_size(communicator, &communicatorSize);
    REQUIRE(communicatorSize == 1);

    SetupTimings setupTimings;
    setupTimings.finiteElementSeconds                  = maximumRankSeconds(finiteElementStart, communicator);

    constexpr double stellarRadius                     = utils::RADIUS;
    constexpr double targetMass                        = utils::MASS;
    const Clock::time_point calibrationStart           = Clock::now();
    const seed::DimensionlessLaneEmdenSolution profile = seed::integrateLaneEmden(3.0, 10.0);
    REQUIRE(profile.firstZeroCoordinate.has_value());
    const double surfaceCoordinate = *profile.firstZeroCoordinate;
    const double surfaceDerivative = profile.thetaDerivative(profile.thetaDerivative.Size() - 1);
    const double dimensionlessMass = -surfaceCoordinate * surfaceCoordinate * surfaceDerivative;
    REQUIRE(dimensionlessMass > 0.0);
    const double massScale          = targetMass / (4.0 * std::numbers::pi_v<double> * dimensionlessMass);
    const double polytropicConstant = std::numbers::pi_v<double> * utils::G * std::pow(massScale, 2.0 / 3.0);
    const double radialScale        = stellarRadius / surfaceCoordinate;
    const double centralDensity =
        std::pow(polytropicConstant / (std::numbers::pi_v<double> * utils::G * radialScale * radialScale), 1.5);
    setupTimings.laneEmdenCalibrationSeconds  = maximumRankSeconds(calibrationStart, communicator);

    const Clock::time_point constructionStart = Clock::now();
    auto model                                = model::StellarModel(
        eos::Polytrope({.n = 3.0, .K = polytropicConstant}),
        surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}),
        integral::FixedTotalMass({.Mtotal = dimensions::MassValue{targetMass}}),
        constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{centralDensity}})
    );
    auto problem                            = equilibrium::discretize(model, std::move(finiteElements));
    setupTimings.problemConstructionSeconds = maximumRankSeconds(constructionStart, communicator);

    announce(communicator, "projecting the n=3 Lane-Emden state");
    const Clock::time_point projectionStart = Clock::now();
    auto projected =
        seed::makeProjectedEquilibriumState(problem, seed::LaneEmden({.radialSampleCount = radialSampleCount}));
    setupTimings.seedProjectionSeconds = maximumRankSeconds(projectionStart, communicator);

    announce(communicator, "preparing the canonical bordered stellar Jacobian");
    const Clock::time_point preparationStart = Clock::now();
    const auto preparation                   = problem.Prepare(projected.values, makeDependencies(), zeroRotation());
    REQUIRE(preparation.assembledResidual);
    setupTimings.operatorPreparationSeconds = maximumRankSeconds(preparationStart, communicator);

    const mfem::Operator &jacobian          = problem.GetLinearizationOperator();
    const mfem::Vector exactCorrection =
        blockBalancedDirection(problem.StateSize(), problem.GetManifest().valueBlocks(), 0.23, communicator);
    mfem::Vector rightHandSide(problem.EquationSize());
    const Clock::time_point rightHandSideStart = Clock::now();
    jacobian.Mult(exactCorrection, rightHandSide);
    setupTimings.manufacturedRightHandSideSeconds = maximumRankSeconds(rightHandSideStart, communicator);
    REQUIRE(std::isfinite(globalNorm(rightHandSide, communicator)));
    const mfem::Vector arnoldiDirection =
        blockBalancedDirection(problem.EquationSize(), problem.GetManifest().residualBlocks(), 0.79, communicator);

    solver::IdentityPreconditioner identity(problem.StateSize());
    measureCandidate(
        {.name                          = "identity",
         .materialFactorization         = "identity",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "identity",
         .gravityFactorization          = "identity",
         .specificationBorder           = "identity"},
        identity, 0.0, problem, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator
    );

    prepareAndMeasureDefault(
        {.name                          = "current_default",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "independent_subsystems",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator
    );

    // P9 found that fitting the surface Riesz scale to A_qq gives the same
    // measured material-surface numerics as fitting the approximate Schur
    // complement while reducing calibration setup by roughly forty percent.
    // Keep both successful material factorizations explicit here: calibrated
    // surface-first is cheaper, while calibrated material LDU gave the best
    // isolated residual.
    constexpr preconditioning::MaterialSurfaceDiagonalOptions surfaceJacobianCalibratedOptions{
        .surfaceCalibration = {
            .target = preconditioning::SurfaceRieszCalibrationTarget::surface_jacobian, .probeCount = 4
        }
    };

    prepareAndMeasureComposed(
        {.name                          = "surface_then_material_aqq_calibrated_independent",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "surface_jacobian_4_probes",
         .stellarStructureFactorization = "independent_subsystems",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::SurfaceThenMaterialTriangular{}, preconditioning::IndependentStellarSubsystems{},
        surfaceJacobianCalibratedOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator
    );
    prepareAndMeasureComposed(
        {.name                          = "material_ldu_aqq_calibrated_independent",
         .materialFactorization         = "approximate_material_surface_ldu",
         .surfaceCalibration            = "surface_jacobian_4_probes",
         .stellarStructureFactorization = "independent_subsystems",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::ApproximateMaterialSurfaceLDU{}, preconditioning::IndependentStellarSubsystems{},
        surfaceJacobianCalibratedOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator
    );
    prepareAndMeasureComposed(
        {.name                          = "surface_then_material_aqq_calibrated_then_gravity",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "surface_jacobian_4_probes",
         .stellarStructureFactorization = "material_then_gravity_triangular",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::SurfaceThenMaterialTriangular{}, preconditioning::MaterialThenGravityTriangular{},
        surfaceJacobianCalibratedOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator
    );
    prepareAndMeasureComposed(
        {.name                          = "gravity_then_surface_then_material_aqq_calibrated",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "surface_jacobian_4_probes",
         .stellarStructureFactorization = "gravity_then_material_triangular",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::SurfaceThenMaterialTriangular{}, preconditioning::GravityThenMaterialTriangular{},
        surfaceJacobianCalibratedOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator
    );
    prepareAndMeasureComposed(
        {.name                          = "surface_then_material_aqq_calibrated_stellar_approximate_ldu",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "surface_jacobian_4_probes",
         .stellarStructureFactorization = "approximate_stellar_block_ldu",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::SurfaceThenMaterialTriangular{}, preconditioning::ApproximateStellarBlockLDU{},
        surfaceJacobianCalibratedOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator
    );
    prepareAndMeasureComposed(
        {.name                          = "material_ldu_aqq_calibrated_stellar_approximate_ldu",
         .materialFactorization         = "approximate_material_surface_ldu",
         .surfaceCalibration            = "surface_jacobian_4_probes",
         .stellarStructureFactorization = "approximate_stellar_block_ldu",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::ApproximateMaterialSurfaceLDU{}, preconditioning::ApproximateStellarBlockLDU{},
        surfaceJacobianCalibratedOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator
    );
}

TEST_CASE(
    "Full Stellar P10 Gravity Backend Finalists",
    "[preconditioning][diagnostics][experiment][p10][full_system][gravity_backend_finalists]"
) {
    using namespace mean_field;

    constexpr int radialSampleCount            = 4096;
    const Clock::time_point finiteElementStart = Clock::now();
    const utils::Args arguments                = test_utils::setup_args();
    fem::FEM finiteElements                    = fem::setup_fem(arguments.mesh_file, arguments, 0);
    REQUIRE(finiteElements.okay());
    const MPI_Comm communicator = finiteElements.mesh->GetComm();
    int communicatorSize        = 0;
    MPI_Comm_size(communicator, &communicatorSize);
    REQUIRE(communicatorSize == 1);

    SetupTimings setupTimings;
    setupTimings.finiteElementSeconds                  = maximumRankSeconds(finiteElementStart, communicator);

    constexpr double stellarRadius                     = utils::RADIUS;
    constexpr double targetMass                        = utils::MASS;
    const Clock::time_point calibrationStart           = Clock::now();
    const seed::DimensionlessLaneEmdenSolution profile = seed::integrateLaneEmden(3.0, 10.0);
    REQUIRE(profile.firstZeroCoordinate.has_value());
    const double surfaceCoordinate = *profile.firstZeroCoordinate;
    const double surfaceDerivative = profile.thetaDerivative(profile.thetaDerivative.Size() - 1);
    const double dimensionlessMass = -surfaceCoordinate * surfaceCoordinate * surfaceDerivative;
    REQUIRE(dimensionlessMass > 0.0);
    const double massScale          = targetMass / (4.0 * std::numbers::pi_v<double> * dimensionlessMass);
    const double polytropicConstant = std::numbers::pi_v<double> * utils::G * std::pow(massScale, 2.0 / 3.0);
    const double radialScale        = stellarRadius / surfaceCoordinate;
    const double centralDensity =
        std::pow(polytropicConstant / (std::numbers::pi_v<double> * utils::G * radialScale * radialScale), 1.5);
    setupTimings.laneEmdenCalibrationSeconds  = maximumRankSeconds(calibrationStart, communicator);

    const Clock::time_point constructionStart = Clock::now();
    auto model                                = model::StellarModel(
        eos::Polytrope({.n = 3.0, .K = polytropicConstant}),
        surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}),
        integral::FixedTotalMass({.Mtotal = dimensions::MassValue{targetMass}}),
        constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{centralDensity}})
    );
    auto problem                            = equilibrium::discretize(model, std::move(finiteElements));
    setupTimings.problemConstructionSeconds = maximumRankSeconds(constructionStart, communicator);

    announce(communicator, "projecting the n=3 Lane-Emden state for gravity backend finalists");
    const Clock::time_point projectionStart = Clock::now();
    auto projected =
        seed::makeProjectedEquilibriumState(problem, seed::LaneEmden({.radialSampleCount = radialSampleCount}));
    setupTimings.seedProjectionSeconds       = maximumRankSeconds(projectionStart, communicator);

    const Clock::time_point preparationStart = Clock::now();
    const auto preparation                   = problem.Prepare(projected.values, makeDependencies(), zeroRotation());
    REQUIRE(preparation.assembledResidual);
    setupTimings.operatorPreparationSeconds = maximumRankSeconds(preparationStart, communicator);

    const mfem::Operator &jacobian          = problem.GetLinearizationOperator();
    const mfem::Vector exactCorrection =
        blockBalancedDirection(problem.StateSize(), problem.GetManifest().valueBlocks(), 0.23, communicator);
    mfem::Vector rightHandSide(problem.EquationSize());
    const Clock::time_point rightHandSideStart = Clock::now();
    jacobian.Mult(exactCorrection, rightHandSide);
    setupTimings.manufacturedRightHandSideSeconds = maximumRankSeconds(rightHandSideStart, communicator);
    const mfem::Vector arnoldiDirection =
        blockBalancedDirection(problem.EquationSize(), problem.GetManifest().residualBlocks(), 0.79, communicator);
    constexpr preconditioning::MaterialSurfaceDiagonalOptions materialOptions{};

    prepareAndMeasureComposed(
        {.name                          = "current_structure_diagonal_mass_amg2",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "independent_subsystems",
         .gravityFactorization          = "approximate_ldu_diagonal_mass_amg2",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::SurfaceThenMaterialTriangular{}, preconditioning::IndependentStellarSubsystems{},
        materialOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator,
        backend::Diagonal{}, 2
    );
    prepareAndMeasureComposed(
        {.name                          = "current_structure_chebyshev3_mass_amg2",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "independent_subsystems",
         .gravityFactorization          = "approximate_ldu_chebyshev3_mass_amg2",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::SurfaceThenMaterialTriangular{}, preconditioning::IndependentStellarSubsystems{},
        materialOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator,
        backend::MatrixFreeChebyshev{.order = 3, .powerIterations = 20}, 2
    );
    prepareAndMeasureComposed(
        {.name                          = "current_structure_chebyshev4_mass_amg3",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "independent_subsystems",
         .gravityFactorization          = "approximate_ldu_chebyshev4_mass_amg3",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::SurfaceThenMaterialTriangular{}, preconditioning::IndependentStellarSubsystems{},
        materialOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator,
        backend::MatrixFreeChebyshev{.order = 4, .powerIterations = 20}, 3
    );
    prepareAndMeasureComposed(
        {.name                          = "current_structure_chebyshev5_mass_amg3",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "independent_subsystems",
         .gravityFactorization          = "approximate_ldu_chebyshev5_mass_amg3",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::SurfaceThenMaterialTriangular{}, preconditioning::IndependentStellarSubsystems{},
        materialOptions, exactCorrection, rightHandSide, arnoldiDirection, setupTimings, communicator,
        backend::MatrixFreeChebyshev{.order = 5, .powerIterations = 20}, 3
    );
}

TEST_CASE(
    "Full Stellar Physical Newton Right Hand Side And Damped Trial States",
    "[preconditioning][diagnostics][experiment][p10][full_system][physical_newton_rhs]"
) {
    using namespace mean_field;

    constexpr int radialSampleCount = 4096;
    const utils::Args arguments     = test_utils::setup_args();
    fem::FEM finiteElements         = fem::setup_fem(arguments.mesh_file, arguments, 0);
    REQUIRE(finiteElements.okay());
    const MPI_Comm communicator = finiteElements.mesh->GetComm();
    int communicatorSize        = 0;
    MPI_Comm_size(communicator, &communicatorSize);
    REQUIRE(communicatorSize == 1);

    constexpr double stellarRadius                     = utils::RADIUS;
    constexpr double targetMass                        = utils::MASS;
    const seed::DimensionlessLaneEmdenSolution profile = seed::integrateLaneEmden(3.0, 10.0);
    REQUIRE(profile.firstZeroCoordinate.has_value());
    const double surfaceCoordinate = *profile.firstZeroCoordinate;
    const double surfaceDerivative = profile.thetaDerivative(profile.thetaDerivative.Size() - 1);
    const double dimensionlessMass = -surfaceCoordinate * surfaceCoordinate * surfaceDerivative;
    REQUIRE(dimensionlessMass > 0.0);
    const double massScale          = targetMass / (4.0 * std::numbers::pi_v<double> * dimensionlessMass);
    const double polytropicConstant = std::numbers::pi_v<double> * utils::G * std::pow(massScale, 2.0 / 3.0);
    const double radialScale        = stellarRadius / surfaceCoordinate;
    const double centralDensity =
        std::pow(polytropicConstant / (std::numbers::pi_v<double> * utils::G * radialScale * radialScale), 1.5);

    auto model = model::StellarModel(
        eos::Polytrope({.n = 3.0, .K = polytropicConstant}),
        surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}),
        integral::FixedTotalMass({.Mtotal = dimensions::MassValue{targetMass}}),
        constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{centralDensity}})
    );
    auto problem = equilibrium::discretize(model, std::move(finiteElements));
    auto projected =
        seed::makeProjectedEquilibriumState(problem, seed::LaneEmden({.radialSampleCount = radialSampleCount}));
    auto dependencies      = makeDependencies();
    const auto rotation    = zeroRotation();
    const auto preparation = problem.Prepare(projected.values, dependencies, rotation);
    REQUIRE(preparation.assembledResidual);

    mfem::Vector baseResidual;
    problem.BuildResidual(baseResidual);
    const double baseResidualNorm = globalNorm(baseResidual, communicator);
    REQUIRE(std::isfinite(baseResidualNorm));
    REQUIRE(baseResidualNorm > 0.0);
    mfem::Vector rightHandSide(baseResidual);
    rightHandSide *= -1.0;

    constexpr preconditioning::MaterialSurfaceDiagonalOptions uncalibratedMaterialOptions{};
    constexpr preconditioning::MaterialSurfaceDiagonalOptions rightCalibratedMaterialOptions{
        .surfaceCalibration = {
            .target     = preconditioning::SurfaceRieszCalibrationTarget::surface_jacobian,
            .probeCount = 4,
            .objective  = preconditioning::SurfaceRieszCalibrationObjective::right_preconditioned_action
        }
    };

    std::vector<NewtonCandidateResult> candidates;
    candidates.reserve(4);
    solver::IdentityPreconditioner identity(problem.StateSize());
    candidates.push_back(solveNewtonCandidate(
        {.name                          = "identity",
         .materialFactorization         = "identity",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "identity",
         .gravityFactorization          = "identity",
         .specificationBorder           = "identity"},
        identity, 0.0, problem, baseResidual, rightHandSide, communicator
    ));
    candidates.push_back(prepareAndSolveNewtonComposed(
        {.name                          = "current_default",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "independent_subsystems",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::IndependentStellarSubsystems{}, uncalibratedMaterialOptions, baseResidual,
        rightHandSide, communicator
    ));
    candidates.push_back(prepareAndSolveNewtonComposed(
        {.name                          = "surface_then_material_uncalibrated_then_gravity",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "none",
         .stellarStructureFactorization = "material_then_gravity_triangular",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::MaterialThenGravityTriangular{}, uncalibratedMaterialOptions, baseResidual,
        rightHandSide, communicator
    ));
    candidates.push_back(prepareAndSolveNewtonComposed(
        {.name                          = "surface_then_material_right_calibrated_then_gravity",
         .materialFactorization         = "surface_then_material_triangular",
         .surfaceCalibration            = "surface_jacobian_right_action_4_probes",
         .stellarStructureFactorization = "material_then_gravity_triangular",
         .gravityFactorization          = "approximate_ldu",
         .specificationBorder           = "compiled_dense_schur"},
        problem, preconditioning::MaterialThenGravityTriangular{}, rightCalibratedMaterialOptions, baseResidual,
        rightHandSide, communicator
    ));

    // All four corrections above were obtained while the problem remained at
    // the identical projected Lane-Emden base point.  Only now do we mutate the
    // prepared state to measure the nonlinear quality of each correction.
    constexpr std::array<double, 7> dampingFactors{1.0, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625};
    const auto valueBlocks        = problem.GetManifest().valueBlocks();
    const auto residualBlocks     = problem.GetManifest().residualBlocks();
    const auto &surfaceBlock      = findBlock(valueBlocks, "surface_deformation");
    const auto &densityBlock      = findBlock(valueBlocks, "density");
    const auto &physical          = problem.GetPreparedOperator().GetPhysicalOperator();
    const auto &domainDeformation = physical.GetDomainDeformation();
    mfem::Vector volumeDisplacement(domainDeformation.volumeDisplacementSize());

    for (const NewtonCandidateResult &candidate : candidates) {
        for (const double alpha : dampingFactors) {
            incrementStateRevisions(dependencies);
            mfem::Vector trialState(projected.values);
            trialState.Add(alpha, candidate.correction);

            int localStateIsFinite = 1;
            for (int index = 0; index < trialState.Size(); ++index) {
                if (!std::isfinite(trialState(index))) {
                    localStateIsFinite = 0;
                }
            }
            int stateIsFinite = 0;
            MPI_Allreduce(&localStateIsFinite, &stateIsFinite, 1, MPI_INT, MPI_MIN, communicator);

            const mfem::Vector density(trialState.GetData() + densityBlock.offset, densityBlock.size);
            double localMinimumDensity = std::numeric_limits<double>::infinity();
            for (int index = 0; index < density.Size(); ++index) {
                localMinimumDensity = std::min(localMinimumDensity, density(index));
            }
            double minimumDensity = 0.0;
            MPI_Allreduce(&localMinimumDensity, &minimumDensity, 1, MPI_DOUBLE, MPI_MIN, communicator);

            const mfem::Vector surfaceParameters(trialState.GetData() + surfaceBlock.offset, surfaceBlock.size);
            domainDeformation.buildVolumeDisplacement(surfaceParameters, volumeDisplacement);
            const deformation::DomainDeformationGeometryReport geometry =
                domainDeformation.inspectMappedGeometry(volumeDisplacement);
            const bool geometryIsValid = geometry.isOrientationPreserving();
            const bool densityIsValid  = std::isfinite(minimumDensity) && minimumDensity >= 0.0;
            const bool trialIsValid    = stateIsFinite != 0 && geometryIsValid && densityIsValid;

            mfem::Vector predictedResidual(baseResidual);
            predictedResidual.Add(alpha, candidate.linearAction);
            const double predictedNorm = globalNorm(predictedResidual, communicator);
            std::map<std::string, double> metrics{
                {"alpha", alpha},
                {"dependency_revision", static_cast<double>(dependencies.density.revision)},
                {"state_is_finite", stateIsFinite != 0 ? 1.0 : 0.0},
                {"density_is_nonnegative", densityIsValid ? 1.0 : 0.0},
                {"minimum_density_dof", minimumDensity},
                {"geometry_is_orientation_preserving", geometryIsValid ? 1.0 : 0.0},
                {"minimum_mapping_jacobian_determinant", geometry.minimumJacobianDeterminant},
                {"trial_is_valid", trialIsValid ? 1.0 : 0.0},
                {"base_residual_norm", baseResidualNorm},
                {"correction_norm", globalNorm(candidate.correction, communicator)},
                {"damped_correction_norm", alpha * globalNorm(candidate.correction, communicator)},
                {"predicted_residual_norm", predictedNorm},
                {"predicted_relative_residual", predictedNorm / baseResidualNorm},
                {"predicted_fractional_reduction", 1.0 - predictedNorm / baseResidualNorm}
            };

            auto parameters =
                newtonParameters(candidate.candidate, "damped_physical_newton_trial", problem.StateSize());
            parameters["trial_status"] = trialIsValid ? "prevalidated" : "rejected_before_residual_evaluation";
            if (!trialIsValid) {
                experiment::record_experiment_result(
                    "stellar_preconditioning_p10_newton", candidate.candidate.name + "_alpha_" + std::to_string(alpha),
                    std::move(parameters), std::move(metrics)
                );
                continue;
            }

            try {
                problem.Prepare(trialState, dependencies, rotation);
                mfem::Vector actualResidual;
                problem.BuildResidual(actualResidual);
                const double actualNorm = globalNorm(actualResidual, communicator);
                mfem::Vector nonlinearRemainder(actualResidual);
                nonlinearRemainder -= predictedResidual;
                const double nonlinearRemainderNorm = globalNorm(nonlinearRemainder, communicator);
                mfem::Vector residualDeparture(actualResidual);
                residualDeparture -= baseResidual;
                const double residualDepartureNorm     = globalNorm(residualDeparture, communicator);

                metrics["residual_evaluated"]          = 1.0;
                metrics["actual_residual_norm"]        = actualNorm;
                metrics["actual_relative_residual"]    = actualNorm / baseResidualNorm;
                metrics["actual_fractional_reduction"] = 1.0 - actualNorm / baseResidualNorm;
                metrics["nonlinear_remainder_norm"]    = nonlinearRemainderNorm;
                metrics["relative_nonlinear_remainder"] =
                    nonlinearRemainderNorm / std::max(residualDepartureNorm, std::numeric_limits<double>::min());
                metrics["actual_to_predicted_norm_ratio"] =
                    actualNorm / std::max(predictedNorm, std::numeric_limits<double>::min());

                for (const auto &block : residualBlocks) {
                    const mfem::Vector baseBlock(baseResidual.GetData() + block.offset, block.size);
                    const mfem::Vector predictedBlock(predictedResidual.GetData() + block.offset, block.size);
                    const mfem::Vector actualBlock(actualResidual.GetData() + block.offset, block.size);
                    const mfem::Vector remainderBlock(nonlinearRemainder.GetData() + block.offset, block.size);
                    const double baseBlockNorm          = globalNorm(baseBlock, communicator);
                    const std::string prefix            = "residual_block." + std::string(block.stableId);
                    metrics[prefix + ".base_norm"]      = baseBlockNorm;
                    metrics[prefix + ".predicted_norm"] = globalNorm(predictedBlock, communicator);
                    metrics[prefix + ".actual_norm"]    = globalNorm(actualBlock, communicator);
                    metrics[prefix + ".actual_ratio"] =
                        metrics[prefix + ".actual_norm"] / std::max(baseBlockNorm, std::numeric_limits<double>::min());
                    metrics[prefix + ".remainder_norm"] = globalNorm(remainderBlock, communicator);
                }
                parameters["trial_status"] = "evaluated";
            } catch (const std::exception &) {
                metrics["residual_evaluated"] = 0.0;
                parameters["trial_status"]    = "residual_evaluation_threw";
            }

            experiment::record_experiment_result(
                "stellar_preconditioning_p10_newton", candidate.candidate.name + "_alpha_" + std::to_string(alpha),
                std::move(parameters), std::move(metrics)
            );
        }
    }
}