#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include <memory>
#include <numbers>
#include <tuple>
#include <utility>
#include <variant>

#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>

#include "mfem/linalg/solvers.hpp"
#include "mfem/general/communication.hpp"

#include "serif/equations/angular_momentum.hpp"
#include "serif/equations/barotropic_density.hpp"
#include "serif/equations/primal_gravity.hpp"
#include "serif/equations/surface_potential.hpp"
#include "serif/equations/mass_normalization.hpp"
#include "serif/backend/mfem/model/equilibrium.hpp"
#include "serif/model/eos_reference.hpp"
#include "serif/solver/preconditioner/default_preconditioners.hpp"
#include "serif/tests/deformation/mesh.hpp"
#include "serif/tests/equations/gravity_fixture.hpp"

// MeanField's prepared_stellar_equilibrium.cpp checks the complete moving root
// against its children and centered differences. Keep both checks: child parity
// verifies composition, while finite differences verify the assembled derivative.
namespace {
    namespace model = serif::model;
    namespace backend = model::backend::mfem;
    namespace execution = serif::operators::backend::mfem;
    namespace dependencies = serif::utils::dependencies;
    namespace equations = serif::equations;
    namespace deformation = serif::geometry::deformation;
    namespace surface = serif::geometry::surface;
    namespace blocks = serif::discretization::blocks;
    namespace tests = serif::tests::equations;
    using Polytrope = serif::eos::models::Polytrope;
    using Density = blocks::density::mass;
    using Potential = blocks::gravity::potential;
    using Constant = blocks::barotropic_constant::mass_normalization;
    using Omega = blocks::rotation::angular_velocity;
    using Shape = blocks::surface_deformation::coefficients;
    using ShapeRow = blocks::surface_deformation::shape_equilibrium;
    using Displacement = blocks::displacement::vector;
    using Kind = blocks::RepresentedValueKind;
    using Vector = mfem::Vector;
    using InteriorRows = dependencies::KeyList<Density, Potential, Constant, Omega>;

    /** @brief Own paired stellar/vacuum meshes and all resources borrowed by the public model binding. */
    struct Fixture {
        stroid::StroidMesh generated{serif::tests::deformation::make_mesh(1)};
        tests::gravity::GravityFixture fields;
        mfem::FiniteElementSpace reference_displacement{&fields.mesh, &fields.enthalpy_collection, 3, mfem::Ordering::byNODES};
        mfem::FiniteElementSpace compactification_space{&fields.mesh, generated.exterior_coordinate->space->FEColl()};
        mfem::GridFunction compactification{&compactification_space};

        /** @brief Preserve the generator's exterior coordinate on the identically numbered mesh copy. */
        explicit Fixture(const mfem::Ordering::Type ordering) : fields(mfem::Mesh(*generated.mesh), ordering) {
            REQUIRE(compactification.Size() == generated.exterior_coordinate->values->Size());
            compactification = *generated.exterior_coordinate->values;
        }
    };

    /** @brief Exercise all stellar equations at a valid, deliberately unconverged rotating state. */
    [[nodiscard]] auto physical_model(const double index = 1.5) {
        const auto definition = [](const auto& physics) {
            const auto eos = model::borrow_eos(physics.equation_of_state());
            const auto& rotation = physics.template specification<model::FixedAngularMomentum>();
            const auto& axis = rotation.axis();
            const auto& center = rotation.center();
            return model::EquilibriumEquations{
                .interior_equations = equations::EquationRegistry{
                    equations::BarotropicDensityBalance{eos, serif::dimensions::PressureValue{0.025}, axis, center}, equations::PrimalPoissonBalance{2.3},
                    equations::MassNormalization{physics.template specification<model::FixedTotalMass>().target().value()},
                    equations::AngularMomentumConstraint{rotation.target().value(), axis, center}
                },
                .shape_force = equations::SurfacePotentialBalance{axis, center},
                .surface_condition = surface::Isobaric{serif::dimensions::PressureValue{0.025}}
            };
        };
        return model::make_stellar_model(Polytrope{index, 0.25}, definition,
            model::FixedTotalMass{serif::dimensions::MassValue{1.1}},
            model::FixedAngularMomentum{{.angular_momentum = serif::dimensions::AngularMomentumValue{0.14},
                .axis = {0.2, -0.3, 1.0}, .center = {0.03, -0.02, 0.01}}});
    }

    /** @brief Use the production radial lift and vacuum extension with nine independent surface harmonics. */
    [[nodiscard]] auto discretize(const std::shared_ptr<Fixture>& fixture) {
        const auto& fields = fixture->fields;
        const auto bindings = dependencies::make_dependency_context(
            dependencies::borrow<Density>(fields.density_binding), dependencies::borrow<Potential>(fields.potential_binding),
            dependencies::borrow<Constant>(fields.constant_binding), dependencies::borrow<Omega>(fields.angular_velocity_binding),
            dependencies::borrow<Displacement>(fields.displacement_binding));
        const auto reference = deformation::backend::mfem::deformation_discretization(
            fields.enthalpy_space, fixture->reference_displacement, *fixture->generated.reference_mesh);
        const deformation::DomainDeformation<deformation::SphericalHarmonicSurface> geometry{
            .surface = {.maximum_degree = 2}, .stellar_interior = {}, .vacuum = {}
        };
        return backend::with_deformation(backend::make_stellar_discretization(fixture, bindings, fixture->compactification, fields.mapper), geometry, reference);
    }

    /** @brief Populate all physical fields, scalar constraints and surface modes in the moving-star fixture. */
    template <typename Layout, typename ShapeCoefficients>
    [[nodiscard]] model::StellarState<Layout> initial_state(const Layout& layout, tests::gravity::GravityFixture& fields, const ShapeCoefficients& shape) {
        model::StellarState state(layout);
        const auto initialize = [&]<typename Term>(const Vector& values) {
            const auto destination = state.template get<Term>();
            REQUIRE(destination.size() == static_cast<std::size_t>(values.Size()));
            std::copy_n(values.GetData(), values.Size(), destination.begin());
        };
        initialize.template operator()<Density>(fields.density());
        initialize.template operator()<Potential>(fields.potential());
        state.template get<Constant>()[0] = 0.72;
        state.template get<Omega>()[0] = 0.31;
        std::ranges::copy(shape, state.template get<Shape>().begin());
        return state;
    }

    /** @brief Make an independent physical candidate without sharing any mutable field storage. */
    template <typename Layout>
    [[nodiscard]] model::StellarState<Layout> perturbed_state(const model::StellarState<Layout>& original, const double step, const Vector& direction) {
        REQUIRE(original.values().size() == static_cast<std::size_t>(direction.Size()));
        model::StellarState result(original);
        for (int index = 0; index < direction.Size(); ++index) {
            result.values()[static_cast<std::size_t>(index)] += step * direction[index];
        }
        return result;
    }

    /** @brief Copy a frozen residual before evaluating another candidate. */
    template <typename System>
    [[nodiscard]] Vector residual(const System& system) {
        Vector result(static_cast<int>(system.residual().size()));
        std::copy(system.residual().begin(), system.residual().end(), result.GetData());
        return result;
    }

    /** @brief Check each equation separately so large gravity rows cannot hide missing scalar or shape rows. */
    template <typename Layout>
    void check_rows(const Vector& actual, const Vector& expected, const Layout& layout, double relative_tolerance, double absolute_tolerance) {
        REQUIRE(actual.Size() == expected.Size());
        int offset = 0;
        for (std::size_t row = 0; row < layout.residual_sizes.size(); ++row) {
            CAPTURE(row);
            double error_squared = 0.0, reference_squared = 0.0;
            for (std::size_t entry = 0; entry < layout.residual_sizes[row]; ++entry) {
                const double difference = actual[offset] - expected[offset];
                error_squared += difference * difference;
                reference_squared += expected[offset] * expected[offset];
                ++offset;
            }
            CHECK(std::sqrt(error_squared) <= absolute_tolerance + relative_tolerance * std::sqrt(reference_squared));
        }
    }

    /** @brief Compare cached scalar functionals with the independently retained moving volume actions. */
    template <typename Row, typename Evaluator>
    void check_scalar_row(const Evaluator& evaluator, const backend::PhysicalSystem& physical) {
        const auto& layout = evaluator.layout();
        const auto& cache = evaluator.cache();
        const auto volume = cache.template share<typename execution::volume_keys::Scoped<Row>::Linearization>();
        const auto& shape = cache.template get<model::computations::ShapeDerivatives>();
        const auto reference = backend::volume_row_action(volume, layout, backend::ColumnContribution{
            static_cast<int>(layout.template offset<Kind::value, Shape>()), shape.tangent});
        const auto& cached = cache.template get<model::computations::ScalarResidualRow<Row>>();
        const Vector direction = tests::deterministic(reference.action()->Width(), 0.47);
        const Vector dual = tests::constant(1, 0.37);
        Vector expected(1), actual(1), expected_transpose, actual_transpose(direction.Size());
        reference.action()->Mult(direction, expected);
        cached.action()->Mult(direction, actual);
        CHECK(actual[0] == Catch::Approx(expected[0]).margin(2.0e-12));
        reference.action()->MultTranspose(dual, expected_transpose);
        cached.action()->MultTranspose(dual, actual_transpose);
        CHECK(tests::relative_error(actual_transpose, expected_transpose) < 2.0e-12);
        const auto check_column = [&]<typename Column>() {
            const int offset = static_cast<int>(layout.template offset<Kind::value, Column>());
            const int size = static_cast<int>(layout.template size<Kind::value, Column>());
            const auto selected = physical.block(static_cast<int>(layout.template offset<Kind::residual, Row>()), 1, offset, size);
            Vector input(size), selected_transpose;
            std::copy_n(direction.GetData() + offset, size, input.GetData());
            Vector complete(direction.Size());
            complete = 0.0;
            std::copy_n(input.GetData(), size, complete.GetData() + offset);
            cached.action()->Mult(complete, expected);
            selected.action()->Mult(input, actual);
            CHECK(actual[0] == Catch::Approx(expected[0]).margin(2.0e-12));
            selected.action()->MultTranspose(dual, selected_transpose);
            Vector expected_column(expected_transpose.GetData() + offset, size);
            CHECK(tests::relative_error(selected_transpose, expected_column) < 2.0e-12);
        };
        check_column.template operator()<Density>();
        check_column.template operator()<Shape>();
        check_column.template operator()<Omega>();
    }

    /** @brief Place the volume and surface children, then independently enforce both physical traces. */
    template <typename Evaluator, typename Layout>
    void check_composition(const Evaluator& evaluator, const model::StellarState<Layout>& state,
        const backend::PhysicalSystem& physical, const serif::discretization::backend::mfem::DofSelection& boundary) {
        const auto& layout = evaluator.layout();
        const auto& cache = evaluator.cache();
        Vector expected(static_cast<int>(layout.template total_size<Kind::residual>()));
        expected = 0.0;
        [&]<typename... Rows>(dependencies::KeyList<Rows...>) {
            const auto place = [&]<typename Row>() {
                const auto& child = cache.template get<typename execution::volume_keys::Scoped<Row>::Linearization>();
                REQUIRE(static_cast<std::size_t>(child.residual.Size()) == layout.template size<Kind::residual, Row>());
                std::copy_n(child.residual.GetData(), child.residual.Size(), expected.GetData() + layout.template offset<Kind::residual, Row>());
            };
            (place.template operator()<Rows>(), ...);
        }(InteriorRows{});
        const auto& surface_child = cache.template get<model::computations::ShapeContribution>();
        REQUIRE(surface_child.residual.size() == 1);
        for (const auto& contribution : surface_child.residual) {
            std::copy(contribution.values.begin(), contribution.values.end(), expected.GetData() + contribution.row_offset);
        }
        const int shape_row = static_cast<int>(layout.template offset<Kind::residual, ShapeRow>());
        const int shape_size = static_cast<int>(layout.template size<Kind::residual, ShapeRow>());
        Vector shape_residual(expected.GetData() + shape_row, shape_size);
        REQUIRE(shape_residual.Norml2() > 1.0e-8);

        // Compare each retained surface child coupling with the corresponding assembled block.
        for (const auto& contribution : surface_child.jacobian) {
            const auto selected = physical.block(contribution.row_offset, contribution.block.action()->Height(),
                contribution.column_offset, contribution.block.action()->Width());
            const Vector direction = tests::deterministic(contribution.block.action()->Width(), 0.47);
            const Vector dual = tests::deterministic(contribution.block.action()->Height(), 0.71);
            Vector actual(contribution.block.action()->Height()), reference(contribution.block.action()->Height());
            Vector actual_transpose(contribution.block.action()->Width()), reference_transpose(contribution.block.action()->Width());
            selected.action()->Mult(direction, actual);
            contribution.block.action()->Mult(direction, reference);
            CHECK(tests::relative_error(actual, reference) < 2.0e-13);
            selected.action()->MultTranspose(dual, actual_transpose);
            contribution.block.action()->MultTranspose(dual, reference_transpose);
            CHECK(tests::relative_error(actual_transpose, reference_transpose) < 2.0e-13);
        }

        const auto density = state.template get<Density>();
        REQUIRE(boundary.size() > 0);
        REQUIRE(boundary.size() < static_cast<int>(density.size()));
        const double surface_density = serif::eos::evaluate<serif::dimensions::MassDensity>(
            evaluator.model().equation_of_state(), serif::dimensions::PressureValue{0.025}).value();
        const int density_row = static_cast<int>(layout.template offset<Kind::residual, Density>());
        for (const int row : boundary.reduced_dofs()) {
            expected[density_row + row] = density[row] - surface_density;
        }
        const auto infinity = serif::discretization::backend::mfem::select_boundary_dofs<
            serif::discretization::domain::InfinitySurfaceBoundary>(
                evaluator.binding().discretization->template binding<Potential>());
        REQUIRE(infinity.size() > 0);
        const int potential_row = static_cast<int>(layout.template offset<Kind::residual, Potential>());
        for (const int row : infinity.reduced_dofs()) {
            expected[potential_row + row] = state.template get<Potential>()[row];
        }
        check_rows(residual(physical), expected, layout, 2.0e-13, 2.0e-13);

        // Each strong trace must reach its own scalar field alone, including in the transpose.
        const auto check_trace = [&]<typename Field>(const auto& selected) {
            Vector dual(expected.Size()), transpose;
            dual = 0.0;
            Vector expected_transpose(physical.jacobian().action()->Width());
            expected_transpose = 0.0;
            const int row_offset = static_cast<int>(layout.template offset<Kind::residual, Field>());
            const int column_offset = static_cast<int>(layout.template offset<Kind::value, Field>());
            for (const int row : selected.reduced_dofs()) {
                dual[row_offset + row] = 0.17 + 0.01 * (row % 7);
                expected_transpose[column_offset + row] = dual[row_offset + row];
            }
            physical.jacobian().action()->MultTranspose(dual, transpose);
            CHECK(tests::relative_error(transpose, expected_transpose) < 2.0e-13);
        };
        check_trace.template operator()<Density>(boundary);
        check_trace.template operator()<Potential>(infinity);
    }

    /** @brief Compare every physical state column and a mixed perturbation through the public evaluation path. */
    template <typename Evaluator, typename Layout>
    void check_derivatives(Evaluator& evaluator, const model::StellarState<Layout>& state) {
        const auto& layout = evaluator.layout();
        const auto evaluation = evaluator.evaluate(state);
        REQUIRE(evaluation);
        const auto frozen = *evaluation;
        const Vector direction_seed = tests::deterministic(static_cast<int>(state.values().size()), 0.43);
        const Vector dual = tests::deterministic(frozen->jacobian().action()->Height(), 1.17);
        Vector transpose, block_transpose;
        frozen->jacobian().action()->MultTranspose(dual, transpose);
        const auto blocks = frozen->block(0, dual.Size(), 0, static_cast<int>(state.values().size()));
        blocks.action()->MultTranspose(dual, block_transpose);
        CHECK(tests::relative_error(transpose, block_transpose) < 2.0e-12);

        // The final iteration moves all physical fields and the surface simultaneously.
        int offset = 0;
        for (std::size_t column = 0; column <= layout.state_sizes.size(); ++column) {
            CAPTURE(column);
            Vector direction(direction_seed);
            if (column < layout.state_sizes.size()) {
                direction = 0.0;
                std::copy_n(direction_seed.GetData() + offset, layout.state_sizes[column], direction.GetData() + offset);
                offset += static_cast<int>(layout.state_sizes[column]);
            }
            Vector action, block_action;
            frozen->jacobian().action()->Mult(direction, action);
            blocks.action()->Mult(direction, block_action);
            REQUIRE(action.Norml2() > 1.0e-10);
            check_rows(action, block_action, layout, 2.0e-12, 2.0e-12);
            CHECK((action * dual) == Catch::Approx(direction * transpose).epsilon(2.0e-12).margin(2.0e-12));

            constexpr double step = 1.0e-5;
            const auto plus = evaluator.evaluate(perturbed_state(state, step, direction));
            REQUIRE(plus);
            Vector difference = residual(**plus);
            const auto minus = evaluator.evaluate(perturbed_state(state, -step, direction));
            REQUIRE(minus);
            difference -= residual(**minus);
            difference /= 2.0 * step;
            check_rows(action, difference, layout, 2.0e-6, 2.0e-9);

            if (column == Layout::template index<Kind::value, Shape>()) {
                // Shape must affect each interior equation and the surface potential, not only its own row.
                int row_offset = 0;
                for (std::size_t row = 0; row < layout.residual_sizes.size(); ++row) {
                    Vector row_action(action.GetData() + row_offset, static_cast<int>(layout.residual_sizes[row]));
                    CHECK(row_action.Norml2() > 1.0e-10);
                    row_offset += row_action.Size();
                }
            }
        }
        const auto restored = evaluator.evaluate(state);
        REQUIRE(restored);
        check_rows(residual(**restored), residual(*frozen), layout, 0.0, 0.0);
    }

    /** @brief Check the assembled solver coordinates, including strong rows, scalar constraints and native surface modes. */
    template <typename Evaluator, typename Layout>
    void check_normalized_system(Evaluator& evaluator, const model::StellarState<Layout>& state,
        const serif::discretization::backend::mfem::DofSelection& boundary) {
        namespace normalization = serif::solver::normalization;
        namespace dimensions = serif::dimensions;
        const auto evaluation = evaluator.evaluate(state);
        REQUIRE(evaluation);
        const auto frozen = *evaluation;
        REQUIRE(frozen->physical() == evaluator.cache().template get<model::keys::PhysicalSystem>());
        const auto& layout = evaluator.layout();
        const auto& coordinates = frozen->coordinates();
        const auto& scales = frozen->reference_scales();
        CHECK(scales.mass.value() == Catch::Approx(1.1));
        CHECK(scales.radius.value() > 0.0);
        CHECK(scales.radius.value() != 1.0);
        REQUIRE(coordinates.state_size() == state.values().size());
        REQUIRE(coordinates.residual_size() == frozen->residual().size());
        CHECK(coordinates.state_factors().front() != 1.0);
        Vector expected_residual(static_cast<int>(frozen->residual().size()));
        coordinates.normalize_residual(frozen->physical()->residual(), {expected_residual.GetData(), frozen->residual().size()});
        check_rows(residual(*frozen), expected_residual, layout, 0.0, 0.0);

        const double radius = scales.radius.value();
        const std::size_t density_row = layout.template offset<Kind::residual, Density>();
        for (const int row : boundary.reduced_dofs()) {
            CHECK(coordinates.residual_factors()[density_row + row] == Catch::Approx(1.0 / scales.density.value()).epsilon(2.0e-13));
        }
        const auto& density_metric = evaluator.cache().template get<normalization::ReferenceVolumeMetric<Density>>();
        const auto& selected_rows = boundary.reduced_dofs();
        int interior = 0;
        while (std::find(selected_rows.begin(), selected_rows.end(), interior) != selected_rows.end()) {
            ++interior;
        }
        REQUIRE(interior < static_cast<int>(state.template get<Density>().size()));
        // The density balance residual is a barotropic pressure potential, a specific energy.
        CHECK(coordinates.residual_factors()[density_row + interior] == Catch::Approx(1.0 /
            (scales.specific_energy.value() * std::sqrt(density_metric[interior] * radius * radius * radius))).epsilon(2.0e-13));
        CHECK(coordinates.state_factors()[layout.template offset<Kind::value, Constant>()] == Catch::Approx(1.0 / scales.specific_energy.value()).epsilon(2.0e-13));
        const auto infinity = serif::discretization::backend::mfem::select_boundary_dofs<
            serif::discretization::domain::InfinitySurfaceBoundary>(
                evaluator.binding().discretization->template binding<Potential>());
        for (const int row : infinity.reduced_dofs()) {
            CHECK(coordinates.residual_factors()[layout.template offset<Kind::residual, Potential>() + row] ==
                Catch::Approx(1.0 / scales.specific_energy.value()).epsilon(2.0e-13));
        }
        CHECK(coordinates.residual_factors()[layout.template offset<Kind::residual, Constant>()] == Catch::Approx(1.0 / scales.mass.value()).epsilon(2.0e-13));
        CHECK(coordinates.residual_factors()[layout.template offset<Kind::residual, Omega>()] == Catch::Approx(1.0 / scales.angular_momentum.value()).epsilon(2.0e-13));
        const auto& surface_metric = evaluator.cache().template get<deformation::keys::SurfaceMetric>();
        for (std::size_t mode = 0; mode < state.template get<Shape>().size(); ++mode) {
            CHECK(coordinates.state_factors()[layout.template offset<Kind::value, Shape>() + mode] == Catch::Approx(std::sqrt(surface_metric[mode]) / (radius * radius)).epsilon(2.0e-13));
            CHECK(coordinates.residual_factors()[layout.template offset<Kind::residual, ShapeRow>() + mode] == Catch::Approx(1.0 / scales.specific_energy.value()).epsilon(2.0e-13));
        }

        const Vector seed = tests::deterministic(static_cast<int>(state.values().size()), 0.43);
        const Vector dual = tests::deterministic(static_cast<int>(coordinates.residual_size()), 1.17);
        Vector transpose;
        frozen->jacobian().MultTranspose(dual, transpose);
        for (const int selected : {-1, 0, 1}) {
            CAPTURE(selected);
            Vector direction(seed);
            if (selected >= 0) {
                direction = 0.0;
                const std::size_t offset = selected == 0 ? layout.template offset<Kind::value, Shape>() : layout.template offset<Kind::value, Constant>();
                const std::size_t count = selected == 0 ? state.template get<Shape>().size() : 1;
                std::copy_n(seed.GetData() + offset, count, direction.GetData() + offset);
            }
            Vector physical_direction(static_cast<int>(state.values().size())), action, physical_action, expected_action(static_cast<int>(coordinates.residual_size()));
            coordinates.denormalize_state({direction.GetData(), coordinates.state_size()}, {physical_direction.GetData(), coordinates.state_size()});
            frozen->jacobian().Mult(direction, action);
            frozen->physical()->jacobian().action()->Mult(physical_direction, physical_action);
            coordinates.normalize_residual({physical_action.GetData(), coordinates.residual_size()}, {expected_action.GetData(), coordinates.residual_size()});
            check_rows(action, expected_action, layout, 2.0e-13, 2.0e-13);
            CHECK((action * dual) == Catch::Approx(direction * transpose).epsilon(3.0e-12).margin(3.0e-12));
            constexpr double step = 2.0e-6;
            const auto plus = evaluator.evaluate(perturbed_state(state, step, physical_direction));
            REQUIRE(plus);
            CHECK(evaluator.cache().report().template status<normalization::keys::Scales>() == dependencies::PreparationStatus::reused);
            CHECK(evaluator.cache().report().template status<normalization::keys::Normalization>() == dependencies::PreparationStatus::reused);
            CHECK(evaluator.cache().report().template status<deformation::keys::SurfaceMetric>() == dependencies::PreparationStatus::reused);
            Vector difference = residual(**plus);
            const auto minus = evaluator.evaluate(perturbed_state(state, -step, physical_direction));
            REQUIRE(minus);
            difference -= residual(**minus);
            difference /= 2.0 * step;
            check_rows(action, difference, layout, 3.0e-6, 3.0e-8);
        }
        const auto restored = evaluator.evaluate(state);
        REQUIRE(restored);
        check_rows(residual(**restored), residual(*frozen), layout, 0.0, 0.0);
        Vector repeated;
        frozen->jacobian().MultTranspose(dual, repeated);
        CHECK(tests::relative_error(repeated, transpose) == 0.0);
    }

    /** @brief Check the dimensional inverse, solver-coordinate adapter and independently measured linear residual. */
    template <typename Prepared, typename Layout>
    void check_stellar_linear_solve(const Prepared& prepared, const Layout& layout,
        const serif::discretization::backend::mfem::DofSelection& boundary) {
        const auto& system = *prepared.system();
        const auto& coordinates = system.coordinates();
        const int size = system.jacobian().Width();
        REQUIRE(size == system.jacobian().Height());
        REQUIRE(coordinates.state_factors().front() != 1.0);
        REQUIRE(coordinates.residual_factors().front() != 1.0);
        CHECK(prepared.contract() == serif::solver::preconditioner::ApplicationContract::fixed_linear);
        const Vector normalized_direction = tests::deterministic(size, 0.37);
        Vector physical_direction(size), physical_rhs, normalized_rhs(size);
        coordinates.denormalize_state({normalized_direction.GetData(), coordinates.state_size()},
            {physical_direction.GetData(), coordinates.state_size()});
        system.physical()->jacobian().action()->Mult(physical_direction, physical_rhs);
        coordinates.normalize_residual({physical_rhs.GetData(), coordinates.residual_size()},
            {normalized_rhs.GetData(), coordinates.residual_size()});

        Vector physical_correction(size), expected_correction(size), normalized_correction(size), preconditioned_action;
        prepared.physical_preconditioner().Mult(physical_rhs, physical_correction);
        coordinates.normalize_state({physical_correction.GetData(), coordinates.state_size()},
            {expected_correction.GetData(), coordinates.state_size()});
        prepared.preconditioner().Mult(normalized_rhs, normalized_correction);
        CHECK(tests::relative_error(normalized_correction, expected_correction) < 3.0e-13);
        system.physical()->jacobian().action()->Mult(physical_correction, preconditioned_action);
        for (const int row : boundary.reduced_dofs()) {
            const int index = static_cast<int>(layout.template offset<Kind::residual, Density>()) + row;
            CHECK(preconditioned_action[index] == Catch::Approx(physical_rhs[index]).epsilon(2.0e-11).margin(2.0e-12));
        }
        // The joint surface/scalar Schur solve must enforce every border equation and mutual coupling.
        for (const std::size_t index : {layout.template offset<Kind::residual, Constant>(), layout.template offset<Kind::residual, Omega>()}) {
            CHECK(preconditioned_action[static_cast<int>(index)] == Catch::Approx(physical_rhs[static_cast<int>(index)]).epsilon(2.0e-11).margin(2.0e-12));
        }

        for (std::size_t mode = 0; mode < layout.template size<Kind::residual, ShapeRow>(); ++mode) {
            const int index = static_cast<int>(layout.template offset<Kind::residual, ShapeRow>() + mode);
            CHECK(preconditioned_action[index] == Catch::Approx(physical_rhs[index]).epsilon(2.0e-11).margin(2.0e-12));
        }

        mfem::FGMRESSolver krylov;
        krylov.SetRelTol(1.0e-10);
        krylov.SetAbsTol(1.0e-13);
        krylov.SetMaxIter(400);
        krylov.SetKDim(100);
        krylov.SetPrintLevel(0);
        krylov.SetPreconditioner(prepared.preconditioner());
        krylov.SetOperator(system.jacobian());
        for (const bool newton_right_hand_side : {false, true}) {
            CAPTURE(newton_right_hand_side);
            Vector solve_rhs(physical_rhs);
            if (newton_right_hand_side) {
                solve_rhs = residual(*system.physical());
                solve_rhs *= -1.0;
            }
            coordinates.normalize_residual({solve_rhs.GetData(), coordinates.residual_size()},
                {normalized_rhs.GetData(), coordinates.residual_size()});
            Vector solution(size);
            solution = 0.0;
            krylov.Mult(normalized_rhs, solution);
            CAPTURE(krylov.GetNumIterations(), krylov.GetFinalNorm());
            CHECK(krylov.GetConverged());
            coordinates.denormalize_state({solution.GetData(), coordinates.state_size()},
                {physical_correction.GetData(), coordinates.state_size()});
            system.physical()->jacobian().action()->Mult(physical_correction, preconditioned_action);
            Vector measured_residual(preconditioned_action);
            measured_residual -= solve_rhs;
            Vector normalized_residual(size);
            coordinates.normalize_residual({measured_residual.GetData(), coordinates.residual_size()},
                {normalized_residual.GetData(), coordinates.residual_size()});
            CHECK(normalized_residual.Norml2() / normalized_rhs.Norml2() < 2.0e-10);
            CHECK(measured_residual.Norml2() / solve_rhs.Norml2() < 2.0e-8);
            check_rows(preconditioned_action, solve_rhs, layout, 2.0e-7, 2.0e-9);
        }
    }

    /** @brief Prepare the production stellar recipe selectively and retain its complete frozen numerical state. */
    void check_moving_preconditioner() {
        namespace normalization = serif::solver::normalization;
        namespace preconditioner = serif::solver::preconditioner;
        using InverseKey = preconditioner::keys::Preconditioner;
        using PreparedKey = model::preconditioning::PreparedSystem;
        using Status = dependencies::PreparationStatus;
#ifdef MFEM_USE_MPI
        // The AMG backend borrows the application's MPI session; the test is its application owner.
        if (!mfem::Mpi::IsInitialized()) {
            mfem::Mpi::Init();
        }
#endif
        auto fixture = std::make_shared<Fixture>(mfem::Ordering::byNODES);
        std::weak_ptr<Fixture> lifetime = fixture;
        Vector right_hand_side, frozen_correction, frozen_action;
        auto retained = [&]() {
            const auto physics = physical_model();
            const auto discretization = discretize(fixture);
            using Description = decltype(model::resolve_equilibrium_equations(physics));
            using Plan = model::EquilibriumPlan<Polytrope, Description>;
            const std::array shape{0.004, -0.002, 0.001, 0.003, -0.001, 0.002, -0.003, 0.001, 0.002};
            const auto layout = backend::equilibrium_layout<Plan>(discretization, shape);
            const auto state = initial_state(layout, fixture->fields, shape);
            const auto recipe = preconditioner::make_primal_preconditioner(layout);
            static_assert(dependencies::ContainsKey<Density, typename decltype(recipe)::Values>);
            static_assert(dependencies::ContainsKey<Potential, typename decltype(recipe)::Values>);
            static_assert(dependencies::ContainsKey<Shape, typename decltype(recipe)::Values>);
            static_assert(!dependencies::ContainsKey<::serif::tests::equations::Enthalpy, typename Plan::Values>);
            const auto configured = ::serif::model::make_discretization(::serif::model::DiscretizationOptions{
                .model = physics, .finite_elements = discretization,
                .normalization = normalization::make_normalization(normalization::NormalizationOptions{
                    .scheme = normalization::LaneEmdenApproximation{.n = 3.0, .central_density = serif::dimensions::DensityValue{0.9}}, .coordinates = normalization::Dimensionless{},
                    .gravitational_constant = serif::dimensions::QuantityValue<serif::dimensions::GravitationalConstant>{2.3 / (4.0 * std::numbers::pi)}
                }),
                .preconditioner = recipe
            });
            serif::discretization::quadrature::RuleSet rules;
            rules.set_fallback_control({.fixed_order = 8});
            const execution::VolumeOptions volume_options{.quadrature_policy = serif::discretization::quadrature::Policy{rules}};
            auto evaluator = backend::bind_stellar_model(configured, state, volume_options);
            const auto physical = evaluator.evaluate(state);
            REQUIRE(physical);
            CHECK(evaluator.cache().report().template status<InverseKey>() == Status::not_required);
            CHECK(evaluator.cache().report().template status<PreparedKey>() == Status::not_required);
            const auto prepared = evaluator.prepare_linear_system(state);
            if (!prepared) {
                UNSCOPED_INFO("Failed computation index: " << prepared.error().computation_index);
                if (std::holds_alternative<preconditioner::Failure>(prepared.error().error)) {
                    UNSCOPED_INFO(std::get<preconditioner::Failure>(prepared.error().error).message);
                }
            }
            REQUIRE(prepared);
            const auto first = *prepared;
            CHECK(first->system() == *physical);
            CHECK(evaluator.cache().report().template status<model::keys::PhysicalSystem>() == Status::reused);
            CHECK(evaluator.cache().report().template status<InverseKey>() == Status::computed);
            const auto repeated = evaluator.prepare_linear_system(state);
            REQUIRE(repeated);
            CHECK(*repeated == first);
            CHECK(evaluator.cache().report().computed_count() == 0);
            CHECK(evaluator.cache().report().template status<InverseKey>() == Status::reused);

            const auto boundary = serif::discretization::backend::mfem::select_boundary_dofs<
                serif::discretization::domain::StellarSurfaceBoundary>(fixture->fields.density_binding);
            check_stellar_linear_solve(*first, layout, boundary);
            right_hand_side = tests::deterministic(static_cast<int>(state.values().size()), 0.82);
            frozen_correction.SetSize(right_hand_side.Size());
            first->preconditioner().Mult(right_hand_side, frozen_correction);
            first->system()->jacobian().Mult(right_hand_side, frozen_action);
            const Vector first_residual = residual(*first->system());

            auto candidate = state;
            candidate.template get<Density>()[0] += 0.01;
            candidate.template get<Shape>()[0] += 0.0002;
            const auto trial = evaluator.evaluate(candidate);
            REQUIRE(trial);
            CHECK(*trial != first->system());
            CHECK(evaluator.cache().report().template status<InverseKey>() == Status::not_required);
            CHECK(evaluator.cache().report().template status<PreparedKey>() == Status::not_required);
            CHECK(tests::relative_error(residual(*first->system()), first_residual) == 0.0);
            const auto updated = evaluator.prepare_linear_system(candidate);
            REQUIRE(updated);
            CHECK(*updated != first);
            CHECK((*updated)->system() == *trial);
            CHECK(evaluator.cache().report().template status<InverseKey>() == Status::computed);
            const auto rebuilt = evaluator.prepare_linear_system(candidate, dependencies::make_dependency_context(), dependencies::PreparationMode::rebuild);
            REQUIRE(rebuilt);
            CHECK(*rebuilt != *updated);
            CHECK(evaluator.cache().report().template status<model::keys::PhysicalSystem>() == Status::computed);
            CHECK(evaluator.cache().report().template status<InverseKey>() == Status::computed);
            Vector repeated_correction(right_hand_side.Size()), repeated_action;
            first->preconditioner().Mult(right_hand_side, repeated_correction);
            first->system()->jacobian().Mult(right_hand_side, repeated_action);
            CHECK(tests::relative_error(repeated_correction, frozen_correction) == 0.0);
            CHECK(tests::relative_error(repeated_action, frozen_action) == 0.0);

            // At n=1.5 the exact density Jacobian is undefined at zero density; reject the EOS trial.
            auto singular = state;
            std::ranges::fill(singular.template get<Density>(), 0.0);
            singular.template get<Omega>()[0] = 0.0;
            const auto preparations = evaluator.cache().preparation_count();
            const auto rejected = evaluator.prepare_linear_system(singular);
            REQUIRE_FALSE(rejected);
            REQUIRE(std::holds_alternative<execution::PreparationRejection>(rejected.error().error));
            CHECK(std::get<execution::PreparationRejection>(rejected.error().error).reason == execution::PreparationRejectionReason::EQUATION);
            CHECK_FALSE(evaluator.cache().ready());
            CHECK_FALSE(evaluator.cache().template contains<PreparedKey>());
            CHECK(evaluator.cache().report().computed_count() == 0);
            CHECK(evaluator.cache().report().reused_count() == 0);
            CHECK(evaluator.cache().preparation_count() == preparations);
            first->preconditioner().Mult(right_hand_side, repeated_correction);
            first->system()->jacobian().Mult(right_hand_side, repeated_action);
            CHECK(tests::relative_error(repeated_correction, frozen_correction) == 0.0);
            CHECK(tests::relative_error(repeated_action, frozen_action) == 0.0);
            CHECK(tests::relative_error(residual(*first->system()), first_residual) == 0.0);
            const auto recovered = evaluator.prepare_linear_system(candidate);
            REQUIRE(recovered);
            CHECK(evaluator.cache().ready());
            return first;
        }();
        fixture.reset();
        CHECK_FALSE(lifetime.expired());
        Vector repeated_correction(right_hand_side.Size()), repeated_action;
        retained->preconditioner().Mult(right_hand_side, repeated_correction);
        retained->system()->jacobian().Mult(right_hand_side, repeated_action);
        CHECK(tests::relative_error(repeated_correction, frozen_correction) == 0.0);
        CHECK(tests::relative_error(repeated_action, frozen_action) == 0.0);
        retained.reset();
        CHECK(lifetime.expired());
    }

    /** @brief Evaluate the full moving star in both native displacement orderings and retain its frozen actions. */
    template <bool Normalize = false>
    void check_moving_equilibrium(const mfem::Ordering::Type ordering) {
        auto fixture = std::make_shared<Fixture>(ordering);
        std::weak_ptr<Fixture> lifetime = fixture;
        Vector direction, dual, action, transpose, normalized_action, normalized_transpose;
        std::shared_ptr<const backend::NormalizedSystem> normalized_retained;
        const auto exercise = [&]() {
            auto& fields = fixture->fields;
            const auto physics = physical_model();
            const auto discretization = discretize(fixture);
            using Description = decltype(model::resolve_equilibrium_equations(physics));
            using Plan = model::EquilibriumPlan<Polytrope, Description>;
            const std::array shape{0.004, -0.002, 0.001, 0.003, -0.001, 0.002, -0.003, 0.001, 0.002};
            const auto layout = backend::equilibrium_layout<Plan>(discretization, shape);
            REQUIRE(layout.state_sizes.size() == 5);
            REQUIRE(layout.residual_sizes.size() == 5);
            REQUIRE(layout.template total_size<Kind::value>() == layout.template total_size<Kind::residual>());
            REQUIRE(shape.size() < static_cast<std::size_t>(fields.displacement_binding.size()));
            const auto state = initial_state(layout, fields, shape);
            const auto boundary = serif::discretization::backend::mfem::select_boundary_dofs<
                serif::discretization::domain::StellarSurfaceBoundary>(fields.density_binding);
            serif::discretization::quadrature::RuleSet rules;
            rules.set_fallback_control({.fixed_order = 8});
            const execution::VolumeOptions volume_options{.quadrature_policy = serif::discretization::quadrature::Policy{rules}};
            auto evaluator = [&]() {
                if constexpr (Normalize) {
                    namespace normalization = serif::solver::normalization;
                    using Reference = normalization::LaneEmdenApproximation<normalization::FromModel, serif::dimensions::DensityValue>;
                    const auto configured = ::serif::model::make_discretization(::serif::model::DiscretizationOptions{
                        .model = physics, .finite_elements = discretization,
                        .normalization = normalization::make_normalization(normalization::NormalizationOptions<Reference>{
                            .scheme = Reference{.n = 3.0, .central_density = serif::dimensions::DensityValue{0.9}},
                            .coordinates = normalization::Dimensionless{},
                            .gravitational_constant = serif::dimensions::QuantityValue<serif::dimensions::GravitationalConstant>{2.3 / (4.0 * std::numbers::pi)}
                        })
                    });
                    return backend::bind_stellar_model(configured, state, volume_options);
                } else {
                    return backend::bind_stellar_model(physics, discretization, state, volume_options);
                }
            }();
            const auto initial = evaluator.evaluate(state);
            REQUIRE(initial);
            const auto first = *initial;
            using MassRow = model::computations::ScalarResidualRow<Constant>;
            using AngularMomentumRow = model::computations::ScalarResidualRow<Omega>;
            const auto mass_row = evaluator.cache().template share<MassRow>();
            const auto angular_momentum_row = evaluator.cache().template share<AngularMomentumRow>();
            CHECK(evaluator.cache().report().template status<MassRow>() == dependencies::PreparationStatus::computed);
            CHECK(evaluator.cache().report().template status<AngularMomentumRow>() == dependencies::PreparationStatus::computed);
            const auto repeated = evaluator.evaluate(state);
            REQUIRE(repeated);
            CHECK(*repeated == first);
            CHECK(evaluator.cache().report().computed_count() == 0);
            CHECK(evaluator.cache().template share<MassRow>() == mass_row);
            CHECK(evaluator.cache().template share<AngularMomentumRow>() == angular_momentum_row);

            // A successful trial is merely another result; the caller retains the previous candidate.
            const auto original = state;
            auto copied = state;
            CHECK(copied.values().data() != state.values().data());
            copied.template get<Density>()[0] += 0.01;
            copied.template get<Shape>()[0] += 0.0002;
            auto moved = std::move(copied);
            CHECK(std::ranges::equal(state.values(), original.values()));
            const Vector first_residual = residual(*first);
            const auto discarded = evaluator.evaluate(moved);
            REQUIRE(discarded);
            CHECK(evaluator.cache().template share<MassRow>() != mass_row);
            CHECK(evaluator.cache().template share<AngularMomentumRow>() != angular_momentum_row);
            CHECK(evaluator.cache().report().template status<MassRow>() == dependencies::PreparationStatus::computed);
            CHECK(evaluator.cache().report().template status<AngularMomentumRow>() == dependencies::PreparationStatus::computed);
            CHECK(*discarded != first);
            CHECK(tests::relative_error(residual(**discarded), first_residual) > 1.0e-10);
            CHECK(tests::relative_error(residual(*first), first_residual) == 0.0);
            CHECK(std::ranges::equal(state.values(), original.values()));
            const auto restored = evaluator.evaluate(state);
            REQUIRE(restored);
            if constexpr (Normalize) {
                check_normalized_system(evaluator, state, boundary);
            } else {
                check_composition(evaluator, state, **restored, boundary);
                check_scalar_row<Constant>(evaluator, **restored);
                check_scalar_row<Omega>(evaluator, **restored);
                check_derivatives(evaluator, state);
            }

            const auto evaluated = evaluator.evaluate(state);
            REQUIRE(evaluated);
            const auto retained = *evaluated;
            const auto physical = [&]() {
                if constexpr (Normalize) {
                    return retained->physical();
                } else {
                    return retained;
                }
            }();
            const auto preparations = evaluator.cache().preparation_count();
            const Vector retained_residual = residual(*retained);
            direction = tests::deterministic(static_cast<int>(state.values().size()), 0.61);
            dual = tests::deterministic(static_cast<int>(state.values().size()), 1.31);
            physical->jacobian().action()->Mult(direction, action);
            physical->jacobian().action()->MultTranspose(dual, transpose);
            if constexpr (Normalize) {
                normalized_retained = retained;
                normalized_retained->jacobian().Mult(direction, normalized_action);
                normalized_retained->jacobian().MultTranspose(dual, normalized_transpose);
            }
            auto invalid = state;
            invalid.template get<Shape>()[0] = std::numeric_limits<double>::quiet_NaN();
            const auto rejected = evaluator.evaluate(invalid);
            REQUIRE_FALSE(rejected);
            CHECK(std::holds_alternative<deformation::DeformationFailure>(rejected.error().error));
            CHECK(tests::relative_error(residual(*retained), retained_residual) == 0.0);
            if constexpr (Normalize) {
                CHECK_FALSE(evaluator.cache().template contains<serif::model::normalization::PreparedSystem>());
            }
            CHECK_FALSE(evaluator.cache().ready());
            CHECK_FALSE(evaluator.cache().template contains<model::keys::PhysicalSystem>());
            CHECK(evaluator.cache().report().computed_count() == 0);
            CHECK(evaluator.cache().report().reused_count() == 0);
            CHECK(evaluator.cache().preparation_count() == preparations);
            const auto recovered = evaluator.evaluate(state);
            REQUIRE(recovered);
            CHECK(*recovered == retained);
            CHECK(evaluator.cache().report().computed_count() == 0);
            return physical;
        };
        auto retained = exercise();
        fixture.reset();
        CHECK_FALSE(lifetime.expired());
        Vector repeated_action, repeated_transpose;
        retained->jacobian().action()->Mult(direction, repeated_action);
        retained->jacobian().action()->MultTranspose(dual, repeated_transpose);
        CHECK(tests::relative_error(repeated_action, action) == 0.0);
        CHECK(tests::relative_error(repeated_transpose, transpose) == 0.0);
        retained.reset();
        if constexpr (Normalize) {
            CHECK_FALSE(lifetime.expired());
            normalized_retained->jacobian().Mult(direction, repeated_action);
            normalized_retained->jacobian().MultTranspose(dual, repeated_transpose);
            CHECK(tests::relative_error(repeated_action, normalized_action) == 0.0);
            CHECK(tests::relative_error(repeated_transpose, normalized_transpose) == 0.0);
            normalized_retained.reset();
        }
        CHECK(lifetime.expired());
    }
}

TEST_CASE("Complete moving stellar assembly preserves every coupled equation and boundary row", "[model][moving_equilibrium][integration]") {
    SECTION("Component-block displacement ordering") {
        check_moving_equilibrium(mfem::Ordering::byNODES);
    }
    SECTION("Interleaved displacement ordering") {
        check_moving_equilibrium(mfem::Ordering::byVDIM);
    }
}

TEST_CASE("Complete moving stellar normalization preserves solver coordinates and frozen ownership", "[model][moving_normalization][integration]") {
    check_moving_equilibrium<true>(mfem::Ordering::byNODES);
}

TEST_CASE("Complete moving stellar preconditioning prepares a frozen solver system and solves its physical Jacobian", "[model][moving_preconditioner][integration]") {
#ifdef MFEM_USE_MPI
    check_moving_preconditioner();
#else
    SKIP("The default stellar recipe uses the one-process BoomerAMG backend, which requires an MPI-enabled MFEM build.");
#endif
}