1366 lines
56 KiB
Markdown
1366 lines
56 KiB
Markdown
# Surface Deformation Prescription and Domain Extension Design
|
|
|
|
**Repository:** `/Users/tboudreaux/Programming/4DSTAR/MeanField`
|
|
**Design date:** 2026-08-30
|
|
**Status:** design only; no production code changes are authorized by this document
|
|
**First implementation target:** `NodalRadialSurface` with one linear stellar-interior extension and one linear vacuum extension
|
|
|
|
## 1. Executive decision
|
|
|
|
The coupled equilibrium system should cease treating every stellar-interior and vacuum displacement degree of freedom as an independent nonlinear unknown. The physical geometric unknown for the present problem is the shape of the pressure-zero stellar surface. Interior and vacuum mesh coordinates should be generated deterministically from that surface.
|
|
|
|
The target nonlinear state is therefore
|
|
|
|
\[
|
|
x = [\rho,\; s,\; g,\; \Phi,\; h,\; C],
|
|
\]
|
|
|
|
where \(s\) contains the parameters of a selected **surface deformation prescription**. For the first implementation, \(s\) contains one radial displacement coefficient per scalar stellar-surface node.
|
|
|
|
The full volume displacement remains an important internal field:
|
|
|
|
\[
|
|
d = \mathcal D(s).
|
|
\]
|
|
|
|
It is still consumed by the existing domain mapper and all geometry-dependent physics operators, but it is no longer stored as an independent block in the root nonlinear vector.
|
|
|
|
The map \(\mathcal D\) is a composition of three independently extensible parts:
|
|
|
|
\[
|
|
s
|
|
\xrightarrow{\mathcal P}
|
|
d_\Gamma
|
|
\xrightarrow{(\mathcal E_\star,\mathcal E_{\rm vac})}
|
|
d.
|
|
\]
|
|
|
|
- \(\mathcal P\) is the selected surface deformation prescription.
|
|
- \(d_\Gamma\) is a vector displacement trace on the stellar surface.
|
|
- \(\mathcal E_\star\) extends that trace through the stellar interior.
|
|
- \(\mathcal E_{\rm vac}\) extends the same trace through the vacuum and fixes the reference infinity boundary.
|
|
|
|
The first implementation should use linear, reference-domain maps for all three operations. This keeps the exact Newton derivative simple, makes transpose application inexpensive, and avoids placing an interior or vacuum mesh solve inside every Krylov matrix-vector product.
|
|
|
|
The physical surface condition remains independently specified as constant pressure, presently \(P=0\). A surface deformation prescription describes how the surface location is parameterized; it does not describe which thermodynamic condition defines the surface.
|
|
|
|
## 2. Goals
|
|
|
|
This phase must provide:
|
|
|
|
1. A public concept named `SurfaceDeformationPrescription`.
|
|
2. A first concrete prescription named `NodalRadialSurface`.
|
|
3. Independently selectable stellar-interior and vacuum deformation extensions.
|
|
4. One initial implementation for each extension domain.
|
|
5. A composed surface-to-volume deformation map with exact forward and transpose actions.
|
|
6. A reduced coupled block system whose geometry block contains surface parameters rather than volume displacement coefficients.
|
|
7. Preservation of the existing element-level gravity, closure, hydrostatic, mass, pressure-force, gravity-force, and rotational-force kernels.
|
|
8. Compile-time contract checks for concrete prescriptions and extension policies.
|
|
9. A runtime bridge so the prepared coupled operator can consume different concrete prescriptions without virtual dispatch in quadrature loops.
|
|
10. Tests that establish topology, geometry, derivative, transpose, block-routing, physical, and parallel invariants.
|
|
11. A migration path that keeps each implementation phase compilable and independently testable.
|
|
|
|
## 3. Non-goals
|
|
|
|
This phase will not:
|
|
|
|
- implement a spherical-harmonic prescription;
|
|
- implement nonlinear elasticity;
|
|
- support non-star-shaped stellar surfaces;
|
|
- support topology changes or multiple disconnected stellar surfaces;
|
|
- change the physical constant-pressure surface condition;
|
|
- change the EOS architecture;
|
|
- implement the nonlinear solver or its final preconditioner;
|
|
- optimize for extreme rotation before the first reduced system is verified;
|
|
- remove the full displacement finite-element space from `FEM`;
|
|
- rewrite existing element kernels to understand surface parameters;
|
|
- introduce Lagrange multipliers for mesh or surface enforcement.
|
|
|
|
## 4. Terminology and naming
|
|
|
|
Clear separation of the following terms is required.
|
|
|
|
### 4.1 Surface condition
|
|
|
|
A **surface condition** specifies the physical relation that must hold on the stellar boundary. The current condition is constant pressure, normally
|
|
|
|
\[
|
|
P=0.
|
|
\]
|
|
|
|
For the polytropic EOS this is converted by the EOS into the appropriate enthalpy-carrier relation. This remains the responsibility of the existing surface-condition and EOS machinery.
|
|
|
|
Recommended code terminology:
|
|
|
|
- SurfaceConditionPrescription
|
|
- `ConstantPressureSurface`
|
|
- `CompiledPressureSurfaceConstraint`
|
|
|
|
The existing `SurfacePrescription` concept in `models/stellar_model.cppm` refers to this physical condition. Because the new geometry system introduces a second legitimate kind of surface prescription, it should be renamed to `SurfaceConditionPrescription` during the semantic-cleanup phase. Likewise, the `surfacePrescription()` accessor should become `surfaceConditionPrescription()`.
|
|
|
|
### 4.2 Surface deformation prescription
|
|
|
|
A **surface deformation prescription** specifies the finite-dimensional coordinates used to represent the surface location.
|
|
|
|
Examples include:
|
|
|
|
- `NodalRadialSurface`;
|
|
- a future `SphericalHarmonicSurface`;
|
|
- a future `AxisymmetricLegendreSurface`;
|
|
- a future `NodalNormalSurface`.
|
|
|
|
It maps parameter coefficients to a vector displacement trace on the stellar surface. It does not know about pressure, enthalpy, the EOS, or the force law.
|
|
|
|
### 4.3 Domain deformation extension
|
|
|
|
A **domain deformation extension** determines the nonphysical placement of mesh coordinates away from the stellar surface.
|
|
|
|
The stellar-interior and vacuum extensions are separate policies because their regularity, anchoring, and mesh-quality requirements differ.
|
|
|
|
Recommended code terminology:
|
|
|
|
- `InteriorDeformationExtension`
|
|
- `VacuumDeformationExtension`
|
|
- `PowerLawRadialInteriorExtension`
|
|
- `FixedInfinityRadialVacuumExtension`
|
|
|
|
### 4.4 Prepared domain deformation
|
|
|
|
The **prepared domain deformation** is the compiled, discretization-specific composition of a surface prescription and the two extensions. It is the only object the root coupled operator needs to understand.
|
|
|
|
Recommended name:
|
|
|
|
- `PreparedDomainDeformation`
|
|
|
|
## 5. Why this change is needed
|
|
|
|
The gravity-completed gauge experiments performed on 2026-08-30 showed:
|
|
|
|
- a non-rigid tangential displacement response of approximately \(1.03\times10^{-7}\) after gravity completion;
|
|
- rigid rotational responses of approximately \(1.32\times10^{-6}\);
|
|
- a vacuum-only displacement response of approximately \(3.81\times10^{-6}\);
|
|
- gravity-completion reductions of approximately \(6.65\times10^4\) for the tangential probe and \(1.20\times10^5\) for the vacuum probe.
|
|
|
|
The raw mesh displacements appeared strongly constrained when gravity was held fixed. Once the gravitational field was permitted to change consistently with the geometry, the same directions became nearly null.
|
|
|
|
This establishes that the current displacement block contains at least three conceptually different classes of variation:
|
|
|
|
1. normal surface motion that changes the physical stellar shape;
|
|
2. tangential or interior coordinate relabeling;
|
|
3. vacuum mesh-coordinate motion.
|
|
|
|
Only the first class should be an independent physical shape unknown. The other two require deterministic gauge choices.
|
|
|
|
The existing center row replacement removes rigid translation at the computational origin but does not control tangential, rotational, or vacuum coordinate freedoms. The near-null families are likely high-dimensional, so deflating a short list of vectors is not a complete solution.
|
|
|
|
## 6. Current implementation state
|
|
|
|
### 6.1 Current root state and residual
|
|
|
|
The root block form in `libmeanfield/interface/utils/blocks.cppm` is
|
|
|
|
\[
|
|
x_{\rm current}=[\rho,d,g,\Phi,h,C],
|
|
\]
|
|
|
|
\[
|
|
R_{\rm current}=[R_g,R_\Phi,R_\rho,R_d,R_h,R_M].
|
|
\]
|
|
|
|
The current Jacobian coupling is:
|
|
|
|
| Residual | \(\rho\) | \(d\) | \(g\) | \(\Phi\) | \(h\) | \(C\) |
|
|
|---|:---:|:---:|:---:|:---:|:---:|:---:|
|
|
| \(R_g\) | | yes | yes | yes | | |
|
|
| \(R_\Phi\) | yes | yes | yes | | | |
|
|
| \(R_\rho\) | yes | yes | | | yes | |
|
|
| \(R_d\) | yes | yes | yes | | yes | |
|
|
| \(R_h\) | | yes | | yes | yes | yes |
|
|
| \(R_M\) | yes | yes | | | | |
|
|
|
|
The displacement value and residual blocks have the same full-volume H1 vector size, making the root operator square.
|
|
|
|
### 6.2 Current displacement field
|
|
|
|
`field::Displacement` is a vector H1 field of order three supported on `Core`, `Envelope`, and `Vacuum`. Its reduced `FieldDofMap` therefore ordinarily contains all displacement true DOFs.
|
|
|
|
The `FEM` object owns:
|
|
|
|
- the displacement finite-element collection;
|
|
- the volume displacement finite-element space;
|
|
- a displacement grid function.
|
|
|
|
These objects remain necessary after the reduction because the domain mapper and existing kernels consume a full volume displacement field.
|
|
|
|
### 6.3 Current coupled orchestration
|
|
|
|
`PreparedStellarEquilibriumOperator` presently:
|
|
|
|
1. extracts the complete displacement block from the root state;
|
|
2. packs it into the gravity state;
|
|
3. supplies it to closure and hydrostatic preparation;
|
|
4. obtains the mechanical displacement residual from `PreparedDisplacementResidualOperator`;
|
|
5. replaces three force rows with the centering constraint;
|
|
6. routes complete displacement directions to every geometry-dependent Jacobian action.
|
|
|
|
This top-level orchestration is the principal integration point for the change.
|
|
|
|
### 6.4 Existing reusable infrastructure
|
|
|
|
The repository already provides:
|
|
|
|
- semantic discovery of the `StellarSurface` through `DomainSchema`;
|
|
- `FieldBoundaryDofMap` for boundary DOFs in reduced field ordering;
|
|
- a full-volume vector H1 displacement space;
|
|
- stateless mapped geometry and first variations;
|
|
- complete displacement-column Jacobian actions for all current child operators;
|
|
- a registered `Displacement::Form::MeshExtension` quadrature form;
|
|
- dependency identity/revision tracking;
|
|
- matrix-free prepared residual and Jacobian composition.
|
|
|
|
The registered mesh-extension form is useful for future harmonic or elastic auxiliary extensions. It is not necessary for the first analytic linear extension.
|
|
|
|
## 7. Target mathematical architecture
|
|
|
|
Let \(q\in\mathbb R^{N_q}\) be the surface deformation parameters. Use \(q\), rather than \(s\), in operator equations below to distinguish parameter coefficients from the geometric surface itself.
|
|
|
|
The surface prescription produces a vector trace:
|
|
|
|
\[
|
|
d_\Gamma=\mathcal P(q),
|
|
\qquad
|
|
d_\Gamma\in V_\Gamma^d.
|
|
\]
|
|
|
|
The two extension policies produce the volume displacement:
|
|
|
|
\[
|
|
d_\star=\mathcal E_\star(d_\Gamma),
|
|
\qquad
|
|
d_{\rm vac}=\mathcal E_{\rm vac}(d_\Gamma).
|
|
\]
|
|
|
|
The combined continuous H1 displacement is
|
|
|
|
\[
|
|
d=\mathcal D(q)
|
|
=\operatorname{merge}(d_\star,d_{\rm vac}).
|
|
\]
|
|
|
|
The merge must assign the shared stellar-surface true DOFs exactly once and must produce identical values from both extension policies on that trace.
|
|
|
|
For the initial implementation, all maps are linear on the fixed reference mesh:
|
|
|
|
\[
|
|
d_\Gamma=Pq,
|
|
\qquad
|
|
d=Ed_\Gamma,
|
|
\qquad
|
|
d=LPq\equiv Lq,
|
|
\]
|
|
|
|
where the composed lift is denoted simply by \(L=EP\).
|
|
|
|
The reduced mechanical residual is obtained from virtual work:
|
|
|
|
\[
|
|
R_q=L^T R_d.
|
|
\]
|
|
|
|
This is not equivalent to selecting force residual entries that happen to lie on the surface. The transpose accumulates the work of every interior and vacuum displacement response induced by each surface parameter.
|
|
|
|
## 8. Surface deformation prescription contract
|
|
|
|
### 8.1 Compile-time contract
|
|
|
|
A concrete prescription should satisfy a `SurfaceDeformationPrescription` concept. The exact spelling may evolve during implementation, but the concept must require the following semantic operations:
|
|
|
|
- validation independent of the mesh where possible;
|
|
- compilation against a concrete surface trace discretization and reference geometry;
|
|
- a parameter count after compilation;
|
|
- evaluation from parameters to surface vector displacement;
|
|
- derivative action;
|
|
- derivative-transpose action;
|
|
- a runtime descriptor of invariant capabilities.
|
|
|
|
The hot loops must remain in concrete compiled implementations. Type erasure, if used, should dispatch once per complete vector operation, not once per surface DOF or quadrature point.
|
|
|
|
### 8.2 Runtime contract
|
|
|
|
The root prepared operator should be able to consume different compiled prescriptions through a common runtime interface. Consumers need reliable answers to questions such as:
|
|
|
|
- parameter dimension;
|
|
- spatial dimension;
|
|
- whether the prescription is linear on the reference geometry;
|
|
- whether it produces radial-only, normal-only, or general vector surface motion;
|
|
- whether exact derivative-transpose application is available;
|
|
- whether the prescription requires a star-shaped reference surface;
|
|
- whether the prescription itself removes translation or orientation modes.
|
|
|
|
The runtime contract should expose capabilities, not a closed enumeration of every concrete prescription type. Adding `SphericalHarmonicSurface` should not require editing a central `switch` statement.
|
|
|
|
### 8.3 Nonlinear future prescriptions
|
|
|
|
A future prescription may use shape-dependent normals or another nonlinear parameterization. The general deformation map is then \(d=\mathcal D(q)\), with Jacobian \(D\mathcal D(q)\).
|
|
|
|
The physical shape residual is
|
|
|
|
\[
|
|
R_q=D\mathcal D(q)^T R_d.
|
|
\]
|
|
|
|
Its exact Newton action contains two terms:
|
|
|
|
\[
|
|
\delta R_q
|
|
=D\mathcal D(q)^T\delta R_d
|
|
+D\!\left[D\mathcal D(q)^T\right][\delta q]R_d.
|
|
\]
|
|
|
|
The second, curvature-like term vanishes for every linear reference-domain prescription and extension planned in the initial implementation. The public contract should nevertheless reserve a clearly named operation for it rather than allowing a future nonlinear implementation to silently return an inexact Newton Jacobian.
|
|
|
|
## 9. First prescription: `NodalRadialSurface`
|
|
|
|
### 9.1 Unknowns
|
|
|
|
`NodalRadialSurface` uses one scalar parameter per scalar H1 trace DOF on `StellarSurface`:
|
|
|
|
\[
|
|
q=(q_1,\ldots,q_{N_q}).
|
|
\]
|
|
|
|
At a surface trace node with fixed reference position \(X_i\), reference center \(X_c\), and radial unit vector
|
|
|
|
\[
|
|
e_{r,i}=\frac{X_i-X_c}{\lVert X_i-X_c\rVert},
|
|
\]
|
|
|
|
the prescribed vector surface displacement is
|
|
|
|
\[
|
|
d_{\Gamma,i}=q_i e_{r,i}.
|
|
\]
|
|
|
|
The new surface point is
|
|
|
|
\[
|
|
x_{\Gamma,i}=X_i+q_i e_{r,i}.
|
|
\]
|
|
|
|
The parameter is a radial displacement with units of length. It is not a radius multiplier and not an arbitrary Cartesian displacement component.
|
|
|
|
### 9.2 Reference assumptions
|
|
|
|
The first implementation requires:
|
|
|
|
- a selected reference center;
|
|
- a star-shaped reference stellar surface relative to that center;
|
|
- strictly positive reference surface radius in every represented direction;
|
|
- one radial intersection of each represented ray with the stellar surface;
|
|
- a conforming surface trace;
|
|
- the current vector H1 displacement ordering and compatible scalar trace order.
|
|
|
|
Construction must fail clearly if these assumptions are not satisfied.
|
|
|
|
### 9.3 Surface trace discretization
|
|
|
|
The surface parameter space should be an actual scalar trace coordinate system, not a list of every Cartesian component in `FieldBoundaryDofMap`.
|
|
|
|
The preferred infrastructure is:
|
|
|
|
- a scalar H1 finite-element space using the same scalar collection/order as the displacement field;
|
|
- a `SurfaceFieldDofMap` or boundary-supported `FieldDofMap` containing only true DOFs on `StellarSurface`;
|
|
- deterministic parallel reduced ordering;
|
|
- a vector surface trace workspace used internally by the compiled prescription.
|
|
|
|
This requires extending field-support machinery to represent boundary-supported fields. It should not reuse the enthalpy trace as the geometry coordinate system merely because both happen to be scalar H1 traces.
|
|
|
|
### 9.4 Linear operator representation
|
|
|
|
The compiled nodal-radial prescription is a linear operator
|
|
|
|
\[
|
|
P:\mathbb R^{N_q}\rightarrow\mathbb R^{N_{\Gamma,d}}.
|
|
\]
|
|
|
|
For each surface scalar coefficient, \(P\) writes the coefficient into the three vector components weighted by the fixed reference radial direction. Its transpose computes radial virtual work:
|
|
|
|
\[
|
|
(P^T f_\Gamma)_i=e_{r,i}\cdot f_{\Gamma,i}.
|
|
\]
|
|
|
|
The implementation should provide exact `Mult` and `MultTranspose` operations. An explicit sparse parallel operator is preferred because it makes the transpose identity auditable and applications cheap.
|
|
|
|
### 9.5 Gauge consequences
|
|
|
|
The prescription cannot represent tangential surface relabeling. Rigid rotation of a spherical reference surface is therefore absent from the shape coordinates.
|
|
|
|
The extension described below fixes the reference center, so exact rigid translation of the complete mesh is also absent. A nodal radial pattern with dominant \(\ell=1\) content can still move the surface approximately like a translation while deforming the anchored interior. The solver must monitor center of mass, and a future surface-space translation projection or center-of-mass condition may still be desirable.
|
|
|
|
## 10. Extensible stellar-interior deformation
|
|
|
|
### 10.1 Policy contract
|
|
|
|
An `InteriorDeformationExtension` must:
|
|
|
|
- accept a complete vector displacement trace on `StellarSurface`;
|
|
- produce volume displacement values on stellar-supported displacement DOFs;
|
|
- reproduce the supplied trace exactly;
|
|
- enforce a documented center behavior;
|
|
- provide derivative and derivative-transpose actions;
|
|
- report whether it is linear and whether it requires a radial foliation;
|
|
- avoid writing vacuum-only DOFs;
|
|
- support parallel true-DOF ownership.
|
|
|
|
The interface should operate on vector surface traces so it can be reused by radial, normal, spherical-harmonic, or future general-vector surface prescriptions.
|
|
|
|
### 10.2 First implementation: `PowerLawRadialInteriorExtension`
|
|
|
|
The initial interior extension uses fixed reference rays. For a reference point
|
|
|
|
\[
|
|
X=X_c+r\omega,
|
|
\]
|
|
|
|
let \(R_\Gamma^0(\omega)\) be the reference surface radius and let \(d_\Gamma(\omega)\) be the surface displacement interpolated where that ray meets the reference surface. Define
|
|
|
|
\[
|
|
d_\star(X)
|
|
=\left(\frac{r}{R_\Gamma^0(\omega)}\right)^p
|
|
d_\Gamma(\omega).
|
|
\]
|
|
|
|
The first configured power should be
|
|
|
|
\[
|
|
p=2.
|
|
\]
|
|
|
|
This choice:
|
|
|
|
- reproduces the supplied displacement at the stellar surface;
|
|
- fixes the center exactly;
|
|
- makes the displacement vanish at least quadratically toward the center;
|
|
- is linear in the surface trace;
|
|
- is inexpensive to apply;
|
|
- agrees with the radial interior extension pattern already exercised in the domain-mapper tests.
|
|
|
|
The power should be a validated construction option, but the first production model should use two rather than expose an untested menu to end users.
|
|
|
|
### 10.3 Radial correspondence construction
|
|
|
|
For every volume displacement interpolation node, the compiled extension must know how to evaluate the surface trace at the node's reference direction. This should be precomputed once from the reference mesh:
|
|
|
|
1. compute the reference direction from the selected center;
|
|
2. intersect the ray with the reference stellar surface;
|
|
3. locate the corresponding surface element;
|
|
4. evaluate its trace basis at the intersection;
|
|
5. store the interpolation row and radial weight.
|
|
|
|
The resulting extension is a sparse linear operator. The center is handled as an explicit zero row because its direction is undefined.
|
|
|
|
The construction must validate that every noncentral stellar node obtains exactly one valid surface correspondence.
|
|
|
|
### 10.4 Future interior policies
|
|
|
|
The same contract should later permit:
|
|
|
|
- harmonic extension;
|
|
- linear-elastic extension;
|
|
- pseudoelastic extension with radius-dependent stiffness;
|
|
- modal analytic extension for spherical harmonics;
|
|
- transfinite interpolation;
|
|
- mesh-quality-optimized nonlinear extension.
|
|
|
|
If an extension requires an elliptic solve, it should not automatically be eliminated into the reduced surface operator. See Section 17.
|
|
|
|
## 11. Extensible vacuum deformation
|
|
|
|
### 11.1 Policy contract
|
|
|
|
A `VacuumDeformationExtension` must:
|
|
|
|
- accept the same surface vector displacement trace as the interior policy;
|
|
- reproduce it exactly at `StellarSurface`;
|
|
- produce volume displacement values on vacuum-supported displacement DOFs;
|
|
- satisfy a documented condition at `InfinitySurface`;
|
|
- provide derivative and derivative-transpose actions;
|
|
- be compatible with the exterior compactification map;
|
|
- avoid writing stellar-interior-only DOFs except for the shared surface trace;
|
|
- support parallel true-DOF ownership.
|
|
|
|
### 11.2 First implementation: `FixedInfinityRadialVacuumExtension`
|
|
|
|
For a reference vacuum point on the ray \(X=X_c+r\omega\), let \(R_\Gamma^0(\omega)\) be the stellar-surface radius and let \(R_\infty^0(\omega)\) be the reference radius of the finite boundary that is mapped to physical infinity. Define
|
|
|
|
\[
|
|
d_{\rm vac}(X)
|
|
=\frac{R_\infty^0(\omega)-r}
|
|
{R_\infty^0(\omega)-R_\Gamma^0(\omega)}
|
|
d_\Gamma(\omega).
|
|
\]
|
|
|
|
This extension:
|
|
|
|
- equals the supplied surface displacement at the stellar boundary;
|
|
- vanishes at the reference infinity boundary;
|
|
- is linear in the surface trace;
|
|
- is compatible with retaining a fixed compactification boundary;
|
|
- agrees with the exterior radial decay pattern already exercised in the domain-mapper tests.
|
|
|
|
As in the interior, surface correspondence and interpolation rows should be precomputed once.
|
|
|
|
### 11.3 Future vacuum policies
|
|
|
|
The same contract should later permit:
|
|
|
|
- harmonic vacuum extension;
|
|
- compactification-coordinate-weighted extension;
|
|
- elastic or Winslow-type mesh smoothing;
|
|
- analytic modal continuation;
|
|
- a moving finite outer boundary, if a future compactification requires it.
|
|
|
|
## 12. Composed prepared deformation operator
|
|
|
|
The prepared geometry component should own or reference three compiled operators:
|
|
|
|
\[
|
|
P:q\mapsto d_\Gamma,
|
|
\]
|
|
|
|
\[
|
|
E_\star:d_\Gamma\mapsto d_\star,
|
|
\]
|
|
|
|
\[
|
|
E_{\rm vac}:d_\Gamma\mapsto d_{\rm vac}.
|
|
\]
|
|
|
|
It should expose the composed map
|
|
|
|
\[
|
|
L:q\mapsto d.
|
|
\]
|
|
|
|
Required operations for the initial linear implementation are:
|
|
|
|
- `BuildVolumeDisplacement(q, d)`;
|
|
- `ApplyJacobian(qDirection, dDirection)`;
|
|
- `ApplyJacobianTranspose(dDual, qDual)`;
|
|
- access to parameter, surface-trace, and full displacement sizes;
|
|
- access to descriptors and reference-geometry validation reports;
|
|
- preparation statistics and dependency stamps where rebuilding can occur.
|
|
|
|
The composed implementation need not form the matrix product \(EP\). It can apply
|
|
|
|
\[
|
|
q\xrightarrow{P}d_\Gamma\xrightarrow{E}d
|
|
\]
|
|
|
|
and
|
|
|
|
\[
|
|
f_d\xrightarrow{E^T}f_\Gamma\xrightarrow{P^T}f_q.
|
|
\]
|
|
|
|
Keeping the factors separate preserves extensibility and avoids unnecessary fill.
|
|
|
|
The merge between interior and vacuum extensions must have explicit ownership rules:
|
|
|
|
- shared stellar-surface DOFs are assigned once;
|
|
- stellar-only DOFs come from the interior extension;
|
|
- vacuum-only DOFs come from the vacuum extension;
|
|
- infinity-surface DOFs satisfy the vacuum policy;
|
|
- no DOF is silently left uninitialized or written inconsistently.
|
|
|
|
## 13. Revised coupled block system
|
|
|
|
### 13.1 Value and residual blocks
|
|
|
|
The root value blocks become
|
|
|
|
```text
|
|
[rho, surface_deformation, gravity_gradient,
|
|
gravity_potential, enthalpy, bernoulli_constant]
|
|
```
|
|
|
|
and the root residual blocks become
|
|
|
|
```text
|
|
[gravity_gradient_relation, poisson,
|
|
barotropic_closure, surface_shape_equilibrium,
|
|
hydrostatic_equilibrium, mass_normalization]
|
|
```
|
|
|
|
In mathematical notation:
|
|
|
|
\[
|
|
x_{\rm new}=[\rho,q,g,\Phi,h,C],
|
|
\]
|
|
|
|
\[
|
|
R_{\rm new}=[R_g,R_\Phi,R_\rho,R_q,R_h,R_M].
|
|
\]
|
|
|
|
The root operator remains square because the full displacement value and residual blocks are both replaced by parameter-sized surface deformation blocks.
|
|
|
|
### 13.2 Compile-time coupling metadata
|
|
|
|
The new Jacobian form is:
|
|
|
|
| Residual | \(\rho\) | \(q\) | \(g\) | \(\Phi\) | \(h\) | \(C\) |
|
|
|---|:---:|:---:|:---:|:---:|:---:|:---:|
|
|
| \(R_g\) | | yes | yes | yes | | |
|
|
| \(R_\Phi\) | yes | yes | yes | | | |
|
|
| \(R_\rho\) | yes | yes | | | yes | |
|
|
| \(R_q\) | yes | yes | yes | | yes | |
|
|
| \(R_h\) | | yes | | yes | yes | yes |
|
|
| \(R_M\) | yes | yes | | | | |
|
|
|
|
The graph is physically unchanged. Every former displacement column becomes a surface-deformation column through the lift \(L\), and the former mechanical residual row becomes a surface-shape row through \(L^T\).
|
|
|
|
Recommended compile-time block names are:
|
|
|
|
```cpp
|
|
surface_deformation::radial_parameters::value
|
|
surface_deformation::shape_equilibrium::residual
|
|
```
|
|
|
|
The value-block type should not be named after `NodalRadialSurface`, because a future spherical-harmonic prescription must use the same root block role. A more representation-neutral final spelling such as
|
|
|
|
```cpp
|
|
surface_deformation::parameters::value
|
|
```
|
|
|
|
is preferable if the nested block system permits it cleanly.
|
|
|
|
### 13.3 Child gravity layout
|
|
|
|
The internal gravity operator currently packs
|
|
|
|
```text
|
|
[rho, d, g, Phi]
|
|
```
|
|
|
|
and should continue to do so. It is an internal child layout, not the root nonlinear layout. `PreparedStellarEquilibriumOperator` will construct \(d=Lq\) before packing the gravity state.
|
|
|
|
This avoids forcing the gravity implementation to know how the surface is represented.
|
|
|
|
### 13.4 Mechanical row pullback
|
|
|
|
`PreparedDisplacementResidualOperator` should continue to assemble the full-volume mechanical dual vector
|
|
|
|
\[
|
|
R_d=R_d^{\rm pressure}+R_d^{\rm gravity}+R_d^{\rm rotation}.
|
|
\]
|
|
|
|
The root composer then constructs
|
|
|
|
\[
|
|
R_q=L^T R_d.
|
|
\]
|
|
|
|
The same rule applies to Jacobian actions. For a complete reduced direction
|
|
|
|
\[
|
|
\delta x=[\delta\rho,\delta q,\delta g,\delta\Phi,\delta h,\delta C],
|
|
\]
|
|
|
|
first construct
|
|
|
|
\[
|
|
\delta d=L\delta q.
|
|
\]
|
|
|
|
Then use the existing child actions to obtain
|
|
|
|
\[
|
|
\delta R_d
|
|
=J_{d\rho}\delta\rho
|
|
+J_{dd}\delta d
|
|
+J_{dg}\delta g
|
|
+J_{dh}\delta h,
|
|
\]
|
|
|
|
and return
|
|
|
|
\[
|
|
\delta R_q=L^T\delta R_d.
|
|
\]
|
|
|
|
For the initial linear deformation map, no second-derivative term is present.
|
|
|
|
### 13.5 Pressure surface row
|
|
|
|
The constant-pressure condition remains a row replacement in the hydrostatic/enthalpy residual block. It is evaluated on the same topological stellar boundary after the geometry has been generated from \(q\).
|
|
|
|
The surface deformation parameter count need not equal the number of enthalpy surface rows. The geometric and thermodynamic surface spaces have different roles and may have different finite-element orders in the future.
|
|
|
|
### 13.6 Centering row
|
|
|
|
The current `PreparedCenteringConstraint` replaces three full displacement-force rows. Those rows no longer exist in the root layout and the constraint cannot be retained unchanged.
|
|
|
|
For the first radial extension, the volume map fixes the reference center exactly. This removes exact rigid translation from the admissible geometry family. The old center row replacement should therefore be removed from the reduced root operator, while center-of-mass displacement remains a required diagnostic.
|
|
|
|
A later surface-space gauge may explicitly remove three translation-like modes or impose center of mass. That decision should be based on the conditioning of the reduced system rather than carried over mechanically from the full displacement formulation.
|
|
|
|
## 14. Residual and Jacobian execution flow
|
|
|
|
### 14.1 Preparation at a nonlinear state
|
|
|
|
At each root `Prepare` call:
|
|
|
|
1. extract \(\rho,q,g,\Phi,h,C\) from the reduced block vector;
|
|
2. apply the prepared deformation map to build the full displacement \(d=Lq\);
|
|
3. validate that the mapped geometry is finite and orientation-preserving;
|
|
4. pack \(\rho,d,g,\Phi\) into the existing gravity child state;
|
|
5. prepare gravity, closure, hydrostatic, mechanical, and mass children at the same generated geometry;
|
|
6. prepare the pressure surface constraint from the carrier field;
|
|
7. build the full mechanical residual \(R_d\);
|
|
8. apply \(L^T\) to obtain \(R_q\);
|
|
9. concatenate the six reduced residual blocks.
|
|
|
|
The generated displacement must be cached for the complete frozen Newton linearization.
|
|
|
|
### 14.2 Matrix-free Jacobian application
|
|
|
|
For each root `Mult` call:
|
|
|
|
1. extract \(\delta\rho,\delta q,\delta g,\delta\Phi,\delta h,\delta C\);
|
|
2. compute \(\delta d=L\delta q\);
|
|
3. call all existing child Jacobian actions with \(\delta d\);
|
|
4. obtain the full mechanical action \(\delta R_d\);
|
|
5. compute \(\delta R_q=L^T\delta R_d\);
|
|
6. concatenate the reduced action blocks.
|
|
|
|
No prescription or extension should reprepare during `Mult`.
|
|
|
|
### 14.3 Performance requirement
|
|
|
|
For the first linear implementation, one reduced root Jacobian application should add only:
|
|
|
|
- one surface-prescription forward application;
|
|
- one domain-extension forward application;
|
|
- one domain-extension transpose application;
|
|
- one surface-prescription transpose application.
|
|
|
|
These operations should be sparse or matrix-free linear maps with no iterative solves and no per-application mesh searches.
|
|
|
|
## 15. Dependency and lifetime design
|
|
|
|
### 15.1 Root dependencies
|
|
|
|
The public `StellarEquilibriumDependencies` should replace its independent `displacement` stamp with a `surfaceDeformation` stamp.
|
|
|
|
The prepared domain deformation owns a derived generated-geometry revision. When the surface parameters change, it rebuilds the cached displacement and advances that geometry revision. Child operators continue receiving a displacement dependency stamp constructed from the generated geometry identity and revision.
|
|
|
|
This preserves the existing child contracts while making the root dependency semantically correct.
|
|
|
|
### 15.2 Discretization changes
|
|
|
|
The compiled ray correspondences and sparse lift/restriction maps depend on:
|
|
|
|
- mesh identity and sequence;
|
|
- displacement finite-element space identity and sequence;
|
|
- surface trace space identity and sequence;
|
|
- reference center;
|
|
- domain schema and boundary attributes;
|
|
- reference compactification boundary.
|
|
|
|
Any change to these inputs invalidates the compiled deformation object. Adaptive mesh changes should require reconstruction initially; incremental update is a deferred feature.
|
|
|
|
### 15.3 Model ownership
|
|
|
|
`StellarModel` should own:
|
|
|
|
- the structure prescription;
|
|
- the physical surface-condition prescription;
|
|
- the surface deformation prescription;
|
|
- the interior deformation-extension prescription;
|
|
- the vacuum deformation-extension prescription.
|
|
|
|
These objects are configuration-level policies and should retain stable addresses across model moves, following the current structure/surface ownership pattern.
|
|
|
|
Compilation against `FEM` and the reference geometry should occur when constructing the prepared coupled operator, because the user-facing model does not own a discretization.
|
|
|
|
### 15.4 User-facing construction
|
|
|
|
The intended construction should read approximately as:
|
|
|
|
```cpp
|
|
models::StellarModel model{
|
|
structure::PolytropicStructure{equationOfState, targetMass},
|
|
surface::ConstantPressureSurface{eos::PressureValue{0.0}},
|
|
deformation::NodalRadialSurface{referenceCenter},
|
|
deformation::PowerLawRadialInteriorExtension{2},
|
|
deformation::FixedInfinityRadialVacuumExtension{}
|
|
};
|
|
```
|
|
|
|
Exact constructor organization may use a named aggregate to avoid a long positional argument list. Astronomer-facing names should communicate physical meaning without requiring knowledge of MFEM.
|
|
|
|
## 16. Compile-time and runtime extensibility
|
|
|
|
The architecture should follow the same broad principle adopted for the EOS:
|
|
|
|
- concrete types satisfy compile-time concepts;
|
|
- compilation produces optimized concrete vector operations;
|
|
- a small runtime bridge permits heterogeneous models to reach common consumers;
|
|
- runtime dispatch occurs once per whole vector operation, outside quadrature loops.
|
|
|
|
Separate concepts should exist for:
|
|
|
|
- `SurfaceDeformationPrescription`;
|
|
- `InteriorDeformationExtension`;
|
|
- `VacuumDeformationExtension`.
|
|
|
|
The concepts should reject incomplete implementations at compile time. Runtime validation should reject incompatible combinations such as:
|
|
|
|
- a radial prescription on a non-star-shaped reference surface;
|
|
- mismatched spatial dimensions;
|
|
- an extension that cannot reproduce the selected surface trace space;
|
|
- a vacuum extension that does not define infinity-boundary behavior;
|
|
- missing derivative-transpose support;
|
|
- a nonlinear prescription that does not provide its pullback-derivative term.
|
|
|
|
## 17. Analytic versus elliptic extensions
|
|
|
|
The reduced surface-only block is the correct first implementation because the selected extensions are explicit linear maps.
|
|
|
|
Future harmonic or elastic extensions require a separate architectural decision. If
|
|
|
|
\[
|
|
K_{ii}d_i=-K_{i\Gamma}d_\Gamma,
|
|
\]
|
|
|
|
then evaluating the eliminated map
|
|
|
|
\[
|
|
d_i=-K_{ii}^{-1}K_{i\Gamma}d_\Gamma
|
|
\]
|
|
|
|
requires a linear solve. Its transpose requires an adjoint solve. Hiding both inside every root Jacobian application may be expensive and may obscure preconditioning.
|
|
|
|
For such policies, two valid realizations should be supported eventually:
|
|
|
|
1. **Eliminated extension:** use a reusable factorization or strongly preconditioned solve inside the deformation map when the cost is demonstrably acceptable.
|
|
2. **Auxiliary mesh block:** retain interior/vacuum displacement as auxiliary unknowns and expose the harmonic or elastic mesh equation as a sparse block.
|
|
|
|
The second realization still treats the surface as the only physical geometric unknown. The additional volume coefficients are deterministic mesh variables, not physical stellar coordinates.
|
|
|
|
The extension policy should therefore describe the mathematical extension, while a separate realization choice determines whether it is explicitly eliminated or retained as an auxiliary block.
|
|
|
|
No auxiliary block is needed for `PowerLawRadialInteriorExtension` or `FixedInfinityRadialVacuumExtension`.
|
|
|
|
## 18. Expected interaction with the preconditioner
|
|
|
|
The reduced geometry block is
|
|
|
|
\[
|
|
J_{qq}=L^TJ_{dd}L
|
|
\]
|
|
|
|
for the first linear deformation map. Other geometry couplings are similarly lifted or pulled back:
|
|
|
|
\[
|
|
J_{gq}=J_{gd}L,
|
|
\qquad
|
|
J_{\Phi q}=J_{\Phi d}L,
|
|
\qquad
|
|
J_{qg}=L^TJ_{dg},
|
|
\]
|
|
|
|
and so forth.
|
|
|
|
The surface block is much smaller than the current full displacement block and excludes the measured tangential and vacuum gauge families by construction. This should improve conditioning before any sophisticated preconditioner is introduced.
|
|
|
|
The first preconditioner need not form \(L^TJ_{dd}L\) exactly. Possible approximations include:
|
|
|
|
- a diagonal or lumped surface response;
|
|
- a surface mass plus surface Laplacian operator;
|
|
- a low-order radial stiffness approximation;
|
|
- a spherical-harmonic diagonal approximation once that prescription exists;
|
|
- a geometry-gravity Schur approximation.
|
|
|
|
The block metadata must describe the reduced surface coupling from the start so later preconditioner code does not depend on obsolete full-displacement sizes.
|
|
|
|
## 19. Code cross section
|
|
|
|
### 19.1 High-impact files
|
|
|
|
| Area | Current file | Required change |
|
|
|---|---|---|
|
|
| Root block form | `libmeanfield/interface/utils/blocks.cppm` | Replace root displacement value/residual roles with surface-deformation parameter/equilibrium roles; preserve internal gravity and pressure-force forms where they still describe full displacement |
|
|
| Root interface | `libmeanfield/interface/operators/prepared_stellar_equilibrium.cppm` | Own/expose prepared deformation, revise dependencies/report/accessors, remove root centering operator |
|
|
| Root implementation | `libmeanfield/impl/operators/prepared_stellar_equilibrium.cpp` | Build full displacement from surface parameters, lift directions, pull back mechanical residual/actions, revise layouts and gravity packing |
|
|
| Stellar model | `libmeanfield/interface/models/stellar_model.cppm` | Own deformation and extension prescriptions; disambiguate physical surface-condition naming |
|
|
| FEM | `libmeanfield/interface/fem.cppm`, `libmeanfield/impl/fem.cpp` | Add scalar surface-trace discretization or the scalar space from which its trace map is built |
|
|
| Build registration | `CMakeLists.txt` | Register new interface, implementation, and test translation units |
|
|
|
|
### 19.2 New focused components
|
|
|
|
Likely new modules include:
|
|
|
|
```text
|
|
libmeanfield/interface/deformation/surface_prescription.cppm
|
|
libmeanfield/interface/deformation/nodal_radial_surface.cppm
|
|
libmeanfield/interface/deformation/interior_extension.cppm
|
|
libmeanfield/interface/deformation/vacuum_extension.cppm
|
|
libmeanfield/interface/deformation/radial_extensions.cppm
|
|
libmeanfield/interface/deformation/prepared_domain_deformation.cppm
|
|
|
|
libmeanfield/impl/deformation/nodal_radial_surface.cpp
|
|
libmeanfield/impl/deformation/radial_extensions.cpp
|
|
libmeanfield/impl/deformation/prepared_domain_deformation.cpp
|
|
```
|
|
|
|
The final directory and module names should remain compact. `deformation` is preferable to placing geometric parameterization under the existing physical `surface` namespace.
|
|
|
|
### 19.3 Moderate-impact infrastructure
|
|
|
|
| Area | Required change |
|
|
|---|---|
|
|
| Field registry | Register a scalar surface-deformation quantity or explicitly document why a dedicated trace map is not a registered volume field |
|
|
| Field-DOF support | Add boundary-supported scalar true-DOF resolution with deterministic MPI ordering |
|
|
| Test tags | Add semantic tags such as `surface_deformation`, `interior_extension`, `vacuum_extension`, and `deformation_pullback` to `test_helpers` |
|
|
| Seed/state helpers | Construct a zero surface state for spherical seeds and project requested radial surface perturbations |
|
|
| Dependency helpers | Translate root surface-deformation revisions into generated displacement revisions for child operators |
|
|
|
|
### 19.4 Components expected to remain unchanged
|
|
|
|
The following production kernels and prepared child operators should not require semantic changes:
|
|
|
|
- `DomainMapper` point/volume/face/variation kernels;
|
|
- mapped H(div) gravity mass and source kernels;
|
|
- gravity divergence and boundary operators;
|
|
- barotropic closure kernels;
|
|
- hydrostatic equilibrium kernels;
|
|
- pressure displacement-force kernels;
|
|
- gravity displacement-force kernels;
|
|
- rotational displacement-force kernels;
|
|
- mass-normalization kernels;
|
|
- EOS pressure-surface conversion;
|
|
- prepared pressure-surface row enforcement.
|
|
|
|
They will continue receiving full displacement states and directions from their existing prepared contexts.
|
|
|
|
### 19.5 Tests with substantial migration impact
|
|
|
|
The largest changes will occur in:
|
|
|
|
- `tests/operators/prepared_stellar_equilibrium.cpp`;
|
|
- coupled-state helpers in `tests/test_helpers.cppm`;
|
|
- `experiments/stellar_null_space.cppm`;
|
|
- `experiments/coupled_gauge_modes.cpp`;
|
|
- other experiments that address the root displacement block directly;
|
|
- block-layout tests in `tests/utils/blocks.cpp`;
|
|
- stellar-model ownership and type-contract tests.
|
|
|
|
Individual child-operator tests should remain full-displacement tests and should not be rewritten merely because the root coordinates changed.
|
|
|
|
## 20. Existing tests that must change
|
|
|
|
### 20.1 Block-layout tests
|
|
|
|
Tests must assert:
|
|
|
|
- the new value and residual ordering;
|
|
- the surface block size equals the selected prescription's parameter size;
|
|
- the root operator remains square;
|
|
- every geometry-dependent residual advertises the surface-parameter coupling;
|
|
- internal gravity forms still advertise full displacement coupling;
|
|
- invalid or duplicate block types fail at compile time as before.
|
|
|
|
### 20.2 Prepared stellar-equilibrium composition tests
|
|
|
|
Tests that currently compare the root Jacobian with restricted child Jacobians must change from direct displacement restriction to the exact reduced relation
|
|
|
|
\[
|
|
J_{\rm reduced}=Q^TJ_{\rm full}Q,
|
|
\]
|
|
|
|
where \(Q\) is identity on nongeometry blocks and \(L\) on the geometry block.
|
|
|
|
In particular, the mechanical row comparison must use \(L^T\), not surface-row extraction.
|
|
|
|
### 20.3 Complete centered-difference tests
|
|
|
|
Finite differences must perturb surface parameters, regenerate the volume displacement for the plus and minus states, and compare against the matrix-free reduced Jacobian action.
|
|
|
|
The test must include simultaneous variations in surface parameters and other coupled fields so routing errors cannot cancel in isolated columns.
|
|
|
|
### 20.4 Centering tests
|
|
|
|
Tests asserting replacement of center force rows should be retired from the reduced root operator. They should be replaced with tests that verify:
|
|
|
|
- the interior extension fixes the center;
|
|
- rigid translation is not representable by the composed map;
|
|
- center of mass is reported for representative reduced states;
|
|
- translation-like nodal radial patterns remain finite and do not create an exact algebraic null direction.
|
|
|
|
### 20.5 Null-space experiments
|
|
|
|
The tangential and vacuum-only root directions should no longer be constructible. The existing experiment should be retained as historical evidence or adapted to compare the old full-displacement operator with the reduced lift range.
|
|
|
|
New reduced experiments should probe:
|
|
|
|
- uniform radial scale/homology direction;
|
|
- low-order \(\ell=1\)-like nodal radial patterns;
|
|
- quadrupolar/oblate surface deformation;
|
|
- high-frequency nodal radial deformation;
|
|
- conditioning as surface resolution changes.
|
|
|
|
Rigid rotation of a spherical surface should be absent from the radial parameter space rather than merely producing a small Jacobian response.
|
|
|
|
## 21. New tests
|
|
|
|
All test names and tags must describe physical or computer-science meaning, not implementation phases.
|
|
|
|
### 21.1 Compile-time contract tests
|
|
|
|
Add tests that:
|
|
|
|
- accept `NodalRadialSurface` as a `SurfaceDeformationPrescription`;
|
|
- reject types missing validation, compilation, forward, derivative, or transpose operations;
|
|
- accept both initial extension policies under their respective concepts;
|
|
- reject an interior policy used as a vacuum policy and vice versa;
|
|
- verify const-correct references and noncopyable/movable ownership decisions;
|
|
- verify that physical surface-condition and geometric deformation types cannot be interchanged.
|
|
|
|
Suggested semantic tags:
|
|
|
|
```text
|
|
[surface_deformation][type_contract][unit]
|
|
[interior_extension][type_contract][unit]
|
|
[vacuum_extension][type_contract][unit]
|
|
```
|
|
|
|
### 21.2 Surface trace topology tests
|
|
|
|
Add tests that:
|
|
|
|
- identify exactly the scalar H1 trace DOFs on `StellarSurface`;
|
|
- exclude interior-only and vacuum-only scalar DOFs;
|
|
- preserve shared parallel surface DOFs;
|
|
- produce sorted, unique, deterministic reduced ordering;
|
|
- distinguish `StellarSurface` from `InfinitySurface` through `DomainSchema` rather than raw attributes;
|
|
- work with alternate schema attribute bindings;
|
|
- fail clearly for nonconforming spaces until explicitly supported.
|
|
|
|
### 21.3 `NodalRadialSurface` tests
|
|
|
|
Add tests that:
|
|
|
|
- map each unit parameter to radial vector displacement at the corresponding surface node;
|
|
- produce no tangential component relative to the reference ray;
|
|
- reproduce uniform radial expansion;
|
|
- reproduce a known quadrupolar radial pattern;
|
|
- preserve physical units and sign;
|
|
- reject a surface node at the reference center;
|
|
- reject a non-star-shaped or multiply intersected reference surface;
|
|
- match centered differences;
|
|
- satisfy the adjoint identity
|
|
|
|
\[
|
|
\langle P\delta q,f_\Gamma\rangle
|
|
=\langle\delta q,P^Tf_\Gamma\rangle.
|
|
\]
|
|
|
|
### 21.4 Interior extension tests
|
|
|
|
Add tests that verify:
|
|
|
|
- exact reproduction of arbitrary vector trace data at `StellarSurface`;
|
|
- exact zero displacement at the center;
|
|
- quadratic radial weight at known radii for the initial power two;
|
|
- no writes to vacuum-only DOFs;
|
|
- continuous shared interface values;
|
|
- linearity;
|
|
- centered-difference agreement;
|
|
- transpose identity;
|
|
- successful construction on the production spherical reference mesh;
|
|
- failure when a radial correspondence cannot be established.
|
|
|
|
### 21.5 Vacuum extension tests
|
|
|
|
Add tests that verify:
|
|
|
|
- exact reproduction of arbitrary vector trace data at `StellarSurface`;
|
|
- exact zero displacement at `InfinitySurface`;
|
|
- the expected linear reference-radius weight;
|
|
- no writes to stellar-interior-only DOFs;
|
|
- continuous shared interface values;
|
|
- compatibility with the compactified mapper;
|
|
- linearity;
|
|
- centered-difference agreement;
|
|
- transpose identity;
|
|
- failure for a missing or incompatible infinity boundary.
|
|
|
|
### 21.6 Composed deformation tests
|
|
|
|
Add tests that verify:
|
|
|
|
- every full displacement true DOF is assigned exactly once under the merge rules;
|
|
- interior and vacuum traces agree exactly;
|
|
- `BuildVolumeDisplacement` maps zero surface parameters to zero displacement;
|
|
- the composed map equals separate prescription and extension applications;
|
|
- the composed transpose reverses the operation order correctly;
|
|
- the adjoint identity
|
|
|
|
\[
|
|
\langle L\delta q,f_d\rangle
|
|
=\langle\delta q,L^Tf_d\rangle
|
|
\]
|
|
|
|
holds to a tolerance appropriate for parallel floating-point reductions;
|
|
- mapped Jacobian determinants remain positive for bounded representative surface amplitudes;
|
|
- invalid amplitudes that invert elements are detected before physics assembly.
|
|
|
|
### 21.7 Reduced root-operator tests
|
|
|
|
Add tests that verify:
|
|
|
|
- the root geometry state has parameter size rather than displacement size;
|
|
- child preparation receives exactly \(d=Lq\);
|
|
- every former displacement-column action equals the corresponding full action applied to \(L\delta q\);
|
|
- the reduced mechanical residual equals \(L^TR_d\);
|
|
- the reduced complete Jacobian equals the lifted/restricted full child Jacobian;
|
|
- simultaneous centered differences match every coupled block;
|
|
- pressure-surface row replacement remains exact;
|
|
- mass normalization sees the generated geometry;
|
|
- selective invalidation rebuilds geometry only when surface parameters or discretization change;
|
|
- repeated `Mult` calls do not rebuild the deformation map or reprepare children.
|
|
|
|
### 21.8 Physical and convergence tests
|
|
|
|
Add tests or experiments that verify:
|
|
|
|
- the spherical Lane-Emden seed is unchanged at zero deformation;
|
|
- a uniform positive surface parameter produces expansion everywhere with the specified weights;
|
|
- a quadrupolar radial deformation produces the intended oblate/prolate sign pattern;
|
|
- the reduced analytic equilibrium accuracy remains at the established projection floor;
|
|
- the result converges under surface and volume refinement;
|
|
- changing only the nonphysical extension policy in a future comparison does not change the converged physical surface beyond discretization error.
|
|
|
|
### 21.9 Performance tests
|
|
|
|
Record, without fragile wall-clock regression thresholds:
|
|
|
|
- surface parameter count;
|
|
- volume displacement count;
|
|
- nonzeros or work estimates in \(P\), \(E\), and their transpose paths;
|
|
- number of geometry rebuilds per nonlinear preparation;
|
|
- number of lift and pullback applications per Jacobian action;
|
|
- time fractions for lift, physics action, and pullback in experiments.
|
|
|
|
The initial implementation must perform no ray search, sparse assembly, or iterative solve during a Krylov `Mult`.
|
|
|
|
## 22. Migration plan
|
|
|
|
Each implementation phase must run the repository `./format` script before compilation. Tests and tags must use semantic physical or software meanings, never phase numbers.
|
|
|
|
### Phase A: terminology and contracts
|
|
|
|
- Rename the existing physical `SurfacePrescription` concept and accessors to explicitly mention `SurfaceCondition`.
|
|
- Add compile-time concepts and descriptors for surface deformation, interior extension, and vacuum extension.
|
|
- Add negative compile-time contract tests.
|
|
- Make no root block change yet.
|
|
|
|
### Phase B: scalar surface trace coordinates
|
|
|
|
- Add boundary-supported scalar DOF resolution.
|
|
- Construct the stellar-surface scalar trace map with deterministic parallel ordering.
|
|
- Add topology, schema, ordering, and MPI tests.
|
|
- Keep the current root displacement system unchanged.
|
|
|
|
### Phase C: `NodalRadialSurface`
|
|
|
|
- Compile fixed reference radial directions.
|
|
- Implement parameter-to-vector-trace forward and transpose actions.
|
|
- Add analytic, finite-difference, and adjoint tests.
|
|
- Keep the current root displacement system unchanged.
|
|
|
|
### Phase D: radial interior and vacuum extensions
|
|
|
|
- Build and validate reference radial correspondences.
|
|
- Implement `PowerLawRadialInteriorExtension` with power two.
|
|
- Implement `FixedInfinityRadialVacuumExtension`.
|
|
- Implement exact transpose paths.
|
|
- Add support, boundary, weight, continuity, determinant, and adjoint tests.
|
|
- Keep the current root displacement system unchanged.
|
|
|
|
### Phase E: composed prepared domain deformation
|
|
|
|
- Compose prescription, interior extension, vacuum extension, and merge ownership.
|
|
- Add cached full displacement and direction workspaces.
|
|
- Add dependency and statistics reporting.
|
|
- Test the complete lift and pullback independently of stellar physics.
|
|
|
|
### Phase F: stellar-model ownership
|
|
|
|
- Extend `StellarModel` to own the three geometry policy objects.
|
|
- Compile them against `FEM` during prepared-operator construction.
|
|
- Preserve stable addresses and move semantics.
|
|
- Add model ownership and runtime-view tests.
|
|
|
|
### Phase G: root block migration
|
|
|
|
- Replace the root displacement block with the surface-deformation parameter block.
|
|
- Generate full displacement during `Prepare`.
|
|
- Lift surface directions during `Mult`.
|
|
- Pull back the full mechanical residual and Jacobian action.
|
|
- Replace the root displacement dependency with surface deformation plus derived geometry revision.
|
|
- Remove the root centering-row replacement.
|
|
- Keep internal child displacement layouts unchanged.
|
|
|
|
This should be treated as one coherent migration rather than leaving two production root forms indefinitely.
|
|
|
|
### Phase H: coupled-test and experiment migration
|
|
|
|
- Update root state builders and block-layout tests.
|
|
- Update exact restricted-Jacobian and centered-difference tests.
|
|
- Replace obsolete root tangential/vacuum probes with reduced-space reachability and conditioning probes.
|
|
- Re-run analytic equilibrium and rotating-response coverage.
|
|
|
|
### Phase I: solver-readiness audit
|
|
|
|
- Measure reduced block sizes and action costs.
|
|
- Recheck gauge directions and homology.
|
|
- verify minimum mapped determinant under representative oblate deformations;
|
|
- define the first surface-block preconditioner approximation;
|
|
- confirm that the reduced operator is square and suitable for nonlinear solver wiring.
|
|
|
|
## 23. Risk register
|
|
|
|
### 23.1 Incorrect residual restriction
|
|
|
|
**Risk:** selecting force entries on the surface instead of applying \(L^T\) loses interior virtual work and produces an inconsistent reduced Jacobian.
|
|
|
|
**Mitigation:** require exact adjoint identities and compare the reduced root action with an independently lifted/pulled-back full action.
|
|
|
|
### 23.2 Missing nonlinear pullback term
|
|
|
|
**Risk:** a future shape-dependent prescription supplies \(D\mathcal D(q)\) but omits the derivative of its transpose, producing an inexact Newton operator.
|
|
|
|
**Mitigation:** reserve and test the pullback-derivative contract now; mark the initial policies explicitly linear.
|
|
|
|
### 23.3 Radial correspondence failure
|
|
|
|
**Risk:** rays fail to intersect the surface uniquely, particularly after mesh changes or for non-star-shaped reference geometry.
|
|
|
|
**Mitigation:** compile correspondences on the reference mesh, validate uniqueness, and fail construction with geometric diagnostics.
|
|
|
|
### 23.4 Mesh inversion
|
|
|
|
**Risk:** high-amplitude or high-frequency nodal deformation produces negative mapped determinants.
|
|
|
|
**Mitigation:** validate determinants during nonlinear preparation, globalize nonlinear steps, report the minimum determinant, and later add smoothing or a stronger extension policy.
|
|
|
|
### 23.5 Extension-dependent discrete equilibrium
|
|
|
|
**Risk:** different interior/vacuum extensions yield measurably different discrete surface shapes because the reduced virtual-work space changes.
|
|
|
|
**Mitigation:** retain full physical field equations, use consistent pullback, test refinement convergence, and later compare at least two extension policies.
|
|
|
|
### 23.6 Translation and orientation
|
|
|
|
**Risk:** anchoring the center removes exact translation but leaves weak translation-like or orientation-like surface patterns.
|
|
|
|
**Mitigation:** monitor center of mass, repeat reduced null-space probes, and add a surface-space gauge only if evidence requires it.
|
|
|
|
### 23.7 Parallel trace ownership
|
|
|
|
**Risk:** shared surface DOFs are duplicated, omitted, or inconsistently ordered across ranks.
|
|
|
|
**Mitigation:** reuse FieldDof synchronization principles and add explicit multi-rank topology and adjoint tests.
|
|
|
|
### 23.8 Compiler and runtime overhead
|
|
|
|
**Risk:** extensibility introduces virtual calls or allocation in per-element or per-quadrature loops.
|
|
|
|
**Mitigation:** compile prescriptions into prepared vector operators, dispatch once per complete operation, preallocate workspaces, and prohibit ray searches or assembly during `Mult`.
|
|
|
|
### 23.9 Over-restricting the physical shape family
|
|
|
|
**Risk:** a radial surface representation cannot describe a later non-star-shaped or materially multivalued surface.
|
|
|
|
**Mitigation:** keep the root block semantic (`surface deformation parameters`) independent of the concrete prescription and make the prescription replaceable at model construction.
|
|
|
|
### 23.10 Test churn masking regressions
|
|
|
|
**Risk:** broad root-vector changes cause established physics tests to be weakened unnecessarily.
|
|
|
|
**Mitigation:** leave child-operator full-displacement tests intact and change only tests whose subject is the root coordinate system or root composition.
|
|
|
|
## 24. Acceptance criteria for the phase
|
|
|
|
The phase is complete only when all of the following hold:
|
|
|
|
1. The root nonlinear state contains surface deformation parameters and no independent volume displacement block.
|
|
2. `NodalRadialSurface` satisfies the compile-time and runtime prescription contracts.
|
|
3. The first interior and vacuum extension policies satisfy their separate contracts.
|
|
4. The composed map supplies exact forward, derivative, and transpose operations.
|
|
5. No ray search, matrix assembly, allocation proportional to mesh size, or iterative extension solve occurs during a root Jacobian application.
|
|
6. The reduced mechanical residual is \(L^TR_d\), not a selected subset of force rows.
|
|
7. All existing child physics kernels remain reusable with full generated displacement vectors.
|
|
8. The physical constant-pressure surface condition remains EOS-mediated and independent of geometry representation.
|
|
9. Adjoint identities pass for the prescription, both extensions, and the complete map.
|
|
10. Complete reduced Jacobian actions match centered differences.
|
|
11. The spherical equilibrium accuracy remains at the established discrete projection floor.
|
|
12. Representative radial and quadrupolar deformations preserve positive mapped determinants within documented amplitude bounds.
|
|
13. Tangential and vacuum-only mesh directions are absent from the root nonlinear coordinate space.
|
|
14. Test tags and names are semantic and registered through `test_helpers`.
|
|
15. The full project test suite passes after migration.
|
|
|
|
## 25. Deferred prescriptions and extensions
|
|
|
|
### 25.1 `SphericalHarmonicSurface`
|
|
|
|
A future spherical-harmonic prescription can use
|
|
|
|
\[
|
|
R(\theta,\phi)
|
|
=R_0(\theta,\phi)
|
|
+\sum_{\ell,m}a_{\ell m}Y_{\ell m}(\theta,\phi).
|
|
\]
|
|
|
|
Its parameter vector consists of \(a_{\ell m}\), but it compiles to the same vector surface-trace contract. It may encode axisymmetry, reflection symmetry, translation removal, and angular truncation directly through its selected modes.
|
|
|
|
No root block or physics operator should need to change when it is added.
|
|
|
|
### 25.2 Modal interior and vacuum extensions
|
|
|
|
Spherical harmonics make analytic mode-dependent extensions possible. Regular interior radial factors and compactification-compatible exterior factors could replace the nodal ray interpolation while preserving the same prepared domain-deformation contract.
|
|
|
|
### 25.3 Elastic extensions
|
|
|
|
An elastic extension can use the already registered mesh-extension form and may be realized either as an eliminated map or an auxiliary sparse block. That choice should be based on measured solve cost and conditioning, not fixed by the prescription interface.
|
|
|
|
## 26. Final architectural boundary
|
|
|
|
The durable boundary should be:
|
|
|
|
```text
|
|
physical surface condition
|
|
ConstantPressureSurface + EOS
|
|
|
|
|
| independent of
|
|
v
|
|
surface deformation prescription
|
|
q -> vector displacement on StellarSurface
|
|
|
|
|
v
|
|
interior and vacuum extension policies
|
|
surface displacement -> full volume displacement
|
|
|
|
|
v
|
|
existing mapped physics operators
|
|
full displacement -> residuals and Jacobian actions
|
|
|
|
|
v
|
|
transpose pullback
|
|
full mechanical dual -> surface-parameter residual
|
|
```
|
|
|
|
This boundary permits future researchers to change how a stellar surface is represented and how the computational mesh follows it without changing the EOS, the physical surface condition, or the verified element-level physics.
|