MeanField/docs/stellar_physics_extension_guide.md

# Extending stellar-model physics

MeanField treats a stellar model as a compile-time collection of independent physics specifications. An EOS, a surface condition, an integral constraint, and a phase condition each describe only their own contribution. The library then infers the model type, generated scalar blocks, residual blocks, Jacobian edges, normalization, and manifest entries from that collection.

The important rule is: **add one declaration and, where needed, one implementation per physics specification—never a specialization for a whole combination of constraints.** A new integral must compose with every already-supported EOS, surface, and other constraint without naming those combinations.

Every downstream layer folds the same canonical specification pack. The symbolic compiler appends generated blocks and Jacobian edges, the manifest and normalization layers project metadata for those exact blocks, the prepared equilibrium operator owns one statically typed runtime slot per specification, and the default preconditioner folds the corresponding border actions. “Variadic” therefore means one inferred C++ type and one inferred block system for each pack, not a runtime list and not a library of hand-written pack combinations.

## The small physics-facing contract

A specification is an ordinary physics class with:

- a `Parameters` aggregate;
- a constructor from `Parameters`; and
- one nested `ModelDefinition` made from a role-specific alias.

The declaration is the single source used by `SpecificationTraits` and `SpecificationContribution`; an extension author does not specialize either trait.

### EOS

The EOS still uses the existing relation vocabulary. `ModelDefinition` only tells the stellar-model compiler which role the object has.

```cpp
class LinearBarotrope final {
public:
    struct Parameters final {
        double densityPerEnthalpy;
    };

    using ModelDefinition = mean_field::eos::ConstitutiveLaw<
        LinearBarotrope,
        "LinearBarotrope">;

    using Relations = mean_field::eos::RelationCatalog<
        mean_field::eos::DensityFromSpecificEnthalpy>;

    explicit constexpr LinearBarotrope(Parameters p) noexcept
        : m_densityPerEnthalpy(p.densityPerEnthalpy) {
    }

    [[nodiscard]] constexpr mean_field::dimensions::DensityValue evaluate(
        mean_field::eos::DensityFromSpecificEnthalpy,
        mean_field::dimensions::SpecificEnthalpyValue h
    ) const noexcept {
        return mean_field::dimensions::DensityValue{m_densityPerEnthalpy * h.value()};
    }

    [[nodiscard]] constexpr auto partialDerivative(
        mean_field::eos::DensityFromSpecificEnthalpy,
        mean_field::eos::WithRespectTo<mean_field::dimensions::quantity::SpecificEnthalpy>,
        mean_field::dimensions::SpecificEnthalpyValue
    ) const noexcept {
        return mean_field::eos::PartialDerivative<
            mean_field::dimensions::quantity::Density,
            mean_field::dimensions::quantity::SpecificEnthalpy>{m_densityPerEnthalpy};
    }

private:
    double m_densityPerEnthalpy;
};

static_assert(mean_field::eos::EquationOfStateModel<LinearBarotrope>);
static_assert(mean_field::models::ModelSpecification<LinearBarotrope>);
```

The explicit quantity type is what prevents accidental pressure/density/enthalpy interchange.

### Surface condition

A surface condition declares the physical quantity it fixes. Describing a new condition is intentionally separate from implementing its finite-element carrier equation.

```cpp
class IsothermalSurface final {
public:
    struct Parameters final {
        mean_field::dimensions::TemperatureValue Tsurface;
    };

    using ModelDefinition = mean_field::surface::BoundaryCondition<
        IsothermalSurface,
        "IsothermalSurface">;
    using PhysicalQuantity = mean_field::dimensions::quantity::Temperature;
    using TargetValue = mean_field::dimensions::TemperatureValue;

    explicit constexpr IsothermalSurface(Parameters p) noexcept
        : m_target(p.Tsurface) {
    }

    [[nodiscard]] constexpr TargetValue targetTemperature() const noexcept {
        return m_target;
    }

private:
    TargetValue m_target;
};
```

### Integral constraint

An integral says which stellar state variables it reads (`stellar::Reads`), which equations its generated coordinate changes (`stellar::Changes`), how its generated value and residual are normalized, and how they appear in diagnostics. The ordinary declaration uses physics names; the equilibrium compiler translates them to numerical blocks.

```cpp
class PreparedFixedBaryonMass;

class FixedBaryonMass final {
public:
    struct Parameters final {
        mean_field::dimensions::MassValue Mbaryon;
    };

    using ScalarDescription = mean_field::stellar::ScalarConstraint<
        mean_field::dimensions::quantity::Mass,
        mean_field::dimensions::quantity::SpecificEnergy,
        mean_field::dimensions::quantity::Mass,
        "fixed_baryon_mass.multiplier", "C_b",
        "fixed_baryon_mass.residual", "R_Mb">;
    using TargetValue = typename ScalarDescription::TargetValue;

    using ModelDefinition = mean_field::integral::FixedScalarWithMultiplier<
        FixedBaryonMass,
        "FixedBaryonMass",
        mean_field::stellar::Reads<
            mean_field::stellar::state::Density,
            mean_field::stellar::state::SurfaceShape>,
        mean_field::stellar::Changes<
            mean_field::stellar::equation::HydrostaticBalance>,
        ScalarDescription>;

    using EquilibriumPhysics =
        mean_field::operators::LocalSpecificationEquilibriumPhysics<
            PreparedFixedBaryonMass>;

    explicit FixedBaryonMass(Parameters p) : m_target(p.Mbaryon) {
    }

    // This conventional accessor lets the generic manifest read the target.
    [[nodiscard]] TargetValue target() const noexcept {
        return m_target;
    }

private:
    TargetValue m_target;
};
```

`FixedIntegralWithPhysicalCoordinate` has the same declaration shape, but its generated scalar is a physical state coordinate rather than a Lagrange multiplier. `FixedAngularMomentum`, whose generated coordinate is angular velocity, is the current example.

`stellar::ScalarConstraint` names three independent physical quantities: the
target, the generated Newton coordinate, and the appended constraint residual.
The wrapper derives numerical scale laws and diagnostic unit labels from those
strong quantity types, so a total energy cannot accidentally be normalized as
a specific energy merely because two free-form strings agree. The quantities
need not be the same: `FixedCentralDensity` has a density target and a
specific-enthalpy phase residual. The wrapper contains no assembly logic. The
lower-level `FixedIntegralWithMultiplier` spelling remains available for
unusual compatibility declarations; opting either its normalization or
manifest into typed dimensional metadata requires the complete typed pair. The
current generated-constraint
compiler supports exactly one scalar value and one scalar residual per
generated specification; a distributed or multi-component generated quantity
requires an explicit extension of that compiler boundary and is not obtained by
switching to the lower-level spelling.

If the target and appended residual have different quantity types, the
specification also supplies the physical reference used to report and scale
that residual:

```cpp
[[nodiscard]] mean_field::dimensions::SpecificEnergyValue
residualReference() const noexcept {
    return m_enthalpyReference;
}
```

The return type must be exactly the `QuantityValue` for
`ScalarDescription::ConstraintResidualQuantity`. Omitting this accessor, or
returning the target's quantity instead, makes
`CompilableRootManifestFor<Model, Form>` false. Thus a density target cannot
silently provide the scale for an energy-valued equation. The built-in
`FixedCentralDensity` phase condition is the deliberate special case: its
specific-enthalpy reference is derived from the selected EOS when the
constraint is compiled.

When an integral residual also depends on the coordinate generated by that same specification, use `stellar::state::OwnGeneratedCoordinate`. The compiler resolves that relative name to the specification's exact generated type. Prefer the exported `mean_field::stellar::{state,equation}` names and the `stellar::Reads`, `stellar::Changes`, `stellar::ScalarConstraint`, and `stellar::Derivative` wrappers in ordinary physics code. These wrappers are aliases over the compiler's canonical types: they improve readability without introducing a parallel implementation or weakening any compile-time check.

### Phase condition

A phase condition uses the same vocabulary but generates a solver-border coordinate. It fixes a degeneracy; it is not automatically a new physical invariant.

```cpp
class FixedCentralEnthalpy final {
public:
    struct Parameters final {
        mean_field::dimensions::SpecificEnthalpyValue hcenter;
    };

    using ScalarDescription = mean_field::stellar::ScalarConstraint<
        mean_field::dimensions::quantity::SpecificEnthalpy,
        mean_field::dimensions::quantity::SpecificEnthalpy,
        mean_field::dimensions::quantity::SpecificEnthalpy,
        "fixed_central_enthalpy.border", "lambda_hc",
        "fixed_central_enthalpy.residual", "R_hc">;
    using TargetValue = typename ScalarDescription::TargetValue;

    using ModelDefinition = mean_field::constraint::ScalarPhaseCondition<
        FixedCentralEnthalpy,
        "FixedCentralEnthalpy",
        mean_field::stellar::Reads<
            mean_field::stellar::state::SpecificEnthalpy>,
        mean_field::stellar::Changes<
            mean_field::stellar::equation::HydrostaticBalance>,
        ScalarDescription>;

    explicit FixedCentralEnthalpy(Parameters p) : m_target(p.hcenter) {
    }

    [[nodiscard]] auto target() const noexcept {
        return m_target;
    }

private:
    TargetValue m_target;
};
```

Generated arity is inferred by the role alias: a constitutive law or boundary condition generates no scalar pair, while the three scalar constraint aliases each generate one balanced scalar value/residual pair. This zero-or-one-pair rule is the deliberate current generated-constraint boundary, not a claim that arbitrary-rank new fields are already supported. A malformed normalization or manifest remains inspectable but fails its corresponding completeness query.

### Exact residual and Jacobian providers

A declaration determines the complete provider topology. The nested `EquilibriumPhysics` alias above connects that declaration to its numerical formulas.

Most prepared physics classes keep the simple `Physics(const Specification&)`
constructor. A specification which declares both `state::Density` and
`state::SurfaceShape` in `Reads` may instead request the current physical
density-volume integral through
`stellar::DensityVolumeIntegralContext<Specification>`:

- `integrateDensity(rho)` evaluates the distributed, mapped-volume finite-element
  integral and returns `dimensions::MassValue`;
- `linearizeDensityIntegral(deltaRho)` returns its exact density action; and
- `linearizeSurfaceShapeIntegral(deltaShape)` returns its exact geometry action.

Both declarations are mandatory because `integral rho dV` depends on the
deformed volume. A density-only declaration is rejected before the context type
can be formed, so it cannot hide a missing geometry column in the inferred
Jacobian. A future fixed/reference-geometry integral would be a separate,
explicitly named service rather than a weakened form of this one.
The context is a small, copyable, non-owning handle. It exposes neither the
physical core nor FEM objects, uses no heap type erasure, and delegates global
reduction and reusable scratch storage to the physical core. Its numerical
methods are valid in `PrepareAfterPhysical` and subsequent residual, Jacobian,
or border-action work, after the core has prepared the current state; a physics
constructor should only retain the handle. This service is deliberately scoped
to the current core's density-volume integral rather than being an untyped
general-purpose backend escape hatch.

`PreparedFixedBaryonMass` is still an ordinary concrete class. It does not need
a model template, block list, offset table, or MFEM access. A representative
context-using shape is:

```cpp
class PreparedFixedBaryonMass final {
public:
    using Report = BaryonMassPreparationReport;
    using IntegralContext =
        mean_field::stellar::DensityVolumeIntegralContext<FixedBaryonMass>;

    PreparedFixedBaryonMass(
        const FixedBaryonMass &specification,
        IntegralContext integrals
    ) : m_target(specification.target()), m_integrals(integrals) {
    }

    template <typename State>
    [[nodiscard]] Report PrepareAfterPhysical(const State &state) {
        // Only names from stellar::Reads plus this specification's generated
        // coordinate are visible here.
        const auto density = state.density();
        const auto multiplier = state.generatedCoordinate();
        m_currentMass = m_integrals.integrateDensity(density);
        m_massResidual = m_currentMass.value() - m_target.value();
        prepareHydrostaticTerm(multiplier);
        m_prepared = true;
        return {};
    }

    template <typename Row>
    [[nodiscard]] auto AddResidual(
        mean_field::stellar::equation::OwnConstraint,
        Row &row
    ) const {
        return row.add(m_massResidual);
    }

    template <typename Row>
    [[nodiscard]] auto AddResidual(
        mean_field::stellar::equation::HydrostaticBalance,
        Row &row
    ) const {
        return row.add(m_hydrostaticTerm);
    }

    template <typename Direction, typename Row>
    [[nodiscard]] auto AddJacobianAction(
        mean_field::stellar::Derivative<
            mean_field::stellar::equation::OwnConstraint,
            mean_field::stellar::state::Density>,
        const Direction &direction,
        Row &row
    ) const {
        return row.add(
            m_integrals.linearizeDensityIntegral(direction.density()).value());
    }

    template <typename Direction, typename Row>
    [[nodiscard]] auto AddJacobianAction(
        mean_field::stellar::Derivative<
            mean_field::stellar::equation::OwnConstraint,
            mean_field::stellar::state::SurfaceShape>,
        const Direction &direction,
        Row &row
    ) const {
        return row.add(m_integrals
            .linearizeSurfaceShapeIntegral(direction.surfaceShape()).value());
    }

    // Reads is conservative for all changed equations. This particular
    // derivative is mathematically absent, so say that explicitly.
    template <typename Direction, typename Row>
    [[nodiscard]] mean_field::stellar::StructuralZero AddJacobianAction(
        mean_field::stellar::Derivative<
            mean_field::stellar::equation::HydrostaticBalance,
            mean_field::stellar::state::Density>,
        const Direction &,
        Row &
    ) const noexcept {
        return mean_field::stellar::zeroDerivative;
    }

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

    template <typename Direction, typename Row>
    [[nodiscard]] auto AddJacobianAction(
        mean_field::stellar::Derivative<
            mean_field::stellar::equation::HydrostaticBalance,
            mean_field::stellar::state::OwnGeneratedCoordinate>,
        const Direction &direction,
        Row &row
    ) const {
        return row.add(applyMultiplierDerivative(
            direction.generatedCoordinate()));
    }

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

private:
    // Physics data and helper functions; no solver block bookkeeping.
    FixedBaryonMass::TargetValue m_target;
    IntegralContext m_integrals;
    mean_field::dimensions::MassValue m_currentMass{0.0};
    double m_massResidual{0.0};
    double m_hydrostaticTerm{0.0};
    bool m_prepared{false};
};
```

The hydrostatic helper formula above is a placeholder for the extension's actual physics; the integral and its two directional actions are the real typed service contract. The compiler derives these calls from `Reads` and `Changes`:

- it requests one `AddResidual(EquationTag, Row&)` provider for the specification's own constraint equation and for each distinct equation in `Changes`;
- it requests `OwnConstraint` derivatives with respect to every state in `Reads`; and
- it requests each `Changes` equation's derivatives with respect to every state in `Reads` plus `OwnGeneratedCoordinate`.

The compiler, not the physics class, enumerates and invokes this set. Every Jacobian call receives a direction containing exactly one declared state and a row that can modify exactly one declared equation. Named accessors such as `density()`, `specificEnthalpy()`, and `generatedCoordinate()` are available only when that source is the selected derivative. `value()` is the backend-independent accessor for a one-source extension tag that does not yet have a named convenience accessor.

Every provider returns one of two statements:

- return the token from `row.add(scalar_or_vector)` or, for a residual entry, `row.addEntry(index, value)`; or
- return `stellar::structuralZero` for an intentionally zero residual contribution, and conventionally `stellar::zeroDerivative` for an intentionally zero derivative.

This is stronger than a broad callback that is merely callable. Omitting one inferred overload makes `CompleteStellarEquilibriumPhysicsProvider<Physics, Specification, Model>` false at compile time. Within a callback, adding twice is rejected at runtime. Returning `ContributionAdded` without one addition, or returning `StructuralZero` after an addition, is also rejected immediately. Separate specifications may still add to the same assembled equation; additive composition across specification slots is intentional.

The exact row API exposes no assignment, mutable root vector, block offset, or MFEM alias synchronization. Consequently one specification cannot redirect a derivative to another legal row or erase the physical core or another specification's contribution. The lower-level `SpecificationEquilibriumPhysics` spelling remains available when a concrete implementation genuinely depends on the complete `Model` type, but it uses this same exact adapter and guarantees.

For the current isobaric surface, the pressure condition owns selected rows of the specific-enthalpy equation. Exact residual and Jacobian rows automatically suppress contributions on those replaced rows. A constraint computes its ordinary full hydrostatic contribution; it does not reproduce the surface-row mask.

## Composition and model access

Construction order does not affect the inferred type:

```cpp
auto model = mean_field::model::StellarModel(
    LinearBarotrope({.densityPerEnthalpy = 2.0}),
    IsothermalSurface({.Tsurface = mean_field::dimensions::TemperatureValue{1.0e6}}),
    FixedBaryonMass({.Mbaryon = mean_field::dimensions::MassValue{1.0}}),
    FixedCentralEnthalpy({.hcenter = mean_field::dimensions::SpecificEnthalpyValue{0.25}})
);

using Model = std::remove_cvref_t<decltype(model)>;

static_assert(Model::template specificationRoleCount<
              mean_field::models::SpecificationRole::invariant> == 1);

const auto &eos = model.equationOfState();
const auto &surface = model.surfaceCondition();
const auto &mass = model.specification<FixedBaryonMass>();
```

`equationOfState()` is always available because a valid model has exactly one constitutive law. `surfaceCondition()` and `specificationForRole<Role>()` are available only when that role is unique. Roles such as integral constraints may occur many times; use `SpecificationsForRole<Role>` at compile time and `specification<ExactType>()` at runtime.

Stable names are identities within a role. Two independently developed invariants therefore choose descriptive stable names; they do not coordinate global numeric ordinals. Reusing a stable name in the same role is rejected at compile time, while the same readable name may intentionally appear in two different roles. Manifest lookup always compares the full `(role, stable name)` identity.

### Coupling constraints to one another

Global constraints are not required to be mutually independent. A specification can refer to another specification in the same model using physics-level owner tags:

```cpp
class CoupledVirialCondition final {
public:
    using ModelDefinition = mean_field::integral::FixedScalarWithMultiplier<
        CoupledVirialCondition,
        "CoupledVirialCondition",
        mean_field::stellar::Reads<
            mean_field::stellar::state::Density,
            mean_field::stellar::state::GeneratedCoordinateOf<FixedBaryonMass>>,
        mean_field::stellar::Changes<
            mean_field::stellar::equation::HydrostaticBalance,
            mean_field::stellar::equation::ConstraintOf<FixedBaryonMass>>,
        VirialScalarDescription>;

    // Parameters, target(), and the two local physics packages follow the
    // same pattern as an uncoupled scalar constraint.
};
```

`GeneratedCoordinateOf<FixedBaryonMass>` means “read the scalar coordinate owned by this exact specification type.” `ConstraintOf<FixedBaryonMass>` means “my generated coordinate changes that specification's scalar equation.” `OwnGeneratedCoordinate` and `OwnConstraint` are the corresponding relative names for the declaring specification. These declarations generate ordinary `stellar::Derivative<Equation, State>` calls, so cross-constraint formulas use the same exact residual/Jacobian and preconditioner protocols as every other edge. For a cross-generated direction, the equilibrium view's generic one-source accessor is `value()`; a border-preconditioner direction may use `values()` or `generatedCoordinate<Owner>()`.

The compiler resolves owner references after collecting the complete canonical specification pack, so declaration order is irrelevant. A referenced owner that is absent, does not generate the requested coordinate/residual, or cannot be mapped into the final form makes symbolic compilation fail. It is never converted into a runtime lookup or silently omitted.

## Capability ladder

Description is deliberately easier than executable solver support. Each rung is
an ordinary compile-time capability query, so an unsupported physics component
is rejected at the public overload boundary instead of producing an error deep
inside a constructor or solver body.

For a model declaration, the first two rungs are:

```cpp
static_assert(mean_field::models::SelfDescribingModelSpecification<FixedBaryonMass>);
static_assert(mean_field::models::ModelSpecification<FixedBaryonMass>);

static_assert(mean_field::operators::StellarEquilibriumSpecificationCompilable<FixedBaryonMass>);
static_assert(mean_field::operators::StellarEquilibriumSymbolicallyCompilable<Model>);

static_assert(mean_field::models::CompleteGeneratedScalarDimensionsFor<FixedBaryonMass>);
static_assert(mean_field::models::CompleteGeneratedNormalizationFor<FixedBaryonMass>);
static_assert(mean_field::normalization::CompleteGeneratedPhysicalRieszNormalizationFor<FixedBaryonMass>);
static_assert(mean_field::models::CompleteGeneratedManifestFor<FixedBaryonMass>);

using Form = mean_field::operators::CompiledStellarEquilibriumForm<Model>;
// This example's new surface has no manifest/runtime adapter yet.
static_assert(!mean_field::operators::CompilableRootManifestFor<Model, Form>);

// These become true after the per-specification runtime, EOS core, and surface
// implementations are added. They are intentionally false for this example.
static_assert(!mean_field::operators::StellarEquilibriumPhysicsAvailableFor<
              FixedBaryonMass, Model>);
static_assert(!mean_field::operators::hasCompleteStellarEquilibriumRuntime<Model>);
static_assert(!mean_field::operators::hasStellarEquilibriumCoreRuntime<Model>);
```

`Model::hasCompleteEquilibriumDeclaration` reports only that the model's physics declarations are structurally complete and balanced. It deliberately makes no claim about backend support. The clearly named query for block and Jacobian compilation is `StellarEquilibriumSymbolicallyCompilable<Model>`; `StellarEquilibriumSystemCompilable<Model>` remains its compatibility spelling. Only `CompiledStellarEquilibriumSystem<Model>` exposes an equilibrium-system compilation classification. This distinction lets an astronomer write and inspect a new physical declaration even when one of its physics names has not yet been mapped to a numerical block.

Once a model has an executable core, surface runtime, manifest, and all of its
per-specification residual/Jacobian actions, it satisfies
`equilibrium::StellarEquilibriumModel` and can be discretized into a `Problem`.
The remaining solver-facing capabilities are then checked on that inferred
problem type:

```cpp
// The selected discretization prescription has a valid preparation path.
static_assert(mean_field::normalization::NormalizableStellarEquilibriumProblem<Problem>);

using PhysicalRieszDiscretization = mean_field::equilibrium::StellarDiscretizationFor<
    mean_field::normalization::PhysicalRieszDiagonal<>>;
using PhysicalRieszProblem = mean_field::equilibrium::StellarEquilibriumProblem<
    Model, PhysicalRieszDiscretization>;

// These stronger queries are specifically for the Physical Riesz path.
static_assert(mean_field::normalization::PhysicalRieszStellarEquilibriumProblem<
              PhysicalRieszProblem>);
static_assert(mean_field::normalization::NormalizableStellarEquilibriumProblem<
              PhysicalRieszProblem>);
static_assert(mean_field::normalization::PhysicalRieszStellarEquilibriumCore<
              typename PhysicalRieszProblem::PhysicalCoreType>);

// The EOS-specific material/surface surrogate supports this exact core.
static_assert(mean_field::preconditioning::MaterialSurfacePreconditionerProblem<Problem>);
static_assert(mean_field::preconditioning::StellarStructurePreconditionerProblem<Problem>);

// Every generated border has an action, and the structure block is supported.
static_assert(mean_field::preconditioning::CompleteSpecificationBorderActionsFor<Problem>);
static_assert(mean_field::preconditioning::DefaultStellarPreconditionerAvailableFor<Problem>);
```

The complete ladder is:

1. **Model description:** the object is constructible and has one valid, self-owned definition.
2. **Blocks and Jacobian:** its block lists are valid, generated rows and columns balance, and every declared edge endpoint exists in the final model. The compiler performs two passes, so one constraint may name a block generated by another constraint without depending on declaration order.
3. **Executable equilibrium root:** `StellarEquilibriumCoreRuntime<EOS>` selects the concrete physical `CoreType` for the constitutive law. Surface compilation/runtime, root-manifest construction, and every specification's nested `EquilibriumPhysics` package (or core-authorized library backend entry) are checked separately. For an ordinary nested package, `CompleteStellarEquilibriumPhysicsProvider` proves that every compiler-enumerated residual row and derivative has an exact provider. Each call receives one read-only source and one additive-only row; the current isobaric row-replacement mask is applied by that row. Passing these checks makes one `PreparedVariadicStellarEquilibriumOperator<Model>`; it does not select a hand-written class for the particular constraint pack.
4. **Normalization:** `NormalizableStellarEquilibriumProblem<Problem>` means the prescription stored in the discretization has a valid `prepareNormalization(problem)` operation. `Unnormalized` is the identity path and therefore remains available to a lean alternate core. A Physical Riesz discretization additionally requires `PhysicalRieszStellarEquilibriumCore<Problem::PhysicalCoreType>`, a structural contract for the density, enthalpy, gravity-gradient, gravity-potential, and surface-deformation maps used to assemble the physical Gram diagonals. It also requires complete per-specification physical-scale contributions. A new core can opt in by exposing that physical interface; it does not inherit from or name the Polytrope core.
5. **Material and stellar-structure preconditioner:** `MaterialSurfacePreconditionerProblem<Problem>` checks that the EOS-specific material/surface backend is registered for the problem's exact selected core type. `StellarStructurePreconditionerProblem<Problem>` is the corresponding public capability for the default coupled material/surface-plus-gravity structure block; it currently inherits the same support boundary.
6. **Generated border and full default preconditioner:** `CompleteSpecificationBorderPhysicsProvider` proves the exact border-incident derivatives for one physics action. `CompleteSpecificationBorderActionsFor<Problem>` folds that requirement over the model's specifications. `DefaultStellarPreconditionerAvailableFor<Problem>` is the final conjunction of a discretized equilibrium problem, complete per-specification border actions, and a supported stellar-structure block. `makePreconditioner(problem)` is available only at this rung.

No rung is specialized for a complete EOS/surface/constraint combination. The
model type determines the packs, the compiler folds their declarations, and
each numerical layer folds the corresponding per-component hooks. If a query
fails, implement the missing hook for the individual EOS, surface, integral, or
phase condition; do not weaken the query and do not add a model-pack
specialization.

### Adding another normalization family

A normalization prescription is a compile-time discretization choice, not a
model constraint. Adding one does not add cases to `discretize`, the stellar
operator, or the solver-facing adapter. Its public
`NormalizationCompilation<Prescription, Form>` specialization owns both the
complete coordinate plan and the early runtime-compatibility declaration for
the selected physical core and specification pack. The library independently
validates that the advertised plan covers every inferred value and residual
block exactly once.

Most new policies use only the public problem/manifest interface when they
prepare their diagonal factors. They can reuse
`RuntimePreparedNormalizationCompilation` rather than spelling core or
specification concepts:

```cpp
namespace my_normalization {
class UniformDiagonal final
    : public mean_field::normalization::NormalizationPrescriptionTag {
public:
    UniformDiagonal(double stateFactor, double residualFactor);

    [[nodiscard]] double stateFactor() const noexcept;
    [[nodiscard]] double residualFactor() const noexcept;
};
} // namespace my_normalization

namespace mean_field::normalization {
template <typename Form>
    requires mean_field::utils::blocks::block_form_is_valid_v<Form>
struct NormalizationCompilation<my_normalization::UniformDiagonal, Form>
    : RuntimePreparedNormalizationCompilation<
          my_normalization::UniformDiagonal,
          Form> { };
} // namespace mean_field::normalization
```

The policy's namespace then provides
`prepareStellarNormalization(const UniformDiagonal&, const Problem&) ->
DiagonalNormalization`. It is found from the exact policy argument carried by
`Problem::NormalizationPrescriptionType`; no central registry or built-in-policy
switch is edited. The generated plan labels every block with
`RuntimePreparedCoordinate<UniformDiagonal>` rather than falsely claiming an
identity or Physical Riesz formula, and the returned diagonal is dimension-checked
against the complete inferred problem. `discretize` rejects a policy
whose compilation omits `runtimeAvailableFor`, and
also rejects a third-party runtime policy whose plan does not name that exact
policy on every coordinate. Thus merely providing a preparation function cannot
bypass or contradict the compile-time plan. The same `discretize` boundary
rejects a policy whose preparation operation is absent or returns the wrong
type, before a solver-facing problem is constructed. A policy that requires private facilities
of a particular physical core writes its own `runtimeAvailableFor<Core,
Specifications>` predicate in the same compilation specialization; Physical
Riesz is the reference example. Ordinary EOS, surface, and constraint authors
do not use this extension point—their dimensional scalar wrappers continue to
generate the appropriate per-block normalization metadata automatically.

## Radial seeds and border preconditioning

The optional numerical hooks also compose per specification. A specification that intentionally leaves a spherical seed unchanged can say so without touching the seed compiler:

```cpp
using RadialProjection = mean_field::seed::projection::Use<
    mean_field::seed::projection::NoStateChange>;
```

A specification with a generated coordinate supplies a similarly wrapped physics rule. The rule receives `RadialProjectionScales`, `RadialProjectionState`, and only its own generated coordinate. It does not calculate block offsets or inspect the other constraints. `RadialProfileProjectableModel<Model>` is true only when every specification has a complete rule and exactly one rule supplies the radial mass scale. Thus a newly introduced EOS, surface condition, or constraint cannot be advertised as projectable and then fail deep inside the function body.

For an ordinary generated constraint, no second Jacobian implementation is
needed for border preconditioning. The default adapter projects the same exact
`AddJacobianAction` providers from `EquilibriumPhysics` onto the subset whose
row or source is a generated border. This is why the `FixedBaryonMass`
declaration above needs only its `EquilibriumPhysics` alias: the residual,
Jacobian, and default border preconditioner have one source of numerical
formulas.

An advanced implementation may deliberately provide a different border
approximation. In that case it opts in explicitly:

```cpp
using SpecificationBorderPhysics =
    mean_field::preconditioning::LocalSpecificationBorderPhysics<MyPreparedBorderAction>;
```

The action uses the same astronomy-facing `stellar::Derivative<Equation, State>` vocabulary as the equilibrium Jacobian. For example, a phase condition declared with `Reads<SpecificEnthalpy>` and `Changes<HydrostaticBalance>` has two border-incident derivatives:

```cpp
class MyPreparedBorderAction {
public:
    explicit MyPreparedBorderAction(const PreparedMyPhaseCondition &physics);

    template <typename Direction, typename Row>
    [[nodiscard]] auto ApplyJacobianAction(
        mean_field::stellar::Derivative<
            mean_field::stellar::equation::OwnConstraint,
            mean_field::stellar::state::SpecificEnthalpy>,
        const Direction &direction,
        Row &row
    ) const {
        return row.add(
            evaluatePhaseDerivative(direction.specificEnthalpy()));
    }

    template <typename Direction, typename Row>
    [[nodiscard]] auto ApplyJacobianAction(
        mean_field::stellar::Derivative<
            mean_field::stellar::equation::HydrostaticBalance,
            mean_field::stellar::state::OwnGeneratedCoordinate>,
        const Direction &direction,
        Row &row
    ) const {
        return row.add(
            evaluateMultiplierAction(direction.generatedCoordinate()(0)));
    }
};
```

The framework derives the full derivative list from `Reads` and `Changes`, then partitions it into structure-to-border, border-to-structure, and border-to-border edges. It invokes only the border-incident edges. Cross-constraint edges are partitioned in exactly the same way, including edges between two generated borders.

A purely structure-to-structure derivative is checked separately. The default
structure preconditioner accepts it only when the selected physical core names
that exact `(specification, coupling)` as a backend capability or the nested
equilibrium provider returns `stellar::StructuralZero`. A custom nonzero term
on even an otherwise familiar base edge therefore makes
`DefaultStellarPreconditionerAvailableFor<Problem>` false instead of being
silently dropped.

Each direction contains exactly one source and each row can update exactly one destination. For this example, asking the first direction for `density()` or the second row to write a gravity-potential equation is ill-formed. An exact-zero border derivative returns `stellar::zeroDerivative`. Missing an incident overload makes `CompleteSpecificationBorderPhysicsProvider<Action, Specification, Problem>` false at compile time; `PreparedSpecificationBorderPhysicsActionFor` additionally checks construction and move requirements. The single-use `row.add(...)` token has the same duplicate and result-consistency checks as equilibrium assembly.

For an explicit approximation, the action constructor receives only the
prepared equilibrium physics for its own specification. It can snapshot
state-dependent coefficients without access to the enclosing problem, physical
core, or other constraints. `LocalSpecificationBorderPhysics<Action>` is the
common wrapper for a concrete local class. `SpecificationBorderPhysics<ActionTemplate>`
is available when the action genuinely needs the inferred `Problem` type; both
names use the same adapter and neither reimplements pack traversal.

The library handles problem binding, operation partitioning, offsets, alias synchronization, surface-row ownership, and aggregation. The preconditioner lifecycle includes the prepared-operator generation, so changing state-dependent contribution data invalidates and rebuilds the border Schur complement even when the caller reuses the same dependency stamps.

These contracts establish structural completeness and safe composition; they do not establish preconditioner quality. The current default structure solve and generated-border Schur approximation are provisional numerical choices. A problem satisfying `DefaultStellarPreconditionerAvailableFor` has every required action and a registered backend, but it may still be poorly conditioned or ineffective for a new regime. Iteration counts, scaled spectra, parameter continuation behavior, and nonlinear convergence must be measured separately for each important physical family.

## What works end to end today

As of the current checkpoint:

- arbitrary well-formed self-described specifications can form distinct canonical `StellarModel` types, subject to the model-pack role and balance contracts;
- generated scalar blocks, residual rows, and Jacobian edges are inferred from `stellar::Reads` and `stellar::Changes`, including named edges between constraints;
- the residual/Jacobian runtime and generated-border preconditioner enumerate those inferred rows and edges and require exact providers or explicit structural zeros;
- generated Physical Riesz declarations and generated manifest metadata are projected generically, while Physical Riesz field diagonals are available to any selected core satisfying `PhysicalRieszStellarEquilibriumCore`;
- the variadic outer layers currently wrap a fixed five-field barotropic physical core: density, surface deformation, gravity gradient, gravity potential, and specific enthalpy;
- each generated constraint currently adds exactly one scalar value block and one scalar residual block around that core;
- the executable equilibrium kernels are registered for `eos::Polytrope` with `surface::Isobaric`;
- the default material/surface and stellar-structure preconditioner kernels are likewise registered for that Polytrope/Isobaric backend and its exact physical core;
- per-specification runtime and prepared border actions exist for `FixedTotalMass`, `FixedAngularMomentum`, and `FixedCentralDensity`; and
- `makePreconditioner(problem)` is exposed only when `DefaultStellarPreconditionerAvailableFor<Problem>` proves both the structure backend and every generated border action are present.

The five-field statement is a real boundary, not merely a missing convenience alias. A scalar constraint over those existing fields can compose without changing the core. A genuinely new distributed field—magnetic induction, temperature, composition, or another PDE state—requires extending the shared physical-field vocabulary and mapping, the selected executable core and discretization, its residual/Jacobian implementation, normalization metadata, manifest, and the applicable structure preconditioner. Supplying a backend block type in `Reads` does not by itself make that field part of the physical core. This work is still per new field/equation family, not per complete combination of constraints.

The nonbarotropic EOS and isothermal-surface mocks in the compile-contract tests are symbolic coverage. They prove that EOS choice, surface choice, and constraint packs produce distinct canonical types and correctly inferred block/Jacobian structures; they are not executable stellar models and do not claim a compatible five-field core or default preconditioner. Their rejection by `equilibrium::StellarEquilibriumModel` and `discretize` is intentional.

A new EOS or surface can therefore be described and type-checked today, but `discretize` will reject it until its per-EOS core or per-surface compilation/runtime support (including carrier/manifest handling) is registered. A new integral or phase condition can obtain its scalar block/Jacobian pair, normalization, and ordinary `target()` manifest structure from `ModelDefinition`, but still needs its residual/Jacobian runtime. The ordinary constraint places that concrete implementation on itself as `EquilibriumPhysics = operators::LocalSpecificationEquilibriumPhysics<MyPreparedPhysics>`; its exact providers are reused automatically by the default generated-border preconditioner. Only a deliberately distinct border approximation adds `SpecificationBorderPhysics`. Neither path specializes a combination-level registry.

The aggregate `StellarEquilibriumRuntimeContribution` registry is not an alternative astronomy-facing API. It can access the complete model, FEM objects, mapper, dependency stamps, and physical core, so an entry is selected only when the exact specification also appears in the selected concrete core type's `BackendSpecifications` list. A raw specialization—or a custom EOS that merely reuses the built-in core type—cannot add itself to that core's list and is therefore ignored. A custom core author may authorize its own backend entries in that single list; if an authorized entry and a nested `EquilibriumPhysics` package both exist for the same specification, the model is rejected as genuinely ambiguous. Ordinary EOS, surface, and constraint authors should use the restricted nested packages and never manage this list.

An alternate executable core may choose the identity `Unnormalized` path without implementing Physical Riesz maps. Selecting `PhysicalRieszDiagonal` for that same core fails `PhysicalRieszStellarEquilibriumProblem` and `NormalizableStellarEquilibriumProblem` cleanly. Likewise, an executable problem whose EOS/core pair has no material-surface backend fails `MaterialSurfacePreconditionerProblem`, while a problem with one missing constraint action fails `CompleteSpecificationBorderActionsFor`; neither failure requires nor invites a combination-specific fallback. These are availability and structural-correctness claims, not evidence that the current preconditioner will deliver acceptable Krylov or nonlinear convergence.

Those remaining hooks are per physics component. Adding `FixedMagneticEnergy` over existing state fields should require work for `FixedMagneticEnergy` itself. Adding a genuinely new magnetic field/equation additionally requires the shared field-family work listed above—not new implementations for `Polytrope + Isobaric + FixedTotalMass + FixedMagneticEnergy`, then another for every other constraint combination. A converged rotating solution and its nonlinear solver remain downstream goals; the compile-time architecture and provisional preconditioner do not, by themselves, establish that convergence result.