MeanField_SERiF/docs/static/normalization.md

# Stellar normalization

Normalization changes the coordinates used by the numerical solver. It does not
change the descriptive stellar model, its EOS, or the physical values supplied to
equation assembly. The default configuration estimates a reference star with the
shared dimensionless Lane–Emden solver and prepares dimensionless Riesz coordinates.

## Public configuration

Include `serif/backend/mfem/model/equilibrium.hpp`. Given a descriptive
`stellar_model`, existing `finite_elements` from `make_stellar_discretization`
(optionally `with_deformation`), and the usual physical `coefficients` and `geometry`:

```cpp
namespace stellar = serif::model::backend::mfem;
namespace numerical = serif::solver::normalization;

auto discretization = stellar::make_discretization(stellar::DiscretizationOptions{
    .model = stellar_model,
    .finite_elements = finite_elements,
    .normalization = numerical::make_normalization()
});
auto evaluator = stellar::bind_stellar_model(discretization, coefficients, geometry);
auto result = evaluator.evaluate(coefficients, geometry);
if (result) {
    auto system = *result;
    // system->residual(), system->jacobian(), system->coordinates()
}
```

`StellarModel` describes the physics. `StellarState` owns one candidate's physical
coefficients. `EquilibriumEvaluator` evaluates candidates and returns owning,
immutable results. A successful evaluation does not accept a nonlinear step;
the solver makes that decision and retains whichever state and result it needs.

For a moving star, one state contains both fields and independent surface
coefficients. Given its physical `layout` and initialized packed coefficients:

```cpp
serif::model::StellarState state(layout, physical_coefficients);
auto evaluator = stellar::bind_stellar_model(discretization, state);
auto result = evaluator.evaluate(state);
```

Use `state.get<Term>()` for typed field access or `state.values()` for a packed
physical update. Copying the state creates an independent trial candidate.
Prescribed fields remain separate inputs, with no solved coordinate. Fixed-geometry
evaluation additionally takes the prescribed displacement. External field views
remain supported through the first example; both paths use the same computations.

`DiscretizationOptions{...}` and `NormalizationOptions{...}` are deduced aggregate
types. The explicit aggregate names let C++ infer their template arguments. A bare
`{...}` passed to a function template cannot do that deduction. The model,
reference strategy, coordinate convention, and explicit plan retain their types;
there is no runtime strategy registry.

The default is `n = 3`, the model's `FixedTotalMass` target, and `G = 1` in the
existing code-unit convention. No model prescribes central density, so supply its
reference value explicitly; the mass may also be given explicitly when the model has no
`FixedTotalMass`. Missing values are a compile-time error. To override the reference without changing the physical model:

```cpp
using namespace serif::dimensions;
auto normalization = numerical::make_normalization(numerical::NormalizationOptions{
    .scheme = numerical::LaneEmdenApproximation{
        .n = 3.0,
        .mass = MassValue{2.0},
        .central_density = DensityValue{0.9}
    },
    .coordinates = numerical::Dimensionless{},
    .gravitational_constant = QuantityValue<GravitationalConstant>{1.0}
});
```

Either dimensional input can be omitted independently when the corresponding
model target exists. Integration controls belong to `scheme.integration`.
The coordinate choice is a type rather than a boolean so the selected plan is checked
at compilation. `IdentityCoordinate` remains available in explicit plans.

## Reference star and units

The Lane–Emden solution supplies its first zero `xi1` and the dimensionless mass
`mu1 = -xi1^2 theta'(xi1)`. With a reference mass `M` and central density `rho_c`,

\[
a = \left(\frac{M}{4\pi\rho_c\mu_1}\right)^{1/3},\qquad R = a\xi_1.
\]

The index and mass alone do not determine the radius. In particular, an `n = 3`
polytrope still needs a density or length scale. The normalization reference is
an estimate; it need not have the same polytropic index as the model's EOS.

The resulting scales include

\[
\rho_0=M/R^3,\quad h_0=GM/R,\quad P_0=GM^2/R^4,\quad
g_0=GM/R^2,\quad \Omega_0=\sqrt{GM/R^3},\quad F_0=GM^2/R^2.
\]

Here `rho_0` is a characteristic density, not the supplied central density.
The implied time scale is `T0 = 1 / Omega0`. Physical `G` must be supplied in the
same units as the model, mesh, EOS, and constraints. It is not inferred from mass
and radius. In SI or CGS, supply the corresponding numerical value of `G` and
consistent values everywhere else. These quantity types enforce dimensions, not
automatic SI/CGS conversion. The reference `G` must agree with the physical
Poisson equation to interpret its dimensionless value as one; generic equation
registries are not inspected to infer or replace their gravitational constant.

Invalid reference inputs, integration limits that do not reach a finite surface,
and nonrepresentable scales return preparation failures. For example, `n = 5`
has a valid dimensionless profile but cannot furnish a finite stellar radius.

## From equations to solver coordinates

The automatic plan reads state quantities from field metadata and residual
quantities from the equations. Their units need not be the same. Each reference
Gram diagonal is assembled once using the field's existing FE representation,
support, and ordering. RT flux coordinates retain the native Piola transform.

For a volume Gram diagonal entry `D`, a field scale `q`, and reference radius `R`,
the dimensionless factors are

\[
S_i=\frac{\sqrt{D_i}}{q R^{3/2}},\qquad
L_i=\frac{1}{q R^{3/2}\sqrt{D_i}}.
\]

The residual `q` describes the equation density before weak integration. Global
scalars use their own quantity scale without a volume measure. Strong surface
replacement rows use the selected boundary formulation's point quantity, which
may be density or pressure. Surface displacement and generalized force use the
reference surface metric and their appropriate primal and dual geometric factors.

With these diagonal transformations, the solver receives

\[
\widehat u=Su,\qquad \widehat r=Lr,\qquad
\widehat J=LJS^{-1}.
\]

`coordinates().denormalize_state()` also converts a normalized Newton correction
back to a physical correction. All transformations are linear and have no offset.
The existing scaled preconditioner adapter uses the same frozen factors and
passes the dimensional Jacobian to the physical inverse.

## Evaluation, ownership, and extension

The reference configuration is resolved once when constructing the configured
discretization. Numerical evaluation is deferred until `evaluate`. Reference
scales and metrics share the ordinary dependency graph and remain independent of
candidate field coefficients and surface displacement. Reference mesh revisions
invalidate the metrics through existing MFEM revision tracking.

A successful `evaluate()` returns an expected value containing a shared pointer to
the matching physical system, coordinate factors, normalized residual, and
normalized Jacobian. Keep that pointer through a linear solve. It survives later
evaluations, rejection, and destruction of the evaluator. A failed physical or
normalization computation leaves previously returned results intact and makes
the current cache unavailable. `system->physical()` exposes the dimensional
system for diagnostics and preconditioning. Evaluation never changes the state.
Operator scratch storage is reused, so one operator instance is not intended for
concurrent applications.

A configured discretization can also select an ordinary preconditioner recipe:

```cpp
auto configuration = make_discretization(DiscretizationOptions{
    .model = stellar_model,
    .finite_elements = finite_elements,
    .normalization = normalization,
    .preconditioner = preconditioner::make_primal_preconditioner(state.layout())
});
auto evaluator = bind_stellar_model(std::move(configuration), state);
auto evaluated = evaluator.evaluate(state);
auto linear_system = evaluator.prepare_linear_system(state);
```

The readable composed defaults live in
`solver/preconditioner/default_preconditioners.hpp`. The primal recipe covers
density, potential, optional surface motion, and all generated scalar constraints. It uses the existing composition factories; custom recipes remain
available through the same configuration member. Omitting the member keeps ordinary
evaluation available without adding inverse computations.

`evaluate()` requests only the evaluated system. `prepare_linear_system()` requests
its matching inverse in the same dependency graph and returns an expected owning
handle. Check that expected result before using `(*linear_system)->system()` and
`(*linear_system)->preconditioner()` with a linear solver. The former exposes the
same normalized residual and Jacobian interface as `evaluate()`. The latter uses
exactly those coordinates, retaining the dimensional inverse and its physical
Jacobian underneath. `contract()` distinguishes fixed linear and variable inverse
applications when selecting the outer Krylov method.

Unchanged inputs reuse prepared results. Candidate changes refresh the required
physical blocks and numerical factors, while reference geometry and normalization
can remain reusable. A trial evaluation does not request the inverse branch.
Numerical inverse setup failures are returned as `preconditioner::Failure` in the
equilibrium failure variant; earlier owning linear systems remain usable after
failure, later evaluation, and evaluator destruction. Successful preparation alone
does not accept a nonlinear step.

To change the reference or coordinate convention, construct a new configuration
and binding. For additional reference strategies, provide
`prepare_reference_scales(strategy, gravitational_constant)` returning
`std::expected<StellarReferenceScales, std::string>` in the strategy's namespace.
Explicit `StellarReferenceScales` are also accepted as the scheme.

Advanced configurations may supply a complete `NormalizationPlan`, plus ordinary
dependency `sources` and `computations`. Compile-time checks require one supported
rule per unknown and equation row. Existing runtime prescriptions and custom
metrics use these same hooks; no separate normalization cache or scheduler is
introduced.