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:

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:

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:

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:

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.