Emily Boudreaux
6b57738330
Note for now much of the docs is LLM written. This was done to save time and focus on code development. The docs have been LOOSLEY reviewed. However, the code itself has been either fully written by hand or fully reviewed. Therefore, the docs should at the moment be seen as a loose resource rather than authoratiative (and this is on top of the general ability for docs to drift)
212 lines
10 KiB
Markdown
212 lines
10 KiB
Markdown
# 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.
|