#include <algorithm>
#include <cmath>
#include <concepts>
#include <limits>
#include <numbers>
#include <type_traits>
#include <utility>

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

import mean_field;
import test_helpers;

namespace outer_manifest_report_test {
    template <mean_field::model::StellarModelType Model> class PreparedEarlierMultiplier;

    class EarlierMultiplier final {
    public:
        struct Parameters final {
            mean_field::dimensions::SpecificEnergyValue target;
        };

        using ModelDefinition = mean_field::integral::FixedWithMultiplier<
            EarlierMultiplier,
            "AardvarkOuterManifestMultiplier",
            mean_field::models::DependsOn<mean_field::models::stellar::state::Density>,
            mean_field::models::Affects<mean_field::models::stellar::equation::HydrostaticBalance>,
            mean_field::models::GlobalScalarNormalization<
                mean_field::models::PhysicalScaleLaw::specific_energy,
                mean_field::models::PhysicalScaleLaw::specific_energy>,
            mean_field::models::GeneratedManifest<
                "aardvark_outer_manifest.value",
                "a",
                "aardvark_outer_manifest.residual",
                "R_a",
                "specific_energy",
                "specific_energy">>;
        using EquilibriumPhysics = mean_field::operators::SpecificationEquilibriumPhysics<PreparedEarlierMultiplier>;

        explicit EarlierMultiplier(const Parameters parameters) noexcept : m_target(parameters.target) {
        }

        [[nodiscard]] mean_field::dimensions::SpecificEnergyValue target() const noexcept {
            return m_target;
        }

    private:
        mean_field::dimensions::SpecificEnergyValue m_target;
    };

    template <mean_field::model::StellarModelType Model> class PreparedEarlierMultiplier final {
    public:
        using Report = mean_field::operators::EmptySpecificationPreparationReport;

        explicit PreparedEarlierMultiplier(const EarlierMultiplier &) noexcept {
        }

        template <typename StateView> [[nodiscard]] Report PrepareAfterPhysical(const StateView &) noexcept {
            return {};
        }

        template <
            typename Equation,
            typename Row>
        [[nodiscard]] mean_field::stellar::StructuralZero AddResidual(
            Equation,
            Row &
        ) const noexcept {
            return mean_field::stellar::structuralZero;
        }

        template <
            typename Equation,
            typename State,
            typename Direction,
            typename Row>
        [[nodiscard]] mean_field::stellar::StructuralZero AddJacobianAction(
            mean_field::stellar::Derivative<
                Equation,
                State>,
            const Direction &,
            Row &
        ) const noexcept {
            return mean_field::stellar::zeroDerivative;
        }

        [[nodiscard]] bool IsPrepared() const noexcept {
            return true;
        }
    };
} // namespace outer_manifest_report_test

namespace {
    using BaseModel                          = mean_field::model::StellarModel<mean_field::models::SpecificationSet<
        mean_field::eos::Polytrope,
        mean_field::surface::Isobaric,
        mean_field::integral::FixedTotalMass>>;

    using CentralDensityModel                = mean_field::model::StellarModel<mean_field::models::SpecificationSet<
        mean_field::eos::Polytrope,
        mean_field::surface::Isobaric,
        mean_field::integral::FixedTotalMass,
        mean_field::constraint::FixedCentralDensity>>;

    using AngularMomentumModel               = mean_field::model::StellarModel<mean_field::models::SpecificationSet<
        mean_field::eos::Polytrope,
        mean_field::surface::Isobaric,
        mean_field::integral::FixedTotalMass,
        mean_field::integral::FixedAngularMomentum>>;

    using AngularMomentumCentralDensityModel = mean_field::model::StellarModel<mean_field::models::SpecificationSet<
        mean_field::eos::Polytrope,
        mean_field::surface::Isobaric,
        mean_field::integral::FixedTotalMass,
        mean_field::integral::FixedAngularMomentum,
        mean_field::constraint::FixedCentralDensity>>;

    using IncompleteModel =
        mean_field::model::StellarModel<mean_field::models::SpecificationSet<mean_field::eos::Polytrope>>;

    using EarlierMultiplierModel = mean_field::model::StellarModel<mean_field::models::SpecificationSet<
        mean_field::eos::Polytrope,
        mean_field::surface::Isobaric,
        outer_manifest_report_test::EarlierMultiplier,
        mean_field::integral::FixedTotalMass>>;

    template <typename Candidate>
    concept HasLegacyNumericalModelAdapter = requires { typename Candidate::NumericalModelAdapter; };

    [[nodiscard]] mean_field::operators::StellarEquilibriumDependencies make_dependencies() {
        return {
            .discretization     = {.identity = 4001, .revision = 1},
            .density            = {.identity = 4003, .revision = 1},
            .surfaceDeformation = {.identity = 4007, .revision = 1},
            .gravityGradient    = {.identity = 4013, .revision = 1},
            .gravityPotential   = {.identity = 4019, .revision = 1},
            .enthalpy           = {.identity = 4021, .revision = 1},
            .bernoulliConstant  = {.identity = 4027, .revision = 1},
            .rotation           = {.identity = 4049, .revision = 1},
            .targetMass         = {.identity = 4051, .revision = 1}
        };
    }

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

    [[nodiscard]] double relative_difference(
        const mfem::Vector &left,
        const mfem::Vector &right
    ) {
        mfem::Vector difference(left);
        difference -= right;
        return difference.Norml2() / std::max({1.0, left.Norml2(), right.Norml2()});
    }

} // namespace

TEST_CASE(
    "Stellar Model Selects A Compile-Time Equilibrium Problem Type",
    tags::stellar_equilibrium_problem_type_contract
) {
    using namespace mean_field;

    using BaseProblem            = equilibrium::StellarEquilibriumProblem<BaseModel>;
    using CentralDensityProblem  = equilibrium::StellarEquilibriumProblem<CentralDensityModel>;
    using AngularMomentumProblem = equilibrium::StellarEquilibriumProblem<AngularMomentumModel>;
    using AngularMomentumCentralDensityProblem =
        equilibrium::StellarEquilibriumProblem<AngularMomentumCentralDensityModel>;

    STATIC_CHECK(equilibrium::StellarEquilibriumModel<BaseModel>);
    STATIC_CHECK(equilibrium::StellarEquilibriumModel<CentralDensityModel>);
    STATIC_CHECK(equilibrium::StellarEquilibriumModel<AngularMomentumModel>);
    STATIC_CHECK(equilibrium::StellarEquilibriumModel<AngularMomentumCentralDensityModel>);
    STATIC_CHECK(equilibrium::StellarEquilibriumModel<EarlierMultiplierModel>);
    STATIC_CHECK_FALSE(equilibrium::StellarEquilibriumModel<IncompleteModel>);
    STATIC_CHECK_FALSE(
        operators::StellarEquilibriumRuntimeContribution<outer_manifest_report_test::EarlierMultiplier>::registered
    );
    STATIC_CHECK_FALSE(
        operators::stellarEquilibriumBackendRuntimeAuthorized<
            outer_manifest_report_test::EarlierMultiplier, EarlierMultiplierModel>
    );
    STATIC_CHECK(
        operators::StellarEquilibriumPhysicsAvailableFor<
            outer_manifest_report_test::EarlierMultiplier, EarlierMultiplierModel>
    );
    STATIC_CHECK_FALSE(std::same_as<BaseProblem, CentralDensityProblem>);
    STATIC_CHECK(BaseProblem::symbolicallySquare);
    STATIC_CHECK(CentralDensityProblem::symbolicallySquare);
    STATIC_CHECK_FALSE(BaseProblem::hasFixedCentralDensity);
    STATIC_CHECK(CentralDensityProblem::hasFixedCentralDensity);
    STATIC_CHECK(AngularMomentumProblem::hasFixedAngularMomentum);
    STATIC_CHECK_FALSE(AngularMomentumProblem::hasFixedCentralDensity);
    STATIC_CHECK(AngularMomentumCentralDensityProblem::hasFixedAngularMomentum);
    STATIC_CHECK(AngularMomentumCentralDensityProblem::hasFixedCentralDensity);
    STATIC_CHECK_FALSE(HasLegacyNumericalModelAdapter<BaseProblem>);
    STATIC_CHECK_FALSE(HasLegacyNumericalModelAdapter<CentralDensityProblem>);
    STATIC_CHECK(std::same_as<BaseProblem, equilibrium::StellarEquilibriumSystem<BaseModel>>);
    STATIC_CHECK(
        std::same_as<
            typename BaseProblem::PreparedOperatorType,
            operators::PreparedVariadicStellarEquilibriumOperator<BaseModel>>
    );
    STATIC_CHECK(
        std::same_as<
            typename CentralDensityProblem::PreparedOperatorType,
            operators::PreparedVariadicStellarEquilibriumOperator<CentralDensityModel>>
    );
    STATIC_CHECK(
        std::same_as<
            typename AngularMomentumProblem::PreparedOperatorType,
            operators::PreparedVariadicStellarEquilibriumOperator<AngularMomentumModel>>
    );
    STATIC_CHECK_FALSE(
        std::same_as<typename BaseProblem::PreparedOperatorType, typename CentralDensityProblem::PreparedOperatorType>
    );
    STATIC_CHECK_FALSE(
        std::same_as<
            typename AngularMomentumProblem::PreparedOperatorType,
            typename AngularMomentumCentralDensityProblem::PreparedOperatorType>
    );
    STATIC_CHECK(AngularMomentumProblem::FormType::value_block_count == 7);
    STATIC_CHECK(AngularMomentumCentralDensityProblem::FormType::value_block_count == 8);
    STATIC_CHECK(
        std::same_as<typename BaseProblem::FormType, utils::blocks::surface_deformed_stellar_equilibrium_form>
    );
    STATIC_CHECK(
        std::same_as<
            typename CentralDensityProblem::FormType, utils::blocks::central_density_bordered_stellar_equilibrium_form>
    );
    STATIC_CHECK(
        std::same_as<
            typename BaseProblem::CompiledSurfaceConstraintType,
            surface::CompiledPressureSurfaceConstraintT<
                typename BaseProblem::ThermodynamicEquationsType::PressureSurfaceFormulation, eos::Polytrope>>
    );
    STATIC_CHECK(material::CompiledThermodynamicEquations<typename BaseProblem::ThermodynamicEquationsType>);
}

TEST_CASE(
    "Fixed Angular Momentum Root Uses Its Generated Angular Velocity In Every Physical Row",
    "[fixed-angular-momentum][stellar-equilibrium][jacobian][integration]"
) {
    using namespace mean_field;

    utils::Args arguments   = test_utils::setup_args();
    fem::FEM finiteElements = fem::setup_fem(arguments.mesh_file, arguments, 0);
    REQUIRE(finiteElements.okay());

    constexpr double radius                = utils::RADIUS;
    constexpr double mass                  = utils::MASS;
    constexpr double targetAngularMomentum = 0.1;
    const double polytropicConstant        = 2.0 * utils::G * radius * radius / std::numbers::pi_v<double>;
    const double seedCentralDensity        = 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{targetAngularMomentum}, .axis = {0.0, 0.0, 3.0}}
        )
    );
    auto problem   = equilibrium::discretize(model, std::move(finiteElements));
    auto projected = seed::makeProjectedEquilibriumState(
        problem,
        seed::LaneEmden({.centralDensity = dimensions::DensityValue{seedCentralDensity}, .radialSampleCount = 1024})
    );
    auto dependencies      = make_dependencies();

    const auto preparation = problem.Prepare(projected.values, dependencies);
    CHECK(preparation.generatedPhysicalControl);
    CHECK(preparation.physical.DidAnyWork());
    CHECK(preparation.template specification<models::FixedAngularMomentum>().constraint.DidAnyWork());
    CHECK(preparation.template specification<models::FixedAngularMomentum>().generatedRotation);
    CHECK(problem.IsPrepared());
    const auto angularReport = problem.GetPreparedOperator().GetAngularMomentumReport();
    CHECK(angularReport.targetAngularMomentum == targetAngularMomentum);
    CHECK(angularReport.angularVelocity > 0.0);
    CHECK(angularReport.momentOfInertia > 0.0);
    CHECK(std::abs(angularReport.scaledResidual) < 7.0e-4);

    mfem::Vector direction(problem.StateSize());
    direction                             = 0.0;
    mfem::Vector angularVelocityDirection = problem.GetManifest().stateView(direction).block(
        utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term
    );
    REQUIRE(angularVelocityDirection.Size() == 1);
    angularVelocityDirection(0) = -0.37;
    angularVelocityDirection.SyncAliasMemory(direction);

    mfem::Vector analyticAction;
    problem.ApplyLinearization(direction, analyticAction);

    constexpr double step = 1.0e-5;
    mfem::Vector plusState(projected.values);
    plusState.Add(step, direction);
    problem.Prepare(plusState, dependencies);
    mfem::Vector plusResidual;
    problem.BuildResidual(plusResidual);
    mfem::Vector minusState(projected.values);
    minusState.Add(-step, direction);
    problem.Prepare(minusState, dependencies);
    mfem::Vector minusResidual;
    problem.BuildResidual(minusResidual);
    plusResidual -= minusResidual;
    plusResidual /= 2.0 * step;

    auto analyticView     = problem.GetManifest().residualView(analyticAction);
    auto differenceView   = problem.GetManifest().residualView(plusResidual);
    const auto blockError = [&](const auto &term) {
        const mfem::Vector analytic   = analyticView.block(term);
        const mfem::Vector difference = differenceView.block(term);
        return relative_difference(analytic, difference);
    };

    const double surfaceError  = blockError(utils::blocks::surface_deformation_field.shape_equilibrium_term);
    const double enthalpyError = blockError(utils::blocks::enthalpy_field.specific_term);
    const double angularMomentumError =
        blockError(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term);
    INFO("Generated-Omega surface-row centered-difference error = " << surfaceError);
    INFO("Generated-Omega hydrostatic-row centered-difference error = " << enthalpyError);
    INFO("Generated-Omega invariant-row centered-difference error = " << angularMomentumError);
    CHECK(surfaceError < 3.0e-7);
    CHECK(enthalpyError < 3.0e-7);
    CHECK(angularMomentumError < 3.0e-10);
    CHECK(analyticView.block(utils::blocks::surface_deformation_field.shape_equilibrium_term).Norml2() > 0.0);
    CHECK(analyticView.block(utils::blocks::enthalpy_field.specific_term).Norml2() > 0.0);
    CHECK(analyticView.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term).Norml2() > 0.0);
    CHECK(analyticView.block(utils::blocks::gravity_field.gradient_term).Norml2() == 0.0);
    CHECK(analyticView.block(utils::blocks::gravity_field.poisson_term).Norml2() == 0.0);
    CHECK(analyticView.block(utils::blocks::density_field.mass_term).Norml2() == 0.0);
    CHECK(analyticView.block(utils::blocks::fixed_total_mass_constraint.mass_normalization_term).Norml2() == 0.0);

    mfem::Vector nonFiniteAngularVelocityState(projected.values);
    auto nonFiniteAngularVelocity = problem.GetManifest()
                                        .stateView(nonFiniteAngularVelocityState)
                                        .block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term);
    REQUIRE(nonFiniteAngularVelocity.Size() == 1);
    nonFiniteAngularVelocity(0) = std::numeric_limits<double>::quiet_NaN();
    nonFiniteAngularVelocity.SyncAliasMemory(nonFiniteAngularVelocityState);
    const auto nonFiniteControl = problem.TryPrepare(nonFiniteAngularVelocityState, dependencies);
    REQUIRE_FALSE(nonFiniteControl.has_value());
    CHECK(
        nonFiniteControl.error().reason == operators::StellarEquilibriumPreparationRejectionReason::non_finite_physics
    );
    CHECK_FALSE(problem.IsPrepared());
    REQUIRE(problem.TryPrepare(projected.values, dependencies).has_value());

    mfem::Vector negativeDensityState(projected.values);
    auto negativeDensity =
        problem.GetManifest().stateView(negativeDensityState).block(utils::blocks::density_field.mass_term);
    negativeDensity *= -1.0;
    negativeDensity.SyncAliasMemory(negativeDensityState);
    auto negativeDensityDependencies = dependencies;
    ++negativeDensityDependencies.density.revision;
    const auto inadmissibleMoment = problem.TryPrepare(negativeDensityState, negativeDensityDependencies);
    REQUIRE_FALSE(inadmissibleMoment.has_value());
    CHECK(
        inadmissibleMoment.error().reason ==
        operators::StellarEquilibriumPreparationRejectionReason::inadmissible_physics
    );
    CHECK_FALSE(problem.IsPrepared());
    ++negativeDensityDependencies.density.revision;
    REQUIRE(problem.TryPrepare(projected.values, negativeDensityDependencies).has_value());
    CHECK(problem.IsPrepared());

    auto zeroModel = 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.0}})
    );
    fem::FEM zeroFiniteElements = fem::setup_fem(arguments.mesh_file, arguments, 0);
    REQUIRE(zeroFiniteElements.okay());
    auto zeroProblem = equilibrium::discretize(zeroModel, std::move(zeroFiniteElements));
    mfem::Vector zeroState(projected.values);
    zeroProblem.GetManifest().stateView(zeroState).block(
        utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term
    ) = 0.0;
    zeroProblem.Prepare(zeroState, dependencies);
    mfem::Vector zeroAction;
    zeroProblem.ApplyLinearization(direction, zeroAction);
    auto zeroView = zeroProblem.GetManifest().residualView(zeroAction);
    CHECK(zeroView.block(utils::blocks::surface_deformation_field.shape_equilibrium_term).Norml2() == 0.0);
    CHECK(zeroView.block(utils::blocks::enthalpy_field.specific_term).Norml2() == 0.0);
    CHECK(zeroView.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term).Norml2() > 0.0);
}

TEST_CASE(
    "Fixed Mass Reports Use The Inferred Outer Manifest Indices",
    "[stellar-equilibrium][manifest][runtime][ordering]"
) {
    using namespace mean_field;
    using Form = operators::CompiledStellarEquilibriumForm<EarlierMultiplierModel>;
    using EarlierValue =
        utils::blocks::generated_value_block<models::MultiplierFor<outer_manifest_report_test::EarlierMultiplier>>;
    using MassValue = utils::blocks::fixed_total_mass::mass_normalization::value;

    STATIC_CHECK(utils::blocks::type_index_v<EarlierValue, typename Form::value_blocks> == 5);
    STATIC_CHECK(utils::blocks::type_index_v<MassValue, typename Form::value_blocks> == 6);

    utils::Args arguments   = test_utils::setup_args();
    fem::FEM finiteElements = fem::setup_fem(arguments.mesh_file, arguments, 0);
    REQUIRE(finiteElements.okay());

    auto model = model::StellarModel(
        eos::Polytrope({.n = 1.0, .K = 0.25}), surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}),
        outer_manifest_report_test::EarlierMultiplier({.target = dimensions::SpecificEnergyValue{0.75}}),
        integral::FixedTotalMass({.Mtotal = dimensions::MassValue{1.25}})
    );
    auto problem = equilibrium::discretize(model, std::move(finiteElements));

    mfem::Vector state(problem.StateSize());
    state                                                        = 0.0;
    const auto stateView                                         = problem.GetManifest().stateView(state);
    stateView.block(utils::blocks::density_field.mass_term)      = 1.0;
    stateView.block(utils::blocks::enthalpy_field.specific_term) = 1.0;
    stateView.block(utils::blocks::fixed_total_mass_constraint.mass_normalization_term) = 0.25;

    const auto preparation = problem.TryPrepare(state, make_dependencies(), make_zero_rotation());
    REQUIRE(preparation.has_value());
    REQUIRE(preparation->physical.DidAnyWork());

    const auto report           = problem.GetPreparedOperator().GetFixedMassReport();
    const auto &outerDescriptor = problem.GetManifest().template specification<models::FixedTotalMass>();
    CHECK(report.descriptor.stableId == outerDescriptor.stableId);
    CHECK(report.descriptor.valueBlock == outerDescriptor.valueBlock);
    CHECK(report.descriptor.residualBlock == outerDescriptor.residualBlock);
    CHECK(report.descriptor.valueBlock == 6);
    CHECK(report.descriptor.residualBlock == 6);
    CHECK(report.descriptor.target == 1.25);
    CHECK(report.dimensionalResidual == report.achieved - report.descriptor.target);
    CHECK(report.scaledResidual == report.dimensionalResidual / report.descriptor.residualScale);
}

TEST_CASE(
    "Discretized Stellar Equilibrium Problem Is Exactly Equivalent To The Legacy Construction Path",
    tags::stellar_equilibrium_problem_integration
) {
    using namespace mean_field;

    utils::Args args              = test_utils::setup_args();
    fem::FEM legacyFiniteElements = fem::setup_fem(args.mesh_file, args, 0);
    REQUIRE(legacyFiniteElements.okay());

    fem::FEM modelDrivenFiniteElements = fem::setup_fem(args.mesh_file, args, 0);
    REQUIRE(modelDrivenFiniteElements.okay());
    const MPI_Comm modelDrivenCommunicator         = modelDrivenFiniteElements.mesh->GetComm();
    const mapping::DomainMapper *modelDrivenMapper = modelDrivenFiniteElements.domainMapperStateless.get();

    models::StellarModel legacyModel{
        models::structure::PolytropicStructure{eos::Polytrope{3.0, 0.25}, 1.25},
        surface::ConstantPressureSurface{eos::PressureValue{0.0}}
    };
    operators::PreparedStellarEquilibriumOperator legacyOperator(
        legacyFiniteElements, *legacyFiniteElements.domainMapperStateless, legacyModel
    );

    auto equilibriumProblem = equilibrium::discretize(
        model::StellarModel(
            integral::FixedTotalMass({.Mtotal = dimensions::MassValue{1.25}}),
            surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}), eos::Polytrope({.n = 3.0, .K = 0.25})
        ),
        std::move(modelDrivenFiniteElements)
    );
    auto &modelDrivenOperator    = equilibriumProblem.GetPreparedOperator();
    const auto &physicalOperator = equilibriumProblem.GetPhysicalOperator();

    CHECK(equilibriumProblem.StateSize() == legacyOperator.Width());
    CHECK(equilibriumProblem.EquationSize() == legacyOperator.Height());
    CHECK(equilibriumProblem.StateSize() == equilibriumProblem.EquationSize());
    CHECK(equilibriumProblem.GetCommunicator() == modelDrivenCommunicator);
    CHECK(&equilibriumProblem.GetDiscretization().domainMapper() == modelDrivenMapper);
    CHECK(equilibriumProblem.GetDiscretization().isCurrent());
    CHECK(physicalOperator.GetTargetMass() == 1.25);
    CHECK(physicalOperator.GetSurfaceConstraintOperator().GetPhysicalCondition().targetPressure == 0.0);
    CHECK(equilibriumProblem.GetCompiledSurfaceConstraint().targetPressure() == dimensions::PressureValue{0.0});
    CHECK(physicalOperator.GetDomainDeformation().matchesCurrentDiscretization());
    CHECK(&equilibriumProblem.GetLinearizationOperator() == &modelDrivenOperator);
    CHECK(equilibriumProblem.GetManifest().template specification<models::FixedTotalMass>().target == 1.25);

    mfem::Vector state(legacyOperator.Width());
    state                                                        = 0.0;
    const auto stateView                                         = legacyOperator.GetRootManifest().stateView(state);
    stateView.block(utils::blocks::density_field.mass_term)      = 1.0;
    stateView.block(utils::blocks::enthalpy_field.specific_term) = 1.0;

    const operators::StellarEquilibriumDependencies dependencies = make_dependencies();
    const physics::RigidRotation rotation                        = make_zero_rotation();
    legacyOperator.Prepare(state, dependencies, rotation);
    equilibriumProblem.Prepare(state, dependencies, rotation);

    mfem::Vector legacyResidual;
    mfem::Vector modelDrivenResidual;
    legacyOperator.BuildResidual(legacyResidual);
    equilibriumProblem.BuildResidual(modelDrivenResidual);
    CHECK(relative_difference(modelDrivenResidual, legacyResidual) < 2.0e-15);

    mfem::Vector direction(state.Size());
    for (int index = 0; index < direction.Size(); ++index) {
        direction(index) = 0.01 * std::sin(0.31 * static_cast<double>(index + 1));
    }
    mfem::Vector legacyAction;
    mfem::Vector modelDrivenAction;
    legacyOperator.Mult(direction, legacyAction);
    equilibriumProblem.ApplyLinearization(direction, modelDrivenAction);
    CHECK(relative_difference(modelDrivenAction, legacyAction) < 2.0e-15);
}

TEST_CASE(
    "Fixed Angular Momentum Composes With The Optional Central Density Phase At Runtime",
    "[fixed-angular-momentum][central-density][stellar-equilibrium][integration]"
) {
    using namespace mean_field;

    utils::Args arguments   = test_utils::setup_args();
    fem::FEM finiteElements = fem::setup_fem(arguments.mesh_file, arguments, 0);
    REQUIRE(finiteElements.okay());
    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}}),
        integral::FixedAngularMomentum({.Jtotal = dimensions::AngularMomentumValue{0.2}}),
        constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{1.0}})
    );
    auto problem  = equilibrium::discretize(model, std::move(finiteElements));
    using Problem = std::remove_cvref_t<decltype(problem)>;
    STATIC_CHECK(Problem::FormType::value_block_count == 8);
    STATIC_CHECK(Problem::FormType::residual_block_count == 8);

    mfem::Vector state(problem.StateSize());
    state                                                        = 0.0;
    const auto stateView                                         = problem.GetManifest().stateView(state);
    stateView.block(utils::blocks::density_field.mass_term)      = 1.0;
    stateView.block(utils::blocks::enthalpy_field.specific_term) = 1.0;
    stateView.block(utils::blocks::fixed_total_mass_constraint.mass_normalization_term)     = 0.25;
    stateView.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term) = 0.4;
    stateView.block(utils::blocks::fixed_central_density_phase.central_value_term)          = 0.03;

    const auto report = problem.Prepare(state, make_dependencies());
    CHECK(report.template specification<models::FixedCentralDensity>().constraint.DidAnyWork());
    CHECK(report.template specification<models::FixedAngularMomentum>().constraint.DidAnyWork());
    CHECK(problem.IsPrepared());
    CHECK(problem.StateSize() == problem.GetPhysicalOperator().Width() + 2);
    REQUIRE(problem.GetManifest().constraints().size() == 4);
    CHECK(
        problem.GetManifest().template specification<models::FixedAngularMomentum>().stableId == "FixedAngularMomentum"
    );
    CHECK(
        problem.GetManifest().template specification<models::FixedCentralDensity>().stableId == "FixedCentralDensity"
    );

    mfem::Vector residual;
    problem.BuildResidual(residual);
    REQUIRE(residual.Size() == problem.EquationSize());
    const auto residualView = problem.GetManifest().residualView(residual);
    CHECK(std::isfinite(residualView.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term)(0)));
    CHECK(std::isfinite(residualView.block(utils::blocks::fixed_central_density_phase.central_value_term)(0)));

    mfem::Vector direction(problem.StateSize());
    for (int index = 0; index < direction.Size(); ++index) {
        direction(index) = 0.01 * std::sin(0.17 * static_cast<double>(index + 1));
    }
    mfem::Vector action;
    problem.ApplyLinearization(direction, action);
    REQUIRE(action.Size() == problem.EquationSize());
    for (int index = 0; index < action.Size(); ++index) {
        CHECK(std::isfinite(action(index)));
    }
}