module;

#include <array>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <expected>
#include <memory>
#include <optional>
#include <stdexcept>
#include <tuple>
#include <type_traits>
#include <utility>

#include <mfem.hpp>
#include <mpi.h>

export module mean_field:operators.prepared_variadic_stellar_equilibrium;

export import :operators.prepared_angular_momentum;
export import :operators.prepared_central_density;
export import :operators.prepared_stellar_equilibrium;
export import :operators.stellar_equilibrium_compiler;

/*
 * A physics-authored residual or derivative provider must make one of two
 * explicit statements for every row/edge inferred from Reads/Changes:
 *
 *   - return the token produced by row.add(...); or
 *   - return structuralZero when the declared edge is identically zero.
 *
 * The outer runtime, rather than the extension, enumerates the compiler's
 * complete incidence set.  These tiny result types let that enumeration
 * distinguish an intentional mathematical zero from an accidentally empty
 * hook without exposing any backend block machinery to a physics author.
 */
export namespace mean_field::stellar {
    struct ContributionAdded final { };
    struct StructuralZero final { };

    inline constexpr StructuralZero structuralZero{};
    inline constexpr StructuralZero zeroDerivative{};

    template <typename Candidate>
    concept ContributionResult = std::same_as<std::remove_cvref_t<Candidate>, ContributionAdded> ||
                                 std::same_as<std::remove_cvref_t<Candidate>, StructuralZero>;
} // namespace mean_field::stellar

export namespace mean_field::operators {
    /**
     * Capability boundary for the coupled finite-element physics core.
     * Describing an EOS is intentionally easier than implementing its
     * finite-element runtime core. Surface equations are compiled and
     * prepared by their own specification contribution, so this backend is
     * selected solely by the constitutive law.
     */
    template <typename EquationOfState> struct StellarEquilibriumCoreRuntime {
        static constexpr bool registered = false;
    };

    template <> struct StellarEquilibriumCoreRuntime<eos::Polytrope> {
        static constexpr bool registered = true;
        using CoreType                   = PreparedStellarEquilibriumOperator;

        [[nodiscard]] static std::unique_ptr<CoreType> Make(
            fem::FEM &finiteElements,
            const mapping::DomainMapper &domainMapper,
            const eos::Polytrope &equationOfState,
            const models::CompiledFixedMass &fixedMass,
            PressureSurfaceConstraintView surfaceConstraint,
            deformation::PreparedDomainDeformationRuntime domainDeformation
        ) {
            return std::make_unique<PreparedStellarEquilibriumOperator>(
                finiteElements, domainMapper, equationOfState, fixedMass, surfaceConstraint,
                std::move(domainDeformation)
            );
        }

        [[nodiscard]] static int SurfaceEquationCount(const CoreType &core) noexcept {
            return static_cast<int>(core.GetSurfaceConstraintOperator().GetSurfaceRows().size());
        }
    };

    /**
     * Minimal common protocol consumed by the variadic outer root.
     *
     * EOS backends may use different concrete core types.  They only need to
     * implement this numerical protocol and expose that type as
     * StellarEquilibriumCoreRuntime<EOS>::CoreType.  Specification runtimes
     * are audited separately against the selected concrete core, so a
     * constraint that needs additional physical facilities is rejected at its
     * own compile-time boundary.
     */
    template <typename Candidate>
    concept PreparedStellarEquilibriumPhysicalCore =
        std::derived_from<std::remove_cvref_t<Candidate>, mfem::Operator> &&
        requires(
            std::remove_cvref_t<Candidate> &core,
            const std::remove_cvref_t<Candidate> &constantCore,
            const mfem::Vector &state,
            mfem::Vector &residual,
            const StellarEquilibriumDependencies &dependencies,
            const physics::RigidRotation &rotation
        ) {
            { constantCore.GetLayout() } -> std::same_as<const StellarEquilibriumLayout &>;
            { core.Prepare(state, dependencies, rotation) } -> std::same_as<PreparedStellarEquilibriumReport>;
            { constantCore.BuildResidual(residual) } -> std::same_as<void>;
            { constantCore.IsPrepared() } -> std::convertible_to<bool>;
            { constantCore.GetFixedMassReport() } -> std::same_as<RootConstraintReport>;
            { constantCore.GetDependencies() } -> std::same_as<const StellarEquilibriumDependencies &>;
            {
                constantCore.GetGeneratedDisplacementDependency()
            } -> std::same_as<const StellarEquilibriumDependencyStamp &>;
            { constantCore.GetSurfaceConstraintOperator() } -> std::same_as<const PreparedPressureSurfaceConstraint &>;
        };

    template <typename Candidate>
    concept FalliblePreparedStellarEquilibriumPhysicalCore =
        PreparedStellarEquilibriumPhysicalCore<Candidate> && requires(
                                                                 std::remove_cvref_t<Candidate> &core,
                                                                 const mfem::Vector &state,
                                                                 const StellarEquilibriumDependencies &dependencies,
                                                                 const physics::RigidRotation &rotation
                                                             ) {
            {
                core.TryPrepare(state, dependencies, rotation)
            } -> std::same_as<StellarEquilibriumPreparationResult<PreparedStellarEquilibriumReport>>;
        };

    namespace detail {
        template <typename Candidate> struct BackendSpecificationListTraits final {
            static constexpr bool valid                        = false;

            template <typename> static constexpr bool contains = false;
        };

        template <typename... Specifications>
        struct BackendSpecificationListTraits<models::ModelTypeList<Specifications...>> final {
            static constexpr bool valid =
                (models::ModelSpecification<Specifications> && ...) &&
                utils::blocks::types_are_unique_v<utils::blocks::type_list<Specifications...>>;

            template <typename Query>
            static constexpr bool contains = (std::same_as<std::remove_cvref_t<Query>, Specifications> || ...);
        };

        template <model::StellarModelType Model, typename = void> struct CoreRuntimeInterfaceAudit {
            using CoreType                 = void;
            static constexpr bool complete = false;
        };

        template <model::StellarModelType Model>
        struct CoreRuntimeInterfaceAudit<
            Model,
            std::void_t<
                typename StellarEquilibriumCoreRuntime<
                    typename std::remove_cvref_t<Model>::EquationOfStateType>::CoreType,
                typename StellarEquilibriumCoreRuntime<
                    typename std::remove_cvref_t<Model>::EquationOfStateType>::CoreType::BackendSpecifications,
                std::bool_constant<static_cast<bool>(StellarEquilibriumCoreRuntime<typename std::remove_cvref_t<
                                                         Model>::EquationOfStateType>::registered)>>> {
        private:
            using ModelType       = std::remove_cvref_t<Model>;
            using EquationOfState = typename ModelType::EquationOfStateType;
            using Runtime         = StellarEquilibriumCoreRuntime<EquationOfState>;

        public:
            using CoreType = typename Runtime::CoreType;

            static constexpr bool complete =
                BackendSpecificationListTraits<typename CoreType::BackendSpecifications>::valid &&
                PreparedStellarEquilibriumPhysicalCore<CoreType> &&
                requires(
                    fem::FEM &finiteElements,
                    const mapping::DomainMapper &domainMapper,
                    const EquationOfState &equationOfState,
                    const models::CompiledFixedMass &fixedMass,
                    PressureSurfaceConstraintView surfaceConstraint,
                    deformation::PreparedDomainDeformationRuntime domainDeformation,
                    const CoreType &core
                ) {
                    requires Runtime::registered;
                    {
                        Runtime::Make(
                            finiteElements, domainMapper, equationOfState, fixedMass, surfaceConstraint,
                            std::move(domainDeformation)
                        )
                    } -> std::same_as<std::unique_ptr<CoreType>>;
                    { Runtime::SurfaceEquationCount(core) } -> std::convertible_to<int>;
                };
        };
    } // namespace detail

    /**
     * Structural contract for the privileged specification list owned by an
     * EOS/core backend.  Ordinary EOS, surface, and constraint authors do not
     * use this facility; their nested EquilibriumPhysics package remains the
     * restricted, astronomy-facing extension path.
     */
    template <typename Candidate>
    concept StellarEquilibriumBackendSpecificationList =
        detail::BackendSpecificationListTraits<std::remove_cvref_t<Candidate>>::valid;

    /**
     * Backend runtime extension point for one physical model specification.
     *
     * The prepared stellar root is assembled by folding this trait over the
     * model's canonical specification list.  A new specification therefore
     * contributes one prepared slot; no specialization for a *combination* of
     * specifications is ever required. New physics-facing specifications
     * should prefer their nested EquilibriumPhysics package below; explicit
     * specializations remain the library/backend registry mechanism, but are
     * selected only when the concrete core owner lists that exact
     * specification in CoreType::BackendSpecifications.
     */
    template <models::ModelSpecification Specification> struct StellarEquilibriumRuntimeContribution {
        static constexpr bool registered               = false;
        static constexpr std::size_t rotationProviders = 0;
    };

    template <template <typename> typename PreparedImplementation, std::size_t RotationProviderCount = 0>
    struct PreparedStellarEquilibriumContribution {
        static constexpr bool registered                        = true;
        static constexpr std::size_t rotationProviders          = RotationProviderCount;

        template <model::StellarModelType Model> using Prepared = PreparedImplementation<Model>;
    };

    /* Physics-facing declaration for a specification's residual/Jacobian
     * runtime. A constraint may expose
     *
     *   using EquilibriumPhysics =
     *       operators::SpecificationEquilibriumPhysics<MyPreparedPhysics>;
     *
     * inside its own class. Unlike an explicit
     * StellarEquilibriumRuntimeContribution specialization, this declaration
     * is adapted through restricted physics-facing views. The implementation
     * receives its exact specification at construction and never receives the
     * full model, FEM backend, domain mapper, dependency set, or physical core.
     * Explicit registry specializations remain a candidate backend extension
     * point for built-in physics which must coordinate core internals. The
     * selected concrete core's non-extendable BackendSpecifications member is what
     * grants that candidate privileged access. */
    template <template <typename> typename PreparedImplementation, std::size_t RotationProviderCount = 0>
    struct SpecificationEquilibriumPhysics final {
        static_assert(
            RotationProviderCount <= 1,
            "One specification runtime can provide at most one rigid-rotation control."
        );

        static constexpr bool registered                       = true;
        static constexpr std::size_t rotationProviders         = RotationProviderCount;

        template <model::StellarModelType Model> using Physics = PreparedImplementation<Model>;
    };

    namespace detail {
        template <typename LocalPhysics> struct BindLocalSpecificationEquilibriumPhysics final {
            template <typename> using Physics = LocalPhysics;
        };
    } // namespace detail

    /**
     * Convenience spelling for ordinary specification-local physics.
     *
     * Most constraint implementations depend only on their exact
     * specification and the restricted block views supplied by the adapter;
     * they do not need the complete Model type. This alias binds such a
     * concrete class into SpecificationEquilibriumPhysics without introducing
     * a second runtime or duplicating any adapter logic.
     */
    template <typename LocalPhysics>
    using LocalSpecificationEquilibriumPhysics = SpecificationEquilibriumPhysics<
        detail::BindLocalSpecificationEquilibriumPhysics<LocalPhysics>::template Physics>;

    namespace detail {
        template <typename Specification>
        concept HasNestedStellarEquilibriumPhysics =
            requires { typename std::remove_cvref_t<Specification>::EquilibriumPhysics; };

        template <typename Specification, typename Model, typename = void>
        struct BackendRuntimeContributionAuthorization : std::false_type { };

        template <models::ModelSpecification Specification, model::StellarModelType Model>
        struct BackendRuntimeContributionAuthorization<
            Specification,
            Model,
            std::void_t<typename StellarEquilibriumCoreRuntime<
                typename std::remove_cvref_t<Model>::EquationOfStateType>::CoreType::BackendSpecifications>>
            final
            : std::bool_constant<
                  CoreRuntimeInterfaceAudit<std::remove_cvref_t<Model>>::complete &&
                  std::remove_cvref_t<Model>::template containsSpecification<std::remove_cvref_t<Specification>> &&
                  BackendSpecificationListTraits<typename StellarEquilibriumCoreRuntime<
                      typename std::remove_cvref_t<Model>::EquationOfStateType>::CoreType::BackendSpecifications>::
                      template contains<std::remove_cvref_t<Specification>>> { };

        template <typename Specification, typename Model, typename = void>
        struct BackendRuntimeContributionCandidate final {
            static constexpr bool available                = false;
            static constexpr bool registered               = false;
            static constexpr std::size_t rotationProviders = 0;
        };

        template <models::ModelSpecification Specification, model::StellarModelType Model>
        struct BackendRuntimeContributionCandidate<
            Specification,
            Model,
            std::void_t<
                std::enable_if_t<BackendRuntimeContributionAuthorization<Specification, Model>::value>,
                typename StellarEquilibriumRuntimeContribution<Specification>::template Prepared<Model>,
                std::bool_constant<static_cast<bool>(StellarEquilibriumRuntimeContribution<Specification>::registered)>,
                std::integral_constant<
                    std::size_t,
                    static_cast<std::size_t>(StellarEquilibriumRuntimeContribution<Specification>::rotationProviders)>>>
            final {
            using Contribution                             = StellarEquilibriumRuntimeContribution<Specification>;
            using Prepared                                 = typename Contribution::template Prepared<Model>;

            static constexpr bool available                = true;
            static constexpr bool registered               = Contribution::registered;
            static constexpr std::size_t rotationProviders = Contribution::rotationProviders;
        };

        template <models::ModelSpecification Specification, model::StellarModelType Model, typename Physics>
        class PhysicsFacingSpecificationRuntime;

        template <
            models::ModelSpecification Specification,
            model::StellarModelType Model,
            bool HasNestedPhysics = HasNestedStellarEquilibriumPhysics<Specification>,
            typename              = void>
        struct RuntimeContributionSelection {
            static constexpr bool available                = false;
            static constexpr bool ambiguous                = false;
            static constexpr std::size_t rotationProviders = 0;
        };

        /* A core-authorized explicit registry specialization is trusted
         * backend code. Its established protocol deliberately retains direct
         * access to the FEM, mapper, physical core, model, and dependency
         * stamps. A specialization alone is ignored, so it cannot confer that
         * privilege on an external specification paired with an existing
         * core. */
        template <models::ModelSpecification Specification, model::StellarModelType Model>
        struct RuntimeContributionSelection<
            Specification,
            Model,
            false,
            std::void_t<typename BackendRuntimeContributionCandidate<Specification, Model>::Prepared>> {
            using Candidate                                = BackendRuntimeContributionCandidate<Specification, Model>;
            using Contribution                             = typename Candidate::Contribution;
            using Prepared                                 = typename Candidate::Prepared;

            static constexpr bool available                = Candidate::available;
            static constexpr bool registered               = Candidate::registered;
            static constexpr bool ambiguous                = false;
            static constexpr std::size_t rotationProviders = Candidate::rotationProviders;
        };

        /* A nested package is the safe physics-author path. The adapter owns
         * all interaction with backend objects and forwards only restricted
         * views to the implementation. Malformed packages remain detection
         * safe through this partial specialization. */
        template <models::ModelSpecification Specification, model::StellarModelType Model>
        struct RuntimeContributionSelection<
            Specification,
            Model,
            true,
            std::void_t<
                typename Specification::EquilibriumPhysics,
                std::bool_constant<static_cast<bool>(Specification::EquilibriumPhysics::registered)>,
                std::integral_constant<
                    std::size_t,
                    static_cast<std::size_t>(Specification::EquilibriumPhysics::rotationProviders)>,
                typename Specification::EquilibriumPhysics::template Physics<Model>>> {
            using Contribution               = typename Specification::EquilibriumPhysics;
            using Physics                    = typename Contribution::template Physics<Model>;
            using Prepared                   = PhysicsFacingSpecificationRuntime<Specification, Model, Physics>;

            static constexpr bool available  = true;
            static constexpr bool registered = Contribution::registered;
            static constexpr bool ambiguous  = BackendRuntimeContributionCandidate<Specification, Model>::registered;
            static constexpr std::size_t rotationProviders = Contribution::rotationProviders;
        };

        template <models::ModelSpecification Specification, model::StellarModelType Model>
        using PreparedRuntimeContribution = typename RuntimeContributionSelection<Specification, Model>::Prepared;
    } // namespace detail

    /**
     * Whether the selected concrete core explicitly permits one specification
     * to use the privileged aggregate runtime registry. This is intentionally
     * false for an otherwise valid external specification paired with the
     * built-in core; such a specification must use its restricted nested
     * EquilibriumPhysics package.
     */
    template <typename Specification, typename Model>
    inline constexpr bool stellarEquilibriumBackendRuntimeAuthorized = detail::
        BackendRuntimeContributionAuthorization<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>::value;

    struct EmptySpecificationPreparationReport final {
        [[nodiscard]] constexpr bool DidAnyWork() const noexcept {
            return false;
        }
    };

    namespace detail {
        struct StellarEquilibriumControlContext final {
            StellarEquilibriumDependencies dependencies;
            std::optional<physics::RigidRotation> rotation;
            std::size_t rotationProviderCount{0};
            bool generatedPhysicalControl{false};
        };

        /*
         * Additional facilities used specifically by FixedAngularMomentum.
         * This is deliberately a constraint-local protocol: an EOS core can
         * support the base root without implementing these operations, and is
         * rejected only when this physical constraint is selected.
         */
        template <typename Candidate>
        concept FixedAngularMomentumPhysicalCore =
            PreparedStellarEquilibriumPhysicalCore<Candidate> && requires(const std::remove_cvref_t<Candidate> &core) {
                {
                    core.GetGravityContext()
                } -> std::same_as<const context::gravity_field::GravityFieldLinearizationContext &>;
                { core.GetDomainDeformation() } -> std::same_as<const deformation::PreparedDomainDeformationRuntime &>;
                { core.GetBarotropicClosureOperator() } -> std::same_as<const PreparedBarotropicClosureOperator &>;
                { core.GetHydrostaticOperator() } -> std::same_as<const PreparedHydrostaticEquilibriumOperator &>;
                { core.GetSurfaceConstraintOperator() } -> std::same_as<const PreparedPressureSurfaceConstraint &>;
                { core.GetDisplacementOperator() } -> std::same_as<const PreparedDisplacementResidualOperator &>;
                { core.GetSurfaceDeformationParameters() } -> std::same_as<const mfem::Vector &>;
                {
                    core.GetGeneratedDisplacementDependency()
                } -> std::same_as<const StellarEquilibriumDependencyStamp &>;
            };

        [[nodiscard]] inline StellarEquilibriumPreparationRejection
        makeAngularMomentumPreparationRejection(const AngularMomentumPreparationRejection &rejection) {
            using ChildReason = AngularMomentumPreparationRejectionReason;
            using RootReason  = StellarEquilibriumPreparationRejectionReason;
            using RootStage   = StellarEquilibriumPreparationStage;

            switch (rejection.reason) {
            case ChildReason::inverted_geometry:
                return {.reason = RootReason::inverted_geometry, .stage = RootStage::model_specification};
            case ChildReason::non_finite_geometry:
                return {.reason = RootReason::non_finite_geometry, .stage = RootStage::model_specification};
            case ChildReason::negative_moment_of_inertia:
                return {.reason = RootReason::inadmissible_physics, .stage = RootStage::model_specification};
            case ChildReason::non_finite_angular_velocity:
            case ChildReason::non_finite_density:
            case ChildReason::non_finite_moment_of_inertia:
            case ChildReason::non_finite_residual:
                return {.reason = RootReason::non_finite_physics, .stage = RootStage::model_specification};
            }
            throw std::logic_error("An unknown angular-momentum trial rejection reached the stellar root.");
        }

        [[nodiscard]] inline StellarEquilibriumPreparationRejection
        markSpecificationRejection(StellarEquilibriumPreparationRejection rejection) noexcept {
            if (rejection.stage == StellarEquilibriumPreparationStage::unspecified) {
                rejection.stage = StellarEquilibriumPreparationStage::model_specification;
            }
            return rejection;
        }

        [[nodiscard]] inline StellarEquilibriumPreparationResult<double> synchronizeReplicatedControl(
            const double localValue,
            const MPI_Comm communicator,
            const char *failureMessage
        ) {
            const int localFinite = std::isfinite(localValue) ? 1 : 0;
            int globallyFinite    = 0;
            if (MPI_Allreduce(&localFinite, &globallyFinite, 1, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS) {
                throw std::runtime_error(failureMessage);
            }
            if (globallyFinite == 0) {
                return std::unexpected(
                    StellarEquilibriumPreparationRejection{
                        .reason = StellarEquilibriumPreparationRejectionReason::non_finite_physics,
                        .stage  = StellarEquilibriumPreparationStage::model_specification
                    }
                );
            }

            double minimumValue     = 0.0;
            double maximumValue     = 0.0;
            const int minimumStatus = MPI_Allreduce(&localValue, &minimumValue, 1, MPI_DOUBLE, MPI_MIN, communicator);
            const int maximumStatus = MPI_Allreduce(&localValue, &maximumValue, 1, MPI_DOUBLE, MPI_MAX, communicator);
            if (minimumStatus != MPI_SUCCESS || maximumStatus != MPI_SUCCESS) {
                throw std::runtime_error(failureMessage);
            }
            if (minimumValue != maximumValue) {
                return std::unexpected(
                    StellarEquilibriumPreparationRejection{
                        .reason = StellarEquilibriumPreparationRejectionReason::inadmissible_physics,
                        .stage  = StellarEquilibriumPreparationStage::model_specification
                    }
                );
            }
            return minimumValue;
        }

        template <typename Specification, model::StellarModelType Model> class EmbeddedSpecificationRuntime final {
        public:
            using Report = EmptySpecificationPreparationReport;

            template <PreparedStellarEquilibriumPhysicalCore PhysicalCore>
            EmbeddedSpecificationRuntime(
                fem::FEM &,
                const mapping::DomainMapper &,
                PhysicalCore &,
                const Model &
            ) noexcept {
            }

            template <typename StateView>
            void ReadPhysicalControls(
                const StateView &,
                StellarEquilibriumControlContext &
            ) noexcept {
            }

            template <typename StateView>
            [[nodiscard]] StellarEquilibriumPreparationResult<void> TryReadPhysicalControls(
                const StateView &,
                StellarEquilibriumControlContext &
            ) noexcept {
                return {};
            }

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

            template <typename ResidualView> void AddResidual(const ResidualView &) const noexcept {
            }

            template <
                typename DirectionView,
                typename ActionView,
                PreparedStellarEquilibriumPhysicalCore PhysicalCore>
            void AddJacobianAction(
                const DirectionView &,
                const ActionView &,
                const PhysicalCore &
            ) const noexcept {
            }

            [[nodiscard]] constexpr bool IsPrepared() const noexcept {
                return true;
            }
        };

        struct FixedAngularMomentumPreparationReport final {
            PreparedAngularMomentumReport constraint;
            bool generatedRotation{false};

            [[nodiscard]] bool DidAnyWork() const noexcept {
                return constraint.DidAnyWork() || generatedRotation;
            }
        };

        template <model::StellarModelType Model> class FixedAngularMomentumRuntime final {
        public:
            using Report            = FixedAngularMomentumPreparationReport;
            using PreparationResult = StellarEquilibriumPreparationResult<Report>;
            using ControlResult     = StellarEquilibriumPreparationResult<void>;

            template <FixedAngularMomentumPhysicalCore PhysicalCore>
            FixedAngularMomentumRuntime(
                fem::FEM &finiteElements,
                const mapping::DomainMapper &domainMapper,
                PhysicalCore &physical,
                const Model &model
            )
                : m_constraint(
                      finiteElements,
                      domainMapper,
                      physical.GetGravityContext(),
                      models::compileConstraint(model.template specification<models::FixedAngularMomentum>())
                  ),
                  m_volumeDisplacementDirection(physical.GetDomainDeformation().volumeDisplacementSize()),
                  m_rotationalAngularVelocityAction(physical.GetDomainDeformation().volumeDisplacementSize()),
                  m_surfaceAngularVelocityAction(physical.GetDomainDeformation().parameterCount()),
                  m_hydrostaticAngularVelocityAction(physical.GetBarotropicClosureOperator().GetEnthalpySize()),
                  m_zeroEnthalpy(physical.GetBarotropicClosureOperator().GetEnthalpySize()),
                  m_communicator(finiteElements.mesh->GetComm()) {
                m_generatedRotationDependency.identity =
                    static_cast<std::uint64_t>(reinterpret_cast<std::uintptr_t>(this));
                m_zeroEnthalpy = 0.0;
            }

            template <typename StateView>
            void ReadPhysicalControls(
                const StateView &state,
                StellarEquilibriumControlContext &context
            ) {
                const auto angularVelocity =
                    state.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term);
                MFEM_VERIFY(
                    angularVelocity.Size() == 1 && std::isfinite(angularVelocity(0)),
                    "FixedAngularMomentum must generate one finite angular-velocity coordinate."
                );

                UpdatePhysicalControls(angularVelocity(0), context);
            }

            template <typename StateView>
            [[nodiscard]] ControlResult TryReadPhysicalControls(
                const StateView &state,
                StellarEquilibriumControlContext &context
            ) {
                const auto angularVelocity =
                    state.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term);
                MFEM_VERIFY(
                    angularVelocity.Size() == 1, "FixedAngularMomentum must generate one angular-velocity coordinate."
                );

                auto synchronizedAngularVelocity = synchronizeReplicatedControl(
                    angularVelocity(0), m_communicator,
                    "FixedAngularMomentum could not synchronize its angular-velocity coordinate."
                );
                if (!synchronizedAngularVelocity.has_value()) {
                    return std::unexpected(synchronizedAngularVelocity.error());
                }

                UpdatePhysicalControls(*synchronizedAngularVelocity, context);
                return {};
            }

            template <
                typename StateView,
                FixedAngularMomentumPhysicalCore PhysicalCore>
            [[nodiscard]] Report PrepareAfterPhysical(
                const StateView &state,
                const StellarEquilibriumDependencies &dependencies,
                const PhysicalCore &physical
            ) {
                auto result = TryPrepareAfterPhysical(state, dependencies, physical);
                if (!result.has_value()) {
                    throwStellarEquilibriumPreparationRejection(result.error());
                }
                return std::move(result).value();
            }

            template <
                typename StateView,
                FixedAngularMomentumPhysicalCore PhysicalCore>
            [[nodiscard]] PreparationResult TryPrepareAfterPhysical(
                const StateView &,
                const StellarEquilibriumDependencies &dependencies,
                const PhysicalCore &physical
            ) {
                m_isPrepared                                          = false;
                const StellarEquilibriumDependencyStamp &displacement = physical.GetGeneratedDisplacementDependency();
                auto constraintResult                                 = m_constraint.TryPrepare(
                    m_angularVelocity,
                    {.discretization =
                         {.identity = dependencies.discretization.identity,
                          .revision = dependencies.discretization.revision},
                     .density = {.identity = dependencies.density.identity, .revision = dependencies.density.revision},
                     .displacement = {.identity = displacement.identity, .revision = displacement.revision},
                     .rotation     = {
                         .identity = dependencies.rotation.identity, .revision = dependencies.rotation.revision
                     }}
                );
                if (!constraintResult.has_value()) {
                    return std::unexpected(makeAngularMomentumPreparationRejection(constraintResult.error()));
                }
                m_isPrepared = true;
                return Report{
                    .constraint = std::move(constraintResult).value(), .generatedRotation = m_generatedRotationChanged
                };
            }

            template <typename ResidualView> void AddResidual(const ResidualView &residual) const {
                mfem::Vector constraintResidual;
                m_constraint.BuildResidual(constraintResidual);
                residual.add(
                    utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term, constraintResidual
                );
            }

            template <
                typename DirectionView,
                typename ActionView,
                FixedAngularMomentumPhysicalCore PhysicalCore>
            void AddJacobianAction(
                const DirectionView &direction,
                const ActionView &action,
                const PhysicalCore &physical
            ) const {
                const auto densityDirection = direction.block(utils::blocks::density_field.mass_term);
                const auto surfaceDirection = direction.block(utils::blocks::surface_deformation_field.parameters_term);
                const auto angularVelocityDirection =
                    direction.block(utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term);
                MFEM_VERIFY(
                    angularVelocityDirection.Size() == 1 && std::isfinite(angularVelocityDirection(0)),
                    "The angular-velocity direction must be finite."
                );

                m_angularMomentumAction.SetSize(1);
                m_angularMomentumAction = 0.0;
                physical.GetDomainDeformation().applyJacobian(
                    physical.GetSurfaceDeformationParameters(), surfaceDirection, m_volumeDisplacementDirection
                );
                m_constraint.ApplyCompleteJacobianAction(
                    densityDirection, m_volumeDisplacementDirection, angularVelocityDirection(0),
                    m_angularMomentumAction
                );
                action.add(
                    utils::blocks::fixed_angular_momentum_constraint.angular_velocity_term, m_angularMomentumAction
                );

                if (m_angularVelocity == 0.0 || angularVelocityDirection(0) == 0.0) {
                    return;
                }

                const double fractionalVariation = angularVelocityDirection(0) / m_angularVelocity;
                physical.GetHydrostaticOperator().ApplyRotationAmplitudeJacobianAction(
                    fractionalVariation, m_hydrostaticAngularVelocityAction
                );
                m_zeroEnthalpy = 0.0;
                physical.GetSurfaceConstraintOperator().ApplyJacobianRows(
                    m_zeroEnthalpy, m_hydrostaticAngularVelocityAction
                );
                action.add(utils::blocks::enthalpy_field.specific_term, m_hydrostaticAngularVelocityAction);

                physical.GetDisplacementOperator().GetRotationalOperator().BuildResidual(
                    m_rotationalAngularVelocityAction
                );
                m_rotationalAngularVelocityAction *= 2.0 * fractionalVariation;
                physical.GetDomainDeformation().applyJacobianTranspose(
                    physical.GetSurfaceDeformationParameters(), m_rotationalAngularVelocityAction,
                    m_surfaceAngularVelocityAction
                );
                action.add(
                    utils::blocks::surface_deformation_field.shape_equilibrium_term, m_surfaceAngularVelocityAction
                );
            }

            [[nodiscard]] bool IsPrepared() const noexcept {
                return m_isPrepared && m_constraint.IsPrepared();
            }

            [[nodiscard]] const PreparedAngularMomentumOperator &constraint() const noexcept {
                return m_constraint;
            }

        private:
            void UpdatePhysicalControls(
                const double angularVelocity,
                StellarEquilibriumControlContext &context
            ) {
                m_generatedRotationChanged = !m_isPrepared || angularVelocity != m_angularVelocity;
                if (m_generatedRotationChanged) {
                    m_angularVelocity = angularVelocity;
                    ++m_generatedRotationDependency.revision;
                }
                context.dependencies.rotation = m_generatedRotationDependency;
                context.rotation              = m_constraint.GetCompiledConstraint().makeRotation(m_angularVelocity);
                ++context.rotationProviderCount;
                context.generatedPhysicalControl = true;
            }
            PreparedAngularMomentumOperator m_constraint;
            StellarEquilibriumDependencyStamp m_generatedRotationDependency;
            double m_angularVelocity{0.0};
            bool m_generatedRotationChanged{false};
            bool m_isPrepared{false};
            mutable mfem::Vector m_volumeDisplacementDirection;
            mutable mfem::Vector m_rotationalAngularVelocityAction;
            mutable mfem::Vector m_surfaceAngularVelocityAction;
            mutable mfem::Vector m_hydrostaticAngularVelocityAction;
            mutable mfem::Vector m_angularMomentumAction;
            mutable mfem::Vector m_zeroEnthalpy;
            MPI_Comm m_communicator{MPI_COMM_NULL};
        };

        struct FixedCentralDensityPreparationReport final {
            PreparedCentralDensityReport constraint;

            [[nodiscard]] bool DidAnyWork() const noexcept {
                return constraint.DidAnyWork();
            }
        };

        template <model::StellarModelType Model> class FixedCentralDensityRuntime final {
        public:
            using Report            = FixedCentralDensityPreparationReport;
            using PreparationResult = StellarEquilibriumPreparationResult<Report>;

            template <PreparedStellarEquilibriumPhysicalCore PhysicalCore>
            FixedCentralDensityRuntime(
                fem::FEM &finiteElements,
                const mapping::DomainMapper &,
                PhysicalCore &,
                const Model &model
            )
                : m_compiled(
                      models::compileConstraint(
                          model.template specification<models::FixedCentralDensity>(),
                          model.equationOfState()
                      )
                  ),
                  m_constraint(
                      MakeCenterDofMap(finiteElements),
                      finiteElements.mesh->GetComm()
                  ),
                  m_enthalpyResidual(m_constraint.GetCenterDof().field_size()),
                  m_phaseResidual(1),
                  m_enthalpyAction(m_constraint.GetCenterDof().field_size()),
                  m_phaseAction(1),
                  m_communicator(finiteElements.mesh->GetComm()) {
            }

            template <typename StateView>
            void ReadPhysicalControls(
                const StateView &,
                StellarEquilibriumControlContext &
            ) noexcept {
            }

            template <
                typename StateView,
                PreparedStellarEquilibriumPhysicalCore PhysicalCore>
            [[nodiscard]] Report PrepareAfterPhysical(
                const StateView &state,
                const StellarEquilibriumDependencies &dependencies,
                const PhysicalCore &physical
            ) {
                auto result = TryPrepareAfterPhysical(state, dependencies, physical);
                if (!result.has_value()) {
                    throwStellarEquilibriumPreparationRejection(result.error());
                }
                return std::move(result).value();
            }

            template <
                typename StateView,
                PreparedStellarEquilibriumPhysicalCore PhysicalCore>
            [[nodiscard]] PreparationResult TryPrepareAfterPhysical(
                const StateView &state,
                const StellarEquilibriumDependencies &dependencies,
                const PhysicalCore &
            ) {
                const auto enthalpy = state.block(utils::blocks::enthalpy_field.specific_term);
                const auto border   = state.block(utils::blocks::fixed_central_density_phase.central_value_term);
                MFEM_VERIFY(border.Size() == 1, "The central-density phase must provide one border coordinate.");

                m_isPrepared            = false;
                auto synchronizedBorder = synchronizeReplicatedControl(
                    border(0), m_communicator, "FixedCentralDensity could not synchronize its phase coordinate."
                );
                if (!synchronizedBorder.has_value()) {
                    return std::unexpected(synchronizedBorder.error());
                }

                auto report = m_constraint.Prepare(
                    m_compiled, enthalpy, *synchronizedBorder,
                    {.enthalpy = {
                         .identity = dependencies.enthalpy.identity, .revision = dependencies.enthalpy.revision
                     }}
                );
                m_isPrepared = true;
                return Report{.constraint = std::move(report)};
            }

            template <typename ResidualView> void AddResidual(const ResidualView &residual) const {
                m_enthalpyResidual = 0.0;
                m_phaseResidual    = 0.0;
                m_constraint.AddResidual(m_enthalpyResidual, m_phaseResidual);
                residual.add(utils::blocks::enthalpy_field.specific_term, m_enthalpyResidual);
                residual.add(utils::blocks::fixed_central_density_phase.central_value_term, m_phaseResidual);
            }

            template <
                typename DirectionView,
                typename ActionView,
                PreparedStellarEquilibriumPhysicalCore PhysicalCore>
            void AddJacobianAction(
                const DirectionView &direction,
                const ActionView &action,
                const PhysicalCore &
            ) const {
                const auto enthalpyDirection = direction.block(utils::blocks::enthalpy_field.specific_term);
                const auto borderDirection =
                    direction.block(utils::blocks::fixed_central_density_phase.central_value_term);
                MFEM_VERIFY(
                    borderDirection.Size() == 1 && std::isfinite(borderDirection(0)),
                    "The central-density phase direction must be finite."
                );
                m_enthalpyAction = 0.0;
                m_phaseAction    = 0.0;
                m_constraint.ApplyJacobian(
                    {.enthalpyVariation = enthalpyDirection, .borderVariation = borderDirection(0)},
                    {.enthalpyAction = m_enthalpyAction, .phaseAction = m_phaseAction}
                );
                action.add(utils::blocks::enthalpy_field.specific_term, m_enthalpyAction);
                action.add(utils::blocks::fixed_central_density_phase.central_value_term, m_phaseAction);
            }

            [[nodiscard]] bool IsPrepared() const noexcept {
                return m_isPrepared && m_constraint.IsPrepared();
            }

            [[nodiscard]] const PreparedCentralDensityConstraint &constraint() const noexcept {
                return m_constraint;
            }

            [[nodiscard]] const models::CompiledFixedCentralDensity &compiled() const noexcept {
                return m_compiled;
            }

        private:
            [[nodiscard]] static field::FieldPointDofMap MakeCenterDofMap(const fem::FEM &finiteElements) {
                using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema;
                MFEM_VERIFY(
                    finiteElements.mesh != nullptr && finiteElements.enthalpyFes != nullptr,
                    "The central-density phase requires a mesh and enthalpy finite-element space."
                );
                const field::FieldDofMap enthalpyMap =
                    field::make_field_dof_map<field::Enthalpy, DomainSchema>(*finiteElements.enthalpyFes);
                mfem::Vector origin(finiteElements.mesh->SpaceDimension());
                origin = 0.0;
                return field::make_field_point_dof_map<field::Enthalpy>(
                    *finiteElements.enthalpyFes, enthalpyMap, origin, 1.0e-12
                );
            }

            models::CompiledFixedCentralDensity m_compiled;
            PreparedCentralDensityConstraint m_constraint;
            mutable mfem::Vector m_enthalpyResidual;
            mutable mfem::Vector m_phaseResidual;
            mutable mfem::Vector m_enthalpyAction;
            mutable mfem::Vector m_phaseAction;
            MPI_Comm m_communicator{MPI_COMM_NULL};
            bool m_isPrepared{false};
        };
    } // namespace detail

    template <> struct StellarEquilibriumRuntimeContribution<eos::Polytrope> {
        static constexpr bool registered               = true;
        static constexpr std::size_t rotationProviders = 0;
        template <model::StellarModelType Model>
        using Prepared = detail::EmbeddedSpecificationRuntime<eos::Polytrope, Model>;
    };

    template <> struct StellarEquilibriumRuntimeContribution<surface::Isobaric> {
        static constexpr bool registered               = true;
        static constexpr std::size_t rotationProviders = 0;
        template <model::StellarModelType Model>
        using Prepared = detail::EmbeddedSpecificationRuntime<surface::Isobaric, Model>;
    };

    template <> struct StellarEquilibriumRuntimeContribution<models::FixedTotalMass> {
        static constexpr bool registered               = true;
        static constexpr std::size_t rotationProviders = 0;
        template <model::StellarModelType Model>
        using Prepared = detail::EmbeddedSpecificationRuntime<models::FixedTotalMass, Model>;
    };

    template <>
    struct StellarEquilibriumRuntimeContribution<models::FixedAngularMomentum>
        : PreparedStellarEquilibriumContribution<detail::FixedAngularMomentumRuntime, 1> { };

    template <>
    struct StellarEquilibriumRuntimeContribution<models::FixedCentralDensity>
        : PreparedStellarEquilibriumContribution<detail::FixedCentralDensityRuntime> { };

    namespace detail {
        template <typename Subset, typename Superset> struct BlockListIsSubset : std::false_type { };

        template <typename... Blocks, typename Superset>
        struct BlockListIsSubset<utils::blocks::type_list<Blocks...>, Superset>
            : std::bool_constant<(utils::blocks::contains_type_v<Blocks, Superset> && ...)> { };

        template <typename Blocks> struct SinglePhysicsBlock;

        template <typename Block> struct SinglePhysicsBlock<utils::blocks::type_list<Block>> final {
            using Type = Block;
        };

        template <typename ValueBlock> struct PhysicsFacingValueTerm final {
            using value = ValueBlock;
        };

        template <typename ResidualBlock> struct PhysicsFacingResidualTerm final {
            using residual = ResidualBlock;
        };

        template <typename Form, typename AllowedValueBlocks, models::ModelSpecification Specification>
        class RestrictedSpecificationStateView final {
        public:
            explicit RestrictedSpecificationStateView(const RootStateView<Form> &state) noexcept : m_state(state) {
            }

            template <typename NarrowedValueBlocks>
                requires BlockListIsSubset<
                    NarrowedValueBlocks,
                    AllowedValueBlocks>::value
            [[nodiscard]] auto narrow() const noexcept {
                return RestrictedSpecificationStateView<Form, NarrowedValueBlocks, Specification>{m_state};
            }

            template <typename Term>
                requires requires { typename std::remove_cvref_t<Term>::value; } &&
                         utils::blocks::contains_type_v<
                             typename std::remove_cvref_t<Term>::value,
                             AllowedValueBlocks>
            [[nodiscard]] auto block(const Term &term) const {
                return m_state.block(term);
            }

            /* Every compiler-enumerated derivative receives a view containing
             * exactly one source.  value() is the backend-agnostic spelling
             * for advanced physics vocabulary that does not yet have a named
             * convenience accessor. */
            [[nodiscard]] auto value() const
                requires(AllowedValueBlocks::size == 1)
            {
                using ValueBlock = typename SinglePhysicsBlock<AllowedValueBlocks>::Type;
                return m_state.block(PhysicsFacingValueTerm<ValueBlock>{});
            }

            [[nodiscard]] auto density() const
                requires StellarDependencyBlock<
                             Specification,
                             models::stellar::state::Density>::mapped
                         && utils::blocks::contains_type_v<
                                typename StellarDependencyBlock<
                                    Specification,
                                    models::stellar::state::Density>::Type,
                                AllowedValueBlocks>
            {
                return physicsBlock<models::stellar::state::Density>();
            }

            [[nodiscard]] auto surfaceShape() const
                requires StellarDependencyBlock<
                             Specification,
                             models::stellar::state::SurfaceShape>::mapped
                         && utils::blocks::contains_type_v<
                                typename StellarDependencyBlock<
                                    Specification,
                                    models::stellar::state::SurfaceShape>::Type,
                                AllowedValueBlocks>
            {
                return physicsBlock<models::stellar::state::SurfaceShape>();
            }

            [[nodiscard]] auto gravityGradient() const
                requires StellarDependencyBlock<
                             Specification,
                             models::stellar::state::GravityGradient>::mapped
                         && utils::blocks::contains_type_v<
                                typename StellarDependencyBlock<
                                    Specification,
                                    models::stellar::state::GravityGradient>::Type,
                                AllowedValueBlocks>
            {
                return physicsBlock<models::stellar::state::GravityGradient>();
            }

            [[nodiscard]] auto gravitationalPotential() const
                requires StellarDependencyBlock<
                             Specification,
                             models::stellar::state::GravitationalPotential>::mapped
                         && utils::blocks::contains_type_v<
                                typename StellarDependencyBlock<
                                    Specification,
                                    models::stellar::state::GravitationalPotential>::Type,
                                AllowedValueBlocks>
            {
                return physicsBlock<models::stellar::state::GravitationalPotential>();
            }

            [[nodiscard]] auto specificEnthalpy() const
                requires StellarDependencyBlock<
                             Specification,
                             models::stellar::state::SpecificEnthalpy>::mapped
                         && utils::blocks::contains_type_v<
                                typename StellarDependencyBlock<
                                    Specification,
                                    models::stellar::state::SpecificEnthalpy>::Type,
                                AllowedValueBlocks>
            {
                return physicsBlock<models::stellar::state::SpecificEnthalpy>();
            }

            [[nodiscard]] auto generatedCoordinate() const
                requires StellarDependencyBlock<
                             Specification,
                             models::stellar::state::OwnGeneratedCoordinate>::mapped
                         && utils::blocks::contains_type_v<
                                typename StellarDependencyBlock<
                                    Specification,
                                    models::stellar::state::OwnGeneratedCoordinate>::Type,
                                AllowedValueBlocks>
            {
                return physicsBlock<models::stellar::state::OwnGeneratedCoordinate>();
            }

            template <models::ModelSpecification Owner>
            [[nodiscard]] auto generatedCoordinate() const
                requires StellarDependencyBlock<
                             Specification,
                             models::stellar::state::GeneratedCoordinateOf<Owner>>::mapped
                         && utils::blocks::contains_type_v<
                                typename StellarDependencyBlock<
                                    Specification,
                                    models::stellar::state::GeneratedCoordinateOf<Owner>>::Type,
                                AllowedValueBlocks>
            {
                return physicsBlock<models::stellar::state::GeneratedCoordinateOf<Owner>>();
            }

        private:
            template <typename PhysicsQuantity> [[nodiscard]] auto physicsBlock() const {
                using ValueBlock = typename StellarDependencyBlock<Specification, PhysicsQuantity>::Type;
                return m_state.block(PhysicsFacingValueTerm<ValueBlock>{});
            }

            RootStateView<Form> m_state;
        };

        template <typename Form, typename AllowedResidualBlocks> class RestrictedSpecificationResidualView final {
        public:
            RestrictedSpecificationResidualView(
                const ResidualView<Form> &residual,
                const PreparedPressureSurfaceConstraint &surfaceConstraint
            ) noexcept
                : m_residual(residual),
                  m_surfaceConstraint(std::addressof(surfaceConstraint)) {
            }

            template <typename OtherAllowedResidualBlocks>
                requires BlockListIsSubset<
                             AllowedResidualBlocks,
                             OtherAllowedResidualBlocks>::value
            explicit RestrictedSpecificationResidualView(
                const RestrictedSpecificationResidualView<
                    Form,
                    OtherAllowedResidualBlocks> &residual
            ) noexcept
                : m_residual(residual.m_residual),
                  m_surfaceConstraint(residual.m_surfaceConstraint) {
            }

            /* Physics-facing assembly is intentionally additive-only.  A
             * specification cannot erase the physical core or an earlier
             * contribution.  Hydrostatic additions are also projected away
             * from rows owned by the surface condition, so extension authors
             * do not need to understand the backend row-replacement policy. */
            template <typename Term>
                requires requires { typename std::remove_cvref_t<Term>::residual; } &&
                         utils::blocks::contains_type_v<
                             typename std::remove_cvref_t<Term>::residual,
                             AllowedResidualBlocks>
            void
            add(const Term &term,
                const double contribution) const {
                mfem::Vector destination = m_residual.block(term);
                destination += contribution;
                RestoreReplacedRows<typename std::remove_cvref_t<Term>::residual>(destination, contribution);
                destination.SyncAliasMemory(m_residual.vector());
            }

            template <typename Term>
                requires requires { typename std::remove_cvref_t<Term>::residual; } &&
                         utils::blocks::contains_type_v<
                             typename std::remove_cvref_t<Term>::residual,
                             AllowedResidualBlocks>
            void
            add(const Term &term,
                const mfem::Vector &contribution) const {
                mfem::Vector destination = m_residual.block(term);
                if (destination.Size() != contribution.Size()) {
                    throw std::invalid_argument(
                        "A specification runtime contribution has the wrong residual block size."
                    );
                }
                destination += contribution;
                RestoreReplacedRows<typename std::remove_cvref_t<Term>::residual>(destination, contribution);
                destination.SyncAliasMemory(m_residual.vector());
            }

            template <typename Term>
                requires requires { typename std::remove_cvref_t<Term>::residual; } &&
                         utils::blocks::contains_type_v<
                             typename std::remove_cvref_t<Term>::residual,
                             AllowedResidualBlocks>
            void addEntry(
                const Term &term,
                const int index,
                const double contribution
            ) const {
                mfem::Vector destination = m_residual.block(term);
                if (index < 0 || index >= destination.Size()) {
                    throw std::out_of_range("A specification runtime contribution selected an invalid residual entry.");
                }
                if (!IsReplacedRow<typename std::remove_cvref_t<Term>::residual>(index)) {
                    destination(index) += contribution;
                }
                destination.SyncAliasMemory(m_residual.vector());
            }

        private:
            template <typename, typename> friend class RestrictedSpecificationResidualView;

            template <typename ResidualBlock> [[nodiscard]] bool IsReplacedRow(const int index) const noexcept {
                if constexpr (!std::same_as<ResidualBlock, utils::blocks::enthalpy::specific::residual>) {
                    return false;
                } else {
                    for (const int row : m_surfaceConstraint->GetSurfaceRows().reduced_dofs()) {
                        if (row == index) {
                            return true;
                        }
                    }
                    return false;
                }
            }

            template <typename ResidualBlock>
            void RestoreReplacedRows(
                mfem::Vector &destination,
                const double contribution
            ) const {
                if constexpr (std::same_as<ResidualBlock, utils::blocks::enthalpy::specific::residual>) {
                    for (const int row : m_surfaceConstraint->GetSurfaceRows().reduced_dofs()) {
                        destination(row) -= contribution;
                    }
                }
            }

            template <typename ResidualBlock>
            void RestoreReplacedRows(
                mfem::Vector &destination,
                const mfem::Vector &contribution
            ) const {
                if constexpr (std::same_as<ResidualBlock, utils::blocks::enthalpy::specific::residual>) {
                    for (const int row : m_surfaceConstraint->GetSurfaceRows().reduced_dofs()) {
                        destination(row) -= contribution(row);
                    }
                }
            }

            ResidualView<Form> m_residual;
            const PreparedPressureSurfaceConstraint *m_surfaceConstraint;
        };

        /*
         * Single-use row handed to one compiler-enumerated physics provider.
         * It deliberately has no row selector: the equation tag selected the
         * row before the extension was called.  A provider therefore cannot
         * redirect a legal source into a different legal residual.  Runtime
         * accounting additionally rejects double assembly or a contribution
         * token inconsistent with what the provider actually did.
         */
        template <typename Form, typename ResidualBlock> class ExactSpecificationResidualRow final {
        public:
            using View = RestrictedSpecificationResidualView<Form, utils::blocks::type_list<ResidualBlock>>;

            explicit ExactSpecificationResidualRow(const View &row) noexcept : m_row(row) {
            }

            ExactSpecificationResidualRow(const ExactSpecificationResidualRow &)            = delete;
            ExactSpecificationResidualRow &operator=(const ExactSpecificationResidualRow &) = delete;
            ExactSpecificationResidualRow(ExactSpecificationResidualRow &&)                 = delete;
            ExactSpecificationResidualRow &operator=(ExactSpecificationResidualRow &&)      = delete;

            [[nodiscard]] stellar::ContributionAdded add(const double contribution) {
                RequireUnused();
                m_row.add(PhysicsFacingResidualTerm<ResidualBlock>{}, contribution);
                m_addCount = 1;
                return {};
            }

            [[nodiscard]] stellar::ContributionAdded add(const mfem::Vector &contribution) {
                RequireUnused();
                m_row.add(PhysicsFacingResidualTerm<ResidualBlock>{}, contribution);
                m_addCount = 1;
                return {};
            }

            [[nodiscard]] stellar::ContributionAdded addEntry(
                const int index,
                const double contribution
            ) {
                RequireUnused();
                m_row.addEntry(PhysicsFacingResidualTerm<ResidualBlock>{}, index, contribution);
                m_addCount = 1;
                return {};
            }

            template <stellar::ContributionResult Result> void Verify(const Result &) const {
                if constexpr (std::same_as<std::remove_cvref_t<Result>, stellar::ContributionAdded>) {
                    if (m_addCount != 1) {
                        throw std::logic_error(
                            "A stellar physics provider returned ContributionAdded without adding exactly once."
                        );
                    }
                } else {
                    if (m_addCount != 0) {
                        throw std::logic_error(
                            "A stellar physics provider returned StructuralZero after adding to its row."
                        );
                    }
                }
            }

        private:
            void RequireUnused() const {
                if (m_addCount != 0) {
                    throw std::logic_error(
                        "A compiler-enumerated stellar residual/Jacobian edge may be assembled only once."
                    );
                }
            }

            View m_row;
            int m_addCount{0};
        };

        template <typename ModelList> struct PhysicsTagList;

        template <typename... Tags> struct PhysicsTagList<models::ModelTypeList<Tags...>> final {
            using Type = utils::blocks::type_list<Tags...>;
        };

        template <typename Equation, typename States> struct DerivativesOfEquation;

        template <typename Equation, typename... States>
        struct DerivativesOfEquation<Equation, utils::blocks::type_list<States...>> final {
            using Type = utils::blocks::type_list<models::stellar::Derivative<Equation, States>...>;
        };

        template <typename Equations, typename States> struct CartesianPhysicsDerivatives;

        template <typename... Equations, typename States>
        struct CartesianPhysicsDerivatives<utils::blocks::type_list<Equations...>, States> final {
            using Type = ConcatenateBlockListsT<typename DerivativesOfEquation<Equations, States>::Type...>;
        };

        /* The topology consumed by physics-facing providers is generated
         * directly from the same Reads/Changes declaration used by the block
         * compiler.  No implementation-owned provider list can get out of
         * sync with the model declaration. */
        template <models::ModelSpecification Specification> struct SpecificationPhysicsTopology final {
        private:
            using Contribution = models::SpecificationContribution<Specification>;

        public:
            using ReadStates        = typename PhysicsTagList<typename Contribution::DependsOn>::Type;
            using ChangedEquations  = typename PhysicsTagList<typename Contribution::Affects>::Type;
            using OwnGeneratedState = std::conditional_t<
                Contribution::generatedValueArity == 0,
                utils::blocks::type_list<>,
                utils::blocks::type_list<models::stellar::state::OwnGeneratedCoordinate>>;
            using OwnConstraintEquation = std::conditional_t<
                Contribution::generatedResidualArity == 0,
                utils::blocks::type_list<>,
                utils::blocks::type_list<models::stellar::equation::OwnConstraint>>;

            using ResidualEquations     = UniqueConcatenateBlockListsT<OwnConstraintEquation, ChangedEquations>;
            using ChangedEquationInputs = UniqueConcatenateBlockListsT<ReadStates, OwnGeneratedState>;
            using ConstraintDerivatives = typename CartesianPhysicsDerivatives<OwnConstraintEquation, ReadStates>::Type;
            using ChangedEquationDerivatives =
                typename CartesianPhysicsDerivatives<ChangedEquations, ChangedEquationInputs>::Type;
            using Derivatives = UniqueConcatenateBlockListsT<ConstraintDerivatives, ChangedEquationDerivatives>;
        };

        template <
            typename Specification,
            typename State,
            bool IsSpecification = models::ModelSpecification<std::remove_cvref_t<Specification>>>
        struct SpecificationReadsState : std::false_type { };

        template <typename Specification, typename State>
        struct SpecificationReadsState<Specification, State, true>
            : std::bool_constant<utils::blocks::contains_type_v<
                  State,
                  typename SpecificationPhysicsTopology<std::remove_cvref_t<Specification>>::ReadStates>> { };
    } // namespace detail

    /**
     * Astronomy-facing access to the current core's physical-volume density
     * integral,
     *
     *     M[rho] = integral_{Omega_star} rho dV.
     *
     * The context is available only to a specification which declares both
     * stellar::state::Density and stellar::state::SurfaceShape in Reads.  That
     * is a mathematical requirement rather than an implementation detail:
     * the mapped physical-volume integral depends on both rho and the domain
     * geometry.  Requiring both declarations prevents a residual from using
     * this service while omitting its geometry column from the inferred
     * Jacobian.  The context is a small, copyable, non-owning service handle:
     * it exposes neither FEM objects nor the physical core, owns no backend
     * object, and uses no allocating type erasure. The selected core remains
     * responsible for quadrature, mapped physical volume, distributed
     * reduction, scratch storage, and exact directional actions.
     */
    template <typename Specification>
    concept DensityVolumeIntegralSpecification =
        detail::SpecificationReadsState<Specification, models::stellar::state::Density>::value &&
        detail::SpecificationReadsState<Specification, models::stellar::state::SurfaceShape>::value;

    template <DensityVolumeIntegralSpecification Specification> class DensityVolumeIntegralContext final {
    public:
        [[nodiscard]] dimensions::MassValue integrateDensity(const mfem::Vector &density) const {
            return dimensions::MassValue{m_densityAction(m_core, density)};
        }

        [[nodiscard]] dimensions::MassValue linearizeDensityIntegral(const mfem::Vector &densityDirection) const {
            return dimensions::MassValue{m_densityAction(m_core, densityDirection)};
        }

        [[nodiscard]] dimensions::MassValue
        linearizeSurfaceShapeIntegral(const mfem::Vector &surfaceShapeDirection) const {
            return dimensions::MassValue{m_surfaceShapeAction(m_core, surfaceShapeDirection)};
        }

    private:
        template <
            models::ModelSpecification OtherSpecification,
            model::StellarModelType OtherModel,
            typename OtherPhysics>
        friend class detail::PhysicsFacingSpecificationRuntime;

        using Action = double (*)(
            const void *,
            const mfem::Vector &
        );

    public:
        /* Public construction is intentionally backend-facing: it accepts a
         * core which already owns the prepared integration service, but the
         * resulting physics-facing handle has no route back to that core. This
         * also permits direct distributed contract tests of the service. */
        template <typename PhysicalCore>
            requires requires(
                         const PhysicalCore &core,
                         const mfem::Vector &direction
                     ) {
                { core.ApplyDensityVolumeIntegralDensityAction(direction) } -> std::convertible_to<double>;
                { core.ApplyDensityVolumeIntegralSurfaceShapeAction(direction) } -> std::convertible_to<double>;
            }
        explicit DensityVolumeIntegralContext(const PhysicalCore &core) noexcept
            : m_core(std::addressof(core)),
              m_densityAction(&ApplyDensityAction<PhysicalCore>),
              m_surfaceShapeAction(&ApplySurfaceShapeAction<PhysicalCore>) {
        }

    private:
        template <typename PhysicalCore>
        [[nodiscard]] static double ApplyDensityAction(
            const void *untypedCore,
            const mfem::Vector &densityDirection
        ) {
            const auto &core = *static_cast<const PhysicalCore *>(untypedCore);
            return core.ApplyDensityVolumeIntegralDensityAction(densityDirection);
        }

        template <typename PhysicalCore>
        [[nodiscard]] static double ApplySurfaceShapeAction(
            const void *untypedCore,
            const mfem::Vector &surfaceShapeDirection
        ) {
            const auto &core = *static_cast<const PhysicalCore *>(untypedCore);
            return core.ApplyDensityVolumeIntegralSurfaceShapeAction(surfaceShapeDirection);
        }

        const void *m_core{nullptr};
        Action m_densityAction{nullptr};
        Action m_surfaceShapeAction{nullptr};
    };

    namespace detail {
        template <
            models::ModelSpecification Specification,
            typename Physics,
            bool ContextAvailable = DensityVolumeIntegralSpecification<Specification>>
        struct DensityVolumeIntegralPhysicsConstruction final {
            using Context                       = void;
            static constexpr bool constructible = false;
        };

        template <models::ModelSpecification Specification, typename Physics>
        struct DensityVolumeIntegralPhysicsConstruction<Specification, Physics, true> final {
            using Context                       = DensityVolumeIntegralContext<Specification>;
            static constexpr bool constructible = std::constructible_from<Physics, const Specification &, Context>;
        };

        template <typename PhysicalCore>
        concept DensityVolumeIntegralCore = requires(const PhysicalCore &core, const mfem::Vector &direction) {
            { core.ApplyDensityVolumeIntegralDensityAction(direction) } -> std::convertible_to<double>;
            { core.ApplyDensityVolumeIntegralSurfaceShapeAction(direction) } -> std::convertible_to<double>;
        };

        /*
         * A Jacobian callback is stricter than residual assembly.  Its
         * row/source pair is checked against the compiler output, then the
         * callback receives a direction view containing only that source and
         * an additive action view containing only that row.  Consequently,
         * independently legal endpoints cannot accidentally be recombined
         * into an undeclared edge inside one callback.  As with any assembly
         * API, this structural contract does not attempt to prove that the
         * callback's arithmetic is the mathematical derivative it claims.
         */
        template <
            typename Form,
            typename AllowedCouplings,
            typename AllowedValueBlocks,
            models::ModelSpecification Specification,
            typename AllowedResidualBlocks>
        class RestrictedSpecificationJacobianView final {
        public:
            RestrictedSpecificationJacobianView(
                const RestrictedSpecificationStateView<
                    Form,
                    AllowedValueBlocks,
                    Specification> &direction,
                const RestrictedSpecificationResidualView<
                    Form,
                    AllowedResidualBlocks> &action
            ) noexcept
                : m_direction(direction),
                  m_action(action) {
            }

            template <
                typename ResidualTerm,
                typename ValueTerm,
                typename Callback>
                requires requires {
                    typename std::remove_cvref_t<ResidualTerm>::residual;
                    typename std::remove_cvref_t<ValueTerm>::value;
                } &&
                         utils::blocks::contains_type_v<
                             typename std::remove_cvref_t<ValueTerm>::value,
                             AllowedValueBlocks> &&
                         utils::blocks::contains_type_v<
                             typename std::remove_cvref_t<ResidualTerm>::residual,
                             AllowedResidualBlocks> &&
                         utils::blocks::contains_type_v<
                             StellarEquilibriumJacobianCoupling<
                                 typename std::remove_cvref_t<ResidualTerm>::residual,
                                 typename std::remove_cvref_t<ValueTerm>::value>,
                             AllowedCouplings> &&
                         requires(
                             Callback &&callback,
                             RestrictedSpecificationStateView<
                                 Form,
                                 utils::blocks::type_list<typename std::remove_cvref_t<ValueTerm>::value>,
                                 Specification> &direction,
                             RestrictedSpecificationResidualView<
                                 Form,
                                 utils::blocks::type_list<typename std::remove_cvref_t<ResidualTerm>::residual>> &action
                         ) {
                             { std::forward<Callback>(callback)(direction, action) } -> std::same_as<void>;
                         }
            void
            add(const ResidualTerm &,
                const ValueTerm &,
                Callback &&callback) const {
                using DirectionView = RestrictedSpecificationStateView<
                    Form, utils::blocks::type_list<typename std::remove_cvref_t<ValueTerm>::value>, Specification>;
                using ActionView = RestrictedSpecificationResidualView<
                    Form, utils::blocks::type_list<typename std::remove_cvref_t<ResidualTerm>::residual>>;

                DirectionView direction =
                    m_direction
                        .template narrow<utils::blocks::type_list<typename std::remove_cvref_t<ValueTerm>::value>>();
                ActionView action{m_action};
                std::forward<Callback>(callback)(direction, action);
            }

        private:
            RestrictedSpecificationStateView<Form, AllowedValueBlocks, Specification> m_direction;
            RestrictedSpecificationResidualView<Form, AllowedResidualBlocks> m_action;
        };

        template <models::ModelSpecification Specification, typename Tags> struct CompilePhysicsTags;

        template <models::ModelSpecification Specification, typename... Tags>
        struct CompilePhysicsTags<Specification, utils::blocks::type_list<Tags...>> final {
            using Type = utils::blocks::type_list<typename StellarDependencyBlock<Specification, Tags>::Type...>;
            static constexpr bool complete = (StellarDependencyBlock<Specification, Tags>::mapped && ...);
        };

        template <models::ModelSpecification Specification, typename Derivatives> struct CompilePhysicsDerivatives;

        template <models::ModelSpecification Specification, typename Derivative> struct CompilePhysicsDerivative;

        template <models::ModelSpecification Specification, typename Equation, typename State>
        struct CompilePhysicsDerivative<Specification, models::stellar::Derivative<Equation, State>> final {
            using Type = StellarEquilibriumJacobianCoupling<
                typename StellarDependencyBlock<Specification, Equation>::Type,
                typename StellarDependencyBlock<Specification, State>::Type>;
            static constexpr bool complete = StellarDependencyBlock<Specification, Equation>::mapped &&
                                             StellarDependencyBlock<Specification, State>::mapped;
        };

        template <models::ModelSpecification Specification, typename... Derivatives>
        struct CompilePhysicsDerivatives<Specification, utils::blocks::type_list<Derivatives...>> final {
            using Type =
                utils::blocks::type_list<typename CompilePhysicsDerivative<Specification, Derivatives>::Type...>;
            static constexpr bool complete = (CompilePhysicsDerivative<Specification, Derivatives>::complete && ...);
        };

        template <models::ModelSpecification Specification, model::StellarModelType Model>
        struct SpecificationRuntimeAccess final {
            using SpecificationType = Specification;
            using Compilation       = StellarEquilibriumSpecificationCompilation<Specification>;
            using Topology          = SpecificationPhysicsTopology<Specification>;
            using Form              = CompiledStellarEquilibriumForm<Model>;
            using ValueBlocks       = UniqueConcatenateBlockListsT<
                typename Compilation::GeneratedValueBlocks,
                typename Compilation::DependsOnValueBlocks>;
            using ResidualBlocks = UniqueConcatenateBlockListsT<
                typename Compilation::GeneratedResidualBlocks,
                typename Compilation::AffectedResidualBlocks>;
            using StateView    = RestrictedSpecificationStateView<Form, ValueBlocks, Specification>;
            using ResidualView = RestrictedSpecificationResidualView<Form, ResidualBlocks>;

            template <typename StateTag>
            using ValueBlockFor = typename StellarDependencyBlock<Specification, StateTag>::Type;

            template <typename EquationTag>
            using ResidualBlockFor = typename StellarDependencyBlock<Specification, EquationTag>::Type;

            template <typename StateTag>
            using DirectionView = RestrictedSpecificationStateView<
                Form,
                utils::blocks::type_list<ValueBlockFor<StateTag>>,
                Specification>;

            template <typename EquationTag>
            using RowView =
                RestrictedSpecificationResidualView<Form, utils::blocks::type_list<ResidualBlockFor<EquationTag>>>;

            template <typename EquationTag>
            using Row          = ExactSpecificationResidualRow<Form, ResidualBlockFor<EquationTag>>;

            using JacobianView = RestrictedSpecificationJacobianView<
                Form,
                typename Compilation::JacobianCouplings,
                ValueBlocks,
                Specification,
                ResidualBlocks>;

            using ProviderResidualBlocks =
                typename CompilePhysicsTags<Specification, typename Topology::ResidualEquations>::Type;
            using ProviderJacobianCouplings =
                typename CompilePhysicsDerivatives<Specification, typename Topology::Derivatives>::Type;

            static_assert(CompilePhysicsTags<
                          Specification,
                          typename Topology::ResidualEquations>::complete);
            static_assert(CompilePhysicsDerivatives<
                          Specification,
                          typename Topology::Derivatives>::complete);
            static_assert(std::same_as<
                          ProviderResidualBlocks,
                          ResidualBlocks>);
            static_assert(std::same_as<
                          ProviderJacobianCouplings,
                          typename Compilation::JacobianCouplings>);
        };

        template <typename Physics, typename Access, typename Equations> struct ExactResidualProviderSet;

        template <typename Physics, typename Access, typename... Equations>
        struct ExactResidualProviderSet<Physics, Access, utils::blocks::type_list<Equations...>> final {
            static constexpr bool complete =
                (requires(const Physics &physics, typename Access::template Row<Equations> &row) {
                    { physics.AddResidual(Equations{}, row) } -> stellar::ContributionResult;
                } && ...);

            static void Apply(
                const Physics &physics,
                const typename Access::ResidualView &residual
            )
                requires complete
            {
                (ApplyOne<Equations>(physics, residual), ...);
            }

        private:
            template <typename Equation>
            static void ApplyOne(
                const Physics &physics,
                const typename Access::ResidualView &residual
            ) {
                typename Access::template RowView<Equation> rowView{residual};
                typename Access::template Row<Equation> row{rowView};
                decltype(auto) result = physics.AddResidual(Equation{}, row);
                row.Verify(result);
            }
        };

        template <typename Physics, typename Access, typename Derivatives> struct ExactJacobianProviderSet;

        template <typename Physics, typename Access, typename... Derivatives>
        struct ExactJacobianProviderSet<Physics, Access, utils::blocks::type_list<Derivatives...>> final {
        private:
            template <typename Derivative> struct Traits;

            template <typename Equation, typename State>
            struct Traits<models::stellar::Derivative<Equation, State>> final {
                using EquationTag = Equation;
                using StateTag    = State;
            };

            template <typename Derivative> [[nodiscard]] static consteval bool ProviderIsComplete() {
                using Equation = typename Traits<Derivative>::EquationTag;
                using State    = typename Traits<Derivative>::StateTag;
                return requires(
                    const Physics &physics, const typename Access::template DirectionView<State> &direction,
                    typename Access::template Row<Equation> &row
                ) {
                    { physics.AddJacobianAction(Derivative{}, direction, row) } -> stellar::ContributionResult;
                };
            }

        public:
            static constexpr bool complete = (ProviderIsComplete<Derivatives>() && ...);

            static void Apply(
                const Physics &physics,
                const typename Access::StateView &direction,
                const typename Access::ResidualView &action
            )
                requires complete
            {
                (ApplyOne<Derivatives>(physics, direction, action), ...);
            }

        private:
            template <typename Derivative>
            static void ApplyOne(
                const Physics &physics,
                const typename Access::StateView &direction,
                const typename Access::ResidualView &action
            ) {
                using Equation = typename Traits<Derivative>::EquationTag;
                using State    = typename Traits<Derivative>::StateTag;
                typename Access::template DirectionView<State> source =
                    direction
                        .template narrow<utils::blocks::type_list<typename Access::template ValueBlockFor<State>>>();
                typename Access::template RowView<Equation> rowView{action};
                typename Access::template Row<Equation> row{rowView};
                decltype(auto) result = physics.AddJacobianAction(Derivative{}, source, row);
                row.Verify(result);
            }
        };

        template <typename Physics, typename Access>
        inline constexpr bool exactSpecificationPhysicsProvidersComplete =
            ExactResidualProviderSet<Physics, Access, typename Access::Topology::ResidualEquations>::complete &&
            ExactJacobianProviderSet<Physics, Access, typename Access::Topology::Derivatives>::complete;

        /* Classify one exact nested derivative provider by its return type.
         * This information is useful outside residual assembly as well: a
         * preconditioner may omit a compiler-declared core-to-core edge only
         * when the physics implementation itself proves that the edge is the
         * identically zero map.  The primary remains well formed so capability
         * queries for malformed providers fail normally rather than producing
         * diagnostics deep in a factory body. */
        template <typename Physics, typename Access, typename Derivative, typename = void>
        struct ExactJacobianProviderResult final {
            using Coupling                       = void;
            using Result                         = void;

            static constexpr bool complete       = false;
            static constexpr bool structuralZero = false;
        };

        template <typename Physics, typename Access, typename Equation, typename State>
        struct ExactJacobianProviderResult<
            Physics,
            Access,
            models::stellar::Derivative<Equation, State>,
            std::void_t<decltype(std::declval<const Physics &>().AddJacobianAction(
                models::stellar::Derivative<Equation, State>{},
                std::declval<const typename Access::template DirectionView<State> &>(),
                std::declval<typename Access::template Row<Equation> &>()
            ))>>
            final {
            using Derivative = models::stellar::Derivative<Equation, State>;
            using Coupling   = typename CompilePhysicsDerivative<typename Access::SpecificationType, Derivative>::Type;
            using Result     = decltype(std::declval<const Physics &>().AddJacobianAction(
                Derivative{},
                std::declval<const typename Access::template DirectionView<State> &>(),
                std::declval<typename Access::template Row<Equation> &>()
            ));

            static constexpr bool complete = stellar::ContributionResult<Result>;
            static constexpr bool structuralZero =
                complete && std::same_as<std::remove_cvref_t<Result>, stellar::StructuralZero>;
        };

        template <typename Physics, typename Access, typename Coupling, typename Derivatives>
        struct ExactCouplingProvidersAreStructuralZero;

        template <typename Physics, typename Access, typename Coupling, typename... Derivatives>
        struct ExactCouplingProvidersAreStructuralZero<
            Physics,
            Access,
            Coupling,
            utils::blocks::type_list<Derivatives...>>
            final {
        private:
            template <typename Derivative> using Provider = ExactJacobianProviderResult<Physics, Access, Derivative>;

            static constexpr bool hasMatchingProvider =
                (false || ... || std::same_as<Coupling, typename Provider<Derivatives>::Coupling>);
            static constexpr bool everyMatchingProviderIsZero =
                (true && ... &&
                 (!std::same_as<Coupling, typename Provider<Derivatives>::Coupling> ||
                  Provider<Derivatives>::structuralZero));

        public:
            static constexpr bool value = hasMatchingProvider && everyMatchingProviderIsZero;
        };

        template <typename Specification, typename Model, typename Coupling, typename = void>
        struct NestedSpecificationCouplingIsStructuralZero : std::false_type { };

        template <models::ModelSpecification Specification, model::StellarModelType Model, typename Coupling>
            requires Model::template
        containsSpecification<Specification> struct NestedSpecificationCouplingIsStructuralZero<
            Specification,
            Model,
            Coupling,
            std::void_t<typename RuntimeContributionSelection<Specification, Model>::Physics>>
            final : std::bool_constant<
                        RuntimeContributionSelection<Specification, Model>::available &&
                        RuntimeContributionSelection<Specification, Model>::registered &&
                        !RuntimeContributionSelection<Specification, Model>::ambiguous &&
                        ExactCouplingProvidersAreStructuralZero<
                            typename RuntimeContributionSelection<Specification, Model>::Physics,
                            SpecificationRuntimeAccess<Specification, Model>,
                            Coupling,
                            typename SpecificationPhysicsTopology<Specification>::Derivatives>::value> { };

        /*
         * Adapter for the physics-facing nested extension protocol.
         *
         * The outer runtime constructs this object through the same internal
         * slot interface as trusted backend contributions.  The authored
         * Physics object behind it sees a deliberately smaller interface:
         *
         *   Physics(const Specification &)
         *   // or, only with Reads<Density, SurfaceShape> and a supporting core:
         *   Physics(const Specification &, DensityVolumeIntegralContext<Specification>)
         *   Report PrepareAfterPhysical(const StateView &)
         *   ContributionResult AddResidual(EquationTag, Row &) const
         *   ContributionResult AddJacobianAction(
         *       stellar::Derivative<EquationTag, StateTag>,
         *       const one-source DirectionView &,
         *       Row &) const
         *
         * The adapter invokes those overloads once for every row/edge inferred
         * from Reads/Changes. row.add(...) returns the required success token;
         * an identically absent term must return stellar::structuralZero.
         * Missing overloads fail the capability query at compile time, while
         * double assembly and inconsistent result tokens fail immediately at
         * runtime. Physics authors never enumerate a backend provider list or
         * see an aggregate direction/action object.
         *
         *   bool IsPrepared() const
         *
         * A nested contribution which owns rotation additionally supplies
         * `RigidRotation GenerateRotation(const StateView &)`.  The adapter,
         * rather than the extension, owns the dependency stamp and compares
         * successive rotations.  In particular, no nested implementation is
         * ever handed the model, FEM/mapper objects, dependency set, control
         * context, or physical core. The optional integral context is a
         * read-only service handle and does not expose any of those objects.
         */
        template <models::ModelSpecification Specification, model::StellarModelType Model, typename Physics>
        class PhysicsFacingSpecificationRuntime final {
        private:
            using Access               = SpecificationRuntimeAccess<Specification, Model>;
            using IntegralConstruction = DensityVolumeIntegralPhysicsConstruction<Specification, Physics>;
            using ResidualProviders =
                ExactResidualProviderSet<Physics, Access, typename Access::Topology::ResidualEquations>;
            using JacobianProviders = ExactJacobianProviderSet<Physics, Access, typename Access::Topology::Derivatives>;
            static constexpr std::size_t rotationProviders =
                RuntimeContributionSelection<Specification, Model>::rotationProviders;

        public:
            using Report = typename Physics::Report;

            template <PreparedStellarEquilibriumPhysicalCore PhysicalCore>
                requires(IntegralConstruction::constructible && DensityVolumeIntegralCore<PhysicalCore>)
            PhysicsFacingSpecificationRuntime(
                fem::FEM &,
                const mapping::DomainMapper &,
                PhysicalCore &physical,
                const Model &model
            )
                : m_physics(
                      model.template specification<Specification>(),
                      typename IntegralConstruction::Context{physical}
                  ) {
                m_generatedRotationDependency.identity =
                    static_cast<std::uint64_t>(reinterpret_cast<std::uintptr_t>(this));
            }

            template <PreparedStellarEquilibriumPhysicalCore PhysicalCore>
                requires(
                    std::constructible_from<
                        Physics,
                        const Specification &> &&
                    (!IntegralConstruction::constructible || !DensityVolumeIntegralCore<PhysicalCore>)
                )
            PhysicsFacingSpecificationRuntime(
                fem::FEM &,
                const mapping::DomainMapper &,
                PhysicalCore &,
                const Model &model
            )
                : m_physics(model.template specification<Specification>()) {
                m_generatedRotationDependency.identity =
                    static_cast<std::uint64_t>(reinterpret_cast<std::uintptr_t>(this));
            }

            template <typename StateView>
                requires(
                    rotationProviders == 0 ||
                    requires(
                        Physics &implementation,
                        const StateView &state
                    ) {
                        { implementation.GenerateRotation(state) } -> std::same_as<physics::RigidRotation>;
                    }
                )
            void ReadPhysicalControls(
                const StateView &state,
                StellarEquilibriumControlContext &context
            ) {
                if constexpr (rotationProviders == 1) {
                    physics::RigidRotation rotation = m_physics.GenerateRotation(state);
                    const bool changed =
                        !m_generatedRotation.has_value() || !SameRotation(*m_generatedRotation, rotation);
                    if (changed) {
                        m_generatedRotation = rotation;
                        ++m_generatedRotationDependency.revision;
                    }
                    context.dependencies.rotation = m_generatedRotationDependency;
                    context.rotation              = std::move(rotation);
                    ++context.rotationProviderCount;
                    context.generatedPhysicalControl = context.generatedPhysicalControl || changed;
                }
            }

            template <
                typename StateView,
                PreparedStellarEquilibriumPhysicalCore PhysicalCore>
                requires requires(
                    Physics &implementation,
                    const StateView &state
                ) {
                    { implementation.PrepareAfterPhysical(state) } -> std::same_as<Report>;
                }
            [[nodiscard]] Report PrepareAfterPhysical(
                const StateView &state,
                const StellarEquilibriumDependencies &,
                const PhysicalCore &
            ) {
                return m_physics.PrepareAfterPhysical(state);
            }

            void AddResidual(const typename Access::ResidualView &residual) const
                requires ResidualProviders::complete
            {
                ResidualProviders::Apply(m_physics, residual);
            }

            template <PreparedStellarEquilibriumPhysicalCore PhysicalCore>
                requires JacobianProviders::complete
            void AddJacobianAction(
                const typename Access::StateView &direction,
                const typename Access::ResidualView &action,
                const PhysicalCore &
            ) const {
                JacobianProviders::Apply(m_physics, direction, action);
            }

            [[nodiscard]] bool IsPrepared() const noexcept
                requires requires(const Physics &implementation) {
                    { implementation.IsPrepared() } -> std::convertible_to<bool>;
                }
            {
                return static_cast<bool>(m_physics.IsPrepared());
            }

            [[nodiscard]] Physics &physics() noexcept {
                return m_physics;
            }

            [[nodiscard]] const Physics &physics() const noexcept {
                return m_physics;
            }

        private:
            [[nodiscard]] static bool SameVector(
                const mfem::Vector &left,
                const mfem::Vector &right
            ) noexcept {
                if (left.Size() != right.Size()) {
                    return false;
                }
                for (int component = 0; component < left.Size(); ++component) {
                    if (left(component) != right(component)) {
                        return false;
                    }
                }
                return true;
            }

            [[nodiscard]] static bool SameRotation(
                const physics::RigidRotation &left,
                const physics::RigidRotation &right
            ) noexcept {
                return SameVector(left.angular_velocity(), right.angular_velocity()) &&
                       SameVector(left.center(), right.center());
            }

            Physics m_physics;
            std::optional<physics::RigidRotation> m_generatedRotation;
            StellarEquilibriumDependencyStamp m_generatedRotationDependency{};
        };

        template <
            model::StellarModelType Model,
            PreparedStellarEquilibriumPhysicalCore PhysicalCore,
            typename SpecificationSet>
        class PreparedSpecificationSet;

        template <model::StellarModelType Model, PreparedStellarEquilibriumPhysicalCore PhysicalCore>
        class PreparedSpecificationSet<Model, PhysicalCore, models::detail::SpecificationSetStorage<>> final {
        public:
            PreparedSpecificationSet(
                fem::FEM &,
                const mapping::DomainMapper &,
                PhysicalCore &,
                const Model &
            ) noexcept {
            }

            template <typename StateView>
            void ReadPhysicalControls(
                const StateView &,
                StellarEquilibriumControlContext &
            ) noexcept {
            }

            template <typename StateView>
            [[nodiscard]] StellarEquilibriumPreparationResult<void> TryReadPhysicalControls(
                const StateView &,
                StellarEquilibriumControlContext &
            ) noexcept {
                return {};
            }

            template <
                std::size_t Index,
                typename StateView,
                typename Reports>
            void PrepareAfterPhysical(
                const StateView &,
                const StellarEquilibriumDependencies &,
                const PhysicalCore &,
                Reports &
            ) noexcept {
            }

            template <
                std::size_t Index,
                typename StateView,
                typename Reports>
            [[nodiscard]] StellarEquilibriumPreparationResult<void> TryPrepareAfterPhysical(
                const StateView &,
                const StellarEquilibriumDependencies &,
                const PhysicalCore &,
                Reports &
            ) noexcept {
                return {};
            }

            template <typename ResidualView>
            void AddResidual(
                const ResidualView &,
                const PhysicalCore &
            ) const noexcept {
            }

            template <
                typename DirectionView,
                typename ActionView>
            void AddJacobianAction(
                const DirectionView &,
                const ActionView &,
                const PhysicalCore &
            ) const noexcept {
            }

            [[nodiscard]] constexpr bool IsPrepared() const noexcept {
                return true;
            }
        };

        template <
            model::StellarModelType Model,
            PreparedStellarEquilibriumPhysicalCore PhysicalCore,
            models::ModelSpecification Head,
            models::ModelSpecification... Tail>
        class PreparedSpecificationSet<Model, PhysicalCore, models::detail::SpecificationSetStorage<Head, Tail...>>
            final {
        private:
            using HeadSlot   = PreparedRuntimeContribution<Head, Model>;
            using HeadAccess = SpecificationRuntimeAccess<Head, Model>;
            using TailSlots =
                PreparedSpecificationSet<Model, PhysicalCore, models::detail::SpecificationSetStorage<Tail...>>;

        public:
            PreparedSpecificationSet(
                fem::FEM &finiteElements,
                const mapping::DomainMapper &domainMapper,
                PhysicalCore &physical,
                const Model &model
            )
                : m_head(
                      finiteElements,
                      domainMapper,
                      physical,
                      model
                  ),
                  m_tail(
                      finiteElements,
                      domainMapper,
                      physical,
                      model
                  ) {
            }

            template <typename StateView>
            void ReadPhysicalControls(
                const StateView &state,
                StellarEquilibriumControlContext &context
            ) {
                m_head.ReadPhysicalControls(typename HeadAccess::StateView{state}, context);
                m_tail.ReadPhysicalControls(state, context);
            }

            template <typename StateView>
            [[nodiscard]] StellarEquilibriumPreparationResult<void> TryReadPhysicalControls(
                const StateView &state,
                StellarEquilibriumControlContext &context
            ) {
                using HeadResult = StellarEquilibriumPreparationResult<void>;
                if constexpr (requires {
                                  {
                                      m_head.TryReadPhysicalControls(typename HeadAccess::StateView{state}, context)
                                  } -> std::same_as<HeadResult>;
                              }) {
                    auto headResult = m_head.TryReadPhysicalControls(typename HeadAccess::StateView{state}, context);
                    if (!headResult.has_value()) {
                        return std::unexpected(markSpecificationRejection(headResult.error()));
                    }
                } else {
                    m_head.ReadPhysicalControls(typename HeadAccess::StateView{state}, context);
                }
                return m_tail.TryReadPhysicalControls(state, context);
            }

            template <
                std::size_t Index,
                typename StateView,
                typename Reports>
            void PrepareAfterPhysical(
                const StateView &state,
                const StellarEquilibriumDependencies &dependencies,
                const PhysicalCore &physical,
                Reports &reports
            ) {
                std::get<Index>(reports) =
                    m_head.PrepareAfterPhysical(typename HeadAccess::StateView{state}, dependencies, physical);
                m_tail.template PrepareAfterPhysical<Index + 1>(state, dependencies, physical, reports);
            }

            template <
                std::size_t Index,
                typename StateView,
                typename Reports>
            [[nodiscard]] StellarEquilibriumPreparationResult<void> TryPrepareAfterPhysical(
                const StateView &state,
                const StellarEquilibriumDependencies &dependencies,
                const PhysicalCore &physical,
                Reports &reports
            ) {
                using HeadResult = StellarEquilibriumPreparationResult<typename HeadSlot::Report>;
                if constexpr (requires {
                                  {
                                      m_head.TryPrepareAfterPhysical(
                                          typename HeadAccess::StateView{state}, dependencies, physical
                                      )
                                  } -> std::same_as<HeadResult>;
                              }) {
                    auto headResult =
                        m_head.TryPrepareAfterPhysical(typename HeadAccess::StateView{state}, dependencies, physical);
                    if (!headResult.has_value()) {
                        return std::unexpected(markSpecificationRejection(headResult.error()));
                    }
                    std::get<Index>(reports) = std::move(headResult).value();
                } else {
                    std::get<Index>(reports) =
                        m_head.PrepareAfterPhysical(typename HeadAccess::StateView{state}, dependencies, physical);
                }
                return m_tail.template TryPrepareAfterPhysical<Index + 1>(state, dependencies, physical, reports);
            }

            template <typename ResidualView>
            void AddResidual(
                const ResidualView &residual,
                const PhysicalCore &physical
            ) const {
                m_head.AddResidual(
                    typename HeadAccess::ResidualView{residual, physical.GetSurfaceConstraintOperator()}
                );
                m_tail.AddResidual(residual, physical);
            }

            template <
                typename DirectionView,
                typename ActionView>
            void AddJacobianAction(
                const DirectionView &direction,
                const ActionView &action,
                const PhysicalCore &physical
            ) const {
                m_head.AddJacobianAction(
                    typename HeadAccess::StateView{direction},
                    typename HeadAccess::ResidualView{action, physical.GetSurfaceConstraintOperator()}, physical
                );
                m_tail.AddJacobianAction(direction, action, physical);
            }

            [[nodiscard]] bool IsPrepared() const noexcept {
                return m_head.IsPrepared() && m_tail.IsPrepared();
            }

            template <models::ModelSpecification Specification> [[nodiscard]] auto &Get() noexcept {
                if constexpr (std::same_as<Specification, Head>) {
                    if constexpr (HasNestedStellarEquilibriumPhysics<Head>) {
                        return m_head.physics();
                    } else {
                        return m_head;
                    }
                } else {
                    return m_tail.template Get<Specification>();
                }
            }

            template <models::ModelSpecification Specification> [[nodiscard]] const auto &Get() const noexcept {
                if constexpr (std::same_as<Specification, Head>) {
                    if constexpr (HasNestedStellarEquilibriumPhysics<Head>) {
                        return m_head.physics();
                    } else {
                        return m_head;
                    }
                } else {
                    return m_tail.template Get<Specification>();
                }
            }

        private:
            HeadSlot m_head;
            TailSlots m_tail;
        };

        template <
            models::ModelSpecification Specification,
            model::StellarModelType Model,
            bool SymbolicallyCompilable = StellarEquilibriumSystemCompilable<Model> &&
                                          CoreRuntimeInterfaceAudit<std::remove_cvref_t<Model>>::complete,
            typename                    = void>
        struct RuntimeContributionInterfaceAudit {
            static constexpr bool complete                 = false;
            static constexpr std::size_t rotationProviders = 0;
        };

        template <models::ModelSpecification Specification, model::StellarModelType Model>
        struct RuntimeContributionInterfaceAudit<
            Specification,
            Model,
            true,
            std::void_t<
                PreparedRuntimeContribution<Specification, Model>,
                typename PreparedRuntimeContribution<Specification, Model>::Report,
                std::bool_constant<static_cast<bool>(RuntimeContributionSelection<Specification, Model>::registered)>,
                std::integral_constant<
                    std::size_t,
                    static_cast<std::size_t>(RuntimeContributionSelection<Specification, Model>::rotationProviders)>>> {
        private:
            using Selection        = RuntimeContributionSelection<Specification, Model>;
            using Prepared         = PreparedRuntimeContribution<Specification, Model>;
            using PhysicalCore     = typename CoreRuntimeInterfaceAudit<Model>::CoreType;
            using Access           = SpecificationRuntimeAccess<Specification, Model>;
            using StateView        = typename Access::StateView;
            using ResidualViewType = typename Access::ResidualView;

        public:
            using Report                   = typename Prepared::Report;

            static constexpr bool complete = requires(
                fem::FEM &finiteElements,
                const mapping::DomainMapper &domainMapper,
                PhysicalCore &physical,
                const PhysicalCore &constantPhysical,
                const Model &model,
                Prepared &prepared,
                const Prepared &constantPrepared,
                const StateView &state,
                const ResidualViewType &residual,
                StellarEquilibriumControlContext &controls,
                const StellarEquilibriumDependencies &dependencies
            ) {
                requires Selection::registered;
                requires !Selection::ambiguous;
                requires std::default_initializable<Report>;
                requires std::assignable_from<Report &, Report>;
                requires std::constructible_from<
                    Prepared, fem::FEM &, const mapping::DomainMapper &, PhysicalCore &, const Model &>;
                { prepared.ReadPhysicalControls(state, controls) } -> std::same_as<void>;
                { prepared.PrepareAfterPhysical(state, dependencies, constantPhysical) } -> std::same_as<Report>;
                { constantPrepared.AddResidual(residual) } -> std::same_as<void>;
                { constantPrepared.AddJacobianAction(state, residual, constantPhysical) } -> std::same_as<void>;
                { constantPrepared.IsPrepared() } -> std::convertible_to<bool>;
            };
            static constexpr std::size_t rotationProviders = complete ? Selection::rotationProviders : 0;
        };

        template <model::StellarModelType Model, typename SpecificationSet, bool Complete>
        struct RuntimeContributionAuditImpl;

        template <model::StellarModelType Model, models::ModelSpecification... Specifications>
        struct RuntimeContributionAuditImpl<Model, models::detail::SpecificationSetStorage<Specifications...>, false> {
            static constexpr bool complete                 = false;
            static constexpr std::size_t rotationProviders = 0;
            using ReportTuple                              = std::tuple<>;
        };

        template <model::StellarModelType Model, models::ModelSpecification... Specifications>
        struct RuntimeContributionAuditImpl<Model, models::detail::SpecificationSetStorage<Specifications...>, true> {
            static constexpr bool complete = true;
            static constexpr std::size_t rotationProviders =
                (std::size_t{0} + ... + RuntimeContributionInterfaceAudit<Specifications, Model>::rotationProviders);
            using ReportTuple =
                std::tuple<typename RuntimeContributionInterfaceAudit<Specifications, Model>::Report...>;
        };

        template <model::StellarModelType Model, typename SpecificationSet> struct RuntimeContributionAudit;

        template <model::StellarModelType Model, models::ModelSpecification... Specifications>
        struct RuntimeContributionAudit<Model, models::detail::SpecificationSetStorage<Specifications...>>
            : RuntimeContributionAuditImpl<
                  Model,
                  models::detail::SpecificationSetStorage<Specifications...>,
                  (RuntimeContributionInterfaceAudit<Specifications, Model>::complete && ...)> { };

        template <typename List> struct MakeValueSizes;

        template <typename... Blocks> struct MakeValueSizes<utils::blocks::type_list<Blocks...>> {
            template <typename PhysicalForm>
            [[nodiscard]] static std::array<
                int,
                sizeof...(Blocks)>
            Apply(const utils::blocks::form_layout<PhysicalForm> &physicalLayout) {
                return {BlockSize<Blocks>(physicalLayout)...};
            }

        private:
            template <
                typename Block,
                typename PhysicalForm>
            [[nodiscard]] static int BlockSize(const utils::blocks::form_layout<PhysicalForm> &physicalLayout) {
                if constexpr (utils::blocks::contains_type_v<Block, typename PhysicalForm::value_blocks>) {
                    constexpr int index = utils::blocks::type_index_v<Block, typename PhysicalForm::value_blocks>;
                    return physicalLayout.value_offsets()[index + 1] - physicalLayout.value_offsets()[index];
                } else {
                    static_assert(
                        Block::static_block_size != utils::blocks::dynamic_block_size,
                        "A generated stellar-equilibrium value block must have a compile-time size."
                    );
                    return Block::static_block_size;
                }
            }
        };

        template <typename List> struct MakeResidualSizes;

        template <typename... Blocks> struct MakeResidualSizes<utils::blocks::type_list<Blocks...>> {
            template <typename PhysicalForm>
            [[nodiscard]] static std::array<
                int,
                sizeof...(Blocks)>
            Apply(const utils::blocks::form_layout<PhysicalForm> &physicalLayout) {
                return {BlockSize<Blocks>(physicalLayout)...};
            }

        private:
            template <
                typename Block,
                typename PhysicalForm>
            [[nodiscard]] static int BlockSize(const utils::blocks::form_layout<PhysicalForm> &physicalLayout) {
                if constexpr (utils::blocks::contains_type_v<Block, typename PhysicalForm::residual_blocks>) {
                    constexpr int index = utils::blocks::type_index_v<Block, typename PhysicalForm::residual_blocks>;
                    return physicalLayout.residual_offsets()[index + 1] - physicalLayout.residual_offsets()[index];
                } else {
                    static_assert(
                        Block::static_block_size != utils::blocks::dynamic_block_size,
                        "A generated stellar-equilibrium residual block must have a compile-time size."
                    );
                    return Block::static_block_size;
                }
            }
        };

        template <typename RequiredBlocks, typename AvailableBlocks> struct AllBlocksBelongToList : std::false_type { };

        template <typename... RequiredBlocks, typename AvailableBlocks>
        struct AllBlocksBelongToList<utils::blocks::type_list<RequiredBlocks...>, AvailableBlocks>
            : std::bool_constant<(utils::blocks::contains_type_v<RequiredBlocks, AvailableBlocks> && ...)> { };

        template <
            model::StellarModelType Model,
            bool HasCompiledPhysicalRoot =
                StellarEquilibriumSystemCompilable<Model> &&
                std::remove_cvref_t<Model>::template containsSpecification<models::FixedTotalMass>>
        struct PhysicalRootCompatibilityAudit {
            static constexpr bool complete = false;
        };

        template <model::StellarModelType Model> struct PhysicalRootCompatibilityAudit<Model, true> {
        private:
            using ModelType    = std::remove_cvref_t<Model>;
            using RootForm     = CompiledStellarEquilibriumForm<ModelType>;
            using PhysicalForm = utils::blocks::surface_deformed_stellar_equilibrium_form;

        public:
            static constexpr bool complete =
                AllBlocksBelongToList<typename PhysicalForm::value_blocks, typename RootForm::value_blocks>::value &&
                AllBlocksBelongToList<typename PhysicalForm::residual_blocks, typename RootForm::residual_blocks>::
                    value;
        };

        template <typename Specification, typename SpecificationSet> struct SpecificationIndex;

        template <typename Specification, models::ModelSpecification... Tail>
        struct SpecificationIndex<Specification, models::detail::SpecificationSetStorage<Specification, Tail...>>
            : std::integral_constant<std::size_t, 0> { };

        template <typename Specification, models::ModelSpecification Head, models::ModelSpecification... Tail>
        struct SpecificationIndex<Specification, models::detail::SpecificationSetStorage<Head, Tail...>>
            : std::integral_constant<
                  std::size_t,
                  1 + SpecificationIndex<Specification, models::detail::SpecificationSetStorage<Tail...>>::value> { };

        template <typename Specification>
        struct SpecificationIndex<Specification, models::detail::SpecificationSetStorage<>>;
    } // namespace detail

    /*
     * Physics-extension views expose exactly the blocks declared by one
     * specification.  Generated coordinates/rows are included automatically;
     * no extension author needs to spell out backend block lists twice.
     */
    template <typename Specification, typename Model>
    concept StellarEquilibriumSpecificationBelongsToModel =
        models::ModelSpecification<std::remove_cvref_t<Specification>> &&
        model::StellarModelType<std::remove_cvref_t<Model>> && requires {
            requires std::remove_cvref_t<Model>::template containsSpecification<std::remove_cvref_t<Specification>>;
        };

    template <typename Specification, typename Model>
        requires StellarEquilibriumSpecificationBelongsToModel<Specification, Model>
    using StellarEquilibriumContributionStateView = typename detail::
        SpecificationRuntimeAccess<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>::StateView;

    template <typename Specification, typename Model>
        requires StellarEquilibriumSpecificationBelongsToModel<Specification, Model>
    using StellarEquilibriumContributionResidualView = typename detail::
        SpecificationRuntimeAccess<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>::ResidualView;

    /* Low-level topology-inspection view retained for compiler tests and
     * backend adapters.  Physics-facing nested runtimes are not handed this
     * imperative object; they use the exact provider protocol below. */
    template <typename Specification, typename Model>
        requires StellarEquilibriumSpecificationBelongsToModel<Specification, Model>
    using StellarEquilibriumContributionJacobianView = typename detail::
        SpecificationRuntimeAccess<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>::JacobianView;

    template <typename Specification, typename Model>
        requires StellarEquilibriumSpecificationBelongsToModel<Specification, Model>
    using StellarEquilibriumContributionTopology = typename detail::
        SpecificationRuntimeAccess<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>::Topology;

    template <typename Specification, typename Model, typename StateTag>
        requires StellarEquilibriumSpecificationBelongsToModel<Specification, Model>
    using StellarEquilibriumContributionDirection =
        typename detail::SpecificationRuntimeAccess<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>::
            template DirectionView<StateTag>;

    template <typename Specification, typename Model, typename EquationTag>
        requires StellarEquilibriumSpecificationBelongsToModel<Specification, Model>
    using StellarEquilibriumContributionRow =
        typename detail::SpecificationRuntimeAccess<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>::
            template Row<EquationTag>;

    /* A focused diagnostic concept for extension authors.  It answers the
     * useful question directly: does this prepared physics class implement
     * every residual row and derivative inferred from my declaration? */
    template <typename Physics, typename Specification, typename Model>
    concept CompleteStellarEquilibriumPhysicsProvider =
        StellarEquilibriumSpecificationBelongsToModel<Specification, Model> &&
        detail::exactSpecificationPhysicsProvidersComplete<
            std::remove_cvref_t<Physics>,
            detail::SpecificationRuntimeAccess<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>>;

    /* True only when the exact nested physics provider corresponding to this
     * compiled coupling returns StructuralZero.  This is intentionally a
     * proof about the provider's type, not a second author-written metadata
     * flag.  Preconditioners use it to distinguish an intentional zero from a
     * nonzero contribution that their selected structure backend must either
     * implement or reject. */
    template <typename Specification, typename Model, typename Coupling>
    inline constexpr bool stellarEquilibriumSpecificationCouplingIsStructuralZero =
        StellarEquilibriumSpecificationBelongsToModel<Specification, Model> &&
        detail::NestedSpecificationCouplingIsStructuralZero<
            std::remove_cvref_t<Specification>,
            std::remove_cvref_t<Model>,
            std::remove_cvref_t<Coupling>>::value;

    template <typename Specification, typename Model>
    concept StellarEquilibriumPhysicsAvailableFor =
        StellarEquilibriumSpecificationBelongsToModel<Specification, Model> &&
        detail::RuntimeContributionInterfaceAudit<std::remove_cvref_t<Specification>, std::remove_cvref_t<Model>>::
            complete;

    template <model::StellarModelType Model>
    inline constexpr bool hasCompleteStellarEquilibriumRuntime = detail::RuntimeContributionAudit<
        std::remove_cvref_t<Model>,
        typename std::remove_cvref_t<Model>::SpecificationTypes>::complete;

    template <model::StellarModelType Model>
    inline constexpr bool hasStellarEquilibriumCoreRuntime =
        detail::CoreRuntimeInterfaceAudit<std::remove_cvref_t<Model>>::complete;

    template <model::StellarModelType Model>
        requires hasStellarEquilibriumCoreRuntime<Model>
    using StellarEquilibriumPhysicalCoreType =
        typename detail::CoreRuntimeInterfaceAudit<std::remove_cvref_t<Model>>::CoreType;

    template <model::StellarModelType Model>
    inline constexpr std::size_t stellarEquilibriumRotationProviderCount = detail::RuntimeContributionAudit<
        std::remove_cvref_t<Model>,
        typename std::remove_cvref_t<Model>::SpecificationTypes>::rotationProviders;

    template <model::StellarModelType Model>
    inline constexpr bool hasCompatibleStellarEquilibriumPhysicalRoot =
        detail::PhysicalRootCompatibilityAudit<std::remove_cvref_t<Model>>::complete;

    template <model::StellarModelType Model>
        requires hasCompleteStellarEquilibriumRuntime<Model>
    struct PreparedVariadicStellarEquilibriumReport final {
        using ModelType          = std::remove_cvref_t<Model>;
        using SpecificationTypes = typename ModelType::SpecificationTypes;
        using SpecificationReports =
            typename detail::RuntimeContributionAudit<ModelType, SpecificationTypes>::ReportTuple;

        PreparedStellarEquilibriumReport physical;
        SpecificationReports specifications;
        bool generatedPhysicalControl{false};
        bool assembledResidual{false};

        template <models::ModelSpecification Specification>
            requires ModelType::template
        containsSpecification<Specification> [[nodiscard]] const auto &specification() const noexcept {
            constexpr std::size_t index = detail::SpecificationIndex<Specification, SpecificationTypes>::value;
            return std::get<index>(specifications);
        }

        [[nodiscard]] bool DidAnyWork() const noexcept {
            return physical.DidAnyWork() || generatedPhysicalControl || assembledResidual;
        }
    };

    /**
     * One runtime root for every fully supported specification pack.
     *
     * The class template itself is the inferred type.  Its slot set is a
     * recursive, statically dispatched fold, so adding a specification never
     * creates a new hand-written combination class or a runtime registry.
     */
    template <model::StellarModelType Model>
        requires hasCompatibleStellarEquilibriumPhysicalRoot<Model> && hasCompleteStellarEquilibriumRuntime<Model> &&
                 hasStellarEquilibriumCoreRuntime<Model> &&
                 CompilableRootManifestFor<
                     std::remove_cvref_t<Model>,
                     CompiledStellarEquilibriumForm<std::remove_cvref_t<Model>>>
    class PreparedVariadicStellarEquilibriumOperator final : public mfem::Operator {
    public:
        using ModelType            = std::remove_cvref_t<Model>;
        using EquationOfStateType  = model::EquationOfStateType<ModelType>;
        using SurfaceConditionType = model::SurfaceConditionType<ModelType>;
        using CoreRuntime          = StellarEquilibriumCoreRuntime<EquationOfStateType>;
        using PhysicalCoreType     = StellarEquilibriumPhysicalCoreType<ModelType>;
        using PhysicalCoreOwner    = std::unique_ptr<PhysicalCoreType>;
        using SpecificationTypes   = typename ModelType::SpecificationTypes;
        using FormType             = CompiledStellarEquilibriumForm<ModelType>;
        using JacobianFormType     = CompiledStellarEquilibriumJacobianForm<ModelType>;
        using Layout               = utils::blocks::form_layout<FormType>;
        using Manifest             = EquilibriumSystemManifest<ModelType, FormType, JacobianFormType>;
        using Report               = PreparedVariadicStellarEquilibriumReport<ModelType>;
        using PreparationResult    = StellarEquilibriumPreparationResult<Report>;

        static constexpr std::size_t rotationProviderCount = stellarEquilibriumRotationProviderCount<ModelType>;

        /*
         * The physical core may retain references to constitutive data, so a
         * raw ModelType reference is intentionally not a construction option.
         * Shared ownership keeps every EOS backend safe for the lifetime of
         * this prepared root.  StellarEquilibriumProblem owns and supplies the
         * same handle on the normal user-facing path.
         */
        PreparedVariadicStellarEquilibriumOperator(
            fem::FEM &finiteElements,
            const mapping::DomainMapper &domainMapper,
            std::shared_ptr<const ModelType> model,
            PressureSurfaceConstraintView surfaceConstraint,
            deformation::PreparedDomainDeformationRuntime domainDeformation
        )
            requires(rotationProviderCount <= 1)
            : PreparedVariadicStellarEquilibriumOperator(
                  finiteElements,
                  domainMapper,
                  model,
                  MakePhysical(
                      finiteElements,
                      domainMapper,
                      RequireModel(model),
                      surfaceConstraint,
                      std::move(domainDeformation)
                  )
              ) {
        }

        PreparedVariadicStellarEquilibriumOperator(const PreparedVariadicStellarEquilibriumOperator &) = delete;
        PreparedVariadicStellarEquilibriumOperator &
        operator=(const PreparedVariadicStellarEquilibriumOperator &)                                        = delete;
        PreparedVariadicStellarEquilibriumOperator(PreparedVariadicStellarEquilibriumOperator &&)            = delete;
        PreparedVariadicStellarEquilibriumOperator &operator=(PreparedVariadicStellarEquilibriumOperator &&) = delete;

        [[nodiscard]] Report Prepare(
            const mfem::Vector &state,
            const StellarEquilibriumDependencies &dependencies,
            const physics::RigidRotation &rotation
        )
            requires(rotationProviderCount == 0)
        {
            detail::StellarEquilibriumControlContext controls{
                .dependencies             = dependencies,
                .rotation                 = rotation,
                .rotationProviderCount    = 1,
                .generatedPhysicalControl = false
            };
            return PrepareWithControls(state, std::move(controls));
        }

        [[nodiscard]] PreparationResult TryPrepare(
            const mfem::Vector &state,
            const StellarEquilibriumDependencies &dependencies,
            const physics::RigidRotation &rotation
        )
            requires(rotationProviderCount == 0 && FalliblePreparedStellarEquilibriumPhysicalCore<PhysicalCoreType>)
        {
            detail::StellarEquilibriumControlContext controls{
                .dependencies             = dependencies,
                .rotation                 = rotation,
                .rotationProviderCount    = 1,
                .generatedPhysicalControl = false
            };
            return TryPrepareWithControls(state, std::move(controls));
        }

        [[nodiscard]] Report Prepare(
            const mfem::Vector &state,
            const StellarEquilibriumDependencies &dependencies
        )
            requires(rotationProviderCount == 1)
        {
            detail::StellarEquilibriumControlContext controls{
                .dependencies             = dependencies,
                .rotation                 = std::nullopt,
                .rotationProviderCount    = 0,
                .generatedPhysicalControl = false
            };
            return PrepareWithControls(state, std::move(controls));
        }

        [[nodiscard]] PreparationResult TryPrepare(
            const mfem::Vector &state,
            const StellarEquilibriumDependencies &dependencies
        )
            requires(rotationProviderCount == 1 && FalliblePreparedStellarEquilibriumPhysicalCore<PhysicalCoreType>)
        {
            detail::StellarEquilibriumControlContext controls{
                .dependencies             = dependencies,
                .rotation                 = std::nullopt,
                .rotationProviderCount    = 0,
                .generatedPhysicalControl = false
            };
            return TryPrepareWithControls(state, std::move(controls));
        }

        void BuildResidual(mfem::Vector &residual) const {
            VerifyPrepared();
            residual = m_cachedResidual;
        }

        void Mult(
            const mfem::Vector &direction,
            mfem::Vector &action
        ) const override {
            VerifyPrepared();
            MFEM_VERIFY(direction.Size() == Width(), "The variadic stellar root received a wrong-sized direction.");

            GatherPhysicalValues(direction, m_physicalDirection, typename PhysicalForm::value_blocks{});
            m_physical->Mult(m_physicalDirection, m_physicalAction);

            action.SetSize(Height());
            action = 0.0;
            ScatterPhysicalResiduals(m_physicalAction, action, typename PhysicalForm::residual_blocks{});

            const auto directionView = m_manifest.directionView(direction);
            const auto actionView    = m_manifest.residualView(action);
            m_specifications.AddJacobianAction(directionView, actionView, *m_physical);
        }

        [[nodiscard]] bool IsPrepared() const noexcept {
            return m_isPrepared && m_physical->IsPrepared() && m_specifications.IsPrepared();
        }

        [[nodiscard]] const Layout &GetLayout() const noexcept {
            return m_manifest.layout();
        }

        [[nodiscard]] const Manifest &GetRootManifest() const noexcept {
            return m_manifest;
        }

        [[nodiscard]] physics::RigidRotation GetRotation() const {
            if (!IsPrepared() || !m_activeRotation.has_value()) {
                throw std::logic_error("The prepared stellar-equilibrium root has no active rotation.");
            }
            return *m_activeRotation;
        }

        [[nodiscard]] const PhysicalCoreType &GetPhysicalOperator() const noexcept {
            return *m_physical;
        }

        void BuildVolumeDisplacementDirection(
            const mfem::Vector &stateDirection,
            mfem::Vector &volumeDisplacementDirection
        ) const {
            if (stateDirection.Size() != Width()) {
                throw std::invalid_argument(
                    "The prepared stellar-equilibrium root received a state direction with the wrong size."
                );
            }
            const auto rootDirection = m_manifest.directionView(stateDirection);
            m_physical->BuildVolumeDisplacementDirection(
                rootDirection.block(utils::blocks::surface_deformation_field.parameters_term),
                volumeDisplacementDirection
            );
        }

        template <models::ModelSpecification Specification>
            requires ModelType::template
        containsSpecification<Specification> [[nodiscard]] const auto &GetPreparedContribution() const noexcept {
            return m_specifications.template Get<Specification>();
        }

        [[nodiscard]] const PreparedAngularMomentumOperator &GetAngularMomentumConstraint() const noexcept
            requires ModelType::template
        containsSpecification<models::FixedAngularMomentum> {
            return GetPreparedContribution<models::FixedAngularMomentum>().constraint();
        }

        [[nodiscard]] const PreparedCentralDensityConstraint &GetCentralDensityConstraint() const noexcept
            requires ModelType::template
        containsSpecification<models::FixedCentralDensity> {
            return GetPreparedContribution<models::FixedCentralDensity>().constraint();
        }

        [[nodiscard]] RootConstraintReport GetFixedMassReport() const {
            VerifyPrepared();
            const RootConstraintReport physicalReport = m_physical->GetFixedMassReport();
            return m_manifest.fixedMassReport(physicalReport.achieved);
        }

        [[nodiscard]] AngularMomentumConstraintReport GetAngularMomentumReport() const
            requires ModelType::template
        containsSpecification<models::FixedAngularMomentum> {
            VerifyPrepared();
            return GetAngularMomentumConstraint().GetConstraintReport();
        }

        [[nodiscard]] CentralDensityConstraintReport GetCentralDensityReport() const
            requires ModelType::template
        containsSpecification<models::FixedCentralDensity> {
            VerifyPrepared();
            return GetCentralDensityConstraint().GetConstraintReport();
        }

    private:
        using PhysicalForm       = utils::blocks::surface_deformed_stellar_equilibrium_form;
        using SpecificationSlots = detail::PreparedSpecificationSet<ModelType, PhysicalCoreType, SpecificationTypes>;

        [[nodiscard]] static const ModelType &RequireModel(const std::shared_ptr<const ModelType> &model) {
            MFEM_VERIFY(
                model != nullptr, "The variadic stellar-equilibrium root requires shared ownership of its model."
            );
            return *model;
        }

        template <typename Block>
        void GatherPhysicalValueBlock(
            const mfem::Vector &root,
            mfem::Vector &physical
        ) const {
            constexpr int rootIndex     = utils::blocks::type_index_v<Block, typename FormType::value_blocks>;
            constexpr int physicalIndex = utils::blocks::type_index_v<Block, typename PhysicalForm::value_blocks>;
            const auto &rootOffsets     = m_manifest.layout().value_offsets();
            const auto &physicalOffsets = m_physical->GetLayout().value_offsets();
            const int rootSize          = rootOffsets[rootIndex + 1] - rootOffsets[rootIndex];
            const int physicalSize      = physicalOffsets[physicalIndex + 1] - physicalOffsets[physicalIndex];
            MFEM_VERIFY(rootSize == physicalSize, "A compiled physical value block changed size in the root layout.");
            const mfem::Vector source(const_cast<mfem::real_t *>(root.GetData()) + rootOffsets[rootIndex], rootSize);
            mfem::Vector destination(physical, physicalOffsets[physicalIndex], physicalSize);
            destination = source;
            destination.SyncAliasMemory(physical);
        }

        template <typename... Blocks>
        void GatherPhysicalValues(
            const mfem::Vector &root,
            mfem::Vector &physical,
            utils::blocks::type_list<Blocks...>
        ) const {
            MFEM_VERIFY(physical.Size() == m_physical->Width(), "The physical-state workspace has the wrong size.");
            (GatherPhysicalValueBlock<Blocks>(root, physical), ...);
        }

        template <typename Block>
        void ScatterPhysicalResidualBlock(
            const mfem::Vector &physical,
            mfem::Vector &root
        ) const {
            constexpr int physicalIndex = utils::blocks::type_index_v<Block, typename PhysicalForm::residual_blocks>;
            constexpr int rootIndex     = utils::blocks::type_index_v<Block, typename FormType::residual_blocks>;
            const auto &physicalOffsets = m_physical->GetLayout().residual_offsets();
            const auto &rootOffsets     = m_manifest.layout().residual_offsets();
            const int physicalSize      = physicalOffsets[physicalIndex + 1] - physicalOffsets[physicalIndex];
            const int rootSize          = rootOffsets[rootIndex + 1] - rootOffsets[rootIndex];
            MFEM_VERIFY(
                rootSize == physicalSize, "A compiled physical residual block changed size in the root layout."
            );
            const mfem::Vector source(
                const_cast<mfem::real_t *>(physical.GetData()) + physicalOffsets[physicalIndex], physicalSize
            );
            mfem::Vector destination(root, rootOffsets[rootIndex], rootSize);
            destination = source;
            destination.SyncAliasMemory(root);
        }

        template <typename... Blocks>
        void ScatterPhysicalResiduals(
            const mfem::Vector &physical,
            mfem::Vector &root,
            utils::blocks::type_list<Blocks...>
        ) const {
            MFEM_VERIFY(physical.Size() == m_physical->Height(), "The physical-action workspace has the wrong size.");
            (ScatterPhysicalResidualBlock<Blocks>(physical, root), ...);
        }

        [[nodiscard]] static PhysicalCoreOwner MakePhysical(
            fem::FEM &finiteElements,
            const mapping::DomainMapper &domainMapper,
            const ModelType &model,
            PressureSurfaceConstraintView surfaceConstraint,
            deformation::PreparedDomainDeformationRuntime domainDeformation
        ) {
            return CoreRuntime::Make(
                finiteElements, domainMapper, model.equationOfState(),
                models::compileConstraint(model.template specification<models::FixedTotalMass>()), surfaceConstraint,
                std::move(domainDeformation)
            );
        }

        [[nodiscard]] static std::array<
            int,
            FormType::value_block_count>
        MakeValueSizes(const StellarEquilibriumLayout &physicalLayout) {
            return detail::MakeValueSizes<typename FormType::value_blocks>::Apply(physicalLayout);
        }

        [[nodiscard]] static std::array<
            int,
            FormType::residual_block_count>
        MakeResidualSizes(const StellarEquilibriumLayout &physicalLayout) {
            return detail::MakeResidualSizes<typename FormType::residual_blocks>::Apply(physicalLayout);
        }

        PreparedVariadicStellarEquilibriumOperator(
            fem::FEM &finiteElements,
            const mapping::DomainMapper &domainMapper,
            std::shared_ptr<const ModelType> model,
            PhysicalCoreOwner physical
        )
            : mfem::Operator(
                  Layout(
                      MakeValueSizes(physical->GetLayout()),
                      MakeResidualSizes(physical->GetLayout())
                  )
                      .residual_offsets()
                      .Last(),
                  Layout(
                      MakeValueSizes(physical->GetLayout()),
                      MakeResidualSizes(physical->GetLayout())
                  )
                      .value_offsets()
                      .Last()
              ),
              m_model(std::move(model)),
              m_physical(std::move(physical)),
              m_specifications(
                  finiteElements,
                  domainMapper,
                  *m_physical,
                  *m_model
              ),
              m_manifest(
                  MakeValueSizes(m_physical->GetLayout()),
                  MakeResidualSizes(m_physical->GetLayout()),
                  *m_model,
                  CoreRuntime::SurfaceEquationCount(*m_physical)
              ),
              m_physicalState(m_physical->Width()),
              m_physicalDirection(m_physical->Width()),
              m_physicalAction(m_physical->Height()) {
            static_assert(rotationProviderCount <= 1, "A stellar root cannot have two rigid-rotation providers.");
            MFEM_VERIFY(
                Width() == m_manifest.layout().value_offsets().Last() &&
                    Height() == m_manifest.layout().residual_offsets().Last(),
                "The variadic stellar root has inconsistent compiled dimensions."
            );
        }

        [[nodiscard]] Report PrepareWithControls(
            const mfem::Vector &state,
            detail::StellarEquilibriumControlContext controls
        ) {
            MFEM_VERIFY(state.Size() == Width(), "The variadic stellar root received a wrong-sized state.");
            const auto stateView = m_manifest.stateView(state);

            m_isPrepared         = false;
            m_specifications.ReadPhysicalControls(stateView, controls);
            MFEM_VERIFY(
                controls.rotationProviderCount == 1 && controls.rotation.has_value(),
                "Exactly one rigid-rotation value must be supplied to the stellar physics core."
            );

            GatherPhysicalValues(state, m_physicalState, typename PhysicalForm::value_blocks{});
            Report report;
            report.generatedPhysicalControl = controls.generatedPhysicalControl;
            report.physical = m_physical->Prepare(m_physicalState, controls.dependencies, *controls.rotation);
            m_specifications.template PrepareAfterPhysical<0>(
                stateView, controls.dependencies, *m_physical, report.specifications
            );
            AssembleResidual();
            m_activeRotation         = *controls.rotation;
            report.assembledResidual = true;
            m_isPrepared             = true;
            return report;
        }

        [[nodiscard]] PreparationResult TryPrepareWithControls(
            const mfem::Vector &state,
            detail::StellarEquilibriumControlContext controls
        )
            requires FalliblePreparedStellarEquilibriumPhysicalCore<PhysicalCoreType>
        {
            MFEM_VERIFY(state.Size() == Width(), "The variadic stellar root received a wrong-sized state.");
            const auto stateView = m_manifest.stateView(state);

            m_isPrepared         = false;
            auto controlResult   = m_specifications.TryReadPhysicalControls(stateView, controls);
            if (!controlResult.has_value()) {
                return std::unexpected(controlResult.error());
            }
            MFEM_VERIFY(
                controls.rotationProviderCount == 1 && controls.rotation.has_value(),
                "Exactly one rigid-rotation value must be supplied to the stellar physics core."
            );

            GatherPhysicalValues(state, m_physicalState, typename PhysicalForm::value_blocks{});
            Report report;
            report.generatedPhysicalControl = controls.generatedPhysicalControl;
            auto physicalResult = m_physical->TryPrepare(m_physicalState, controls.dependencies, *controls.rotation);
            if (!physicalResult.has_value()) {
                return std::unexpected(physicalResult.error());
            }
            report.physical          = std::move(physicalResult).value();
            auto specificationResult = m_specifications.template TryPrepareAfterPhysical<0>(
                stateView, controls.dependencies, *m_physical, report.specifications
            );
            if (!specificationResult.has_value()) {
                return std::unexpected(specificationResult.error());
            }
            AssembleResidual();
            m_activeRotation         = *controls.rotation;
            report.assembledResidual = true;
            m_isPrepared             = true;
            return report;
        }

        void AssembleResidual() {
            mfem::Vector physicalResidual;
            m_physical->BuildResidual(physicalResidual);
            m_cachedResidual.SetSize(Height());
            m_cachedResidual = 0.0;
            ScatterPhysicalResiduals(physicalResidual, m_cachedResidual, typename PhysicalForm::residual_blocks{});
            const auto residualView = m_manifest.residualView(m_cachedResidual);
            m_specifications.AddResidual(residualView, *m_physical);
        }

        void VerifyPrepared() const {
            MFEM_VERIFY(IsPrepared(), "The variadic stellar-equilibrium root must be prepared before application.");
        }

        std::shared_ptr<const ModelType> m_model;
        PhysicalCoreOwner m_physical;
        SpecificationSlots m_specifications;
        Manifest m_manifest;
        mfem::Vector m_physicalState;
        mutable mfem::Vector m_physicalDirection;
        mutable mfem::Vector m_physicalAction;
        mfem::Vector m_cachedResidual;
        std::optional<physics::RigidRotation> m_activeRotation;
        bool m_isPrepared{false};
    };
} // namespace mean_field::operators

export namespace mean_field::stellar {
    template <operators::DensityVolumeIntegralSpecification Specification>
    using DensityVolumeIntegralContext = operators::DensityVolumeIntegralContext<Specification>;
} // namespace mean_field::stellar