From 48c6bac1cc15be997d8b1b7e8524664100fc94b6 Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <thomas@boudreauxmail.com>
Date: Mon, 9 Nov 2026 13:20:46 -0500
Subject: [PATCH] feat(model): added model sub system

SERiF defines a so-called "model" subsystem. A model is a fully compile-time description of a problem. Note that this is a description not an evaluation tool. The basic idea is that since we have this robust description we can more easily serialize and be confident that others who load a model are getting the exact same thing. Further, this provides a convient place to do things like validate that all options are compatible, do some limited initial construction of the needed computational graph, allocate needed memory, etc... Some of the stuff that might initially seem to belong to a model will be in the discritization system though (e.g. FEM object lifetimes)
---
 .../serif/model/configured_discretization.hpp |  35 ++++
 src/include/serif/model/description.hpp       |   7 +
 .../model/detail/specification_endpoints.hpp  | 111 +++++++++++
 .../serif/model/detail/specification_list.hpp |  97 +++++++++
 .../serif/model/detail/surface_plan.hpp       |  80 ++++++++
 src/include/serif/model/eos_reference.hpp     |  96 +++++++++
 .../serif/model/equilibrium_equations.hpp     | 101 ++++++++++
 src/include/serif/model/equilibrium_plan.hpp  |  66 +++++++
 .../serif/model/execution_validation.hpp      | 144 ++++++++++++++
 src/include/serif/model/finite_elements.hpp   |  13 ++
 src/include/serif/model/keys.hpp              |   9 +
 src/include/serif/model/mesh_options.hpp      |  34 ++++
 src/include/serif/model/model.hpp             |   7 +
 .../serif/model/normalization/plan.hpp        | 131 ++++++++++++
 .../serif/model/normalization/reference.hpp   |  63 ++++++
 src/include/serif/model/plan.hpp              |  68 +++++++
 src/include/serif/model/specification.hpp     | 155 +++++++++++++++
 .../specifications/fixed_angular_momentum.hpp |  58 ++++++
 .../model/specifications/fixed_total_mass.hpp |  41 ++++
 .../model/specifications/physical_names.hpp   |  36 ++++
 src/include/serif/model/state.hpp             |  80 ++++++++
 src/include/serif/model/stellar_model.hpp     | 186 ++++++++++++++++++
 src/include/serif/model/validation.hpp        | 138 +++++++++++++
 src/lib/serif/model/specifications.cpp        |  60 ++++++
 24 files changed, 1816 insertions(+)
 create mode 100644 src/include/serif/model/configured_discretization.hpp
 create mode 100644 src/include/serif/model/description.hpp
 create mode 100644 src/include/serif/model/detail/specification_endpoints.hpp
 create mode 100644 src/include/serif/model/detail/specification_list.hpp
 create mode 100644 src/include/serif/model/detail/surface_plan.hpp
 create mode 100644 src/include/serif/model/eos_reference.hpp
 create mode 100644 src/include/serif/model/equilibrium_equations.hpp
 create mode 100644 src/include/serif/model/equilibrium_plan.hpp
 create mode 100644 src/include/serif/model/execution_validation.hpp
 create mode 100644 src/include/serif/model/finite_elements.hpp
 create mode 100644 src/include/serif/model/keys.hpp
 create mode 100644 src/include/serif/model/mesh_options.hpp
 create mode 100644 src/include/serif/model/model.hpp
 create mode 100644 src/include/serif/model/normalization/plan.hpp
 create mode 100644 src/include/serif/model/normalization/reference.hpp
 create mode 100644 src/include/serif/model/plan.hpp
 create mode 100644 src/include/serif/model/specification.hpp
 create mode 100644 src/include/serif/model/specifications/fixed_angular_momentum.hpp
 create mode 100644 src/include/serif/model/specifications/fixed_total_mass.hpp
 create mode 100644 src/include/serif/model/specifications/physical_names.hpp
 create mode 100644 src/include/serif/model/state.hpp
 create mode 100644 src/include/serif/model/stellar_model.hpp
 create mode 100644 src/include/serif/model/validation.hpp
 create mode 100644 src/lib/serif/model/specifications.cpp

diff --git a/src/include/serif/model/configured_discretization.hpp b/src/include/serif/model/configured_discretization.hpp
new file mode 100644
index 00000000..0a23c097
--- /dev/null
+++ b/src/include/serif/model/configured_discretization.hpp
@@ -0,0 +1,35 @@
+#pragma once
+
+#include <utility>
+
+#include "serif/model/normalization/reference.hpp"
+
+namespace serif::model::preconditioning {
+    /** @brief Evaluate physical equations without preparing an inverse until a recipe is selected. */
+    struct None final {};
+}
+
+namespace serif::model {
+    /** @brief Select a descriptive star, its existing finite-element resources, and its numerical coordinates. */
+    template <StellarModelType Model, typename FiniteElements, typename Normalization = solver::normalization::NormalizationOptions<>, typename Preconditioner = preconditioning::None>
+    struct DiscretizationOptions final {
+        Model model;
+        FiniteElements finite_elements;
+        Normalization normalization;
+        Preconditioner preconditioner;
+    };
+
+    /**
+     * @brief Resolve model defaults while preserving compile-time field, scheme, and coordinate choices.
+     * @note Finite-element resources remain explicit: normalization does not choose a mesh or approximation order.
+     */
+    template <StellarModelType Model, typename FiniteElements, typename Normalization, typename Preconditioner>
+        requires normalization::ReferenceResolvable<Normalization, Model>
+    [[nodiscard]] auto make_discretization(DiscretizationOptions<Model, FiniteElements, Normalization, Preconditioner> options) {
+        auto reference = normalization::resolve_configuration(std::move(options.normalization), options.model);
+        return DiscretizationOptions<Model, FiniteElements, decltype(reference), Preconditioner>{
+            .model = std::move(options.model), .finite_elements = std::move(options.finite_elements), .normalization = std::move(reference),
+            .preconditioner = std::move(options.preconditioner)
+        };
+    }
+}
diff --git a/src/include/serif/model/description.hpp b/src/include/serif/model/description.hpp
new file mode 100644
index 00000000..d6c80b42
--- /dev/null
+++ b/src/include/serif/model/description.hpp
@@ -0,0 +1,7 @@
+#pragma once
+
+#include "serif/model/equilibrium_equations.hpp"
+#include "serif/model/specification.hpp"
+#include "serif/model/specifications/fixed_angular_momentum.hpp"
+#include "serif/model/specifications/fixed_total_mass.hpp"
+#include "serif/model/stellar_model.hpp"
diff --git a/src/include/serif/model/detail/specification_endpoints.hpp b/src/include/serif/model/detail/specification_endpoints.hpp
new file mode 100644
index 00000000..ae0c6401
--- /dev/null
+++ b/src/include/serif/model/detail/specification_endpoints.hpp
@@ -0,0 +1,111 @@
+#pragma once
+
+#include "serif/model/specification.hpp"
+#include "serif/model/specifications/physical_names.hpp"
+#include "serif/utils/dependencies/binding.hpp"
+#include "serif/utils/dependencies/key_list_operations.hpp"
+#include "serif/utils/misc/concepts/filter.hpp"
+
+namespace serif::model::detail
+{
+    template <typename Specification>
+    using GeneratesState = std::bool_constant<Specification::ModelDefinition::generated_state_kind != GeneratedStateKind::none>;
+
+    /** @brief Numerical checks operate only on prescriptions that introduce a coordinate and equation. */
+    template <typename Model>
+    using GeneratedSpecifications = utils::misc::concepts::FilterT<typename Model::SpecificationTypes, GeneratesState>;
+
+    /** @brief Recognize generated endpoint identities without interpreting unrelated physical names. */
+    template <typename Endpoint>
+    struct EndpointOwner { using Type = void; };
+    template <typename Specification>
+    struct EndpointOwner<GeneratedCoordinateOf<Specification>> { using Type = Specification; };
+    template <typename Specification>
+    struct EndpointOwner<ConstraintOf<Specification>> { using Type = Specification; };
+
+    /** @brief Resolve owned aliases and explicit bindings without merging same-unit identities. */
+    template <typename Owner, typename Endpoint, typename Bindings>
+    [[nodiscard]] consteval auto resolved_endpoint()
+    {
+        using namespace utils::dependencies;
+        if constexpr (std::same_as<Endpoint, state::OwnGeneratedCoordinate>)
+        {
+            return resolved_endpoint<Owner, GeneratedCoordinateOf<Owner>, Bindings>();
+        }
+        else if constexpr (std::same_as<Endpoint, equation::OwnConstraint>)
+        {
+            return resolved_endpoint<Owner, ConstraintOf<Owner>, Bindings>();
+        }
+        else if constexpr (!std::same_as<typename EndpointOwner<Endpoint>::Type, void> &&
+            BindingsWithBinding<Bindings, EnforcementOf<typename EndpointOwner<Endpoint>::Type>>)
+        {
+            using Specification = typename EndpointOwner<Endpoint>::Type;
+            using Enforcement = bound_type_t<Bindings, EnforcementOf<Specification>>;
+            if constexpr (std::same_as<Endpoint, GeneratedCoordinateOf<Specification>>)
+            {
+                return std::type_identity<typename Enforcement::Coordinate>{};
+            }
+            else
+            {
+                return std::type_identity<typename Enforcement::Constraint>{};
+            }
+        }
+        else if constexpr (BindingsWithBinding<Bindings, Endpoint>)
+        {
+            return std::type_identity<bound_type_t<Bindings, Endpoint>>{};
+        }
+        else
+        {
+            return std::type_identity<Endpoint>{};
+        }
+    }
+
+    template <typename Owner, typename Endpoint, typename Bindings>
+    using ResolvedEndpoint = typename decltype(resolved_endpoint<Owner, Endpoint, Bindings>())::type;
+
+    /** @brief A coordinate binding cannot supply an absent or non-generating physical prescription. */
+    template <typename Model, typename Endpoint>
+    [[nodiscard]] consteval bool endpoint_owner_present()
+    {
+        if constexpr (requires { typename Endpoint::SpecificationType; })
+        {
+            return Model::template contains_specification<typename Endpoint::SpecificationType> && DimensionalScalarSpecification<typename Endpoint::SpecificationType>;
+        }
+        return true;
+    }
+
+    template <typename Owner, typename Endpoints, typename Bindings>
+    struct ResolvedEndpoints;
+
+    template <typename Owner, typename... Endpoints, typename Bindings>
+    struct ResolvedEndpoints<Owner, utils::types::TypeList<Endpoints...>, Bindings>
+    {
+        using Keys = utils::dependencies::CombinedKeyList<utils::dependencies::KeyList<ResolvedEndpoint<Owner, Endpoints, Bindings>>...>;
+
+        template <typename Model>
+        static constexpr bool owners_present = (endpoint_owner_present<Model, Endpoints>() && ...);
+    };
+
+    /** @brief Legacy prescriptions retain their original dependencies unless an equation explicitly registers another enforcement. */
+    template <typename Specification, typename Bindings, bool = utils::dependencies::BindingsWithBinding<Bindings, EnforcementOf<Specification>>>
+    struct EnforcementDependencies
+    {
+        using Type = typename Specification::ModelDefinition;
+    };
+    template <typename Specification, typename Bindings>
+    struct EnforcementDependencies<Specification, Bindings, true>
+    {
+        using Type = utils::dependencies::bound_type_t<Bindings, EnforcementOf<Specification>>;
+    };
+
+    /** @brief Share resolved physical identities between endpoint and numerical-provider checks. */
+    template <DimensionalScalarSpecification Specification, typename Bindings>
+    struct ResolvedSpecification
+    {
+        using Coordinate = ResolvedEndpoint<Specification, GeneratedCoordinateOf<Specification>, Bindings>;
+        using Constraint = ResolvedEndpoint<Specification, ConstraintOf<Specification>, Bindings>;
+        using Enforcement = typename EnforcementDependencies<Specification, Bindings>::Type;
+        using Reads = ResolvedEndpoints<Specification, typename Enforcement::DependsOn, Bindings>;
+        using Changes = ResolvedEndpoints<Specification, typename Enforcement::Affects, Bindings>;
+    };
+}
diff --git a/src/include/serif/model/detail/specification_list.hpp b/src/include/serif/model/detail/specification_list.hpp
new file mode 100644
index 00000000..38bde844
--- /dev/null
+++ b/src/include/serif/model/detail/specification_list.hpp
@@ -0,0 +1,97 @@
+#pragma once
+
+#include <tuple>
+#include <utility>
+
+#include "serif/model/specification.hpp"
+#include "serif/utils/misc/concepts/filter.hpp"
+
+namespace serif::model::detail
+{
+    /** @brief Compare physical identities in the historical role, coordinate-kind, stable-name order. */
+    template <ModelSpecification Left, ModelSpecification Right>
+    [[nodiscard]] consteval bool specification_precedes()
+    {
+        using L = typename Left::ModelDefinition;
+        using R = typename Right::ModelDefinition;
+        return std::tuple{L::role, L::generated_state_kind, L::name} < std::tuple{R::role, R::generated_state_kind, R::name};
+    }
+
+    template <typename Specification, typename List>
+    struct InsertSpecification;
+
+    template <typename Specification>
+    struct InsertSpecification<Specification, utils::types::TypeList<>>
+    {
+        using Type = utils::types::TypeList<Specification>;
+    };
+
+    template <typename Specification, typename Head, typename... Tail>
+    struct InsertSpecification<Specification, utils::types::TypeList<Head, Tail...>>
+    {
+        using Type = std::conditional_t<specification_precedes<Specification, Head>(), utils::types::TypeList<Specification, Head, Tail...>,
+                                       utils::misc::concepts::PrependT<Head, typename InsertSpecification<Specification, utils::types::TypeList<Tail...>>::Type>>;
+    };
+
+    template <typename... Specifications>
+    struct CanonicalSpecifications;
+
+    template <>
+    struct CanonicalSpecifications<>
+    {
+        using Type = utils::types::TypeList<>;
+    };
+
+    template <typename Head, typename... Tail>
+    struct CanonicalSpecifications<Head, Tail...>
+    {
+        using Type = typename InsertSpecification<Head, typename CanonicalSpecifications<Tail...>::Type>::Type;
+    };
+
+    /** @brief Reject semantic collisions even when two prescriptions have different C++ types or scalar kinds. */
+    template <ModelSpecification... Specifications>
+    [[nodiscard]] consteval bool unique_specification_names()
+    {
+        constexpr std::array<SpecificationDescriptor, sizeof...(Specifications)> descriptions{{
+            {Specifications::ModelDefinition::name, Specifications::ModelDefinition::role, Specifications::ModelDefinition::generated_state_kind}...
+        }};
+        for (std::size_t left = 0; left < descriptions.size(); ++left)
+        {
+            for (std::size_t right = left + 1; right < descriptions.size(); ++right)
+            {
+                if (descriptions[left].role == descriptions[right].role && descriptions[left].name == descriptions[right].name)
+                {
+                    return false;
+                }
+            }
+        }
+        return true;
+    }
+
+    template <SpecificationRole Role, typename... Specifications>
+    struct SelectRole
+    {
+        template <typename Specification>
+        using HasRole = std::bool_constant<Specification::ModelDefinition::role == Role>;
+        using Type = utils::misc::concepts::FilterT<utils::types::TypeList<Specifications...>, HasRole>;
+    };
+
+    template <typename List, typename... Arguments>
+    struct ArgumentsMatch;
+
+    template <typename... Specifications, typename... Arguments>
+    struct ArgumentsMatch<utils::types::TypeList<Specifications...>, Arguments...>
+    {
+        template <typename Specification>
+        static constexpr std::size_t count = (std::size_t{0} + ... + std::size_t{std::same_as<Specification, std::remove_cvref_t<Arguments>>});
+        static constexpr bool value = sizeof...(Specifications) == sizeof...(Arguments) && ((count<Specifications> == 1) && ...);
+    };
+
+    /** @brief Consume the argument tuple once, moving each specification into its canonical slot once. */
+    template <typename... Specifications, typename... Arguments>
+    [[nodiscard]] std::tuple<Specifications...> canonical_specifications(utils::types::TypeList<Specifications...>, Arguments&&... arguments)
+    {
+        std::tuple<std::remove_cvref_t<Arguments>...> values{std::forward<Arguments>(arguments)...};
+        return {std::get<Specifications>(std::move(values))...};
+    }
+} // namespace serif::model::detail
diff --git a/src/include/serif/model/detail/surface_plan.hpp b/src/include/serif/model/detail/surface_plan.hpp
new file mode 100644
index 00000000..4abfb3c5
--- /dev/null
+++ b/src/include/serif/model/detail/surface_plan.hpp
@@ -0,0 +1,80 @@
+#pragma once
+
+#include "serif/geometry/surface/automatic.hpp"
+#include "serif/model/equilibrium_equations.hpp"
+
+namespace serif::model::detail
+{
+    template <typename Equation, typename = void>
+    struct SurfaceFieldHint
+    {
+        using Type = utils::dependencies::KeyList<>;
+    };
+
+    /** @brief An equation author may disambiguate surface enforcement without changing the physical condition. */
+    template <typename Equation>
+    struct SurfaceFieldHint<Equation, std::void_t<typename Equation::IsobaricSurfaceField>>
+    {
+        using Type = utils::dependencies::KeyList<typename Equation::IsobaricSurfaceField>;
+    };
+
+    template <typename Fields>
+    struct UniqueSurfaceField;
+
+    template <>
+    struct UniqueSurfaceField<utils::dependencies::KeyList<>>
+    {
+        using Type = void;
+    };
+
+    template <typename Field, typename... Others>
+    struct UniqueSurfaceField<utils::dependencies::KeyList<Field, Others...>>
+    {
+        static_assert(sizeof...(Others) == 0, "The selected equations name conflicting IsobaricSurfaceField types. Choose one physical field for the isobaric condition.");
+        using Type = Field;
+    };
+
+    template <typename Equations>
+    struct PreferredSurfaceField;
+
+    template <typename... Equations>
+    struct PreferredSurfaceField<utils::types::TypeList<Equations...>> : UniqueSurfaceField<utils::dependencies::CombinedKeyList<typename SurfaceFieldHint<Equations>::Type...>>
+    {
+    };
+
+    /** @brief A disabled or unavailable boundary has no compiled field or row requirements. */
+    template <typename Selection, bool Supported>
+    struct SelectedSurface
+    {
+        using Formulation = void;
+        using Coefficients = utils::dependencies::KeyList<>;
+        using Rows = utils::dependencies::KeyList<>;
+        static constexpr bool supported = Supported;
+    };
+
+    template <typename Selection>
+    struct SelectedSurface<Selection, true>
+    {
+        using Formulation = typename Selection::Formulation;
+        using Coefficients = typename Formulation::StateTerms;
+        using Rows = utils::dependencies::KeyList<typename Formulation::Test>;
+        static constexpr bool supported = true;
+    };
+
+    template <typename Condition, typename EOS, typename Fields, typename Rows, typename Solved, typename Equations>
+    struct SurfacePlan;
+
+    template <typename EOS, typename Fields, typename Rows, typename Solved, typename Equations>
+    struct SurfacePlan<NoSurfaceCondition, EOS, Fields, Rows, Solved, Equations> : SelectedSurface<void, false>
+    {
+        static constexpr bool supported = true;
+    };
+
+    /** @brief Compile the condition once; planning, graph assembly and validation share the selected formulation. */
+    template <typename EOS, typename Fields, typename Rows, typename Solved, typename Equations>
+    struct SurfacePlan<geometry::surface::Isobaric, EOS, Fields, Rows, Solved, Equations> : SelectedSurface<
+        geometry::surface::AutomaticPressureSurface<EOS, Fields, Rows, Solved, typename PreferredSurfaceField<Equations>::Type>,
+        geometry::surface::AutomaticPressureSurfaceCompilable<EOS, Fields, Rows, Solved, typename PreferredSurfaceField<Equations>::Type>>
+    {
+    };
+}
diff --git a/src/include/serif/model/eos_reference.hpp b/src/include/serif/model/eos_reference.hpp
new file mode 100644
index 00000000..8665001c
--- /dev/null
+++ b/src/include/serif/model/eos_reference.hpp
@@ -0,0 +1,96 @@
+#pragma once
+
+#include <memory>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+
+#include "serif/utils/error/errors.hpp"
+
+namespace serif::model
+{
+    /** @brief Reuse one immutable EOS in equation objects, optionally retaining shared ownership. */
+    template <typename EOS>
+        requires requires { typename EOS::Relations; }
+    class EosReference final
+    {
+    public:
+        using Relations = typename EOS::Relations;
+
+        /** @brief Retain an owned EOS, including an aliasing handle into an immutable model. */
+        explicit EosReference(std::shared_ptr<const EOS> equation_of_state) : m_owner(std::move(equation_of_state)), m_eos(m_owner.get())
+        {
+            using namespace utils::errors;
+            serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(m_eos != nullptr, "An EOS reference requires an equation of state; supply a non-null shared owner.");
+        }
+
+        /** @brief Borrow an EOS whose model owner must outlive all equations and prepared actions using it. */
+        explicit EosReference(const EOS& equation_of_state) noexcept : m_eos(&equation_of_state)
+        {
+        }
+
+        /** @brief Never retain a pointer into a temporary EOS. */
+        EosReference(EOS&&) = delete;
+        /** @brief Const temporaries have the same borrowing lifetime restriction. */
+        EosReference(const EOS&&) = delete;
+
+        /** @brief Forward only the relations actually implemented by the underlying EOS. */
+        template <typename Relation, typename... Values>
+            requires requires(const EOS& eos, Relation relation, Values... values) { eos.evaluate(relation, values...); }
+        [[nodiscard]] decltype(auto) evaluate(Relation relation, Values... values) const noexcept(noexcept(m_eos->evaluate(relation, values...)))
+        {
+            return m_eos->evaluate(relation, values...);
+        }
+
+        /** @brief Preserve exact derivative availability and physical return types. */
+        template <typename Relation, typename Input, typename... Values>
+            requires requires(const EOS& eos, Relation relation, Input input, Values... values) { eos.partial_derivative(relation, input, values...); }
+        [[nodiscard]] decltype(auto) partial_derivative(Relation relation, Input input, Values... values) const noexcept(noexcept(m_eos->partial_derivative(relation, input, values...)))
+        {
+            return m_eos->partial_derivative(relation, input, values...);
+        }
+
+        /** @brief Forward second partials only when the underlying EOS provides them. */
+        template <typename Relation, typename FirstInput, typename SecondInput, typename... Values>
+            requires requires(const EOS& eos, Relation relation, FirstInput first, SecondInput second, Values... values) {
+                eos.second_partial_derivative(relation, first, second, values...);
+            }
+        [[nodiscard]] decltype(auto) second_partial_derivative(Relation relation, FirstInput first, SecondInput second, Values... values) const
+            noexcept(noexcept(m_eos->second_partial_derivative(relation, first, second, values...))) {
+            return m_eos->second_partial_derivative(relation, first, second, values...);
+        }
+
+        /** @brief Inspect EOS-specific physical parameters without storing another copy. */
+        [[nodiscard]] const EOS& get() const noexcept
+        {
+            return *m_eos;
+        }
+
+    private:
+        std::shared_ptr<const EOS> m_owner;
+        const EOS* m_eos;
+    };
+
+    /** @brief Share the EOS subobject while retaining the model that owns it. */
+    template <typename Model>
+    [[nodiscard]] auto share_equation_of_state(std::shared_ptr<const Model> model)
+    {
+        using namespace utils::errors;
+        serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(model != nullptr, "Sharing an equation of state requires a non-null stellar model owner.");
+        using EOS = std::remove_cvref_t<decltype(model->equation_of_state())>;
+        const EOS* eos = &model->equation_of_state();
+        return EosReference<EOS>{std::shared_ptr<const EOS>{std::move(model), eos}};
+    }
+
+    /** @brief Borrow the model-owned EOS when the prepared graph already retains its immutable model owner. */
+    template <typename EOS>
+    [[nodiscard]] EosReference<EOS> borrow_eos(const EOS& eos) noexcept
+    {
+        return EosReference<EOS>{eos};
+    }
+
+    /** @brief Borrowing never extends a temporary EOS lifetime. */
+    template <typename EOS>
+        requires (!std::is_lvalue_reference_v<EOS>)
+    void borrow_eos(EOS&&) = delete;
+} // namespace serif::model
diff --git a/src/include/serif/model/equilibrium_equations.hpp b/src/include/serif/model/equilibrium_equations.hpp
new file mode 100644
index 00000000..3c9f57ca
--- /dev/null
+++ b/src/include/serif/model/equilibrium_equations.hpp
@@ -0,0 +1,101 @@
+#pragma once
+
+#include <functional>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+namespace serif::model
+{
+    /** @brief Explicitly omit a shape equation from a declarative formulation. */
+    struct NoShapeForce
+    {
+    };
+
+    /** @brief Explicitly omit a surface prescription from a declarative formulation. */
+    struct NoSurfaceCondition
+    {
+    };
+
+    /** @brief State the chosen equations explicitly, without selecting physics from the EOS. */
+    template <typename Interior, typename Shape = NoShapeForce, typename Surface = NoSurfaceCondition>
+    struct EquilibriumEquations
+    {
+        using InteriorEquations = std::remove_cvref_t<Interior>;
+        using ShapeForce = std::remove_cvref_t<Shape>;
+        using SurfaceCondition = std::remove_cvref_t<Surface>;
+
+        Interior interior_equations;
+        Shape shape_force{};
+        Surface surface_condition{};
+    };
+
+    /**
+     * @brief Preserve every authored member when a declaration passes through a templated model factory.
+     * @note Explicit deduction avoids compiler substitution dropping defaulted aggregate members. Deduced
+     * declarations supply preceding slots: use NoShapeForce before a surface condition when the shape force is
+     * absent. Explicitly typed aggregates may still skip members.
+     */
+    template <typename Interior, typename Shape = NoShapeForce, typename Surface = NoSurfaceCondition>
+    EquilibriumEquations(Interior, Shape = {}, Surface = {}) -> EquilibriumEquations<Interior, Shape, Surface>;
+
+    namespace detail
+    {
+        template <typename Candidate>
+        inline constexpr bool is_equilibrium_equations = false;
+
+        template <typename Interior, typename Shape, typename Surface>
+        inline constexpr bool is_equilibrium_equations<EquilibriumEquations<Interior, Shape, Surface>> = true;
+
+        /** @brief Evaluate a physical definition or retain its immutable object without requiring a copy. */
+        template <typename Definition, typename Model>
+        [[nodiscard]] decltype(auto) resolve_equation_member(Definition&& definition, const Model& model)
+        {
+            if constexpr (std::invocable<Definition, const Model&>)
+            {
+                return std::invoke(std::forward<Definition>(definition), model);
+            }
+            else if constexpr (std::is_lvalue_reference_v<Definition&&>)
+            {
+                return std::forward<Definition>(definition);
+            }
+            else
+            {
+                return std::remove_cvref_t<Definition>{std::forward<Definition>(definition)};
+            }
+        }
+
+        /** @brief Resolve each explicit part, owning temporary results and borrowing unchanged model members. */
+        template <typename Definition, typename Model>
+        [[nodiscard]] auto resolve_equation_members(Definition&& definition, const Model& model)
+        {
+            using Interior = decltype(resolve_equation_member(std::forward<Definition>(definition).interior_equations, model));
+            using Shape = decltype(resolve_equation_member(std::forward<Definition>(definition).shape_force, model));
+            using Surface = decltype(resolve_equation_member(std::forward<Definition>(definition).surface_condition, model));
+            return EquilibriumEquations<Interior, Shape, Surface>
+            {
+                .interior_equations = resolve_equation_member(std::forward<Definition>(definition).interior_equations, model),
+                .shape_force = resolve_equation_member(std::forward<Definition>(definition).shape_force, model),
+                .surface_condition = resolve_equation_member(std::forward<Definition>(definition).surface_condition, model)
+            };
+        }
+    } // namespace detail
+
+    /** @brief Recognize composition without instantiating any numerical equation contract. */
+    template <typename Candidate>
+    concept EquilibriumEquationDescription = detail::is_equilibrium_equations<std::remove_cvref_t<Candidate>>;
+
+    /** @brief Resolve an explicit composition at binding; retain the model while the result borrows its members. */
+    template <typename Model>
+        requires requires(const Model& model) { model.equations(); } &&
+                 EquilibriumEquationDescription<decltype(detail::resolve_equation_member(std::declval<const Model&>().equations(), std::declval<const Model&>()))>
+    [[nodiscard]] auto resolve_equilibrium_equations(const Model& model)
+    {
+        return detail::resolve_equation_members(detail::resolve_equation_member(model.equations(), model), model);
+    }
+
+    /** @brief A temporary model cannot retain the equation members borrowed by the resolved description. */
+    template <typename Model>
+        requires (!std::is_lvalue_reference_v<Model>)
+    void resolve_equilibrium_equations(Model&&) = delete;
+} // namespace serif::model
diff --git a/src/include/serif/model/equilibrium_plan.hpp b/src/include/serif/model/equilibrium_plan.hpp
new file mode 100644
index 00000000..11fcc19a
--- /dev/null
+++ b/src/include/serif/model/equilibrium_plan.hpp
@@ -0,0 +1,66 @@
+#pragma once
+
+#include "serif/model/equilibrium_equations.hpp"
+#include "serif/model/plan.hpp"
+#include "serif/model/detail/surface_plan.hpp"
+#include "serif/operators/concepts.hpp"
+#include "serif/equations/surface_potential.hpp"
+#include "serif/equations/primal_gravity.hpp"
+
+namespace serif::model
+{
+    namespace detail
+    {
+        template <typename Shape>
+        struct ShapePlan;
+
+        template <>
+        struct ShapePlan<NoShapeForce>
+        {
+            using Coefficients = utils::dependencies::KeyList<>;
+            using Values = utils::dependencies::KeyList<>;
+            using Rows = utils::dependencies::KeyList<>;
+        };
+
+        /** @brief The isobar level couples only existing scalar potential, rotation and mass coordinates. */
+        template <>
+        struct ShapePlan<equations::SurfacePotentialBalance>
+        {
+            using Coefficients = typename OperandTerms<equations::SurfacePotentialBalance::Operands>::Type;
+            using Values = utils::dependencies::KeyList<discretization::blocks::surface_deformation::coefficients>;
+            using Rows = utils::dependencies::KeyList<discretization::blocks::surface_deformation::shape_equilibrium>;
+        };
+
+    }
+
+    /**
+     * @brief Collect the explicitly selected physical coefficients and supported constraint coordinates.
+     * @note This plan verifies representation, not mathematical closure. Equal units never merge distinct
+     * physical identities, and prescribed coefficients remain available as graph inputs without becoming unknowns.
+     */
+    template <typename EOS, EquilibriumEquationDescription Description, utils::dependencies::IsKeyList Prescribed = utils::dependencies::KeyList<>,
+              utils::dependencies::IsKeyList AdditionalFields = utils::dependencies::KeyList<>>
+    struct EquilibriumPlan final
+    {
+        static constexpr bool has_shape = !std::same_as<typename Description::ShapeForce, NoShapeForce>;
+        static constexpr bool has_surface = !std::same_as<typename Description::SurfaceCondition, NoSurfaceCondition>;
+        static constexpr bool has_primal_gravity = []<typename... Equations>(utils::types::TypeList<Equations...>) {
+            return (std::same_as<Equations, equations::PrimalPoissonBalance> || ...);
+        }(typename Description::InteriorEquations::Equations{});
+        using Interior = detail::InteriorPlan<typename Description::InteriorEquations>;
+        using Shape = detail::ShapePlan<typename Description::ShapeForce>;
+        using ResidualEquations = typename Interior::EquationsList;
+        using AvailableFields = utils::dependencies::CombinedKeyList<typename Interior::Coefficients, typename Shape::Coefficients, AdditionalFields>;
+        using SolvedFields = typename detail::RemoveKeys<AvailableFields, Prescribed>::Type;
+        using Residuals = utils::dependencies::CombinedKeyList<typename Interior::Rows, typename Shape::Rows>;
+        using Surface = detail::SurfacePlan<typename Description::SurfaceCondition, EOS, AvailableFields, Residuals, SolvedFields, typename Interior::EquationsList>;
+        using Coefficients = utils::dependencies::CombinedKeyList<typename Interior::Coefficients, typename Shape::Coefficients, typename Surface::Coefficients>;
+        static_assert(!has_shape || !utils::dependencies::ContainsKey<discretization::blocks::surface_deformation::coefficients, Prescribed>,
+            "Prescribe geometry by selecting NoShapeForce and supplying a displacement field. An active shape-force equation uses independent solved surface coordinates.");
+        using Values = utils::dependencies::CombinedKeyList<typename detail::RemoveKeys<Coefficients, Prescribed>::Type, typename Shape::Values>;
+        // Geometry is supplied separately even when no equation solves for the shape.
+        using PhysicalInputs = utils::dependencies::CombinedKeyList<Coefficients, Values, utils::dependencies::KeyList<discretization::blocks::surface_deformation::coefficients>>;
+        using InactiveRows = std::conditional_t<has_shape, utils::dependencies::KeyList<>, utils::dependencies::KeyList<discretization::blocks::surface_deformation::shape_equilibrium>>;
+        using Layout = discretization::blocks::SystemLayout<Values, Residuals>;
+    };
+}
diff --git a/src/include/serif/model/execution_validation.hpp b/src/include/serif/model/execution_validation.hpp
new file mode 100644
index 00000000..71124362
--- /dev/null
+++ b/src/include/serif/model/execution_validation.hpp
@@ -0,0 +1,144 @@
+#pragma once
+
+#include "serif/model/equilibrium_plan.hpp"
+#include "serif/model/validation.hpp"
+
+namespace serif::model
+{
+    namespace detail
+    {
+        template <typename Rows, typename Fields>
+        struct CartesianCouplings;
+
+        template <typename... Rows, typename... Fields>
+        struct CartesianCouplings<utils::dependencies::KeyList<Rows...>, utils::dependencies::KeyList<Fields...>>
+        {
+            template <typename Row>
+            using ForRow = utils::dependencies::KeyList<Derivative<Row, Fields>...>;
+            using Type = utils::dependencies::CombinedKeyList<ForRow<Rows>...>;
+        };
+
+        template <typename Plan, typename = void>
+        struct ExecutionSurface : SelectedSurface<void, false>
+        {
+        };
+
+        template <typename Plan>
+        struct ExecutionSurface<Plan, std::void_t<typename Plan::Surface>> : Plan::Surface
+        {
+        };
+
+        template <typename Plan, typename = void>
+        struct ExecutionInputs
+        {
+            using PhysicalInputs = typename Plan::Values;
+            using InactiveRows = utils::dependencies::KeyList<>;
+        };
+
+        template <typename Plan>
+        struct ExecutionInputs<Plan, std::void_t<typename Plan::PhysicalInputs, typename Plan::InactiveRows>>
+        {
+            using PhysicalInputs = typename Plan::PhysicalInputs;
+            using InactiveRows = typename Plan::InactiveRows;
+        };
+
+        /** @brief Demand an explicitly supplied generated equation and effect before inferring closed-row zeros. */
+        template <DimensionalScalarSpecification Specification, typename Rows, typename InactiveRows, typename Actual, typename Bindings>
+        [[nodiscard]] consteval bool generated_physics_is_present()
+        {
+            using namespace utils::dependencies;
+            using Endpoints = ResolvedSpecification<Specification, Bindings>;
+            return ContainsKey<typename Endpoints::Constraint, Rows> && []<typename... Changed>(KeyList<Changed...>)
+            {
+                return (((ContainsKey<Changed, InactiveRows> && !ContainsKey<Changed, Rows>) || ContainsKey<Derivative<Changed, typename Endpoints::Coordinate>, Actual>) && ...);
+            }(typename Endpoints::Changes::Keys{});
+        }
+
+        /** @brief Check only numerical-provider admission; physical endpoint ownership is checked separately. */
+        template <typename Rows, typename InactiveRows, typename Actual, typename Bindings, typename... Specifications>
+        [[nodiscard]] consteval bool generated_physics_are_present(utils::types::TypeList<Specifications...>)
+        {
+            return (generated_physics_is_present<Specifications, Rows, InactiveRows, Actual, Bindings>() && ...);
+        }
+
+        /** @brief Every surface carrier must have an equation row in the selected layout. */
+        template <typename Rows, typename... Required>
+        [[nodiscard]] consteval bool equation_rows_exist(utils::dependencies::KeyList<Required...>)
+        {
+            return (utils::dependencies::ContainsKey<Required, Rows> && ...);
+        }
+    } // namespace detail
+
+    /** @brief Recognize current equation backends without instantiating an invalid complete plan. */
+    template <typename Description>
+    concept ExecutableEquilibriumDescription = EquilibriumEquationDescription<Description> && ExecutableInteriorEquations<typename Description::InteriorEquations> && requires
+    {
+        typename detail::ShapePlan<typename Description::ShapeForce>::Rows;
+        requires std::same_as<typename Description::SurfaceCondition, NoSurfaceCondition> || std::same_as<typename Description::SurfaceCondition, geometry::surface::Isobaric>;
+        requires std::same_as<typename Description::ShapeForce, NoShapeForce> ||
+            std::same_as<typename Description::ShapeForce, equations::SurfacePotentialBalance>;
+    };
+
+    namespace detail
+    {
+        template <ExecutableEquilibriumDescription Description, typename Surface>
+        struct EquilibriumCouplings
+        {
+            using Interior = InteriorPlan<typename Description::InteriorEquations>;
+            using Shape = ShapePlan<typename Description::ShapeForce>;
+            using GeometryRows = std::conditional_t<std::same_as<typename Description::ShapeForce, NoShapeForce>, utils::dependencies::KeyList<>,
+                                                   utils::dependencies::CombinedKeyList<typename Interior::Rows, typename Shape::Rows>>;
+            using Geometry = typename CartesianCouplings<GeometryRows, utils::dependencies::KeyList<discretization::blocks::surface_deformation::coefficients>>::Type;
+            using Force = typename CartesianCouplings<typename Shape::Rows, typename Shape::Coefficients>::Type;
+            using Type = utils::dependencies::CombinedKeyList<DeclaredEquationCouplings<typename Description::InteriorEquations>, Force, Geometry,
+                typename CartesianCouplings<typename Surface::Rows, typename Surface::Coefficients>::Type>;
+            using AuthoredRows = utils::dependencies::CombinedKeyList<typename Interior::Rows, typename Shape::Rows>;
+            using RequiredRows = utils::dependencies::CombinedKeyList<AuthoredRows, typename Surface::Rows>;
+        };
+    } // namespace detail
+
+    /** @brief Inspect the actual derivative producers admitted by the explicitly supplied physical formulation. */
+    template <ExecutableEquilibriumDescription Description, typename Plan>
+    using ActualEquilibriumCouplings = typename detail::EquilibriumCouplings<Description, detail::ExecutionSurface<Plan>>::Type;
+
+    /**
+     * @brief Validate a complete selected formulation without inferring solvability or selecting physics.
+     * @note Closed equations define zero derivatives for other coordinates, but that fact cannot admit a
+     * missing custom multiplier force. ExplicitZeroCouplings is an explicit author registration, not a fallback.
+     */
+    template <typename Model, typename Description, typename Plan, typename Bindings = DefaultSpecificationBindings,
+              typename ExplicitZeroCouplings = utils::dependencies::KeyList<>>
+    [[nodiscard]] consteval bool executable_model_topology()
+    {
+        using namespace utils::dependencies;
+        if constexpr (!StellarModelType<Model> || !ExecutableEquilibriumDescription<Description> || !IsBindings<Bindings> || !IsKeyList<ExplicitZeroCouplings>)
+        {
+            return false;
+        }
+        else if constexpr (!requires { typename Plan::Values; typename Plan::Residuals; })
+        {
+            return false;
+        }
+        else if constexpr (!IsKeyList<typename Plan::Values> || !IsKeyList<typename Plan::Residuals>)
+        {
+            return false;
+        }
+        else
+        {
+            using Surface = detail::ExecutionSurface<Plan>;
+            using Inputs = detail::ExecutionInputs<Plan>;
+            using Couplings = detail::EquilibriumCouplings<Description, Surface>;
+            if constexpr (!SpecificationEndpointsValid<Model, typename Plan::Values, typename Inputs::PhysicalInputs, typename Plan::Residuals, typename Inputs::InactiveRows, Bindings>)
+            {
+                return false;
+            }
+            else
+            {
+                using Actual = CombinedKeyList<typename Couplings::Type, ExplicitZeroCouplings>;
+                return (std::same_as<typename Description::SurfaceCondition, NoSurfaceCondition> || Surface::supported) &&
+                       detail::equation_rows_exist<typename Plan::Residuals>(typename Couplings::RequiredRows{}) &&
+                       detail::generated_physics_are_present<typename Couplings::AuthoredRows, typename Inputs::InactiveRows, Actual, Bindings>(detail::GeneratedSpecifications<std::remove_cvref_t<Model>>{});
+            }
+        }
+    }
+} // namespace serif::model
diff --git a/src/include/serif/model/finite_elements.hpp b/src/include/serif/model/finite_elements.hpp
new file mode 100644
index 00000000..4d78d3ed
--- /dev/null
+++ b/src/include/serif/model/finite_elements.hpp
@@ -0,0 +1,13 @@
+#pragma once
+
+#include "serif/backend/mfem/model/finite_elements.hpp"
+
+namespace serif::model {
+    // Public construction names retain the backend's native resource ownership without forwarding calls.
+    using backend::mfem::make_stellar_discretization;
+    using backend::mfem::with_deformation;
+}
+
+namespace serif::geometry::deformation {
+    using backend::mfem::deformation_discretization;
+}
diff --git a/src/include/serif/model/keys.hpp b/src/include/serif/model/keys.hpp
new file mode 100644
index 00000000..e3fe140e
--- /dev/null
+++ b/src/include/serif/model/keys.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+namespace serif::model::keys
+{
+    /** @brief The complete accepted physical residual and its matching Jacobian. */
+    struct PhysicalSystem
+    {
+    };
+}
diff --git a/src/include/serif/model/mesh_options.hpp b/src/include/serif/model/mesh_options.hpp
new file mode 100644
index 00000000..248b647f
--- /dev/null
+++ b/src/include/serif/model/mesh_options.hpp
@@ -0,0 +1,34 @@
+#pragma once
+
+#include <array>
+#include "serif/dimensions/dimensions.hpp"
+
+namespace serif::model {
+    /** @brief Add mesh resolution at the stellar interface before creating field spaces. */
+    struct StellarSurfaceRefinementOptions {
+        /// Additional isotropic levels on both sides of the entire stellar surface.
+        int levels{0};
+        /// Additional latitude splits in closed rings around the rotation axis.
+        int equatorial_levels{0};
+        /// Minimum half-width of the resolved equatorial band, in degrees; 0 < angle <= 35.
+        double equatorial_half_angle_degrees{15.0};
+        /// Additional radial halvings of the stellar-side surface shell only; angular surface faces are unchanged.
+        int radial_levels{0};
+        /// Align the mesh's polar direction with this finite nonzero axis through the stellar center.
+        std::array<double, 3> axis{0.0, 0.0, 1.0};
+    };
+
+    /**
+     * @brief Resolve a spherical stellar interior and a compactified exterior on a curved reference mesh.
+     * @note Field approximation orders come from the model's term types. Geometry order controls only
+     * the mesh; its default resolves the coarse spherical boundary for radial seed projection.
+     */
+    struct StellarMeshOptions {
+        dimensions::LengthValue radius{1.0};
+        double core_radius_fraction{0.25};
+        double infinity_radius_fraction{4.0};
+        int refinement{0};
+        int geometry_order{6};
+        StellarSurfaceRefinementOptions surface_refinement{};
+    };
+}
diff --git a/src/include/serif/model/model.hpp b/src/include/serif/model/model.hpp
new file mode 100644
index 00000000..a8aade37
--- /dev/null
+++ b/src/include/serif/model/model.hpp
@@ -0,0 +1,7 @@
+#pragma once
+
+#include "serif/model/description.hpp"
+#include "serif/model/eos_reference.hpp"
+#include "serif/model/equilibrium_plan.hpp"
+#include "serif/model/validation.hpp"
+#include "serif/model/state.hpp"
diff --git a/src/include/serif/model/normalization/plan.hpp b/src/include/serif/model/normalization/plan.hpp
new file mode 100644
index 00000000..c42c6794
--- /dev/null
+++ b/src/include/serif/model/normalization/plan.hpp
@@ -0,0 +1,131 @@
+#pragma once
+
+#include "serif/model/equilibrium_plan.hpp"
+#include "serif/solver/normalization/configuration.hpp"
+#include "serif/solver/normalization/physical_coordinates.hpp"
+
+namespace serif::model::normalization {
+    /** @brief Reduced equation rows whose volume residual is replaced by the selected surface condition. */
+    struct SurfaceRows final {};
+    /** @brief Potential rows replaced by the zero trace at compactified infinity. */
+    struct InfinityRows final {};
+
+    namespace detail {
+        namespace numerical = solver::normalization;
+        namespace dependencies = utils::dependencies;
+        namespace blocks = discretization::blocks;
+
+        template <typename Coordinates, typename Quantity, typename Metric>
+        using Coordinate = numerical::DimensionlessRieszCoordinate<Quantity, Metric>;
+
+        template <typename Term, bool Spatial = blocks::SpatialTerm<Term>>
+        struct Metric {
+            using Type = numerical::UnitMetric;
+        };
+
+        template <typename Term>
+        struct Metric<Term, true> {
+            using Type = numerical::ReferenceVolumeMetric<Term>;
+        };
+
+        template <typename Coordinates, typename Term>
+        struct ValueRule {
+            using Type = dependencies::Binding<Term, Coordinate<Coordinates, typename Term::Quantity, typename Metric<Term>::Type>>;
+        };
+
+        template <typename Coordinates>
+        struct ValueRule<Coordinates, blocks::surface_deformation::coefficients> {
+            using Type = dependencies::Binding<blocks::surface_deformation::coefficients,
+                Coordinate<Coordinates, dimensions::Length, geometry::deformation::keys::SurfaceMetric>>;
+        };
+
+        template <typename Method, typename Row, typename Plan, bool Surface = Plan::has_surface>
+        struct SurfaceRule {
+            using Type = Method;
+        };
+
+        template <typename Method, typename Row, typename Plan>
+        struct SurfaceRule<Method, Row, Plan, true> {
+            using Formulation = typename Plan::Surface::Formulation;
+            using Type = std::conditional_t<std::same_as<Row, typename Formulation::Test>,
+                numerical::WithPointRows<Method, typename Formulation::ResidualQuantity, SurfaceRows>, Method>;
+        };
+
+        /** @brief Pointwise infinity values have potential units rather than integrated Poisson units. */
+        template <typename Method, typename Row, typename Plan>
+        struct InfinityRule {
+            using Type = std::conditional_t<Plan::has_primal_gravity && std::same_as<Row, blocks::gravity::potential>,
+                numerical::WithPointRows<Method, dimensions::GravitationalPotential, InfinityRows>, Method>;
+        };
+
+        template <typename Interior, typename Shape>
+        struct ResidualRules;
+
+        template <typename... Interior, typename... Shape>
+        struct ResidualRules<dependencies::Bindings<Interior...>, dependencies::Bindings<Shape...>> {
+            using Type = dependencies::Bindings<Interior..., Shape...>;
+        };
+
+        template <typename Coordinates, typename Plan, typename Values, typename Equations>
+        struct AutomaticRules;
+
+        /** @brief Equation residual quantities remain independent from the unknowns occupying their rows. */
+        template <typename Coordinates, typename Plan, typename... Values, typename... Equations>
+        struct AutomaticRules<Coordinates, Plan, dependencies::KeyList<Values...>, utils::types::TypeList<Equations...>> {
+            static_assert((requires { typename Equations::ResidualQuantity; } && ...),
+                "Automatic normalization requires a residual quantity for each equation; otherwise supply an explicit normalization plan.");
+            using ValueRules = dependencies::Bindings<typename ValueRule<Coordinates, Values>::Type...>;
+            using InteriorRules = dependencies::Bindings<dependencies::Binding<typename Equations::Test,
+                typename InfinityRule<typename SurfaceRule<Coordinate<Coordinates, typename Equations::ResidualQuantity, typename Metric<typename Equations::Test>::Type>, typename Equations::Test, Plan>::Type,
+                    typename Equations::Test, Plan>::Type>...>;
+            using ShapeRules = std::conditional_t<Plan::has_shape,
+                dependencies::Bindings<dependencies::Binding<blocks::surface_deformation::shape_equilibrium,
+                    Coordinate<Coordinates, dimensions::GravitationalPotential, numerical::UnitMetric>>>, dependencies::Bindings<>>;
+            using Type = numerical::NormalizationPlan<ValueRules, typename ResidualRules<InteriorRules, ShapeRules>::Type>;
+        };
+
+        template <typename Selection, typename Coordinates, typename Bound>
+        struct SelectPlan {
+            using Type = Selection;
+        };
+
+        template <typename Coordinates, typename Bound>
+        struct SelectPlan<numerical::AutomaticPlan, Coordinates, Bound> : AutomaticRules<Coordinates, typename Bound::Plan,
+            typename Bound::Plan::Values, typename Bound::Plan::ResidualEquations> {};
+
+        template <typename Method>
+        struct MethodMetrics {
+            using Type = dependencies::KeyList<>;
+        };
+
+        template <typename Quantity, typename Metric>
+        struct MethodMetrics<numerical::PhysicalRieszCoordinate<Quantity, Metric>> {
+            using Type = dependencies::KeyList<Metric>;
+        };
+
+        template <typename Quantity, typename Metric, numerical::ReferenceGeometry Geometry>
+        struct MethodMetrics<numerical::DimensionlessRieszCoordinate<Quantity, Metric, Geometry>> {
+            using Type = dependencies::KeyList<Metric>;
+        };
+
+        template <typename Method, typename Quantity, typename Selection>
+        struct MethodMetrics<numerical::WithPointRows<Method, Quantity, Selection>> : MethodMetrics<Method> {};
+
+        template <typename Rules>
+        struct RuleMetrics;
+
+        template <typename... Rules>
+        struct RuleMetrics<dependencies::Bindings<Rules...>> {
+            using Type = dependencies::CombinedKeyList<typename MethodMetrics<typename Rules::Target>::Type...>;
+        };
+    }
+
+    /** @brief Compile automatic physics-aware rules or retain the user's complete explicit plan. */
+    template <typename Selection, typename Coordinates, typename Bound>
+    using PlanFor = typename detail::SelectPlan<Selection, Coordinates, Bound>::Type;
+
+    /** @brief Assemble each metric once even when several state and residual rules share it. */
+    template <typename Plan>
+    using MetricsFor = utils::dependencies::CombinedKeyList<typename detail::RuleMetrics<typename Plan::Values>::Type,
+        typename detail::RuleMetrics<typename Plan::Residuals>::Type>;
+}
diff --git a/src/include/serif/model/normalization/reference.hpp b/src/include/serif/model/normalization/reference.hpp
new file mode 100644
index 00000000..e7e94b97
--- /dev/null
+++ b/src/include/serif/model/normalization/reference.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <utility>
+
+#include "serif/model/description.hpp"
+#include "serif/solver/normalization/configuration.hpp"
+
+namespace serif::model::normalization {
+    namespace detail {
+        /** @brief Explicit reference values take precedence over model defaults without changing physical prescriptions. */
+        template <typename Source, typename Model, typename Specification>
+        concept ReferenceAvailable = !std::same_as<Source, solver::normalization::FromModel> || Model::template contains_specification<Specification>;
+
+        /** @brief Resolve a typed reference source only when its model default is available. */
+        template <typename Specification, typename Source, typename Model>
+            requires ReferenceAvailable<Source, Model, Specification>
+        [[nodiscard]] auto reference_value(const Source& source, const Model& model) {
+            if constexpr (std::same_as<Source, solver::normalization::FromModel>) {
+                return model.template specification<Specification>().target();
+            } else {
+                return source;
+            }
+        }
+    }
+
+    /** @brief Resolve the Lane–Emden reference from explicit inputs or the model's mass target; no model prescribes central density. */
+    template <typename Mass, typename Density, StellarModelType Model>
+        requires detail::ReferenceAvailable<Mass, Model, FixedTotalMass> && std::same_as<Density, dimensions::DensityValue>
+    [[nodiscard]] auto resolve_reference(const solver::normalization::LaneEmdenApproximation<Mass, Density>& scheme, const Model& model) {
+        return solver::normalization::LaneEmdenApproximation<dimensions::MassValue, dimensions::DensityValue>{
+            .n = scheme.n,
+            .mass = detail::reference_value<FixedTotalMass>(scheme.mass, model),
+            .central_density = scheme.central_density,
+            .integration = scheme.integration
+        };
+    }
+
+    /** @brief Explicit or custom reference schemes need no model-specific defaults. */
+    template <solver::normalization::ReferenceScheme Scheme, StellarModelType Model>
+    [[nodiscard]] Scheme resolve_reference(const Scheme& scheme, const Model&) {
+        return scheme;
+    }
+
+    /** @brief Check model-provided reference requirements before instantiating numerical preparation. */
+    template <typename Options, typename Model>
+    concept ReferenceResolvable = requires(const Options& options, const Model& model) {
+        resolve_reference(options.scheme, model);
+    };
+
+    /** @brief Freeze resolved reference choices while retaining the deduced coordinate and plan types. */
+    template <typename Options, StellarModelType Model>
+        requires ReferenceResolvable<Options, Model>
+    [[nodiscard]] auto resolve_configuration(Options options, const Model& model) {
+        return solver::normalization::NormalizationOptions{
+            .scheme = resolve_reference(options.scheme, model),
+            .coordinates = options.coordinates,
+            .plan = options.plan,
+            .gravitational_constant = options.gravitational_constant,
+            .sources = std::move(options.sources),
+            .computations = std::move(options.computations)
+        };
+    }
+}
diff --git a/src/include/serif/model/plan.hpp b/src/include/serif/model/plan.hpp
new file mode 100644
index 00000000..ebf6a245
--- /dev/null
+++ b/src/include/serif/model/plan.hpp
@@ -0,0 +1,68 @@
+#pragma once
+
+#include <tuple>
+#include <type_traits>
+
+#include "serif/discretization/blocks/system_layout.hpp"
+#include "serif/equations/registry.hpp"
+#include "serif/operators/plan.hpp"
+
+namespace serif::model
+{
+    namespace detail
+    {
+        template <typename Operands>
+        struct OperandTerms;
+
+        template <typename... Operands>
+        struct OperandTerms<utils::types::TypeList<Operands...>>
+        {
+            using Type = utils::dependencies::CombinedKeyList<utils::dependencies::KeyList<typename Operands::term_type>...>;
+        };
+
+        template <typename Registry>
+        struct InteriorPlan;
+
+        template <operators::VolumeEquation... Equations>
+        struct InteriorPlan<equations::EquationRegistry<Equations...>>
+        {
+            using EquationsList = utils::types::TypeList<Equations...>;
+            using Rows = utils::dependencies::KeyList<typename Equations::Test...>;
+            using Coefficients = utils::dependencies::CombinedKeyList<typename OperandTerms<typename operators::EquationPlan<Equations>::Columns>::Type...>;
+        };
+
+        template <typename Keys, typename Excluded>
+        struct RemoveKeys;
+
+        template <typename... Keys, typename Excluded>
+        struct RemoveKeys<utils::dependencies::KeyList<Keys...>, Excluded>
+        {
+            using Type = utils::dependencies::CombinedKeyList<std::conditional_t<utils::dependencies::ContainsKey<Keys, Excluded>, utils::dependencies::KeyList<>, utils::dependencies::KeyList<Keys>>...>;
+        };
+    }
+
+    /** @brief The current interior backend can evaluate every explicitly listed equation. */
+    template <typename Registry>
+    concept ExecutableInteriorEquations = requires
+    {
+        typename detail::InteriorPlan<std::remove_cvref_t<Registry>>::Rows;
+    };
+
+    /**
+     * @brief Derive physical coordinates from explicit equations, independently of numerical solvability.
+     * @note Prescribed coefficients remain equation inputs but are excluded from solved coordinates.
+     * Additional coordinates and rows describe explicit shape or constraint contributions. Neither
+     * equation counting nor this plan claims absence of null spaces or existence of an equilibrium.
+     */
+    template <ExecutableInteriorEquations Registry, utils::dependencies::IsKeyList Prescribed = utils::dependencies::KeyList<>,
+              utils::dependencies::IsKeyList AdditionalValues = utils::dependencies::KeyList<>,
+              utils::dependencies::IsKeyList AdditionalRows = utils::dependencies::KeyList<>>
+    struct ModelPlan final : detail::InteriorPlan<std::remove_cvref_t<Registry>>
+    {
+        using Interior = detail::InteriorPlan<std::remove_cvref_t<Registry>>;
+        using PrescribedCoefficients = Prescribed;
+        using Values = utils::dependencies::CombinedKeyList<typename detail::RemoveKeys<typename Interior::Coefficients, Prescribed>::Type, AdditionalValues>;
+        using Residuals = utils::dependencies::CombinedKeyList<typename Interior::Rows, AdditionalRows>;
+        using Layout = discretization::blocks::SystemLayout<Values, Residuals>;
+    };
+}
diff --git a/src/include/serif/model/specification.hpp b/src/include/serif/model/specification.hpp
new file mode 100644
index 00000000..1a80254a
--- /dev/null
+++ b/src/include/serif/model/specification.hpp
@@ -0,0 +1,155 @@
+#pragma once
+
+#include <array>
+#include <concepts>
+#include <limits>
+#include <string_view>
+#include <type_traits>
+
+#include "serif/dimensions/quantities.hpp"
+#include "serif/utils/types/strings.hpp"
+#include "serif/utils/types/type_list.hpp"
+
+namespace serif::model
+{
+    /** @brief Physical purpose of a prescription, independently of its numerical implementation. */
+    enum class SpecificationRole
+    {
+        constitutive_law, boundary_condition, invariant, phase_condition, gauge_choice, rotation_law
+    };
+
+    /** @brief Distinguish prescribed physics from the physical meaning of a generated scalar. */
+    enum class GeneratedStateKind
+    {
+        none, multiplier, physical_coordinate, solver_border
+    };
+
+    /** @brief Declare independent target, unknown and residual units and stable descriptive identities. */
+    template <dimensions::PhysicalQuantityType Target, dimensions::PhysicalQuantityType Coordinate, dimensions::PhysicalQuantityType Residual,
+              utils::types::FixedString CoordinateId, utils::types::FixedString CoordinateSymbol, utils::types::FixedString ResidualId, utils::types::FixedString ResidualSymbol>
+    struct ScalarConstraint
+    {
+        using TargetQuantity = Target;
+        using CoordinateQuantity = Coordinate;
+        using ResidualQuantity = Residual;
+        using TargetValue = dimensions::QuantityValue<Target>;
+        static constexpr std::string_view coordinate_id{CoordinateId.data(), CoordinateId.size()};
+        static constexpr std::string_view coordinate_symbol{CoordinateSymbol.data(), CoordinateSymbol.size()};
+        static constexpr std::string_view residual_id{ResidualId.data(), ResidualId.size()};
+        static constexpr std::string_view residual_symbol{ResidualSymbol.data(), ResidualSymbol.size()};
+    };
+
+    template <typename... Fields>
+    using Reads = utils::types::TypeList<Fields...>;
+
+    template <typename... Equations>
+    using Changes = utils::types::TypeList<Equations...>;
+
+    /** @brief Declare a prescription once; numerical providers and normalization consume this metadata. */
+    template <typename Owner, utils::types::FixedString Name, SpecificationRole Role, GeneratedStateKind Kind = GeneratedStateKind::none,
+              typename Scalar = void, typename ReadFields = Reads<>, typename ChangedEquations = Changes<>>
+    struct SpecificationDefinition
+    {
+        using SpecificationType = Owner;
+        using ScalarDescription = Scalar;
+        using DependsOn = ReadFields;
+        using Affects = ChangedEquations;
+        static constexpr std::string_view name{Name.data(), Name.size()};
+        static constexpr SpecificationRole role = Role;
+        static constexpr GeneratedStateKind generated_state_kind = Kind;
+    };
+
+    namespace detail
+    {
+        /** @brief Reject unknown enums and mismatched physical roles before inspecting optional scalar metadata. */
+        [[nodiscard]] consteval bool compatible_generated_state(SpecificationRole role, GeneratedStateKind kind)
+        {
+            switch (role)
+            {
+                case SpecificationRole::constitutive_law:
+                case SpecificationRole::boundary_condition:
+                case SpecificationRole::rotation_law:
+                    return kind == GeneratedStateKind::none;
+                case SpecificationRole::invariant:
+                    return kind == GeneratedStateKind::multiplier || kind == GeneratedStateKind::physical_coordinate;
+                case SpecificationRole::phase_condition:
+                case SpecificationRole::gauge_choice:
+                    return kind == GeneratedStateKind::solver_border;
+            }
+            return false;
+        }
+
+        template <typename Candidate>
+        inline constexpr bool is_type_list = false;
+
+        template <typename... Types>
+        inline constexpr bool is_type_list<utils::types::TypeList<Types...>> = true;
+    } // namespace detail
+
+    /** @brief Recognize a self-describing prescription without demanding executable physics or optional metadata. */
+    template <typename Candidate>
+    concept ModelSpecification = requires
+    {
+        typename std::remove_cvref_t<Candidate>::ModelDefinition;
+        requires std::same_as<typename std::remove_cvref_t<Candidate>::ModelDefinition::SpecificationType, std::remove_cvref_t<Candidate>>;
+        { std::remove_cvref_t<Candidate>::ModelDefinition::name } -> std::convertible_to<std::string_view>;
+        requires (!std::remove_cvref_t<Candidate>::ModelDefinition::name.empty());
+        requires detail::compatible_generated_state(std::remove_cvref_t<Candidate>::ModelDefinition::role, std::remove_cvref_t<Candidate>::ModelDefinition::generated_state_kind);
+        requires detail::is_type_list<typename std::remove_cvref_t<Candidate>::ModelDefinition::DependsOn>;
+        requires detail::is_type_list<typename std::remove_cvref_t<Candidate>::ModelDefinition::Affects>;
+    };
+
+    /** @brief Require exact dimensional agreement only when scalar numerical metadata is needed. */
+    template <typename Candidate>
+    concept DimensionalScalarSpecification = ModelSpecification<Candidate> && requires(const std::remove_cvref_t<Candidate>& specification)
+    {
+        requires (std::remove_cvref_t<Candidate>::ModelDefinition::generated_state_kind != GeneratedStateKind::none);
+        typename std::remove_cvref_t<Candidate>::ModelDefinition::ScalarDescription::TargetQuantity;
+        typename std::remove_cvref_t<Candidate>::ModelDefinition::ScalarDescription::CoordinateQuantity;
+        typename std::remove_cvref_t<Candidate>::ModelDefinition::ScalarDescription::ResidualQuantity;
+        requires dimensions::PhysicalQuantityType<typename std::remove_cvref_t<Candidate>::ModelDefinition::ScalarDescription::CoordinateQuantity>;
+        requires dimensions::PhysicalQuantityType<typename std::remove_cvref_t<Candidate>::ModelDefinition::ScalarDescription::ResidualQuantity>;
+        { specification.target() } -> std::same_as<dimensions::QuantityValue<typename std::remove_cvref_t<Candidate>::ModelDefinition::ScalarDescription::TargetQuantity>>;
+    };
+
+    /** @brief Name another prescription, including a forward declaration, without erasing its identity. */
+    template <typename Specification>
+    struct GeneratedCoordinateOf
+    {
+        using SpecificationType = Specification;
+        static constexpr std::size_t scalar_arity = 1;
+        static constexpr bool global_scalar = true;
+    };
+
+    /** @brief Identify a generated residual independently from its target and unknown units. */
+    template <typename Specification>
+    struct ConstraintOf
+    {
+        using SpecificationType = Specification;
+        static constexpr std::size_t scalar_arity = 1;
+        static constexpr bool global_scalar = true;
+    };
+
+    /** @brief Select an explicit numerical enforcement while preserving the public physical prescription. */
+    template <typename Specification>
+    struct EnforcementOf
+    {
+        using SpecificationType = Specification;
+    };
+
+    /** @brief Query scalar units only once the owning prescription is complete. */
+    template <DimensionalScalarSpecification Specification>
+    using GeneratedCoordinateQuantity = typename Specification::ModelDefinition::ScalarDescription::CoordinateQuantity;
+
+    /** @brief Residual units remain independent from the target and generated coordinate. */
+    template <DimensionalScalarSpecification Specification>
+    using ConstraintQuantity = typename Specification::ModelDefinition::ScalarDescription::ResidualQuantity;
+
+    /** @brief Inspect canonical prescription identity without claiming numerical availability. */
+    struct SpecificationDescriptor
+    {
+        std::string_view name;
+        SpecificationRole role;
+        GeneratedStateKind generated_state_kind;
+    };
+} // namespace serif::model
diff --git a/src/include/serif/model/specifications/fixed_angular_momentum.hpp b/src/include/serif/model/specifications/fixed_angular_momentum.hpp
new file mode 100644
index 00000000..00f795fb
--- /dev/null
+++ b/src/include/serif/model/specifications/fixed_angular_momentum.hpp
@@ -0,0 +1,58 @@
+#pragma once
+
+#include <array>
+
+#include "serif/dimensions/type_alias.hpp"
+#include "serif/model/specification.hpp"
+#include "serif/model/specifications/physical_names.hpp"
+
+namespace serif::model
+{
+    /** @brief Prescribe angular momentum while angular velocity remains an independent solved coordinate. */
+    class FixedAngularMomentum final
+    {
+    public:
+        struct Parameters
+        {
+            dimensions::AngularMomentumValue angular_momentum;
+            std::array<double, 3> axis{0.0, 0.0, 1.0};
+            std::array<double, 3> center{};
+        };
+
+        using ScalarDescription = ScalarConstraint<dimensions::AngularMomentum, dimensions::AngularVelocity, dimensions::AngularMomentum,
+                                                   "fixed_angular_momentum.angular_velocity", "Omega", "fixed_angular_momentum.residual", "R_J">;
+        using TargetValue = ScalarDescription::TargetValue;
+        using ModelDefinition = SpecificationDefinition<FixedAngularMomentum, "FixedAngularMomentum", SpecificationRole::invariant, GeneratedStateKind::physical_coordinate,
+                                                        ScalarDescription, Reads<state::Density, state::SurfaceShape, state::OwnGeneratedCoordinate>,
+                                                        Changes<equation::SurfaceShapeBalance, equation::HydrostaticBalance>>;
+
+        /** @brief Require finite nonnegative angular momentum about the default z axis. */
+        explicit FixedAngularMomentum(TargetValue angular_momentum);
+
+        /** @brief Validate the physical target and normalize a finite nonzero rotation axis. */
+        explicit FixedAngularMomentum(Parameters parameters);
+
+        /** @brief Inspect the prescribed total angular momentum. */
+        [[nodiscard]] TargetValue target() const noexcept
+        {
+            return m_target;
+        }
+
+        /** @brief Read the normalized rotation axis. */
+        [[nodiscard]] const std::array<double, 3>& axis() const noexcept
+        {
+            return m_axis;
+        }
+
+        /** @brief Read the point through which the rotation axis passes. */
+        [[nodiscard]] const std::array<double, 3>& center() const noexcept
+        {
+            return m_center;
+        }
+
+    private:
+        TargetValue m_target;
+        std::array<double, 3> m_axis;
+        std::array<double, 3> m_center;
+    };
+} // namespace serif::model
diff --git a/src/include/serif/model/specifications/fixed_total_mass.hpp b/src/include/serif/model/specifications/fixed_total_mass.hpp
new file mode 100644
index 00000000..c886a3c9
--- /dev/null
+++ b/src/include/serif/model/specifications/fixed_total_mass.hpp
@@ -0,0 +1,41 @@
+#pragma once
+
+#include "serif/dimensions/type_alias.hpp"
+#include "serif/model/specification.hpp"
+#include "serif/model/specifications/physical_names.hpp"
+
+namespace serif::model
+{
+    /** @brief Prescribe total stellar mass while the effective surface-potential level remains an independent solved coordinate. */
+    class FixedTotalMass final
+    {
+    public:
+        struct Parameters
+        {
+            dimensions::MassValue total_mass;
+        };
+
+        // Preserve the serialized coordinate identity while describing its present physical meaning.
+        using ScalarDescription = ScalarConstraint<dimensions::Mass, dimensions::GravitationalPotential, dimensions::Mass,
+                                                   "fixed_total_mass.multiplier", "Us", "fixed_total_mass.residual", "R_M">;
+        using TargetValue = ScalarDescription::TargetValue;
+        // Us retains the mass-constraint multiplier role; changing its physical units must not reorder public prescriptions.
+        using ModelDefinition = SpecificationDefinition<FixedTotalMass, "FixedTotalMass", SpecificationRole::invariant, GeneratedStateKind::multiplier,
+                                                        ScalarDescription, Reads<state::Density, state::SurfaceShape>, Changes<equation::SurfaceShapeBalance>>;
+
+        /** @brief Require a positive finite target mass. */
+        explicit FixedTotalMass(TargetValue mass);
+
+        /** @brief Construct from named physical parameters. */
+        explicit FixedTotalMass(Parameters parameters);
+
+        /** @brief Inspect the prescribed total mass. */
+        [[nodiscard]] TargetValue target() const noexcept
+        {
+            return m_target;
+        }
+
+    private:
+        TargetValue m_target;
+    };
+} // namespace serif::model
diff --git a/src/include/serif/model/specifications/physical_names.hpp b/src/include/serif/model/specifications/physical_names.hpp
new file mode 100644
index 00000000..56c2ec47
--- /dev/null
+++ b/src/include/serif/model/specifications/physical_names.hpp
@@ -0,0 +1,36 @@
+#pragma once
+
+namespace serif::model
+{
+    namespace state
+    {
+        struct Density
+        {
+        };
+        struct SurfaceShape
+        {
+        };
+        struct GravitationalPotential
+        {
+        };
+        struct OwnGeneratedCoordinate
+        {
+        };
+    } // namespace state
+
+    namespace equation
+    {
+        struct PoissonBalance
+        {
+        };
+        struct SurfaceShapeBalance
+        {
+        };
+        struct HydrostaticBalance
+        {
+        };
+        struct OwnConstraint
+        {
+        };
+    } // namespace equation
+} // namespace serif::model
diff --git a/src/include/serif/model/state.hpp b/src/include/serif/model/state.hpp
new file mode 100644
index 00000000..183e7ef7
--- /dev/null
+++ b/src/include/serif/model/state.hpp
@@ -0,0 +1,80 @@
+#pragma once
+
+#include <span>
+#include <utility>
+#include <vector>
+
+#include "serif/discretization/blocks/vector_view.hpp"
+
+namespace serif::model {
+    /**
+     * @brief Own one stellar candidate in physical coordinates, including any solved surface shape.
+     * @note Copying a state makes an independent candidate. Its layout fixes the field identities and
+     * sizes; normalization changes solver coordinates, never these physical coefficients. Prescribed
+     * inputs are supplied separately to evaluation because they have no solved coordinate in this layout.
+     */
+    template <discretization::blocks::SystemLayoutType LayoutT>
+    struct StellarState final {
+        using Layout = std::remove_cvref_t<LayoutT>;
+
+        /** @brief Allocate one zero-initialized coefficient vector in the declared physical ordering. */
+        explicit StellarState(Layout layout) : m_layout(std::move(layout)),
+            m_values(m_layout.template total_size<discretization::blocks::RepresentedValueKind::value>()) {
+        }
+
+        /** @brief Adopt existing physical coefficients after checking their complete coordinate count. */
+        StellarState(Layout layout, std::vector<double> values) : m_layout(std::move(layout)), m_values(std::move(values)) {
+            (void)discretization::blocks::StateView<Layout>(m_layout, std::span<double>{m_values});
+        }
+
+        /** @brief Inspect the fixed field and equation ordering associated with this state. */
+        [[nodiscard]] const Layout& layout() const & noexcept {
+            return m_layout;
+        }
+
+        /** @brief Layout references must not outlive a temporary state. */
+        const Layout& layout() const && = delete;
+
+        /** @brief Edit one physical field or scalar without maintaining a second coefficient copy. */
+        template <typename Term> requires utils::dependencies::ContainsKey<Term, typename Layout::Values>
+        [[nodiscard]] std::span<double> get() & {
+            return discretization::blocks::StateView<Layout>(m_layout, std::span<double>{m_values}).template get<Term>();
+        }
+
+        /** @brief Inspect one physical field while preserving the constness of its owning state. */
+        template <typename Term> requires utils::dependencies::ContainsKey<Term, typename Layout::Values>
+        [[nodiscard]] std::span<const double> get() const & {
+            return discretization::blocks::StateView<Layout, const double>(m_layout, std::span<const double>{m_values}).template get<Term>();
+        }
+
+        /** @brief A field view cannot outlive a temporary owning state. */
+        template <typename Term>
+        std::span<const double> get() const && = delete;
+
+        /** @brief Access packed physical coordinates for an explicitly chosen numerical update. */
+        [[nodiscard]] std::span<double> values() & noexcept {
+            return m_values;
+        }
+
+        /** @brief Read packed physical coordinates without permitting a change in their count. */
+        [[nodiscard]] std::span<const double> values() const & noexcept {
+            return m_values;
+        }
+
+        /** @brief Do not return a span into a temporary owning state. */
+        std::span<const double> values() const && = delete;
+
+    private:
+        Layout m_layout;
+        std::vector<double> m_values;
+    };
+
+    /** @brief Infer a state's layout type from its physical coordinate description. */
+    template <discretization::blocks::SystemLayoutType Layout>
+    StellarState(Layout) -> StellarState<std::remove_cvref_t<Layout>>;
+
+    /** @brief Infer a state's layout type when adopting physical coefficients. */
+    template <discretization::blocks::SystemLayoutType Layout>
+    StellarState(Layout, std::vector<double>) -> StellarState<std::remove_cvref_t<Layout>>;
+
+}
diff --git a/src/include/serif/model/stellar_model.hpp b/src/include/serif/model/stellar_model.hpp
new file mode 100644
index 00000000..400ec6f4
--- /dev/null
+++ b/src/include/serif/model/stellar_model.hpp
@@ -0,0 +1,186 @@
+#pragma once
+
+#include <functional>
+#include <span>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+#include "serif/model/detail/specification_list.hpp"
+#include "serif/model/equilibrium_equations.hpp"
+
+namespace serif::model
+{
+    /** @brief Prescriptions accompany one separately owned EOS; no numerical capability is assumed. */
+    template <typename... Specifications>
+    concept ValidSpecificationPack = (ModelSpecification<Specifications> && ...) &&
+        ((Specifications::ModelDefinition::role != SpecificationRole::constitutive_law) && ...) && detail::unique_specification_names<Specifications...>();
+
+    template <typename EquationOfStateT, typename EquationDefinitionT, typename SpecificationList>
+    class StellarModel;
+
+    /** @brief Own immutable physics and an explicit formulation, independently of meshes, caches and solvers. */
+    template <typename EquationOfStateT, typename EquationDefinitionT, typename... Specifications>
+        requires ValidSpecificationPack<Specifications...> && std::is_object_v<EquationOfStateT> && std::is_object_v<EquationDefinitionT>
+    class StellarModel<EquationOfStateT, EquationDefinitionT, utils::types::TypeList<Specifications...>> final
+    {
+    public:
+        using EquationOfState = EquationOfStateT;
+        using Equations = EquationDefinitionT;
+        using SpecificationTypes = utils::types::TypeList<Specifications...>;
+        static constexpr std::size_t specification_count = sizeof...(Specifications) + 1;
+
+        template <SpecificationRole Role>
+        using SpecificationsForRole = std::conditional_t<Role == SpecificationRole::constitutive_law, utils::types::TypeList<EquationOfState>, typename detail::SelectRole<Role, Specifications...>::Type>;
+
+        template <SpecificationRole Role>
+        static constexpr std::size_t specification_role_count = std::size_t{Role == SpecificationRole::constitutive_law} + (std::size_t{0} + ... + std::size_t{Specifications::ModelDefinition::role == Role});
+
+        template <typename Specification>
+        static constexpr bool contains_specification = std::same_as<std::remove_cvref_t<Specification>, EquationOfState> || (std::same_as<std::remove_cvref_t<Specification>, Specifications> || ...);
+
+        /** @brief Capture each physical input once; constructor order does not change specification identity. */
+        template <typename... Arguments>
+            requires detail::ArgumentsMatch<SpecificationTypes, Arguments...>::value && (std::constructible_from<std::remove_cvref_t<Arguments>, Arguments> && ...)
+        explicit StellarModel(EquationOfState equation_of_state, Equations equations, Arguments&&... specifications) :
+            m_equation_of_state(std::move(equation_of_state)), m_equations(std::move(equations)),
+            m_specifications(detail::canonical_specifications(SpecificationTypes{}, std::forward<Arguments>(specifications)...))
+        {
+        }
+
+        /** @brief Inspect the authoritative EOS without requiring its numerical relations. */
+        [[nodiscard]] const EquationOfState& equation_of_state() const noexcept
+        {
+            return m_equation_of_state;
+        }
+
+        /** @brief Inspect the supplied equation composition or its deferred definition. */
+        [[nodiscard]] const Equations& equations() const noexcept
+        {
+            return m_equations;
+        }
+
+        /** @brief Access a physical prescription by exact type, never by a coincidentally shared unit. */
+        template <typename Specification>
+            requires contains_specification<Specification>
+        [[nodiscard]] const std::remove_cvref_t<Specification>& specification() const noexcept
+        {
+            if constexpr (std::same_as<std::remove_cvref_t<Specification>, EquationOfState>)
+            {
+                return m_equation_of_state;
+            }
+            else
+            {
+                return std::get<std::remove_cvref_t<Specification>>(m_specifications);
+            }
+        }
+
+        /** @brief Return immutable references to every prescription of one physical role in canonical order. */
+        template <SpecificationRole Role>
+        [[nodiscard]] auto specifications_for_role() const noexcept
+        {
+            return [this]<typename... Selected>(utils::types::TypeList<Selected...>)
+            {
+                return std::tuple<const Selected&...>{specification<Selected>()...};
+            }(SpecificationsForRole<Role>{});
+        }
+
+        /** @brief Offer a singular physical accessor only when that role is unambiguous. */
+        template <SpecificationRole Role>
+            requires (specification_role_count<Role> == 1)
+        [[nodiscard]] const auto& specification_for_role() const noexcept
+        {
+            return std::get<0>(specifications_for_role<Role>());
+        }
+
+        /** @brief Visit the EOS first, then physical prescriptions in canonical order. */
+        template <typename Visitor>
+            requires std::invocable<Visitor&, const EquationOfState&> && (std::invocable<Visitor&, const Specifications&> && ...)
+        void visit_specifications(Visitor&& visitor) const
+        {
+            std::invoke(visitor, m_equation_of_state);
+            (std::invoke(visitor, specification<Specifications>()), ...);
+        }
+
+        /** @brief Rebuild same-type physical values; deferred equations will read the new values at binding. */
+        template <typename Transform>
+            requires std::copy_constructible<Equations> && std::invocable<Transform&, const EquationOfState&> &&
+                     std::same_as<std::remove_cvref_t<std::invoke_result_t<Transform&, const EquationOfState&>>, EquationOfState> &&
+                     (std::invocable<Transform&, const Specifications&> && ...) &&
+                     (std::same_as<std::remove_cvref_t<std::invoke_result_t<Transform&, const Specifications&>>, Specifications> && ...)
+        [[nodiscard]] StellarModel transform_specifications(Transform&& transform) const
+        {
+            return transform_specifications(std::forward<Transform>(transform), m_equations);
+        }
+
+        /** @brief Rebuild move-only equation definitions explicitly instead of imposing a cloning policy. */
+        template <typename Transform>
+            requires std::invocable<Transform&, const EquationOfState&> &&
+                     std::same_as<std::remove_cvref_t<std::invoke_result_t<Transform&, const EquationOfState&>>, EquationOfState> &&
+                     (std::invocable<Transform&, const Specifications&> && ...) &&
+                     (std::same_as<std::remove_cvref_t<std::invoke_result_t<Transform&, const Specifications&>>, Specifications> && ...)
+        [[nodiscard]] StellarModel transform_specifications(Transform&& transform, Equations equations) const
+        {
+            return StellarModel{std::invoke(transform, m_equation_of_state), std::move(equations), std::invoke(transform, specification<Specifications>())...};
+        }
+
+        /** @brief Describe canonical non-EOS prescriptions without predicting runtime support. */
+        [[nodiscard]] static constexpr std::span<const SpecificationDescriptor> specification_descriptors() noexcept
+        {
+            static constexpr std::array<SpecificationDescriptor, sizeof...(Specifications)> descriptions{{
+                {Specifications::ModelDefinition::name, Specifications::ModelDefinition::role, Specifications::ModelDefinition::generated_state_kind}...
+            }};
+            return descriptions;
+        }
+
+    private:
+        EquationOfState m_equation_of_state;
+        Equations m_equations;
+        std::tuple<Specifications...> m_specifications;
+    };
+
+    template <typename EOS, typename Definition, typename... Specifications>
+        requires ValidSpecificationPack<std::remove_cvref_t<Specifications>...>
+    StellarModel(EOS&&, Definition&&, Specifications&&...) -> StellarModel<std::remove_cvref_t<EOS>, std::remove_cvref_t<Definition>, typename detail::CanonicalSpecifications<std::remove_cvref_t<Specifications>...>::Type>;
+
+    /** @brief Construct a typed physical description with programmer-selected equations. */
+    template <typename EOS, typename Definition, typename... Specifications>
+        requires ValidSpecificationPack<std::remove_cvref_t<Specifications>...>
+    [[nodiscard]] auto make_stellar_model(EOS&& eos, Definition&& definition, Specifications&&... specifications)
+    {
+        return StellarModel{std::forward<EOS>(eos), std::forward<Definition>(definition), std::forward<Specifications>(specifications)...};
+    }
+
+    namespace detail
+    {
+        template <typename Candidate, typename = void>
+        inline constexpr bool is_stellar_model = false;
+
+        template <typename EOS, typename Definition, typename List>
+        inline constexpr bool is_stellar_model<StellarModel<EOS, Definition, List>, std::void_t<typename StellarModel<EOS, Definition, List>::SpecificationTypes>> = true;
+    } // namespace detail
+
+    /** @brief Safely recognize descriptions, including when presented with an incomplete specialization. */
+    template <typename Candidate>
+    concept StellarModelType = detail::is_stellar_model<std::remove_cvref_t<Candidate>>;
+
+    /** @brief Inspect role counts without instantiating invalid model specializations. */
+    template <SpecificationRole Role, typename Candidate>
+    inline constexpr std::size_t specification_role_count = []
+    {
+        if constexpr (StellarModelType<Candidate>)
+        {
+            return std::remove_cvref_t<Candidate>::template specification_role_count<Role>;
+        }
+        else
+        {
+            return std::size_t{0};
+        }
+    }();
+
+    template <SpecificationRole Role, typename Candidate>
+    concept HasUniqueSpecificationForRole = StellarModelType<Candidate> && specification_role_count<Role, Candidate> == 1;
+
+    template <SpecificationRole Role, StellarModelType Candidate>
+    using SpecificationsForRole = typename std::remove_cvref_t<Candidate>::template SpecificationsForRole<Role>;
+} // namespace serif::model
diff --git a/src/include/serif/model/validation.hpp b/src/include/serif/model/validation.hpp
new file mode 100644
index 00000000..f68257a7
--- /dev/null
+++ b/src/include/serif/model/validation.hpp
@@ -0,0 +1,138 @@
+#pragma once
+
+#include "serif/model/description.hpp"
+#include "serif/model/detail/specification_endpoints.hpp"
+#include "serif/model/plan.hpp"
+#include "serif/operators/linearized_system_view.hpp"
+
+namespace serif::model
+{
+    /** @brief Name a physical derivative using the shared block-source key convention. */
+    template <typename Row, typename Field>
+    using Derivative = operators::PhysicalBlock<utils::dependencies::KeyList<Row>, utils::dependencies::KeyList<Field>>;
+
+    /** @brief Map physical names and built-in generated identities once; custom identities require explicit bindings. */
+    template <utils::dependencies::IsBinding... Additional>
+    using SpecificationBindings = utils::dependencies::Bindings<
+        utils::dependencies::Binding<state::Density, discretization::blocks::density::mass>,
+        utils::dependencies::Binding<state::SurfaceShape, discretization::blocks::surface_deformation::coefficients>,
+        utils::dependencies::Binding<state::GravitationalPotential, discretization::blocks::gravity::potential>,
+        utils::dependencies::Binding<equation::PoissonBalance, discretization::blocks::gravity::potential>,
+        utils::dependencies::Binding<equation::SurfaceShapeBalance, discretization::blocks::surface_deformation::shape_equilibrium>,
+        utils::dependencies::Binding<equation::HydrostaticBalance, discretization::blocks::density::mass>,
+        utils::dependencies::Binding<GeneratedCoordinateOf<FixedTotalMass>, discretization::blocks::barotropic_constant::mass_normalization>,
+        utils::dependencies::Binding<ConstraintOf<FixedTotalMass>, discretization::blocks::barotropic_constant::mass_normalization>,
+        utils::dependencies::Binding<GeneratedCoordinateOf<FixedAngularMomentum>, discretization::blocks::rotation::angular_velocity>,
+        utils::dependencies::Binding<ConstraintOf<FixedAngularMomentum>, discretization::blocks::rotation::angular_velocity>, Additional...>;
+
+    using DefaultSpecificationBindings = SpecificationBindings<>;
+
+    namespace detail
+    {
+        template <typename Equation, typename Columns>
+        struct EquationCouplings;
+
+        template <typename Equation, typename... Columns>
+        struct EquationCouplings<Equation, utils::types::TypeList<Columns...>>
+        {
+            using Type = utils::dependencies::KeyList<Derivative<typename Equation::Test, typename Columns::term_type>...>;
+        };
+
+        template <typename Registry>
+        struct InteriorCouplings;
+
+        template <operators::VolumeEquation... Equations>
+        struct InteriorCouplings<equations::EquationRegistry<Equations...>>
+        {
+            using Type = utils::dependencies::CombinedKeyList<typename EquationCouplings<Equations, typename operators::EquationPlan<Equations>::Columns>::Type...>;
+        };
+
+        /** @brief Check physical inputs separately from unknowns and explicitly inactive equations. */
+        template <DimensionalScalarSpecification Specification, typename Model, typename Values, typename Inputs, typename Rows, typename InactiveRows, typename Bindings>
+        [[nodiscard]] consteval bool valid_specification_endpoints()
+        {
+            using namespace utils::dependencies;
+            using Endpoints = ResolvedSpecification<Specification, Bindings>;
+            if constexpr (!ContainsKey<typename Endpoints::Coordinate, Values> || !ContainsKey<typename Endpoints::Constraint, Rows> ||
+                          !Endpoints::Reads::template owners_present<Model> || !Endpoints::Changes::template owners_present<Model>)
+            {
+                return false;
+            }
+            return []<typename... Read, typename... Changed>(KeyList<Read...>, KeyList<Changed...>)
+            {
+                return sizeof...(Read) != 0 && (ContainsKey<Read, Inputs> && ...) && (ContainsKey<Changed, Rows> || ...) &&
+                       ((ContainsKey<Changed, Rows> || ContainsKey<Changed, InactiveRows>) && ...);
+            }(typename Endpoints::Reads::Keys{}, typename Endpoints::Changes::Keys{});
+        }
+
+        /** @brief Fold independent physical declarations; no equation counting or rank inference is involved. */
+        template <typename Model, typename Values, typename Inputs, typename Rows, typename InactiveRows, typename Bindings, typename... Specifications>
+        [[nodiscard]] consteval bool valid_specifications(utils::types::TypeList<Specifications...>)
+        {
+            using namespace utils::dependencies;
+            constexpr bool declarations = (((!GeneratesState<Specifications>::value && std::same_as<typename Specifications::ModelDefinition::DependsOn, Reads<>> &&
+                std::same_as<typename Specifications::ModelDefinition::Affects, Changes<>>) || DimensionalScalarSpecification<Specifications>) && ...);
+            if constexpr (!declarations)
+            {
+                return false;
+            }
+            else
+            {
+                return []<typename... Generated>(utils::types::TypeList<Generated...>)
+                {
+                    using Coordinates = CombinedKeyList<KeyList<typename ResolvedSpecification<Generated, Bindings>::Coordinate>...>;
+                    using Residuals = CombinedKeyList<KeyList<typename ResolvedSpecification<Generated, Bindings>::Constraint>...>;
+                    return key_count_v<Coordinates> == sizeof...(Generated) && key_count_v<Residuals> == sizeof...(Generated) &&
+                           (valid_specification_endpoints<Generated, Model, Values, Inputs, Rows, InactiveRows, Bindings>() && ...);
+                }(GeneratedSpecifications<Model>{});
+            }
+        }
+
+        /** @brief Validate declared derivative contracts without inferring zeros from missing providers. */
+        template <DimensionalScalarSpecification Specification, typename Couplings, typename Bindings>
+        [[nodiscard]] consteval bool specification_derivatives_available()
+        {
+            using namespace utils::dependencies;
+            using Endpoints = ResolvedSpecification<Specification, Bindings>;
+            return []<typename... Read, typename... Changed>(KeyList<Read...>, KeyList<Changed...>)
+            {
+                constexpr auto changed_row = []<typename Row>()
+                {
+                    return ContainsKey<Derivative<Row, typename Endpoints::Coordinate>, Couplings> && (ContainsKey<Derivative<Row, Read>, Couplings> && ...);
+                };
+                return (ContainsKey<Derivative<typename Endpoints::Constraint, Read>, Couplings> && ...) && (changed_row.template operator()<Changed>() && ...);
+            }(typename Endpoints::Reads::Keys{}, typename Endpoints::Changes::Keys{});
+        }
+
+        /** @brief Check each prescription against independently supplied derivative registrations. */
+        template <typename Couplings, typename Bindings, typename... Specifications>
+        [[nodiscard]] consteval bool specification_derivatives_available(utils::types::TypeList<Specifications...>)
+        {
+            return (specification_derivatives_available<Specifications, Couplings, Bindings>() && ...);
+        }
+    } // namespace detail
+
+    /** @brief Derive declared interior couplings from actual operands rather than a second physical sparsity table. */
+    template <ExecutableInteriorEquations Registry>
+    using DeclaredEquationCouplings = typename detail::InteriorCouplings<std::remove_cvref_t<Registry>>::Type;
+
+    /** @brief Validate physical ownership, available inputs and uniquely owned generated equations. */
+    template <typename Model, typename Values, typename Inputs, typename Rows, typename InactiveRows = utils::dependencies::KeyList<>, typename Bindings = DefaultSpecificationBindings>
+    concept SpecificationEndpointsValid = StellarModelType<Model> && utils::dependencies::IsKeyList<Values> && utils::dependencies::IsKeyList<Inputs> &&
+        utils::dependencies::IsKeyList<Rows> && utils::dependencies::IsKeyList<InactiveRows> && utils::dependencies::IsBindings<Bindings> &&
+        detail::valid_specifications<std::remove_cvref_t<Model>, Values, Inputs, Rows, InactiveRows, Bindings>(typename std::remove_cvref_t<Model>::SpecificationTypes{});
+
+    /** @brief Validate physical endpoints and every derivative in an explicitly declared contract. */
+    template <typename Model, typename Values, typename Rows, typename DeclaredCouplings, typename Bindings = DefaultSpecificationBindings>
+    concept SpecificationTopologyValid = SpecificationEndpointsValid<Model, Values, Values, Rows, utils::dependencies::KeyList<>, Bindings> &&
+        utils::dependencies::IsKeyList<DeclaredCouplings> && detail::specification_derivatives_available<DeclaredCouplings, Bindings>(detail::GeneratedSpecifications<std::remove_cvref_t<Model>>{});
+
+    /**
+     * @brief Require a registered implementation for every declared generated derivative, including explicit zeros.
+     * @note ImplementedCouplings must come from admitted providers. Missing methods never imply zero; this checks
+     * specification contributions only, leaving the equation backend's own execution contract independent.
+     */
+    template <typename Model, typename Values, typename Rows, typename DeclaredCouplings, typename ImplementedCouplings, typename Bindings = DefaultSpecificationBindings>
+    concept SpecificationBackendAvailable = SpecificationTopologyValid<Model, Values, Rows, DeclaredCouplings, Bindings> &&
+        SpecificationTopologyValid<Model, Values, Rows, ImplementedCouplings, Bindings>;
+} // namespace serif::model
diff --git a/src/lib/serif/model/specifications.cpp b/src/lib/serif/model/specifications.cpp
new file mode 100644
index 00000000..80d7af35
--- /dev/null
+++ b/src/lib/serif/model/specifications.cpp
@@ -0,0 +1,60 @@
+#include <algorithm>
+#include <cmath>
+#include <stdexcept>
+
+#include "serif/model/specifications/fixed_angular_momentum.hpp"
+#include "serif/model/specifications/fixed_total_mass.hpp"
+#include "serif/utils/error/errors.hpp"
+
+namespace serif::model
+{
+    /** @brief Validate a dimensional scalar target before storing a physical prescription. */
+    static void validate_target(double value, bool zero_allowed, std::string_view name)
+    {
+        using namespace utils::errors;
+        serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(std::isfinite(value) && (zero_allowed ? value >= 0.0 : value > 0.0),
+            "The {} target must be finite and {}; received {}. Supply a physical target in the permitted range.", name, zero_allowed ? "nonnegative" : "positive", value);
+    }
+
+    /** @brief Store a physically admissible total mass. */
+    FixedTotalMass::FixedTotalMass(TargetValue mass) : m_target(mass)
+    {
+        validate_target(mass.value(), false, "total mass");
+    }
+
+    /** @brief Delegate named mass parameters to the dimensional constructor. */
+    FixedTotalMass::FixedTotalMass(Parameters parameters) : FixedTotalMass(parameters.total_mass)
+    {
+    }
+
+    /** @brief Use the default rotation axis for a dimensional angular momentum target. */
+    FixedAngularMomentum::FixedAngularMomentum(TargetValue angular_momentum) : FixedAngularMomentum(Parameters{.angular_momentum = angular_momentum})
+    {
+    }
+
+    /** @brief Preserve arbitrary finite axes and centers without overflowing the normalization calculation. */
+    FixedAngularMomentum::FixedAngularMomentum(Parameters parameters) : m_target(parameters.angular_momentum), m_axis(parameters.axis), m_center(parameters.center)
+    {
+        using namespace utils::errors;
+        validate_target(m_target.value(), true, "angular momentum");
+
+        double maximum_component = 0.0;
+        for (std::size_t component = 0; component < m_axis.size(); ++component)
+        {
+            serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(std::isfinite(m_axis[component]) && std::isfinite(m_center[component]),
+                "Rotation axis and center components must be finite; component {} contains axis {} and center {}. Supply a finite rotation axis and center.", component, m_axis[component], m_center[component]);
+            maximum_component = std::max(maximum_component, std::abs(m_axis[component]));
+        }
+        serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(maximum_component > 0.0, "A rotation axis must be nonzero; supply a direction such as (0, 0, 1).");
+
+        for (double& component : m_axis)
+        {
+            component /= maximum_component;
+        }
+        const double norm = std::hypot(m_axis[0], m_axis[1], m_axis[2]);
+        for (double& component : m_axis)
+        {
+            component /= norm;
+        }
+    }
+} // namespace serif::model