\documentclass[11pt]{article}

\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage[margin=0.82in]{geometry}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{booktabs}
\usepackage{enumitem}
\usepackage{xcolor}
\usepackage{listings}
\usepackage{tabularx}
\usepackage{longtable}
\usepackage{fancyhdr}
\usepackage{microtype}
\usepackage{hyperref}

\definecolor{MeanFieldBlue}{HTML}{174A67}
\definecolor{MeanFieldLightBlue}{HTML}{EAF3F8}
\definecolor{MeanFieldGold}{HTML}{A66500}
\definecolor{CodeBackground}{HTML}{F5F7F8}
\definecolor{CodeComment}{HTML}{476B45}
\definecolor{CodeKeyword}{HTML}{6B327A}

\hypersetup{
    colorlinks=true,
    linkcolor=MeanFieldBlue,
    urlcolor=MeanFieldBlue,
    pdftitle={MeanField Physics Developer Manual},
    pdfauthor={MeanField extension example}
}

\pagestyle{fancy}
\fancyhf{}
\lhead{MeanField Physics Developer Manual}
\rhead{Extension Example}
\cfoot{\thepage}
\setlength{\headheight}{14pt}

\setlist[itemize]{leftmargin=1.35em,itemsep=0.25em,topsep=0.35em}
\setlist[enumerate]{leftmargin=1.5em,itemsep=0.35em,topsep=0.35em}

\lstdefinestyle{meanfieldcpp}{
    language=C++,
    basicstyle=\ttfamily\small,
    keywordstyle=\color{CodeKeyword}\bfseries,
    commentstyle=\color{CodeComment},
    stringstyle=\color{MeanFieldGold},
    backgroundcolor=\color{CodeBackground},
    frame=single,
    rulecolor=\color{MeanFieldBlue!35},
    numbers=left,
    numberstyle=\tiny\color{black!55},
    numbersep=7pt,
    breaklines=true,
    breakatwhitespace=true,
    columns=fullflexible,
    keepspaces=true,
    showstringspaces=false,
    tabsize=4
}

\lstdefinestyle{meanfieldshell}{
    basicstyle=\ttfamily\small,
    backgroundcolor=\color{CodeBackground},
    frame=single,
    rulecolor=\color{MeanFieldBlue!35},
    breaklines=true,
    columns=fullflexible,
    keepspaces=true,
    showstringspaces=false
}

\newcommand{\MeanField}{\textnormal{\textsc{MeanField}}}
\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\dd}{\mathrm{d}}
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}

\newenvironment{developerbox}
{\begin{center}\begin{minipage}{0.94\textwidth}\color{MeanFieldBlue}\hrule\vspace{0.55em}\color{black}}
{\vspace{0.55em}\color{MeanFieldBlue}\hrule\end{minipage}\end{center}}

\title{\vspace{-1.5em}\color{MeanFieldBlue}\Huge\bfseries
       Developing Physics Extensions for \MeanField\\[0.4em]
       \Large A worked ideal-gas plus radiation equation of state}
\author{Physics developer example module}
\date{Code checkpoint \code{71423d543f1775d3b943d9c0361f49edd56b59c1}\\
      Manual built 2026-09-06}

\begin{document}

\maketitle

\begin{developerbox}
\textbf{What this manual promises.}
You will implement a real two-variable equation of state, make its physical
relations and Jacobian derivatives compile-time contracts, and compose it with
the existing isobaric surface, fixed-total-mass invariant, and
fixed-angular-momentum invariant. You do not need to write block indices or
repeat the model assembly logic.

\medskip
\textbf{What this manual does not claim.}
The current numerical stellar-equilibrium core is barotropic. The worked EOS
depends independently on density and temperature. It therefore forms a valid
typed stellar-model specification, but it cannot yet be passed to the current
\code{discretize} operation. A temperature or entropy field and a governing
thermal equation must be added to the equilibrium formulation first. The
example encodes this boundary as a compile-time assertion.
\end{developerbox}

\tableofcontents
\clearpage

\section{Audience and objective}

This manual is for an astronomer or physicist who knows the equations they want
to add but may not be interested in the template machinery used to assemble a
nonlinear root system. The central design rule is:

\begin{quote}
State the physics once, in physics vocabulary. Let the library infer the
structural type of the model, residual, Jacobian, normalization plan, and
preconditioner wherever a runtime implementation exists. If a required piece
does not exist, reject the operation at compile time.
\end{quote}

The worked extension lives entirely under \code{extension\_example/}. It does
not modify a source file in \code{libmeanfield/}. This separation matters. It
demonstrates the public extension surface rather than relying on privileged
access to the implementation.

At the end, you should be able to answer four questions for a proposed physics
extension:

\begin{enumerate}
    \item What quantities and relations does the new physics provide?
    \item Which analytic derivatives does a Newton Jacobian require?
    \item Which existing specifications can be composed with it immediately?
    \item Is the complete numerical runtime available, or is only the symbolic
          model declaration available?
\end{enumerate}

\section{The mental model: five layers}

It helps to separate five ideas that are often blended together in scientific
software.

\begin{longtable}{p{0.16\textwidth} p{0.31\textwidth} p{0.43\textwidth}}
\toprule
\textbf{Layer} & \textbf{Physics question} & \textbf{MeanField representation} \\
\midrule
\endhead
Quantity & What does this scalar mean? & A type such as
\code{Density}, \code{Temperature}, or \code{Pressure}, wrapped in a
\code{QuantityValue}. \\
\addlinespace
Relation & Which variables determine which result? & A compile-time sentence
such as $P=P(\rho,T)$, represented by \code{eos::Relation<Pressure, Density,
Temperature>}. \\
\addlinespace
Specification & What physical choice did the model make? & A concrete EOS,
surface condition, invariant, phase condition, gauge, or rotation law with a
small \code{ModelDefinition}. \\
\addlinespace
Compiled system & Which unknowns, residuals, and Jacobian couplings follow from
the selected specifications? & A type inferred from the complete, variadic
specification pack. \\
\addlinespace
Runtime provider & How are those residuals and derivatives evaluated on a
mesh? & Numerical assembly, normalization, and preconditioning code that must
exist for the compiled type. \\
\bottomrule
\end{longtable}

A model may be valid at one layer and intentionally unavailable at a later
layer. In this example, $P(\rho,T)$ is a valid EOS relation and the four chosen
specifications form a valid stellar-model type. The present runtime provider
has a barotropic material core, so it rejects the thermal model before a
discretization object can be formed.

This distinction prevents two dangerous outcomes:

\begin{itemize}
    \item silently dropping a supplied physical constraint; and
    \item pretending a multidimensional EOS is barotropic by choosing an
          undocumented path through thermodynamic state space.
\end{itemize}

\section{Map of the example module}

\begin{tabularx}{\textwidth}{p{0.37\textwidth} X}
\toprule
\textbf{File} & \textbf{Purpose} \\
\midrule
\code{ideal\_gas\_radiation.cppm} & Declares relations, implements the EOS,
analytic derivatives, input validation, and compile-time protocol checks. \\
\code{rotating\_stellar\_model.cppm} & Provides the small physics-facing
composition wrapper and documents the current runtime boundary. \\
\code{demo.cpp} & Evaluates one CGS thermodynamic state and constructs the
typed stellar-model specification. \\
\code{tests/ideal\_gas\_radiation.cpp} & Tests dimensions, thermodynamic
identities, scaling laws, derivative accuracy, and domain handling. \\
\code{tests/rotating\_stellar\_}\newline\code{model.cpp} & Tests inferred roles, exact
specification types, values, and capability rejection. \\
\code{manual/physics\_developer\_}\newline\code{manual.tex} & The source of this manual. \\
\bottomrule
\end{tabularx}

The module has its own library, demonstration executable, and test executable.
It links to \MeanField, but its files are not added to the \code{mean\_field}
library target and its tests are not added to the main regression executable.

\section{Worked physics: ideal gas plus radiation}

\subsection{Assumptions}

The example describes a deliberately simple local thermodynamic model:

\begin{itemize}
    \item matter is a classical, nondegenerate, monatomic ideal gas;
    \item the mean molecular weight $\mu$ is fixed;
    \item gas and radiation share a single temperature $T$;
    \item radiation is in local thermodynamic equilibrium;
    \item ionization, composition evolution, pairs, degeneracy, Coulomb
          corrections, and opacity do not enter the EOS.
\end{itemize}

These assumptions are appropriate for demonstrating the extension protocol.
They are not a universal stellar EOS.

Let $k_{\mathrm B}$ be Boltzmann's constant, $m_{\mathrm u}$ the atomic mass
unit, $a$ the radiation energy-density constant, and

\begin{equation}
    \mathcal{R} = \frac{k_{\mathrm B}}{\mu m_{\mathrm u}}
\end{equation}

the gas constant per unit mass. The pressure contributions are

\begin{align}
    P_{\mathrm{gas}} &= \rho \mathcal{R} T, \\
    P_{\mathrm{rad}} &= \frac{aT^4}{3}, \\
    P(\rho,T) &= P_{\mathrm{gas}} + P_{\mathrm{rad}}.
\end{align}

For a monatomic gas, the specific internal energies are

\begin{align}
    u_{\mathrm{gas}} &= \frac{3}{2}\mathcal{R}T, \\
    u_{\mathrm{rad}} &= \frac{aT^4}{\rho}, \\
    u(\rho,T) &= u_{\mathrm{gas}} + u_{\mathrm{rad}}.
\end{align}

Using $h=u+P/\rho$, the specific enthalpy is

\begin{equation}
    h(\rho,T)
    = \frac{5}{2}\mathcal{R}T
    + \frac{4aT^4}{3\rho}.
\end{equation}

\subsection{Analytic derivatives}

A Newton method needs derivatives of the residual with respect to its state.
For this EOS, the useful local derivatives are

\begin{align}
\left.\pd{P}{\rho}\right|_T
    &= \mathcal{R}T, &
\left.\pd{P}{T}\right|_\rho
    &= \rho\mathcal{R}+\frac{4aT^3}{3}, \\
\left.\pd{u}{\rho}\right|_T
    &= -\frac{aT^4}{\rho^2}, &
\left.\pd{u}{T}\right|_\rho
    &= \frac{3}{2}\mathcal{R}+\frac{4aT^3}{\rho}, \\
\left.\pd{h}{\rho}\right|_T
    &= -\frac{4aT^4}{3\rho^2}, &
\left.\pd{h}{T}\right|_\rho
    &= \frac{5}{2}\mathcal{R}+\frac{16aT^3}{3\rho}.
\end{align}

The subscript is not optional notation. It tells the developer what is held
fixed and corresponds directly to the \code{WithRespectTo<...>} tag in the
code. A derivative with the right numerical return type but the wrong held
variable is a physics bug.

\subsection{Domain and units}

The implemented material domain is

\begin{equation}
    \rho > 0, \qquad T \geq 0.
\end{equation}

Positive density is required because the specific radiation energy and
enthalpy contain $1/\rho$. The radiation constant may be set to zero for a
pure-gas verification problem. Other constants and $\mu$ must be finite and
positive.

\MeanField's \code{QuantityValue} gives semantic type safety. It does not
perform dimensional algebra or unit conversion. The defaults in this example
are CGS:

\begin{center}
\begin{tabular}{lll}
\toprule
Constant & Default & CGS unit \\
\midrule
$k_{\mathrm B}$ & $1.380649\times10^{-16}$ & erg K$^{-1}$ \\
$m_{\mathrm u}$ & $1.66053906660\times10^{-24}$ & g \\
$a$ & $7.5657\times10^{-15}$ & erg cm$^{-3}$ K$^{-4}$ \\
\bottomrule
\end{tabular}
\end{center}

If the rest of a model uses nondimensional or code units, all three constants,
density, temperature, and returned thermodynamic values must be transformed
coherently. A typed value prevents passing temperature where density belongs;
it cannot detect a CGS value mixed with a nondimensional one.

\section{Implementing the EOS}

\subsection{Step 1: name the relations}

The first code says what the EOS knows, without saying how it computes it.

\begin{lstlisting}[style=meanfieldcpp]
namespace q = mean_field::dimensions::quantity;

using PressureFromDensityAndTemperature = mean_field::eos::Relation<
    q::Pressure,
    q::Density,
    q::Temperature>;

using SpecificEnthalpyFromDensityAndTemperature = mean_field::eos::Relation<
    q::SpecificEnthalpy,
    q::Density,
    q::Temperature>;
\end{lstlisting}

The first template argument is the output. Remaining arguments are inputs in
call order. Consequently,

\begin{lstlisting}[style=meanfieldcpp]
eos::evaluate<q::Pressure>(eosModel, density, temperature);
\end{lstlisting}

is valid, while the same call with \code{temperature, density} is rejected at
compile time.

Use an existing quantity type when it has the intended physical meaning.
Temperature and specific internal energy already exist in
\code{mean\_field::dimensions::quantity}. If a genuinely new quantity is
needed, it can derive from \code{ThermodynamicQuantity} and provide a stable
identifier, but adding a quantity to numerical normalization also requires an
appropriate scale law and runtime support.

\subsection{Step 2: identify the specification role}

Inside the EOS class, one alias connects the physics type to the stellar-model
front end:

\begin{lstlisting}[style=meanfieldcpp]
using ModelDefinition = mean_field::eos::ConstitutiveLaw<
    IdealGasRadiation,
    "IdealGasRadiation">;
\end{lstlisting}

The stable name is used in compile-time keys and runtime descriptions. Choose
a specific, durable name. Do not encode parameter values in it, and do not
reuse the same role-name pair for unrelated physics.

This alias is the wrapper intended for physics developers. It projects to the
more general model-definition machinery, but the EOS author does not need to
name generated blocks, residual rows, or a normalization topology because a
constitutive law owns no scalar solver border by itself.

\subsection{Step 3: publish the complete relation catalog}

\begin{lstlisting}[style=meanfieldcpp]
using Relations = mean_field::eos::RelationCatalog<
    PressureFromDensityAndTemperature,
    SpecificInternalEnergyFromDensityAndTemperature,
    SpecificEnthalpyFromDensityAndTemperature>;
\end{lstlisting}

Treat this catalog as an API promise. Every listed relation must have an exact
\code{evaluate} overload. A relation that exists as a helper function but is
missing from the catalog is not visible to generic consumers. A relation in
the catalog without a matching evaluator makes \code{EquationOfStateModel}
false.

\subsection{Step 4: implement evaluations in physics notation}

The pressure overload is small because named component functions expose useful
diagnostics:

\begin{lstlisting}[style=meanfieldcpp]
[[nodiscard]] dimensions::PressureValue evaluate(
    PressureFromDensityAndTemperature,
    dimensions::DensityValue density,
    dimensions::TemperatureValue temperature
) const {
    return pressureContributions(density, temperature).total();
}
\end{lstlisting}

The relation tag is an empty compile-time object. It selects the overload but
does not cost storage. Return the exact typed value requested by the relation.
Returning a raw \code{double}, even with the correct number, fails the EOS
concept.

The generic user-facing call is

\begin{lstlisting}[style=meanfieldcpp]
const auto pressure = mean_field::eos::evaluate<
    mean_field::dimensions::quantity::Pressure>(
        equationOfState,
        density,
        temperature
    );
\end{lstlisting}

This call infers the relation from output type and typed input sequence. Client
code does not construct relation tags directly.

\subsection{Step 5: implement Jacobian derivatives}

Each derivative overload repeats the relation and states the differentiation
variable explicitly:

\begin{lstlisting}[style=meanfieldcpp]
[[nodiscard]] eos::PartialDerivative<q::Pressure, q::Temperature>
partialDerivative(
    PressureFromDensityAndTemperature,
    eos::WithRespectTo<q::Temperature>,
    dimensions::DensityValue density,
    dimensions::TemperatureValue temperature
) const {
    const double T = temperature.value();
    return eos::PartialDerivative<q::Pressure, q::Temperature>{
        density.value() * specificGasConstant() +
        4.0 * parameters().radiationConstant * T * T * T / 3.0
    };
}
\end{lstlisting}

The important contracts are:

\begin{itemize}
    \item relation inputs occur in the same order as in \code{evaluate};
    \item \code{WithRespectTo<Temperature>} identifies the independent
          variable;
    \item the return type is exactly
          \code{PartialDerivative<Pressure, Temperature>}; and
    \item the formula is a partial derivative at fixed density.
\end{itemize}

\subsection{Step 6: make the compiler check the public claim}

The source finishes with assertions such as

\begin{lstlisting}[style=meanfieldcpp]
static_assert(eos::EquationOfStateModel<IdealGasRadiation>);

static_assert(eos::SupportsPartialDerivative<
    IdealGasRadiation,
    PressureFromDensityAndTemperature,
    q::Temperature>);
\end{lstlisting}

These assertions are not substitutes for numerical tests. They prove shape:
the relation is declared, the argument order matches, and the exact output type
exists. Numerical tests must still prove that the derivative computes the
right mathematics.

\section{Composing a rotating stellar model}

The model itself is the direct expression of the physical choices:

\begin{lstlisting}[style=meanfieldcpp]
auto model = mean_field::model::StellarModel(
    IdealGasRadiation(eosParameters),
    mean_field::surface::Isobaric({.Psurf = surfacePressure}),
    mean_field::integral::FixedTotalMass({.Mtotal = totalMass}),
    mean_field::integral::FixedAngularMomentum({
        .Jtotal = totalAngularMomentum,
        .axis = rotationAxis,
        .center = rotationCenter
    })
);
\end{lstlisting}

The example wraps this expression in \code{makeRotatingStellarModel}. The
wrapper groups frequently used inputs; it does not duplicate model assembly.
An advanced caller can always use \code{StellarModel(...)} directly.

The compiler infers four distinct specification types and their roles:

\begin{center}
\begin{tabularx}{0.94\textwidth}{p{0.28\textwidth} p{0.21\textwidth} X}
\toprule
Specification & Role & Structural contribution \\
\midrule
\code{IdealGasRadiation} & constitutive law & EOS relations; no generated
scalar coordinate \\
\code{Isobaric} & boundary condition & fixed surface pressure; no generated
scalar coordinate \\
\code{FixedTotalMass} & invariant & mass constraint residual and its scalar
multiplier \\
\code{FixedAngularMomentum} & invariant & angular-momentum constraint residual
and angular velocity as a physical coordinate \\
\bottomrule
\end{tabularx}
\end{center}

The two invariants each contribute one scalar unknown and one scalar residual.
The complete symbolic operator remains square. Reordering the constructor
arguments does not define a new physical combination: the library canonicalizes
the specification set by stable compile-time keys while storing exactly one
object of each inferred type.

\begin{samepage}
The model exposes physics-facing accessors:

\begin{lstlisting}[style=meanfieldcpp]
model.equationOfState();
model.surfaceCondition();
model.specification<integral::FixedTotalMass>();
model.specification<integral::FixedAngularMomentum>();
\end{lstlisting}

These preserve exact types. There is no base-class downcast and no string
lookup in the inner numerical path.
\end{samepage}

\section{What compiles today, and why}

\subsection{The three useful questions}

For a new model, ask these separately:

\begin{enumerate}
    \item \textbf{Is each object a valid specification?} The type has a valid
          role definition, parameters, and contribution metadata.
    \item \textbf{Can the objects form a stellar-model type?} Their stable keys
          are unique, exactly one constitutive law is present, and generated
          scalar arities make a square symbolic system.
    \item \textbf{Can the current numerical backend discretize that exact
          type?} Every core equation, surface conversion, residual assembly,
          Jacobian coupling, normalization action, and runtime provider exists.
\end{enumerate}

For this example, the answers are yes, yes, and no.

\subsection{Why the current core rejects this EOS}

The current stellar material core is barotropic. Its state can close density
through a relation of the form

\begin{equation}
    \rho = \rho(h),
\end{equation}

with the corresponding derivative $\dd\rho/\dd h$. This is sufficient for a
polytropic barotrope.

The example instead supplies

\begin{equation}
    P=P(\rho,T), \qquad u=u(\rho,T), \qquad h=h(\rho,T).
\end{equation}

Hydrostatic balance plus an isobaric surface does not determine both $\rho$
and $T$ throughout the star. A second physical equation is required. Depending
on the intended model, it might be an entropy prescription, radiative energy
transport, convective closure, or a full energy equation with sources and
sinks.

The example therefore includes

\begin{lstlisting}[style=meanfieldcpp]
static_assert(model::StellarModelType<RotatingStellarModel>);
static_assert(!equilibrium::StellarEquilibriumModel<RotatingStellarModel>);
\end{lstlisting}

The second assertion is a capability tripwire. When a thermal equilibrium core
is eventually implemented, the assertion and this manual must be revised
together.

\subsection{What must be added for a coupled thermal solve}

A physically complete extension will need all of the following, designed as a
single compile-time path rather than a special case for this EOS:

\begin{enumerate}
    \item A selected thermal state coordinate, such as temperature or specific
          entropy, with a finite-element field family.
    \item A thermal residual equation whose physical closure is explicit.
    \item EOS relations and partial derivatives required by both hydrostatic
          and thermal residuals.
    \item Surface data appropriate to the thermal equation, such as a fixed
          temperature, luminosity, or atmosphere matching condition.
    \item Generated block and Jacobian couplings inferred from the selected
          thermodynamic equation set.
    \item Normalization metadata for the new field and residual, including a
          physically meaningful Riesz map and scale.
    \item A preconditioner component that approximates the new thermal block
          and its important coupling to structure.
    \item Runtime providers for assembly and linearization of every inferred
          contribution.
\end{enumerate}

The EOS in this module is already expressed in the multi-input relation
vocabulary needed by such a future core. It should not be rewritten as a
bespoke EOS-backend combination.

\section{Testing a physics extension}

A robust extension test suite has four layers. Passing only the first layer is
not enough.

\subsection{Compile-time protocol tests}

These answer structural questions:

\begin{itemize}
    \item Is the class a self-describing constitutive-law specification?
    \item Does every catalog relation have an exact evaluator?
    \item Are required partial derivatives available?
    \item Are input order and output quantity types enforced?
    \item Does composition preserve the exact EOS, surface, and invariant
          types?
    \item Is the inferred scalar root system square?
\end{itemize}

The reversed pressure call $(T,\rho)$ is tested as a non-callable expression.
This is a useful negative contract: it proves that strong quantity types catch
an error that two raw \code{double} arguments would not catch.

\subsection{Exact worked state}

The tests choose artificial constants that make hand calculation easy:

\begin{equation}
    \mu=2, \quad k_{\mathrm B}=12, \quad m_{\mathrm u}=3,
    \quad a=9,
\end{equation}

so $\mathcal{R}=2$. At $\rho=4$ and $T=2$,

\begin{align}
    P_{\mathrm{gas}} &= 16, & P_{\mathrm{rad}} &= 48, & P &= 64,\\
    u_{\mathrm{gas}} &= 6, & u_{\mathrm{rad}} &= 36, & u &= 42,\\
    h &= u + P/\rho = 42 + 16 = 58.
\end{align}

This single state catches wrong factors of $1/3$, $3/2$, $4/3$, and $5/2$.
It also tests the independent thermodynamic identity $h=u+P/\rho$.

\subsection{Scaling and limiting behavior}

The tests verify transformations more diagnostic than a list of isolated
numbers:

\begin{itemize}
    \item doubling $\rho$ doubles $P_{\mathrm{gas}}$;
    \item changing $\rho$ leaves $P_{\mathrm{rad}}$ unchanged;
    \item doubling $T$ doubles $P_{\mathrm{gas}}$;
    \item doubling $T$ multiplies $P_{\mathrm{rad}}$ by 16; and
    \item at
          $T_{\mathrm{cross}}=(3\rho\mathcal{R}/a)^{1/3}$, gas and radiation
          pressures agree.
\end{itemize}

These claims exercise the physical exponents and parameter dependencies. They
would fail even if a few reference values happened to be fitted correctly.

\subsection{Analytic versus numerical derivatives}

Every one of the six analytic partial derivatives is compared with a centered
difference at three materially different states. For a scalar function $f$,

\begin{equation}
    f'(x) \approx \frac{f(x+\epsilon)-f(x-\epsilon)}{2\epsilon}.
\end{equation}

The test uses

\begin{equation}
    \epsilon = \epsilon_{\mathrm{machine}}^{1/3}
               \max(1,|x|).
\end{equation}

For a centered difference, this scale balances truncation and roundoff for
smooth double-precision functions. The comparison tolerance is selected from
the observed error order and is much tighter than a solver tolerance. It tests
the local Jacobian implementation, not convergence of a nonlinear solve.

When extending a tabulated or iterative EOS, use tolerances justified by that
algorithm's interpolation and inner-solve error. Do not copy the analytic-EOS
tolerance mechanically.

\subsection{Domain tests}

Tests require rejection of invalid constants, nonfinite values, zero density,
and negative temperature. Domain tests should check error categories where a
caller can recover differently from a nonfinite input and an out-of-domain
state.

\section{Adding other kinds of physics}

The same pattern extends beyond an EOS: declare a physics-facing wrapper, state
dependencies and effects once, and let the compiler generate structure. The
details below are a map, not a replacement for tests of the actual equations.

\subsection{Surface conditions}

A surface condition uses

\begin{lstlisting}[style=meanfieldcpp]
using ModelDefinition = mean_field::surface::BoundaryCondition<
    MySurfaceCondition,
    "MySurfaceCondition">;
\end{lstlisting}

It should name its physical target quantity and typed target value. A new
surface condition can enter a \code{StellarModel} as a distinct type before a
runtime surface compiler exists. Full discretization is available only if the
chosen EOS can convert the boundary statement into the thermodynamic variable
used by the selected material formulation.

For example, an isothermal surface is physically meaningful for a thermal
model. It does not become meaningful for a barotropic formulation merely
because both types compile as specifications.

\subsection{Integral invariants}

An invariant such as fixed mass or fixed angular momentum owns a global scalar
equation and usually a conjugate scalar coordinate. Physics-facing aliases
include \code{FixedScalarWithMultiplier} and
\code{FixedScalarWithPhysicalCoordinate}. Their declarations identify:

\begin{itemize}
    \item the target physical quantity;
    \item the generated coordinate quantity;
    \item the constraint residual quantity;
    \item which state quantities the integral reads;
    \item which equations the generated coordinate changes;
    \item stable symbols and identifiers; and
    \item normalization metadata.
\end{itemize}

That information is enough to generate the border shape and required Jacobian
incidence. A runtime provider must still implement the integral, its derivative,
and its action on affected equations. If that provider is absent, compilation
must stop at the operation boundary rather than ignore the invariant.

\subsection{Phase conditions}

A phase condition removes a neutral direction or selects one representative of
a family. It is not necessarily a conserved physical quantity. Fixed central
density, for example, can serve as a phase constraint while fixed total mass is
the physically informed invariant.

Use \code{ScalarPhaseCondition} for the physics-facing declaration. Keep the
distinction explicit in naming, documentation, and tests: an invariant changes
the physical problem, while a phase condition selects a solution within a
degenerate representation.

\section{Normalization and preconditioning}

Normalization and preconditioning solve different problems.

The normalization operator gives each field and residual a physically
meaningful inner product and scale. Conceptually, it maps the raw root system
to coordinates where comparisons such as residual norms and line-search
acceptance are not dominated merely by units. A pressure-like row near
$10^{16}$ and a dimensionless row near unity should not be compared as raw
numbers.

The preconditioner approximates the inverse action of the scaled Jacobian. It
is concerned with coupling and spectral difficulty, not just magnitude. A good
normalization can improve the numerical setting in which the preconditioner
operates, but it cannot replace an approximation to elliptic, material,
surface, or global-border couplings.

For a new generated scalar, the specification must provide normalization
metadata that is dimensionally coherent with its target, coordinate, and
residual quantities. For a new distributed thermal field, the discretization
must select a compatible topology and physical scale. If no scale law exists,
the correct outcome is a compile-time rejection of normalized discretization,
not an identity scale inserted without documentation.

Physics tests for a new normalization should include:

\begin{itemize}
    \item invariance under the intended change of physical units;
    \item expected mesh-refinement behavior of the discrete norm;
    \item exact results for constant or low-order fields when available;
    \item consistent scaling of a residual and its Jacobian action; and
    \item quantified interaction with the chosen preconditioner.
\end{itemize}

\section{A practical development workflow}

\begin{enumerate}
    \item \textbf{Write the equations first.} List state variables, parameters,
          outputs, domains, held-fixed variables, and expected limits.
    \item \textbf{Reuse or define quantity types.} Avoid raw scalars at public
          physics boundaries.
    \item \textbf{Declare the smallest honest relation catalog.} Do not promise
          an inverse relation that is multivalued or requires an undocumented
          closure.
    \item \textbf{Implement evaluations and analytic derivatives together.}
          Keeping them adjacent makes sign and factor audits easier.
    \item \textbf{Add compile-time assertions.} Check the EOS concept, each
          required derivative, the exact composed model type, and negative
          capability cases.
    \item \textbf{Add hand-calculated physics tests.} Use states that expose
          every coefficient and term.
    \item \textbf{Add property tests.} Check scaling laws, limits, monotonicity,
          conservation identities, symmetry, or convexity as appropriate.
    \item \textbf{Compare derivatives numerically.} Cover the relevant state
          domain, not just one comfortable point.
    \item \textbf{Compose with several existing specifications.} Confirm that
          the model infers roles and preserves exact types without a bespoke
          combination.
    \item \textbf{Ask the runtime capability question explicitly.} If false,
          document which physical equation or provider is absent.
    \item \textbf{Build only the extension target while iterating.} Run the
          wider regression suite only when shared library code changes.
\end{enumerate}

\section{Build and run}

From the repository root, use the configured build directory and request only
the extension targets:

\begin{lstlisting}[style=meanfieldshell]
cmake --build cmake-build-profile-homebrew-llvm \
  --target extension_example_demo extension_example_tests
\end{lstlisting}

Run only the extension tests:

\begin{lstlisting}[style=meanfieldshell]
./cmake-build-profile-homebrew-llvm/extension_example/extension_example_tests
\end{lstlisting}

Run the demonstration:

\begin{lstlisting}[style=meanfieldshell]
./cmake-build-profile-homebrew-llvm/extension_example/extension_example_demo
\end{lstlisting}

Compile this manual into the build directory:

\begin{lstlisting}[style=meanfieldshell]
cmake --build cmake-build-profile-homebrew-llvm \
  --target extension_example_manual
\end{lstlisting}

The extension tests use a separate Catch2 executable. Running it does not run
the main \MeanField{} test suite.

\section{Gotchas}

\begin{enumerate}
    \item \textbf{A typed scalar is not a unit library.} A
          \code{TemperatureValue} says ``temperature,'' not ``kelvin.'' Keep
          the whole parameter set in one coherent unit system.
    \item \textbf{Relation input order is part of the type.}
          $P(\rho,T)$ and $P(T,\rho)$ are different relation types even if a
          human could reorder the arguments.
    \item \textbf{The catalog is authoritative.} An overload omitted from
          \code{Relations} is invisible to generic dispatch. A catalog entry
          without an exact overload invalidates the EOS concept.
    \item \textbf{State what is held fixed.} A total derivative along an
          isentrope is not the same as a partial derivative at fixed
          temperature.
    \item \textbf{Do not invent an inverse closure.} A two-variable EOS does
          not generally supply $\rho(h)$ without another thermodynamic
          condition.
    \item \textbf{Specification success is not runtime success.} A valid
          \code{StellarModel} type can be rejected by \code{discretize} when a
          required numerical provider is absent.
    \item \textbf{A surface statement must match the formulation.} Fixed
          pressure alone does not set both density and temperature at a
          thermal surface.
    \item \textbf{Stable names are identities.} Keep them unique within a role
          and stable across harmless refactors.
    \item \textbf{An invariant is not a phase condition.} Fixed total mass is
          physically informative; fixed central density may be used to select
          a phase. Their generated coordinates have different meanings.
    \item \textbf{Analytic formulas still need numerical audits.} Templates can
          prove that $\partial P/\partial T$ exists, not that a factor of four
          is correct.
    \item \textbf{Avoid EOS-specific solver branches.} Extend generic relation,
          equation, normalization, and provider mechanisms so future EOS types
          follow the same path.
    \item \textbf{Validate before powers or division.} Nonfinite inputs and
          invalid density should fail with an EOS-domain error rather than
          propagate a NaN through a global residual.
\end{enumerate}

\section{Review checklists}

\subsection{EOS author checklist}

\begin{itemize}
    \item[$\square$] Assumptions and validity domain are written down.
    \item[$\square$] Parameters have physical names and coherent units.
    \item[$\square$] Public inputs and outputs use quantity types.
    \item[$\square$] The relation catalog is minimal and complete.
    \item[$\square$] Every catalog relation has an exact evaluator.
    \item[$\square$] Every required Jacobian partial has the correct held-fixed
          variable and return type.
    \item[$\square$] Nonfinite and out-of-domain states are rejected.
    \item[$\square$] Hand-worked values test all coefficients.
    \item[$\square$] Physical scaling laws and limiting regimes are tested.
    \item[$\square$] Analytic derivatives agree with numerical derivatives.
\end{itemize}

\subsection{Coupled-model checklist}

\begin{itemize}
    \item[$\square$] Exactly one constitutive law is present.
    \item[$\square$] Surface and EOS relations are thermodynamically compatible.
    \item[$\square$] Every supplied invariant and phase condition appears in
          the inferred model type.
    \item[$\square$] Generated unknown and residual arities are square.
    \item[$\square$] All Jacobian incidences required by declared dependencies
          and effects are present.
    \item[$\square$] Every new coordinate and residual has normalization
          metadata and a runtime action.
    \item[$\square$] The preconditioner includes the important new couplings.
    \item[$\square$] Runtime capability is asserted before discretization.
    \item[$\square$] Failure of an unsupported combination occurs at compile
          time with a useful boundary and message.
\end{itemize}

\clearpage
\section{Glossary}

\begin{longtable}{p{0.24\textwidth} p{0.69\textwidth}}
\toprule
\textbf{Term} & \textbf{Meaning in this code base} \\
\midrule
\endhead
Barotropic EOS & An EOS whose thermodynamic closure can be expressed with one
independent material coordinate for the current problem, such as
$\rho=\rho(h)$. \\
\addlinespace
Boundary condition & A physical statement applied at a domain boundary, such
as fixed surface pressure. It is a model specification role. \\
\addlinespace
Canonical specification set & The order-independent compile-time collection of
the exact specification types supplied to \code{StellarModel}. \\
\addlinespace
Compile-time contract & A property checked by C++ type formation or a concept,
before a numerical run begins. \\
\addlinespace
Constitutive law & A specification that relates material state quantities. An
EOS is the principal example. \\
\addlinespace
Discretization & The mesh, finite-element spaces, ordering, normalization
prescription, and runtime data needed to turn a compiled physical system into
a finite-dimensional root problem. \\
\addlinespace
Equation of state (EOS) & A set of thermodynamic relations with a declared
validity domain and parameters. \\
\addlinespace
Generated coordinate & A scalar unknown introduced by a specification, such as
the mass-constraint multiplier or angular velocity associated with fixed
angular momentum. \\
\addlinespace
Invariant & A physically prescribed global quantity enforced by a residual,
such as total mass or total angular momentum. \\
\addlinespace
Jacobian & The derivative of the full residual vector with respect to the full
state vector. Its block structure is inferred from compiled equations and
specification contributions. \\
\addlinespace
Model specification & One exact physical choice carrying a role, stable key,
parameters, and structural contribution metadata. \\
\addlinespace
Normalization & A physically informed mapping that supplies comparable
coordinates and residual norms across fields with different units, topologies,
and magnitudes. \\
\addlinespace
Partial derivative & A derivative with respect to one declared relation input
while its other inputs are held fixed. \\
\addlinespace
Phase condition & A condition that removes a neutral direction or chooses one
representative from a family of equivalent solutions. It is distinct from a
physical invariant. \\
\addlinespace
Physical Riesz map & The discrete map induced by a selected physical inner
product and scale. It connects a field or residual with its dual coordinate in
a topology-aware way. \\
\addlinespace
Preconditioner & An efficient approximation to the inverse scaled Jacobian,
used to make the linearized solve tractable. \\
\addlinespace
Quantity type & A zero-storage type that names physical meaning, such as
\code{Temperature}. A \code{QuantityValue<Temperature>} stores the scalar. \\
\addlinespace
Relation catalog & The complete compile-time list of input-output relations an
EOS claims to implement. \\
\addlinespace
Residual & One equation written as a quantity that should be zero at a
solution. The full nonlinear problem is a vector of residual blocks. \\
\addlinespace
Runtime provider & Numerical code that evaluates or assembles an inferred
physics contribution for a particular formulation. \\
\addlinespace
Solver border & The small set of generated global scalar rows and columns
coupled to the distributed physical core. \\
\addlinespace
Specification role & The category of a physical choice: constitutive law,
boundary condition, invariant, phase condition, gauge choice, or rotation law. \\
\addlinespace
Stable name & A compile-time string used with a role to identify one
specification mechanism consistently across inferred structures and runtime
descriptions. \\
\addlinespace
Stellar model & The strongly typed, canonical composition of the exact physical
specification objects selected by a user. \\
\addlinespace
Symbolically square & The compiled state and residual descriptions have equal
total scalar arity before mesh-dependent sizes are known. \\
\addlinespace
Thermal closure & The additional equation or prescription needed to determine
an independent thermal variable such as temperature or entropy. \\
\addlinespace
Typed value & A scalar wrapper whose C++ type carries its physical meaning and
prevents accidental interchange of unrelated quantities. \\
\bottomrule
\end{longtable}

\section{Final perspective}

The extension mechanism is successful when a physics developer can read the
new source primarily as equations, parameters, domains, and derivatives. The
compiler should then answer structural questions: whether the relations are
complete, whether the selected model is square, which exact constraints are
present, and whether the current numerical backend implements every inferred
piece.

The ideal-gas plus radiation example intentionally reaches that boundary. It
demonstrates a valid nonbarotropic EOS and a valid rotating stellar-model
specification without claiming a nonexistent thermal equilibrium solve. That
is the behavior future extensions should preserve: composable when supported,
precisely rejected when incomplete, and always honest about the physics.

\end{document}