#include <algorithm>
#include <catch2/catch_test_macros.hpp>
#include <cmath>
#include <cstdint>
#include <limits>
#include <memory>
#include <mfem.hpp>
#include <numbers>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>

import mean_field;
import test_helpers;

namespace {
    bool vector_is_finite(const mfem::Vector &vector) {
        for (int index = 0; index < vector.Size(); ++index) {
            if (!std::isfinite(vector(index))) {
                return false;
            }
        }
        return true;
    }

    mfem::Vector make_deterministic_vector(
        const int size,
        const double phase
    ) {
        mfem::Vector vector(size);
        for (int index = 0; index < size; ++index) {
            const double coordinate = static_cast<double>(index + 1);
            vector(index)           = std::sin(phase + 0.017 * coordinate) + 0.25 * std::cos(0.031 * coordinate);
        }
        return vector;
    }

    void require_all_ranks(
        const bool localCondition,
        const MPI_Comm communicator,
        const char *description
    ) {
        int rank = 0;
        int size = 0;
        MPI_Comm_rank(communicator, &rank);
        MPI_Comm_size(communicator, &size);

        const int localFailure = localCondition ? size : rank;
        int firstFailure       = size;
        const int result       = MPI_Allreduce(&localFailure, &firstFailure, 1, MPI_INT, MPI_MIN, communicator);
        REQUIRE(result == MPI_SUCCESS);
        CAPTURE(description, localCondition, firstFailure);
        REQUIRE(firstFailure == size);
    }

    double global_dot(
        const mfem::Vector &left,
        const mfem::Vector &right,
        const MPI_Comm communicator
    ) {
        require_all_ranks(left.Size() == right.Size(), communicator, "global dot-product vector sizes");
        const double local = left * right;
        double global      = 0.0;
        REQUIRE(MPI_Allreduce(&local, &global, 1, MPI_DOUBLE, MPI_SUM, communicator) == MPI_SUCCESS);
        return global;
    }

    double global_norm(
        const mfem::Vector &vector,
        const MPI_Comm communicator
    ) {
        return std::sqrt(global_dot(vector, vector, communicator));
    }

    [[nodiscard]] mean_field::operators::StellarEquilibriumDependencies
    make_stellar_dependencies(const std::uint64_t revision = 1) {
        return {
            .discretization     = {.identity = 16101, .revision = 1},
            .density            = {.identity = 16103, .revision = revision},
            .surfaceDeformation = {.identity = 16111, .revision = revision},
            .gravityGradient    = {.identity = 16127, .revision = revision},
            .gravityPotential   = {.identity = 16139, .revision = revision},
            .enthalpy           = {.identity = 16141, .revision = revision},
            .bernoulliConstant  = {.identity = 16183, .revision = revision},
            .rotation           = {.identity = 16187, .revision = revision},
            .targetMass         = {.identity = 16189, .revision = 1}
        };
    }

    class FixedCentralDensityPhysicalCoreStub final : public mfem::Operator {
    public:
        FixedCentralDensityPhysicalCoreStub() : mfem::Operator(0) {
        }

        void Mult(
            const mfem::Vector &,
            mfem::Vector &
        ) const override {
        }

        const mean_field::operators::StellarEquilibriumLayout &GetLayout() const;
        mean_field::operators::PreparedStellarEquilibriumReport Prepare(
            const mfem::Vector &,
            const mean_field::operators::StellarEquilibriumDependencies &,
            const mean_field::physics::RigidRotation &
        );
        void BuildResidual(mfem::Vector &) const;
        bool IsPrepared() const noexcept;
        mean_field::operators::RootConstraintReport GetFixedMassReport() const;
        const mean_field::operators::StellarEquilibriumDependencies &GetDependencies() const;
        const mean_field::operators::StellarEquilibriumDependencyStamp &GetGeneratedDisplacementDependency() const;
        const mean_field::operators::PreparedPressureSurfaceConstraint &GetSurfaceConstraintOperator() const;
    };

    static_assert(mean_field::operators::PreparedStellarEquilibriumPhysicalCore<FixedCentralDensityPhysicalCoreStub>);

    struct FixedCentralDensityStateView final {
        const mfem::Vector &enthalpy;
        const mfem::Vector &phase;

        const mfem::Vector &block(const mean_field::utils::blocks::enthalpy::specific &) const noexcept {
            return enthalpy;
        }

        const mfem::Vector &
        block(const mean_field::utils::blocks::fixed_central_density::central_value &) const noexcept {
            return phase;
        }
    };

    void check_rank_consistent_scalar(
        const double value,
        const MPI_Comm communicator,
        const double relativeTolerance = 2.0e-13
    ) {
        double minimum = 0.0;
        double maximum = 0.0;
        REQUIRE(MPI_Allreduce(&value, &minimum, 1, MPI_DOUBLE, MPI_MIN, communicator) == MPI_SUCCESS);
        REQUIRE(MPI_Allreduce(&value, &maximum, 1, MPI_DOUBLE, MPI_MAX, communicator) == MPI_SUCCESS);
        CAPTURE(value, minimum, maximum);
        CHECK(std::isfinite(minimum));
        CHECK(std::isfinite(maximum));
        CHECK(std::abs(maximum - minimum) <= relativeTolerance * std::max({1.0, std::abs(minimum), std::abs(maximum)}));
    }

    class DistributedDiagonalOperator final : public mfem::Operator {
    public:
        DistributedDiagonalOperator() : mfem::Operator(2) {
        }

        void Mult(
            const mfem::Vector &input,
            mfem::Vector &output
        ) const override {
            output.SetSize(2);
            output(0) = input(0);
            output(1) = 2.0 * input(1);
        }
    };

    class DistributedIdentityInverse final : public mfem::Solver {
    public:
        static constexpr mean_field::preconditioning::ApplicationContract applicationContract =
            mean_field::preconditioning::ApplicationContract::flexible;

        DistributedIdentityInverse() : mfem::Solver(2) {
        }

        void SetOperator(const mfem::Operator &operation) override {
            if (operation.Height() != 2 || operation.Width() != 2) {
                throw std::invalid_argument("The distributed identity inverse requires a two-by-two operator.");
            }
        }

        void Mult(
            const mfem::Vector &input,
            mfem::Vector &output
        ) const override {
            output = input;
        }
    };

} // namespace

TEST_CASE(
    "MPI MFEM FGMRES Reports Global Warm-Start And Restart Diagnostics Consistently",
    "[mpi][distributed][solver][fgmres]"
) {
    using namespace mean_field;

    const MPI_Comm communicator = MPI_COMM_WORLD;
    int rank                    = 0;
    REQUIRE(MPI_Comm_rank(communicator, &rank) == MPI_SUCCESS);

    DistributedDiagonalOperator operation;
    DistributedIdentityInverse inverse;
    auto backend = solver::linear::prepareLinearBackend(
        solver::linear::FGMRES({.restartLength = 1, .printLevel = -1}), operation, inverse, communicator
    );
    require_all_ranks(backend.IsReady(), communicator, "distributed FGMRES backend readiness");

    mfem::Vector rightHandSide(2);
    rightHandSide(0) = 1.0 + static_cast<double>(rank);
    rightHandSide(1) = 2.0 + 0.5 * static_cast<double>(rank);

    mfem::Vector correction(2);
    correction(0) = 0.125 * static_cast<double>(rank + 1);
    correction(1) = -0.25 * static_cast<double>(rank + 1);

    mfem::Vector initialAction(2);
    operation.Mult(correction, initialAction);
    mfem::Vector initialResidual(rightHandSide);
    initialResidual -= initialAction;
    const double expectedInitialResidualNorm = global_norm(initialResidual, communicator);

    const solver::LinearSolveControl control{
        .relativeTolerance = 1.0e-10, .absoluteTolerance = 1.0e-13, .maximumIterations = 100
    };
    const auto report                 = backend.Solve(rightHandSide, correction, control);

    const auto checkConsistentInteger = [&](const int value) {
        int minimum = 0;
        int maximum = 0;
        REQUIRE(MPI_Allreduce(&value, &minimum, 1, MPI_INT, MPI_MIN, communicator) == MPI_SUCCESS);
        REQUIRE(MPI_Allreduce(&value, &maximum, 1, MPI_INT, MPI_MAX, communicator) == MPI_SUCCESS);
        CHECK(minimum == maximum);
    };
    checkConsistentInteger(static_cast<int>(report.status));
    checkConsistentInteger(report.iterations);
    checkConsistentInteger(report.restarts);
    check_rank_consistent_scalar(report.initialResidualNorm, communicator);
    check_rank_consistent_scalar(report.reportedResidualNorm, communicator);
    check_rank_consistent_scalar(report.trueResidualNorm, communicator);
    check_rank_consistent_scalar(report.relativeTrueResidualNorm, communicator);
    require_all_ranks(report.Converged(), communicator, "distributed FGMRES convergence");

    mfem::Vector finalAction(2);
    operation.Mult(correction, finalAction);
    mfem::Vector trueResidual(rightHandSide);
    trueResidual -= finalAction;
    const double expectedTrueResidualNorm = global_norm(trueResidual, communicator);
    const double rightHandSideNorm        = global_norm(rightHandSide, communicator);
    CHECK(report.status == solver::LinearSolveStatus::converged);
    CHECK(report.iterations > 1);
    CHECK(report.iterations <= control.maximumIterations);
    CHECK(report.restarts == report.iterations - 1);
    CHECK(
        std::abs(report.initialResidualNorm - expectedInitialResidualNorm) <=
        2.0e-13 * std::max(1.0, expectedInitialResidualNorm)
    );
    CHECK(
        std::abs(report.trueResidualNorm - expectedTrueResidualNorm) <=
        2.0e-13 * std::max(1.0, expectedTrueResidualNorm)
    );
    CHECK(report.trueResidualNorm <= control.ConvergenceThreshold(rightHandSideNorm));
    CHECK(
        std::abs(report.relativeTrueResidualNorm - report.trueResidualNorm / rightHandSideNorm) <=
        2.0e-13 * std::max(1.0, report.relativeTrueResidualNorm)
    );

    mfem::Vector invalidRightHandSide(rightHandSide);
    if (rank == 0) {
        invalidRightHandSide(0) = std::numeric_limits<double>::quiet_NaN();
    }
    mfem::Vector rejectedCorrection(correction);
    bool collectivelyRejected = false;
    try {
        (void)backend.Solve(invalidRightHandSide, rejectedCorrection, control);
    } catch (const std::invalid_argument &) {
        collectivelyRejected = true;
    }
    require_all_ranks(collectivelyRejected, communicator, "collective non-finite FGMRES input rejection");

    auto rankDependentBackend = solver::linear::prepareLinearBackend(
        solver::linear::FGMRES({.restartLength = rank + 1, .printLevel = -1}), operation, inverse, communicator
    );
    mfem::Vector rankDependentCorrection(2);
    rankDependentCorrection           = 0.0;
    bool rankDependentBackendRejected = false;
    try {
        (void)rankDependentBackend.Solve(rightHandSide, rankDependentCorrection, control);
    } catch (const std::invalid_argument &) {
        rankDependentBackendRejected = true;
    }
    require_all_ranks(rankDependentBackendRejected, communicator, "rank-dependent FGMRES option rejection");

    auto rankDependentControl              = control;
    rankDependentControl.maximumIterations = control.maximumIterations + rank;
    rankDependentCorrection                = 0.0;
    bool rankDependentControlRejected      = false;
    try {
        (void)backend.Solve(rightHandSide, rankDependentCorrection, rankDependentControl);
    } catch (const std::invalid_argument &) {
        rankDependentControlRejected = true;
    }
    require_all_ranks(rankDependentControlRejected, communicator, "rank-dependent FGMRES solve-control rejection");
}

TEST_CASE(
    "MPI Observer Exceptions Cause Rank-Coordinated Solver Shutdown",
    "[mpi][distributed][solver][newton][observer][exception-safety]"
) {
    using namespace mean_field;

    const MPI_Comm communicator = MPI_COMM_WORLD;
    int rank                    = 0;
    int size                    = 0;
    REQUIRE(MPI_Comm_rank(communicator, &rank) == MPI_SUCCESS);
    REQUIRE(MPI_Comm_size(communicator, &size) == MPI_SUCCESS);
    REQUIRE(size >= 2);

    auto observer = solver::nonlinear::makeObserver(
        [rank](const solver::nonlinear::BeforeIteration &) {
            if (rank == 1) {
                throw std::runtime_error("rank-local before observer failure");
            }
        },
        [rank](const solver::nonlinear::AfterLineSearchTrial &) {
            if (rank == 1) {
                throw std::runtime_error("rank-local trial observer failure");
            }
        },
        [rank](const solver::nonlinear::AfterIteration &) {
            if (rank == 1) {
                throw std::runtime_error("rank-local after observer failure");
            }
        }
    );

    const auto requireCoordinatedEscape = [&](auto &&invoke, const char *description) {
        bool caught = false;
        try {
            invoke();
        } catch (const std::runtime_error &) {
            caught = true;
        }
        require_all_ranks(caught, communicator, description);

        const int localHandshake = rank + 1;
        int globalHandshake      = 0;
        REQUIRE(MPI_Allreduce(&localHandshake, &globalHandshake, 1, MPI_INT, MPI_SUM, communicator) == MPI_SUCCESS);
        CHECK(globalHandshake == size * (size + 1) / 2);
    };

    requireCoordinatedEscape(
        [&] {
            solver::nonlinear::detail::InvokeBeforeIteration(
                observer, solver::nonlinear::BeforeIteration{.communicator = communicator}
            );
        },
        "before observer exception escape"
    );
    requireCoordinatedEscape(
        [&] {
            solver::nonlinear::detail::InvokeAfterLineSearchTrial(
                observer, solver::nonlinear::AfterLineSearchTrial{.communicator = communicator}
            );
        },
        "line-search observer exception escape"
    );
    requireCoordinatedEscape(
        [&] {
            solver::nonlinear::detail::InvokeAfterIteration(
                observer, solver::nonlinear::AfterIteration{.communicator = communicator}
            );
        },
        "after observer exception escape"
    );
}

TEST_CASE(
    "MPI Runtime Preserves World And Split Communicator Membership",
    "[mpi][distributed][unit]"
) {
    int rank = 0;
    int size = 1;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    std::vector<int> ranks(static_cast<std::size_t>(size), -1);
    MPI_Allgather(&rank, 1, MPI_INT, ranks.data(), 1, MPI_INT, MPI_COMM_WORLD);

    CHECK(size >= 2);
    for (int expected = 0; expected < size; ++expected) {
        CHECK(ranks[expected] == expected);
    }

    MPI_Comm parity_communicator = MPI_COMM_NULL;
    MPI_Comm_split(MPI_COMM_WORLD, rank % 2, rank, &parity_communicator);

    int parity_size = 0;
    MPI_Comm_size(parity_communicator, &parity_size);
    const int expected_parity_size = (size + 1 - rank % 2) / 2;
    CHECK(parity_size == expected_parity_size);

    MPI_Comm_free(&parity_communicator);
}

TEST_CASE(
    "MPI FEM Setup Partitions Every Element Exactly Once",
    "[mpi][distributed][mesh][integration]"
) {
    const mean_field::utils::Args args = test_utils::setup_args();
    const mean_field::fem::FEM f       = mean_field::fem::setup_fem(args.mesh_file, args, 0);

    const long long local_elements     = f.mesh->GetNE();
    long long global_elements          = 0;
    long long minimum_elements         = 0;

    MPI_Allreduce(&local_elements, &global_elements, 1, MPI_LONG_LONG, MPI_SUM, f.mesh->GetComm());
    MPI_Allreduce(&local_elements, &minimum_elements, 1, MPI_LONG_LONG, MPI_MIN, f.mesh->GetComm());

    CHECK(global_elements == f.smesh.mesh->GetNE());
    CHECK(minimum_elements > 0);
    CHECK(f.logicalReferenceMesh->GetNE() == f.mesh->GetNE());
}

TEST_CASE(
    "MPI Fixed Angular Momentum Produces One Consistent Global Invariant Row",
    "[mpi][distributed][fixed-angular-momentum][physics][jacobian]"
) {
    using namespace mean_field;

    const auto args         = test_utils::setup_args();
    fem::FEM finiteElements = fem::setup_fem(args.mesh_file, args, 0);
    require_all_ranks(finiteElements.okay(), MPI_COMM_WORLD, "fixed-angular-momentum FEM setup");
    const MPI_Comm communicator = finiteElements.mesh->GetComm();

    mfem::ParGridFunction densityField(finiteElements.densityFes.get());
    mfem::ConstantCoefficient densityCoefficient(1.23);
    densityField.ProjectCoefficient(densityCoefficient);
    mfem::Vector densityTrue;
    densityField.GetTrueDofs(densityTrue);

    mfem::Vector displacementTrue(finiteElements.displacementFes->GetTrueVSize());
    displacementTrue = 0.0;
    finiteElements.displacement->SetFromTrueDofs(displacementTrue);
    mfem::Vector gravityGradientTrue(finiteElements.gravityFluxFes->GetTrueVSize());
    mfem::Vector gravityPotentialTrue(finiteElements.gravityPotentialFes->GetTrueVSize());
    gravityGradientTrue  = 0.0;
    gravityPotentialTrue = 0.0;

    operators::context::gravity_field::GravityFieldLinearizationContext gravityContext(
        finiteElements, *finiteElements.domainMapperStateless
    );
    gravityContext.Prepare(
        {.density           = gravityContext.GetDensityMap().gather(densityTrue),
         .displacement      = gravityContext.GetDisplacementMap().gather(displacementTrue),
         .gravity_gradient  = gravityContext.GetGravityGradientMap().gather(gravityGradientTrue),
         .gravity_potential = gravityContext.GetGravityPotentialMap().gather(gravityPotentialTrue)},
        {.discretization    = {.value = 2},
         .displacement      = {.value = 3},
         .density           = {.value = 5},
         .gravity_gradient  = {.value = 7},
         .gravity_potential = {.value = 11}}
    );

    constexpr double targetAngularMomentum = 0.37;
    constexpr double angularVelocity       = 0.61;
    operators::PreparedAngularMomentumOperator invariant(
        finiteElements, *finiteElements.domainMapperStateless, gravityContext,
        models::compileConstraint(
            integral::FixedAngularMomentum({.Jtotal = dimensions::AngularMomentumValue{targetAngularMomentum}})
        )
    );
    const operators::AngularMomentumDependencies dependencies{
        .discretization = {.identity = 101, .revision = 2},
        .density        = {.identity = 103, .revision = 5},
        .displacement   = {.identity = 107, .revision = 3},
        .rotation       = {.identity = 109, .revision = 13}
    };
    const auto preparation = invariant.Prepare(angularVelocity, dependencies);
    CHECK(preparation.rebuiltStaticPlan);
    CHECK(preparation.refreshedGeometry);
    CHECK(preparation.refreshedDensity);
    CHECK(preparation.updatedAngularVelocity);
    CHECK(preparation.assembledResidual);

    const double independentMoment = analysis::get_moment_of_inertia(finiteElements, densityField);
    const double preparedMoment    = invariant.GetMomentOfInertia();
    const double comparisonScale   = std::max({std::abs(independentMoment), std::abs(preparedMoment), 1.0e-300});
    CHECK(std::abs(preparedMoment - independentMoment) / comparisonScale <= 3.0e-13);

    double minimumMoment = 0.0;
    double maximumMoment = 0.0;
    MPI_Allreduce(&preparedMoment, &minimumMoment, 1, MPI_DOUBLE, MPI_MIN, finiteElements.mesh->GetComm());
    MPI_Allreduce(&preparedMoment, &maximumMoment, 1, MPI_DOUBLE, MPI_MAX, finiteElements.mesh->GetComm());
    CHECK(std::abs(maximumMoment - minimumMoment) / comparisonScale <= 2.0e-15);

    mfem::Vector residual;
    invariant.BuildResidual(residual);
    require_all_ranks(residual.Size() == 1, communicator, "fixed-angular-momentum residual size");
    CHECK(
        std::abs(residual(0) - (angularVelocity * independentMoment - targetAngularMomentum)) /
            std::max({std::abs(residual(0)), std::abs(angularVelocity * independentMoment), 1.0}) <=
        3.0e-13
    );

    mfem::Vector densityAction;
    invariant.ApplyDensityJacobianAction(gravityContext.GetDensityMap().gather(densityTrue), densityAction);
    require_all_ranks(densityAction.Size() == 1, communicator, "fixed-angular-momentum density action size");
    CHECK(std::abs(densityAction(0) - angularVelocity * preparedMoment) / comparisonScale <= 3.0e-13);

    constexpr double angularVelocityVariation = -0.29;
    mfem::Vector angularVelocityAction;
    invariant.ApplyAngularVelocityJacobianAction(angularVelocityVariation, angularVelocityAction);
    require_all_ranks(angularVelocityAction.Size() == 1, communicator, "fixed-angular-momentum rotation action size");
    CHECK(std::abs(angularVelocityAction(0) - angularVelocityVariation * preparedMoment) / comparisonScale <= 2.0e-15);

    int rank = 0;
    REQUIRE(MPI_Comm_rank(communicator, &rank) == MPI_SUCCESS);
    auto rejectedDependencies = dependencies;
    ++rejectedDependencies.rotation.revision;
    const double rankLocalAngularVelocity = rank == 0 ? std::numeric_limits<double>::quiet_NaN() : angularVelocity;
    const auto synchronizedRejection      = invariant.TryPrepare(rankLocalAngularVelocity, rejectedDependencies);
    require_all_ranks(
        !synchronizedRejection.has_value(), communicator, "fixed-angular-momentum synchronized trial rejection"
    );
    CHECK(
        synchronizedRejection.error().reason ==
        operators::AngularMomentumPreparationRejectionReason::non_finite_angular_velocity
    );
    CHECK_FALSE(invariant.IsPrepared());

    ++rejectedDependencies.rotation.revision;
    const auto recovered = invariant.TryPrepare(angularVelocity, rejectedDependencies);
    require_all_ranks(recovered.has_value(), communicator, "fixed-angular-momentum recovery after rejection");
    CHECK(invariant.IsPrepared());
}

TEST_CASE(
    "MPI Fixed Central Density Coordinates Reject Collectively Without Unwinding",
    "[mpi][distributed][fixed-central-density][candidate-rejection]"
) {
    using namespace mean_field;

    const utils::Args arguments = test_utils::setup_args();
    fem::FEM finiteElements     = fem::setup_fem(arguments.mesh_file, arguments, 0);
    require_all_ranks(finiteElements.okay(), MPI_COMM_WORLD, "fixed-central-density FEM setup");
    const MPI_Comm communicator = finiteElements.mesh->GetComm();
    require_all_ranks(
        finiteElements.domainMapperStateless != nullptr, communicator, "fixed-central-density mapper setup"
    );

    auto model = model::StellarModel(
        eos::Polytrope({.n = 1.0, .K = 0.25}), surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}),
        integral::FixedTotalMass({.Mtotal = dimensions::MassValue{1.0}}),
        constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{1.0}})
    );
    FixedCentralDensityPhysicalCoreStub physicalCore;
    using Runtime = operators::detail::FixedCentralDensityRuntime<std::remove_cvref_t<decltype(model)>>;
    Runtime runtime(finiteElements, *finiteElements.domainMapperStateless, physicalCore, model);

    mfem::Vector enthalpy(runtime.constraint().GetCenterDof().field_size());
    enthalpy = 1.0;
    mfem::Vector phase(1);
    phase = 0.03;
    const FixedCentralDensityStateView state{enthalpy, phase};

    const auto initial = runtime.TryPrepareAfterPhysical(state, make_stellar_dependencies(1), physicalCore);
    require_all_ranks(initial.has_value(), communicator, "initial fixed-central-density preparation");

    int rank = 0;
    REQUIRE(MPI_Comm_rank(communicator, &rank) == MPI_SUCCESS);
    phase(0)                = rank == 0 ? 0.04 : 0.03;
    const auto inconsistent = runtime.TryPrepareAfterPhysical(state, make_stellar_dependencies(2), physicalCore);
    require_all_ranks(!inconsistent.has_value(), communicator, "fixed-central-density inconsistent scalar rejection");
    CHECK(inconsistent.error().reason == operators::StellarEquilibriumPreparationRejectionReason::inadmissible_physics);
    CHECK(inconsistent.error().stage == operators::StellarEquilibriumPreparationStage::model_specification);
    require_all_ranks(!runtime.IsPrepared(), communicator, "fixed-central-density inconsistent preparation state");

    phase(0)            = rank == 0 ? std::numeric_limits<double>::quiet_NaN() : 0.03;
    const auto rejected = runtime.TryPrepareAfterPhysical(state, make_stellar_dependencies(3), physicalCore);
    require_all_ranks(!rejected.has_value(), communicator, "fixed-central-density synchronized trial rejection");
    CHECK(rejected.error().reason == operators::StellarEquilibriumPreparationRejectionReason::non_finite_physics);
    CHECK(rejected.error().stage == operators::StellarEquilibriumPreparationStage::model_specification);
    require_all_ranks(!runtime.IsPrepared(), communicator, "fixed-central-density rejected preparation state");

    phase                = 0.03;
    const auto recovered = runtime.TryPrepareAfterPhysical(state, make_stellar_dependencies(4), physicalCore);
    require_all_ranks(recovered.has_value(), communicator, "fixed-central-density recovery after rejection");
    require_all_ranks(runtime.IsPrepared(), communicator, "fixed-central-density recovered preparation state");
}

TEST_CASE(
    "MPI Material Surface Regularization Rejects Rank-Local Non-Finite Data Collectively",
    "[mpi][distributed][preconditioner][material-surface][exception-safety]"
) {
    using namespace mean_field;

    const MPI_Comm communicator = MPI_COMM_WORLD;
    int rank                    = 0;
    REQUIRE(MPI_Comm_rank(communicator, &rank) == MPI_SUCCESS);

    mfem::Vector diagonal(1);
    diagonal(0)               = rank == 0 ? std::numeric_limits<double>::quiet_NaN() : 1.0;
    bool rejectedCollectively = false;
    try {
        (void)preconditioning::detail::regularizeMaterialSurfaceDiagonal(
            diagonal, preconditioning::MaterialSurfaceDiagonalOptions{}, communicator
        );
    } catch (const std::invalid_argument &) {
        rejectedCollectively = true;
    }
    require_all_ranks(rejectedCollectively, communicator, "material-surface non-finite diagonal collective rejection");

    diagonal             = 0.0;
    const auto recovered = preconditioning::detail::regularizeMaterialSurfaceDiagonal(
        diagonal, preconditioning::MaterialSurfaceDiagonalOptions{}, communicator
    );
    CHECK(std::isfinite(diagonal(0)));
    CHECK(diagonal(0) == recovered.appliedFloor);

    std::uint64_t communicatorSize = 0;
    int size                       = 0;
    REQUIRE(MPI_Comm_size(communicator, &size) == MPI_SUCCESS);
    communicatorSize = static_cast<std::uint64_t>(size);
    CHECK(recovered.regularizedEntries == communicatorSize);
}

TEST_CASE(
    "MPI Density Volume Context Forwards The Prepared Global Mass Integral",
    "[mpi][distributed][integral-context][physics-extension]"
) {
    using namespace mean_field;

    const utils::Args arguments = test_utils::setup_args();
    fem::FEM finiteElements     = fem::setup_fem(arguments.mesh_file, arguments, 0);
    require_all_ranks(finiteElements.okay(), MPI_COMM_WORLD, "density-volume context FEM setup");
    const MPI_Comm communicator   = finiteElements.mesh->GetComm();

    constexpr double densityValue = 1.23;
    mfem::ParGridFunction densityField(finiteElements.densityFes.get());
    mfem::ConstantCoefficient densityCoefficient(densityValue);
    densityField.ProjectCoefficient(densityCoefficient);
    mfem::Vector densityTrue;
    densityField.GetTrueDofs(densityTrue);

    mfem::Vector displacementTrue(finiteElements.displacementFes->GetTrueVSize());
    displacementTrue = 0.0;
    finiteElements.displacement->SetFromTrueDofs(displacementTrue);
    mfem::Vector gravityGradientTrue(finiteElements.gravityFluxFes->GetTrueVSize());
    mfem::Vector gravityPotentialTrue(finiteElements.gravityPotentialFes->GetTrueVSize());
    gravityGradientTrue  = 0.0;
    gravityPotentialTrue = 0.0;

    operators::context::gravity_field::GravityFieldLinearizationContext gravityContext(
        finiteElements, *finiteElements.domainMapperStateless
    );
    const mfem::Vector reducedDensity = gravityContext.GetDensityMap().gather(densityTrue);
    gravityContext.Prepare(
        {.density           = reducedDensity,
         .displacement      = gravityContext.GetDisplacementMap().gather(displacementTrue),
         .gravity_gradient  = gravityContext.GetGravityGradientMap().gather(gravityGradientTrue),
         .gravity_potential = gravityContext.GetGravityPotentialMap().gather(gravityPotentialTrue)},
        {.discretization    = {.value = 2},
         .displacement      = {.value = 3},
         .density           = {.value = 5},
         .gravity_gradient  = {.value = 7},
         .gravity_potential = {.value = 11}}
    );

    operators::PreparedMassNormalizationOperator massIntegral(
        finiteElements, *finiteElements.domainMapperStateless, gravityContext
    );
    massIntegral.Prepare(
        models::compileConstraint(integral::FixedTotalMass({.Mtotal = dimensions::MassValue{1.0}})),
        {.discretization = {.identity = 101, .revision = 2},
         .density        = {.identity = 103, .revision = 5},
         .displacement   = {.identity = 107, .revision = 3},
         .targetMass     = {.identity = 109, .revision = 1}}
    );

    class DistributedMassIntegralCore final {
    public:
        explicit DistributedMassIntegralCore(const operators::PreparedMassNormalizationOperator &mass) noexcept
            : m_mass(&mass) {
        }

        [[nodiscard]] double ApplyDensityVolumeIntegralDensityAction(const mfem::Vector &direction) const {
            mfem::Vector action;
            m_mass->ApplyDensityJacobianAction(direction, action);
            return action(0);
        }

        [[nodiscard]] double ApplyDensityVolumeIntegralSurfaceShapeAction(const mfem::Vector &) const noexcept {
            return 0.0;
        }

    private:
        const operators::PreparedMassNormalizationOperator *m_mass;
    };

    const DistributedMassIntegralCore testCore{massIntegral};

    const stellar::DensityVolumeIntegralContext<integral::FixedTotalMass> densityIntegral{testCore};
    const double integratedMass        = densityIntegral.integrateDensity(reducedDensity).value();
    const double linearizedDensityMass = densityIntegral.linearizeDensityIntegral(reducedDensity).value();
    const double preparedMass          = massIntegral.GetCurrentMass();
    const double independentMass       = densityValue * analysis::get_mesh_volume(finiteElements);

    check_rank_consistent_scalar(integratedMass, communicator);
    check_rank_consistent_scalar(linearizedDensityMass, communicator);
    check_rank_consistent_scalar(preparedMass, communicator);
    check_rank_consistent_scalar(independentMass, communicator);
    const double massScale = std::max({1.0, std::abs(integratedMass), std::abs(independentMass)});
    CHECK(std::abs(integratedMass - independentMass) <= 3.0e-12 * massScale);
    CHECK(std::abs(integratedMass - preparedMass) <= 3.0e-13 * massScale);
    CHECK(std::abs(linearizedDensityMass - integratedMass) <= 3.0e-13 * massScale);
}

TEST_CASE(
    "Single-Rank MPI Assembled Variadic Stellar Root Normalizes And Applies Its Inferred Preconditioner",
    "[mpi][single-rank][stellar-equilibrium][normalization][preconditioning][integration]"
) {
    using namespace mean_field;

    const utils::Args arguments = test_utils::setup_args();
    fem::FEM finiteElements     = fem::setup_fem(arguments.mesh_file, arguments, 0);
    require_all_ranks(finiteElements.okay(), MPI_COMM_WORLD, "variadic stellar-root FEM setup");
    const MPI_Comm communicator     = finiteElements.mesh->GetComm();

    constexpr double radius         = utils::RADIUS;
    constexpr double mass           = utils::MASS;
    const double polytropicConstant = 2.0 * utils::G * radius * radius / std::numbers::pi_v<double>;
    const double centralDensity     = std::numbers::pi_v<double> * mass / (4.0 * radius * radius * radius);

    auto model                      = model::StellarModel(
        eos::Polytrope({.n = 1.0, .K = polytropicConstant}),
        surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}),
        integral::FixedTotalMass({.Mtotal = dimensions::MassValue{mass}}),
        integral::FixedAngularMomentum({.Jtotal = dimensions::AngularMomentumValue{0.05}, .axis = {0.0, 0.0, 1.0}}),
        constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{centralDensity}})
    );
    auto problem = equilibrium::discretize(
        model,
        equilibrium::makeStellarDiscretization(
            std::move(finiteElements), normalization::PhysicalRieszDiagonal{dimensions::LengthValue{radius}, utils::G}
        )
    );
    auto projected = seed::makeProjectedEquilibriumState(
        problem, seed::LaneEmden({.centralDensity = dimensions::DensityValue{centralDensity}, .radialSampleCount = 512})
    );
    auto normalized = normalization::makeNormalizedStellarEquilibriumOperator(problem);
    using Problem   = std::remove_cvref_t<decltype(problem)>;
    using Form      = typename Problem::FormType;
    constexpr auto massResidualBlock =
        utils::blocks::get_residual_block<Form>(utils::blocks::fixed_total_mass_constraint.mass_normalization_term);
    constexpr auto angularMomentumResidualBlock =
        utils::blocks::get_residual_block<Form>(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term);
    constexpr auto centralDensityResidualBlock =
        utils::blocks::get_residual_block<Form>(utils::blocks::fixed_central_density_phase.central_value_term);
    require_all_ranks(
        problem.StateSize() == problem.EquationSize(), communicator, "variadic stellar-root square layout"
    );

    mfem::Vector normalizedState;
    normalized.NormalizeState(projected.values, normalizedState);
    const auto preparation = normalized.Prepare(normalizedState, make_stellar_dependencies(1));
    require_all_ranks(
        preparation.DidAnyWork() && preparation.generatedPhysicalControl &&
            preparation.template specification<models::FixedAngularMomentum>().generatedRotation &&
            problem.IsPrepared() && normalized.IsPrepared(),
        communicator, "initial variadic stellar-root preparation"
    );

    // The astronomy-facing integral handle must delegate to the same mapped
    // physical-volume quadrature and collective reduction as the prepared
    // mass equation, without exposing the FEM/core objects to extension
    // physics.  Linearity in density gives an independent check of both
    // public operations on every rank.
    const auto physicalState   = problem.GetManifest().stateView(projected.values);
    const auto physicalDensity = physicalState.block(utils::blocks::density_field.mass_term);
    const stellar::DensityVolumeIntegralContext<integral::FixedTotalMass> densityIntegral{
        problem.GetPhysicalOperator()
    };
    const double integratedMass = densityIntegral.integrateDensity(physicalDensity).value();
    const double linearizedMass = densityIntegral.linearizeDensityIntegral(physicalDensity).value();
    const double preparedMass   = problem.GetPhysicalOperator().GetFixedMassReport().achieved;
    check_rank_consistent_scalar(integratedMass, communicator);
    check_rank_consistent_scalar(linearizedMass, communicator);
    const double massComparisonScale = std::max({std::abs(integratedMass), std::abs(preparedMass), 1.0e-300});
    CHECK(std::abs(integratedMass - preparedMass) / massComparisonScale <= 3.0e-13);
    CHECK(std::abs(linearizedMass - preparedMass) / massComparisonScale <= 3.0e-13);

    mfem::Vector normalizedResidual;
    normalized.BuildResidual(normalizedResidual);
    require_all_ranks(
        normalizedResidual.Size() == problem.EquationSize(), communicator, "normalized variadic residual size"
    );

    auto residualView                      = problem.GetManifest().residualView(normalizedResidual);
    const auto &layout                     = problem.GetManifest().layout();
    const auto &residualFactors            = normalized.GetNormalization().ResidualFactors();
    const auto massReport                  = problem.GetPreparedOperator().GetFixedMassReport();
    const auto angularMomentumReport       = problem.GetPreparedOperator().GetAngularMomentumReport();
    const auto centralDensityReport        = problem.GetPreparedOperator().GetCentralDensityReport();
    const auto checkReportedScalarResidual = [&](const mfem::Vector &block, const double expected) {
        require_all_ranks(block.Size() == 1, communicator, "reported scalar residual block size");
        const double actual = block(0);
        CAPTURE(actual, expected);
        CHECK(std::isfinite(expected));
        CHECK(std::abs(actual - expected) <= 5.0e-13 * std::max({1.0, std::abs(actual), std::abs(expected)}));
        check_rank_consistent_scalar(actual, communicator);
    };
    checkReportedScalarResidual(
        residualView.block(utils::blocks::fixed_total_mass_constraint.mass_normalization_term),
        massReport.dimensionalResidual * residualFactors(layout.offset(massResidualBlock))
    );
    checkReportedScalarResidual(
        residualView.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term),
        angularMomentumReport.dimensionalResidual * residualFactors(layout.offset(angularMomentumResidualBlock))
    );
    checkReportedScalarResidual(
        residualView.block(utils::blocks::fixed_central_density_phase.central_value_term),
        centralDensityReport.enthalpyResidual * residualFactors(layout.offset(centralDensityResidualBlock))
    );

    mfem::Vector normalizedDirection = make_deterministic_vector(problem.StateSize(), 0.37);
    const auto directionState        = problem.GetManifest().stateView(normalizedDirection);
    mfem::Vector massDirection =
        directionState.block(utils::blocks::fixed_total_mass_constraint.mass_normalization_term);
    mfem::Vector angularVelocityDirection =
        directionState.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term);
    mfem::Vector phaseDirection = directionState.block(utils::blocks::fixed_central_density_phase.central_value_term);
    require_all_ranks(
        massDirection.Size() == 1 && angularVelocityDirection.Size() == 1 && phaseDirection.Size() == 1, communicator,
        "generated scalar direction block sizes"
    );
    massDirection(0)            = 0.17;
    angularVelocityDirection(0) = -0.23;
    phaseDirection(0)           = 0.31;
    massDirection.SyncAliasMemory(normalizedDirection);
    angularVelocityDirection.SyncAliasMemory(normalizedDirection);
    phaseDirection.SyncAliasMemory(normalizedDirection);
    normalizedDirection /= global_norm(normalizedDirection, communicator);

    mfem::Vector normalizedAction;
    normalized.Mult(normalizedDirection, normalizedAction);
    require_all_ranks(
        normalizedAction.Size() == problem.EquationSize(), communicator, "normalized variadic Jacobian-action size"
    );

    /* Different dependency revisions are intentional: the state changes in
     * each difference evaluation, so the distributed physical contexts must
     * be rebuilt even though all persistent identities remain the same. */
    constexpr double differenceStep = 1.0e-5;
    mfem::Vector plusState(normalizedState);
    plusState.Add(differenceStep, normalizedDirection);
    const auto plusPreparation = normalized.Prepare(plusState, make_stellar_dependencies(2));
    require_all_ranks(plusPreparation.DidAnyWork(), communicator, "positive finite-difference preparation");
    mfem::Vector plusResidual;
    normalized.BuildResidual(plusResidual);

    mfem::Vector minusState(normalizedState);
    minusState.Add(-differenceStep, normalizedDirection);
    const auto minusPreparation = normalized.Prepare(minusState, make_stellar_dependencies(3));
    require_all_ranks(minusPreparation.DidAnyWork(), communicator, "negative finite-difference preparation");
    mfem::Vector minusResidual;
    normalized.BuildResidual(minusResidual);

    mfem::Vector finiteDifference(plusResidual);
    finiteDifference -= minusResidual;
    finiteDifference /= 2.0 * differenceStep;
    const auto restoredPreparation = normalized.Prepare(normalizedState, make_stellar_dependencies(4));
    require_all_ranks(
        restoredPreparation.DidAnyWork() && normalized.IsPrepared(), communicator,
        "restored finite-difference preparation"
    );
    mfem::Vector restoredResidual;
    normalized.BuildResidual(restoredResidual);
    mfem::Vector restoredResidualDifference(restoredResidual);
    restoredResidualDifference -= normalizedResidual;
    const double restoredResidualError =
        global_norm(restoredResidualDifference, communicator) /
        std::max(
            {global_norm(restoredResidual, communicator), global_norm(normalizedResidual, communicator),
             std::numeric_limits<double>::epsilon()}
        );
    CAPTURE(restoredResidualError);
    CHECK(restoredResidualError <= 2.0e-12);

    mfem::Vector finiteDifferenceError(normalizedAction);
    finiteDifferenceError -= finiteDifference;
    const double actionNorm           = global_norm(normalizedAction, communicator);
    const double finiteDifferenceNorm = global_norm(finiteDifference, communicator);
    const double completeDifferenceError =
        global_norm(finiteDifferenceError, communicator) /
        std::max({actionNorm, finiteDifferenceNorm, std::numeric_limits<double>::epsilon()});
    CAPTURE(actionNorm, finiteDifferenceNorm, completeDifferenceError);
    CHECK(actionNorm > std::numeric_limits<double>::min());
    CHECK(finiteDifferenceNorm > std::numeric_limits<double>::min());
    CHECK(completeDifferenceError <= 8.0e-5);

    const int locallyFinite =
        vector_is_finite(normalizedResidual) && vector_is_finite(normalizedAction) && vector_is_finite(finiteDifference)
            ? 1
            : 0;
    int globallyFinite = 0;
    REQUIRE(MPI_Allreduce(&locallyFinite, &globallyFinite, 1, MPI_INT, MPI_MIN, communicator) == MPI_SUCCESS);
    CHECK(globallyFinite == 1);
    CHECK(global_norm(normalizedResidual, communicator) > 0.0);
    CHECK(global_norm(normalizedAction, communicator) > 0.0);

    auto actionView           = problem.GetManifest().residualView(normalizedAction);
    auto finiteDifferenceView = problem.GetManifest().residualView(finiteDifference);
    const auto checkGlobalRow = [&](const auto &term, const char *rowName) {
        const mfem::Vector residualBlock   = residualView.block(term);
        const mfem::Vector actionBlock     = actionView.block(term);
        const mfem::Vector differenceBlock = finiteDifferenceView.block(term);
        require_all_ranks(
            residualBlock.Size() == 1 && actionBlock.Size() == 1 && differenceBlock.Size() == 1, communicator,
            "global scalar residual/Jacobian block sizes"
        );
        check_rank_consistent_scalar(residualBlock(0), communicator);
        check_rank_consistent_scalar(actionBlock(0), communicator);
        check_rank_consistent_scalar(differenceBlock(0), communicator);
        const double rowMagnitude       = std::max(std::abs(actionBlock(0)), std::abs(differenceBlock(0)));
        const double rowDifferenceError = std::abs(actionBlock(0) - differenceBlock(0)) /
                                          std::max(rowMagnitude, std::numeric_limits<double>::epsilon());
        CAPTURE(rowName, actionBlock(0), differenceBlock(0), rowMagnitude, rowDifferenceError);
        CHECK(rowMagnitude > 1.0e-10);
        CHECK(rowDifferenceError <= 8.0e-5);
    };
    checkGlobalRow(utils::blocks::fixed_total_mass_constraint.mass_normalization_term, "fixed-total-mass");
    checkGlobalRow(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term, "fixed-angular-momentum");
    checkGlobalRow(utils::blocks::fixed_central_density_phase.central_value_term, "fixed-central-density phase");

    auto component = preconditioning::makePreconditioner(problem);
    STATIC_CHECK(decltype(component)::borderValueArity == 3);
    STATIC_CHECK(decltype(component)::borderResidualArity == 3);
    auto physicalInverse = preconditioning::prepare(problem, component);
    require_all_ranks(physicalInverse.IsCurrent(), communicator, "prepared physical preconditioner currentness");
    auto scaledInverse = normalized.MakeScaledPreconditioner(physicalInverse);
    require_all_ranks(scaledInverse.IsCurrent(), communicator, "prepared normalized preconditioner currentness");

    mfem::Vector normalizedRightHandSide = make_deterministic_vector(problem.EquationSize(), 0.73);
    auto rightHandSideView               = problem.GetManifest().residualView(normalizedRightHandSide);
    mfem::Vector massRightHandSide =
        rightHandSideView.block(utils::blocks::fixed_total_mass_constraint.mass_normalization_term);
    mfem::Vector angularMomentumRightHandSide =
        rightHandSideView.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term);
    mfem::Vector phaseRightHandSide =
        rightHandSideView.block(utils::blocks::fixed_central_density_phase.central_value_term);
    massRightHandSide(0)            = 0.11;
    angularMomentumRightHandSide(0) = -0.19;
    phaseRightHandSide(0)           = 0.29;
    massRightHandSide.SyncAliasMemory(normalizedRightHandSide);
    angularMomentumRightHandSide.SyncAliasMemory(normalizedRightHandSide);
    phaseRightHandSide.SyncAliasMemory(normalizedRightHandSide);
    normalizedRightHandSide /= global_norm(normalizedRightHandSide, communicator);

    mfem::Vector firstCorrection(scaledInverse.Height());
    mfem::Vector repeatedCorrection(scaledInverse.Height());
    firstCorrection    = 0.0;
    repeatedCorrection = 0.0;
    scaledInverse.Mult(normalizedRightHandSide, firstCorrection);
    scaledInverse.Mult(normalizedRightHandSide, repeatedCorrection);

    mfem::Vector repeatDifference(repeatedCorrection);
    repeatDifference -= firstCorrection;
    const double correctionNorm = global_norm(firstCorrection, communicator);
    const double repeatError =
        global_norm(repeatDifference, communicator) / std::max(correctionNorm, std::numeric_limits<double>::epsilon());
    CAPTURE(correctionNorm, repeatError);
    const int locallyFiniteCorrections =
        vector_is_finite(firstCorrection) && vector_is_finite(repeatedCorrection) ? 1 : 0;
    int globallyFiniteCorrections = 0;
    REQUIRE(
        MPI_Allreduce(&locallyFiniteCorrections, &globallyFiniteCorrections, 1, MPI_INT, MPI_MIN, communicator) ==
        MPI_SUCCESS
    );
    CHECK(globallyFiniteCorrections == 1);
    CHECK(std::isfinite(correctionNorm));
    CHECK(correctionNorm > 0.0);
    CHECK(repeatError <= 2.0e-12);
    require_all_ranks(
        physicalInverse.IsCurrent() && scaledInverse.IsCurrent(), communicator,
        "preconditioner currentness after repeated application"
    );

    const mfem::Vector &readOnlyCorrection = firstCorrection;
    const auto correctionView              = problem.GetManifest().stateView(readOnlyCorrection);
    const double massCorrection =
        correctionView.block(utils::blocks::fixed_total_mass_constraint.mass_normalization_term)(0);
    const double angularVelocityCorrection =
        correctionView.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term)(0);
    const double phaseCorrection =
        correctionView.block(utils::blocks::fixed_central_density_phase.central_value_term)(0);
    CAPTURE(massCorrection, angularVelocityCorrection, phaseCorrection);
    check_rank_consistent_scalar(massCorrection, communicator);
    check_rank_consistent_scalar(angularVelocityCorrection, communicator);
    check_rank_consistent_scalar(phaseCorrection, communicator);

    /* A Newton iteration reparses the same normalized state under fresh
     * dependency revisions.  Every rank must observe the stale inverse, and
     * refresh must reconstruct the inferred border actions and Schur data. */
    const auto secondPreparation = normalized.Prepare(normalizedState, make_stellar_dependencies(5));
    require_all_ranks(
        secondPreparation.DidAnyWork() && normalized.IsPrepared(), communicator,
        "second variadic stellar-root preparation"
    );
    mfem::Vector secondPreparedResidual;
    normalized.BuildResidual(secondPreparedResidual);
    mfem::Vector secondPreparedResidualDifference(secondPreparedResidual);
    secondPreparedResidualDifference -= normalizedResidual;
    const double secondPreparedResidualError =
        global_norm(secondPreparedResidualDifference, communicator) /
        std::max(
            {global_norm(secondPreparedResidual, communicator), global_norm(normalizedResidual, communicator),
             std::numeric_limits<double>::epsilon()}
        );
    CAPTURE(secondPreparedResidualError);
    CHECK(secondPreparedResidualError <= 2.0e-12);
    require_all_ranks(
        !physicalInverse.IsCurrent() && !scaledInverse.IsCurrent(), communicator,
        "preconditioners become stale together"
    );
    CHECK_THROWS_AS(scaledInverse.Mult(normalizedRightHandSide, repeatedCorrection), std::logic_error);

    const auto refresh = physicalInverse.Refresh();
    CHECK(refresh.specificationActionsRefreshed);
    CHECK(refresh.rebuiltSchurComplement);
    CHECK(refresh.DidAnyWork());
    require_all_ranks(
        physicalInverse.IsCurrent() && scaledInverse.IsCurrent(), communicator, "refreshed preconditioner currentness"
    );
    const auto noOpRefresh = physicalInverse.Refresh();
    CHECK_FALSE(noOpRefresh.DidAnyWork());

    mfem::Vector refreshedCorrection(scaledInverse.Height());
    refreshedCorrection = 0.0;
    scaledInverse.Mult(normalizedRightHandSide, refreshedCorrection);
    mfem::Vector refreshDifference(refreshedCorrection);
    refreshDifference -= firstCorrection;
    const double refreshError =
        global_norm(refreshDifference, communicator) / std::max(correctionNorm, std::numeric_limits<double>::epsilon());
    const int locallyFiniteRefresh = vector_is_finite(refreshedCorrection) ? 1 : 0;
    int globallyFiniteRefresh      = 0;
    REQUIRE(
        MPI_Allreduce(&locallyFiniteRefresh, &globallyFiniteRefresh, 1, MPI_INT, MPI_MIN, communicator) == MPI_SUCCESS
    );
    CAPTURE(refreshError);
    CHECK(globallyFiniteRefresh == 1);
    CHECK(refreshError <= 2.0e-10);
}

TEST_CASE(
    "MPI Prepared Gravity Operators Preserve Global Algebraic Identities",
    "[mpi][distributed][gravity][operators][unit]"
) {
    const auto args        = test_utils::setup_args();
    mean_field::fem::FEM f = mean_field::fem::setup_fem(args.mesh_file, args, 0);

    using GeometryContext  = mean_field::operators::context::gravity_field::GravityFieldGeometryContext;
    GeometryContext geometry_context(f, *f.domainMapperStateless);

    mfem::Vector displacement_true(f.displacementFes->GetTrueVSize());
    displacement_true               = 0.0;
    const mfem::Vector displacement = geometry_context.GetDisplacementMap().gather(displacement_true);
    geometry_context.PreparePrimal(displacement, {0}, {0});

    const mfem::Operator &mass = geometry_context.GetMassOperator();
    const mfem::Vector first   = make_deterministic_vector(mass.Width(), 0.17);
    const mfem::Vector second  = make_deterministic_vector(mass.Width(), 0.83);
    mfem::Vector combination(first);
    combination *= 1.7;
    combination.Add(-0.4, second);

    mfem::Vector first_action;
    mfem::Vector second_action;
    mfem::Vector combination_action;
    mass.Mult(first, first_action);
    mass.Mult(second, second_action);
    mass.Mult(combination, combination_action);

    mfem::Vector expected_combination(first_action);
    expected_combination *= 1.7;
    expected_combination.Add(-0.4, second_action);
    mfem::Vector linearity_difference(combination_action);
    linearity_difference -= expected_combination;

    const MPI_Comm communicator = f.mesh->GetComm();
    const double symmetry_scale = std::max(
        {std::abs(global_dot(first, second_action, communicator)),
         std::abs(global_dot(second, first_action, communicator)), std::numeric_limits<double>::epsilon()}
    );
    const double symmetry_error =
        std::abs(global_dot(first, second_action, communicator) - global_dot(second, first_action, communicator)) /
        symmetry_scale;
    const double linearity_error =
        global_norm(linearity_difference, communicator) /
        std::max(global_norm(expected_combination, communicator), std::numeric_limits<double>::epsilon());

    CHECK(symmetry_error <= 2.0e-12);
    CHECK(linearity_error <= 2.0e-12);

    const mfem::Operator &divergence           = geometry_context.GetDivergenceOperator();
    const mfem::Operator &transpose_divergence = geometry_context.GetTransposeDivergenceOperator();
    const mfem::Vector flux                    = make_deterministic_vector(divergence.Width(), 0.41);
    const mfem::Vector potential               = make_deterministic_vector(divergence.Height(), 0.67);
    mfem::Vector divergence_action;
    mfem::Vector transpose_action;
    divergence.Mult(flux, divergence_action);
    transpose_divergence.Mult(potential, transpose_action);

    const double forward_product   = global_dot(potential, divergence_action, communicator);
    const double transpose_product = global_dot(flux, transpose_action, communicator);
    const double adjoint_scale =
        std::max({std::abs(forward_product), std::abs(transpose_product), std::numeric_limits<double>::epsilon()});
    const double adjoint_error = std::abs(forward_product - transpose_product) / adjoint_scale;

    CHECK(adjoint_error <= 2.0e-12);
}

TEST_CASE(
    "MPI Coupled Gravity LDU Is Stationary Linear And Does Not Reprepare Geometry",
    "[mpi][distributed][gravity][preconditioning][integration]"
) {
    const auto args        = test_utils::setup_args();
    mean_field::fem::FEM f = mean_field::fem::setup_fem(args.mesh_file, args, 0);

    using GeometryContext  = mean_field::operators::context::gravity_field::GravityFieldGeometryContext;
    GeometryContext geometryContext(f, *f.domainMapperStateless);

    mfem::Vector displacementTrue(f.displacementFes->GetTrueVSize());
    displacementTrue                = 0.0;
    const mfem::Vector displacement = geometryContext.GetDisplacementMap().gather(displacementTrue);
    geometryContext.PreparePrimal(displacement, {.value = 1}, {.value = 1});

    namespace backend         = mean_field::preconditioning::backend;
    namespace preconditioning = mean_field::preconditioning;
    const auto block          = preconditioning::GravityFieldBlock(
        backend::Diagonal{}, backend::HypreBoomerAMG{backend::FixedCycles{.cycles = 1}},
        preconditioning::GravityApproximateLDU{}
    );
    auto prepared             = preconditioning::prepare(f, geometryContext, block);

    const mfem::Vector first  = make_deterministic_vector(prepared.Width(), 0.23);
    const mfem::Vector second = make_deterministic_vector(prepared.Width(), 0.79);
    mfem::Vector combined(first);
    combined *= 1.3;
    combined.Add(-0.45, second);

    mfem::Vector firstAction(prepared.Height());
    mfem::Vector secondAction(prepared.Height());
    mfem::Vector combinedAction(prepared.Height());
    mfem::Vector repeatedAction(prepared.Height());
    firstAction                            = 0.0;
    secondAction                           = 0.0;
    combinedAction                         = 0.0;
    repeatedAction                         = 0.0;

    const std::uint64_t massPreparations   = geometryContext.GetMassOperator().GetPreparationCount();
    const std::uint64_t sourcePreparations = geometryContext.GetSourceOperator().GetPreparationCount();
    double *const combinedStorage          = combinedAction.GetData();

    prepared.Mult(first, firstAction);
    prepared.Mult(second, secondAction);
    prepared.Mult(combined, combinedAction);
    prepared.Mult(first, repeatedAction);

    mfem::Vector expectedCombined(firstAction);
    expectedCombined *= 1.3;
    expectedCombined.Add(-0.45, secondAction);

    const MPI_Comm communicator = f.mesh->GetComm();
    mfem::Vector linearityDifference(combinedAction);
    linearityDifference -= expectedCombined;
    mfem::Vector determinismDifference(repeatedAction);
    determinismDifference -= firstAction;
    const double linearityError =
        global_norm(linearityDifference, communicator) /
        std::max(global_norm(expectedCombined, communicator), std::numeric_limits<double>::epsilon());

    CHECK(vector_is_finite(combinedAction));
    CHECK(linearityError <= 5.0e-12);
    CHECK(global_norm(determinismDifference, communicator) <= 5.0e-14);
    CHECK(combinedAction.GetData() == combinedStorage);
    CHECK(geometryContext.GetMassOperator().GetPreparationCount() == massPreparations);
    CHECK(geometryContext.GetSourceOperator().GetPreparationCount() == sourcePreparations);
    CHECK(prepared.GetFactorization().GetStatistics().applications == 4);
}

TEST_CASE(
    "MPI Gravity Analysis And Solve Produce Finite Distributed Fields",
    "[mpi][distributed][gravity][integration]"
) {
    auto args              = test_utils::setup_args();
    mean_field::fem::FEM f = mean_field::fem::setup_fem(args.mesh_file, args, 0);

    *f.displacement        = 0.0;

    using DomainSchema     = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema;
    mfem::Vector attribute_density(f.smesh.mesh->attributes.Max());
    attribute_density = 0.0;
    for (int index = 0; index < f.smesh.mesh->attributes.Size(); ++index) {
        const int attribute = f.smesh.mesh->attributes[index];
        if (DomainSchema::template attribute_belongs_to<mean_field::utils::domain::Stellar>(attribute)) {
            attribute_density(attribute - 1) = 1.0;
        }
    }

    mfem::PWConstCoefficient density_coefficient(attribute_density);
    mfem::ParGridFunction density(f.densityFes.get());
    density.ProjectCoefficient(density_coefficient);
    mean_field::analysis::conserve_mass(f, density, mean_field::utils::MASS);

    const double integrated_mass = mean_field::analysis::domain_integrate_grid_function(
        f, density, mean_field::utils::DOMAINS::STELLAR, mean_field::mapping::COORDINATE_SPACE::PHYSICAL
    );
    f.com = mean_field::analysis::get_com(f, density);
    f.Q   = mean_field::physics::compute_quadrupole_moment_tensor(f, density, f.com);

    const mean_field::physics::GravitySolution solution = mean_field::physics::solve_gravity_field(
        f,
        mean_field::physics::GravitySolveOptions{
            .relativeTolerance = 1.0e-12, .absoluteTolerance = 1.0e-15, .maximumIterations = 1000
        },
        density, *f.displacement
    );

    mfem::Vector flux_true;
    mfem::Vector potential_true;
    solution.gradPhi.GetTrueDofs(flux_true);
    solution.phi.GetTrueDofs(potential_true);

    const int local_finite = vector_is_finite(flux_true) && vector_is_finite(potential_true) ? 1 : 0;
    int globally_finite    = 0;
    MPI_Allreduce(&local_finite, &globally_finite, 1, MPI_INT, MPI_MIN, f.mesh->GetComm());

    const double local_norms[2]{flux_true * flux_true, potential_true * potential_true};
    double global_norms[2]{};
    MPI_Allreduce(local_norms, global_norms, 2, MPI_DOUBLE, MPI_SUM, f.mesh->GetComm());

    CHECK(globally_finite == 1);
    CHECK(std::abs(integrated_mass - mean_field::utils::MASS) <= 1.0e-12 * mean_field::utils::MASS);
    CHECK(global_norms[0] > std::numeric_limits<double>::min());
    CHECK(global_norms[1] > std::numeric_limits<double>::min());
}