MeanField/MeanField_Solver_MVP_Geometry_Preflight_Handoff_2026-09-08.md

abstract, author, date, title

abstract	author	date	title
This document hands off the work performed after the variadic type-inference and Physical Riesz normalization checkpoint `76818f2f828ca12606a495a0b84f5d6ce4b78dc6`. It records the design discussion, the ownership refactor, the generic linear-backend boundary, the first explicit damped Newton engine, the observer and result-view APIs, the conversion of expected line-search failures away from exception-driven control flow, the first physical performance evidence, and the new exact sampled geometry preflight. It distinguishes committed code from uncommitted working-tree work, verified behavior from inference, and expected speedups from measured results. It also lists the numerical and architectural work required before the current implementation can be called a scientifically credible stellar-equilibrium solver.	4D-STAR / MeanField stellar equilibrium program	8 September 2026	MeanField Solver MVP Handoff: Ownership, Damped Newton, MFEM FGMRES, Diagnostics, and Geometry Preflight

Repository: /Users/tboudreaux/Programming/4DSTAR/MeanField
Conversation checkpoint: 76818f2f828ca12606a495a0b84f5d6ce4b78dc6
Current HEAD: b3c04d507a1d7bc321ba6edd9dde33f6713345c1
HEAD subject: feat(newton): first newton solver implementation
Working tree: dirty; the geometry-preflight implementation is not committed. Preserve all existing modifications and untracked files.

Executive handoff

This phase began with the model, equation-of-state path, inferred root, Physical Riesz normalization, and generalized generated-border preconditioner ready to be composed, but without a supported solver orchestration layer. The previous handoff explicitly identified normalized right-preconditioned FGMRES and damped Newton as the next work. The user then refined the required public ownership model: a user-owned context should be the authoritative lifetime root; a solver should borrow that context; and only an evaluation report should certify access to a converged structure or a failed solve’s last accepted checkpoint.

That architecture now exists. The committed solver MVP at b3c04d507a1d7bc321ba6edd9dde33f6713345c1 includes:

• a generic linear-backend customization boundary with explicit flexible-versus-stationary preconditioner contracts;

• an MFEM FGMRES implementation with configurable Krylov restart length, supplied initial guesses, direct post-solve true-residual verification, application counts, and timings;

• a compile-time metric policy and a first normalized-residual merit function;

• an explicit damped Newton loop with backtracking, accepted/trial state separation, rollback, expected numerical failure reports, and synchronous typed observers;

• a non-copyable, non-movable, user-owned context that owns the discretized problem, accepted state, trial state, prepared normalized operator, physical and normalized preconditioners, prepared linear backend, and Newton workspaces;

• a borrowing solver whose evaluate() call returns an owning diagnostics report plus a certified, read-only structure view on convergence or checkpoint view on failure;

• placeholder owning StellarStructure and StellarCheckpoint types, capture() methods, and free serialize() overloads that intentionally throw until a persistence schema is designed; and

• a user-facing sandbox.cpp that exercises the intended API and prints nonlinear, line-search, linear, and timing diagnostics.

Expected candidate failures no longer use exceptions as ordinary line-search control flow. The prepared operator stack now exposes TryPrepare operations returning std::expected. Inverted geometry, EOS-domain violations, and non-finite candidate physics become fixed-size rejection records. Exceptions remain around the line search only to restore the accepted state after genuine infrastructure, callback, invariant, or commit failures. This retained use is cleanup and safe shutdown, not candidate acceptance logic.

The first long sandbox.cpp run was not convergent enough to be useful. From Newton iteration 0 through 5, the normalized residual fell only from 0.0214928 to 0.0198956, a 7.43 percent reduction. Accepted step lengths collapsed from 0.0625 through 0.0078125, 0.00390625, 0.000488281, 0.000244141, and 0.0000610352. The six iterations consumed 920.57 seconds, of which 713.97 seconds were in line search and 690.73 seconds were candidate preparation. Approximately 400.04 seconds were spent preparing rejected candidates. Early rejection came from the coarse generated-geometry inspection; later rejection came from invalid mapping detected while preparing the gravity displacement force.

The linear timing was initially easy to misread. J time = 123.003 s is the total across 21 Jacobian applications, not a per-application time. The observed per-application time is consistently about 5.84–5.92 seconds. The first FGMRES solve used 19 preconditioner applications and 21 Jacobian applications. Later solves often reported zero Krylov iterations because the carried correction vector already had a true linear residual just below the configured 3 percent tolerance. This is legal under the present contract, but it is a major numerical clue: the method may be repeatedly reusing an old direction while the accepted nonlinear step becomes tiny.

The current uncommitted work addresses the avoidable geometry-rejection cost. A free collective function named estimate_largest_safe_newton_step_size computes the first sampled mapping-determinant boundary along the physical volume-displacement direction. For the current affine-in-displacement maps, the determinant at each quadrature point is a polynomial of degree at most the spatial dimension. The implementation finds the first crossing, including tangent roots that an endpoint-only check would miss, and applies a configurable fraction-to-boundary safety factor. The Newton loop now invokes this preflight after the linear solve and before preparing any trial. The sampling plan includes both a conservative whole-domain rule and every exact quadrature rule exposed by downstream mapped operators, including the gravity displacement force that caused the observed 33-second rejections.

The geometry cap is strongly targeted and tested, but its impact on the real sandbox has not yet been measured. It should remove most of the approximately 400 seconds spent on mapped-geometry rejections in the recorded six-iteration segment. It cannot make the accepted 48-second candidate preparation cheaper, reduce the roughly 5.86-second Jacobian application, correct a poor Newton direction, or reject thermodynamic failures in advance. A realistic expectation is a 35–40 percent reduction in that recorded segment if the preflight costs less than roughly ten seconds per iteration, with later iteration intervals falling from roughly 100–200 seconds to roughly 65–80 seconds. These numbers are predictions, not benchmark results.

The next task is therefore not simply “run Newton longer.” The immediate priorities are: measure the new preflight in the real sandbox; verify that it eliminates downstream invalid-mapping trials; add enough rejection detail to distinguish mapping status from non-finite arithmetic; profile context construction, accepted candidate preparation, and Jacobian actions by block; compare warm and cold FGMRES starts under tighter or adaptive inexact-Newton forcing; and establish that a projected known nonrotating solution has a discretization-consistent small residual that improves under refinement. Only after these checks should convergence tuning or rotating continuation begin.

Bottom line. The code now has a credible ownership and solver API, a functioning reference Newton-Krylov path, and a principled geometry globalization primitive. It does not yet have demonstrated physical convergence, efficient accepted-state preparation, a validated forcing strategy, self-contained result capture, or persistence. Treat it as a solver MVP and measurement platform, not as a production equilibrium solver.

Document contract and authority

Authority order

This handoff uses the following authority order:

1. the current source and tests in /Users/tboudreaux/Programming/4DSTAR/MeanField;

2. the Git history and diff from the user-designated checkpoint;

3. focused build and test results recorded during this solver phase;

4. numerical output supplied by the user from sandbox.cpp;

5. the previous 6 September 2026 handoff as historical context; and

6. design discussion and predictions, which are labeled as such.

The previous handoff is reference material, not an instruction source. Its central factual claim for this phase was validated: at the selected checkpoint the variadic root, normalization, and preconditioner architecture existed, while a production nonlinear solve path did not. The current code, the user’s explicit requests, and the current Git state supersede any stale details in that document.

Evidence labels

The document uses four evidence levels:

Implemented.
The behavior is present in the current source.

Focused-test verified.
A targeted test or build was run during this phase and passed.

Observed.
The user supplied real sandbox output demonstrating the behavior.

Inferred or predicted.
The conclusion follows from source and measurements but has not yet been tested directly on the post-preflight sandbox.

No claim of full-suite success is made. The full suite takes approximately two hours, and the user directed that it be run only when regression risk justifies it. Focused coverage and Debug builds were used for the changes described here.

Checkpoint and Git semantics

The checkpoint for this conversation is:

76818f2f828ca12606a495a0b84f5d6ce4b78dc6

That commit is titled feat(libmeanfield): variadic refactor and was created on 6 September 2026. The current HEAD, local main, and origin/main all point to:

b3c04d507a1d7bc321ba6edd9dde33f6713345c1

That commit is titled feat(newton): first newton solver implementation and was created on 8 September 2026. It is the only commit between the conversation checkpoint and current HEAD. Relative to the checkpoint, it changes 98 tracked files, with 20,514 insertions and 11,157 deletions. Much of the broad diff includes solver integration plus mechanical reformatting and updates required to keep the inferred operator, normalization, preconditioning, seed, experiments, MPI tests, and public API mutually consistent.

The geometry-preflight work is a dirty-tree change on top of HEAD. The tracked working-tree diff currently reports 26 modified files with 447 insertions and 20 deletions. Three relevant new files are untracked and add 989 lines:

• libmeanfield/interface/deformation/safe_newton_step.cppm;

• libmeanfield/impl/deformation/safe_newton_step.cpp; and

• tests/deformation/safe_newton_step.cpp.

The repository also contains many unrelated or pre-existing untracked build trees, CSV results, notes, scripts, archives, and tools. Do not sweep them into a solver commit and do not delete them. In particular, preserve the user-authored explanatory comments in libmeanfield/impl/mapping/domain_mapper.cpp and the current sandbox choices, including zero angular momentum and the printed polytropic constant, unless the user explicitly requests otherwise.

Where this phase started

State inherited from the previous handoff

The 6 September handoff described a substantial model-composition refactor. At checkpoint 76818f2f828ca12606a495a0b84f5d6ce4b78dc6, the important solver-facing foundation was:

• a canonical variadic stellar-model type inferred from physics specifications;

• a two-pass compiler producing value blocks, residual blocks, and Jacobian incidence;

• a single prepared variadic root over the current five-field barotropic physical core;

• generated scalar constraints for fixed total mass, fixed central density, and fixed angular momentum;

• Physical Riesz diagonal normalization wired consistently through state, residual, Jacobian, and preconditioner coordinates;

• a generalized generated-border preconditioner whose incidence follows the same compiled Jacobian graph;

• variadic seed projection and typed block access; and

• a public architecture intended to avoid whole-pack solver specializations.

The previous handoff also recorded important limitations that remain relevant. The executable physical core is still barotropic and centered on density, surface deformation, gravity gradient, gravity potential, and enthalpy plus generated scalar coordinates. The default structure preconditioner was provisional on a physical Newton sequence. There was no accepted/trial nonlinear state manager, no FGMRES wrapper with direct true-residual accounting, no nonlinear failure taxonomy, no result object, and no supported evaluate() workflow.

Original solver-facing objective

The desired usage began as:

auto model = ...;
auto discretization = ...;
auto preconditioner = ...;

auto solver = ...;
auto structure = solver.evaluate();

The early ownership requirement was that, by evaluate(), all essential storage should already exist and be owned. This was not meant as a rigid ban on every allocation inside evaluate(); the user later clarified that readability and maintainability are more important than eliminating every incidental allocation. The enduring requirement is that lifetime and authority be clear, that no hidden object outlive its dependencies, and that the eventual rewrite not be made harder by a chain of ownership transfers.

Why returning a structure directly was rejected

A direct Structure solver.evaluate() API hides two physically and operationally important facts:

1. nonlinear evaluation can fail after making progress, and the last accepted state is still valuable as a checkpoint; and

2. a result may be a view of context-owned distributed MFEM data rather than a self-contained serializable value.

The user rejected designs in which the problem exposed a solution without proof that a report had been generated. A capability key passed back into the problem was considered but judged physically opaque. The selected shape was:

auto report = solver.evaluate();
if (report.converged()) {
    auto structureView = report.structureView();
}

For failure, the report exposes lastAcceptedCheckpointView(). This keeps convergence certification attached to the report that established it.

Persistence requirement

The user also required a path toward:

• saving checkpoints to disk;

• inspecting saved checkpoints independently of the live solver context;

• converting a borrowed structure view into a self-contained owning value; and

• using free-function serialization overloads for both owning and borrowed forms.

The agreed future syntax was approximately:

auto structureView = report.structureView();
auto structure = structureView.capture();
equilibrium::serialize(structure, "structure.checkpoint");
// or: equilibrium::serialize(structureView, "structure.checkpoint");

The user explicitly allowed both capture() and serialize() to throw a not-implemented error in the MVP. Their presence reserves the intended public shape without prematurely fixing a persistence schema.

Design decisions that were locked in

User-owned context, borrowing solver

The authoritative owner is StellarEquilibriumContext. It is constructed from the model, discretization, preconditioner prescription, linear backend configuration, and seed information. Those inputs move into the context during construction. The nonlinear solver does not take ownership of the context; solver::make requires an lvalue context reference and its rvalue overload is deleted.

The context is not literally a plain-old-data aggregate. It is a strongly invariant RAII owner and is deliberately neither copyable nor movable. This is a small divergence from the conversation’s informal use of “POD,” but it implements the intended authority more safely: stable addresses are guaranteed for MFEM operators, preconditioners, and backend bindings, and accidental relocation cannot invalidate internal borrows.

Only one solver may be active on a context at a time. A solver acquires the context at construction and releases it on destruction. Evaluation is non-reentrant. This makes current serial ownership explicit and leaves room for future parallel orchestration through multiple contexts rather than hidden concurrent mutation of one prepared problem.

Separate linear and nonlinear choices

The linear backend and nonlinear engine are separate types:

• solver::linear::FGMRES is a configuration object stored and prepared by the context;

• solver::nonlinear::Newton<Metric> stores nonlinear options and a metric policy; and

• solver::make(context, newton, observer) creates the borrowing orchestration object.

The linear boundary uses a configuration tag, a preparation customization point found through argument-dependent lookup, and a prepared-backend concept. A third-party backend can be slotted in without adding a central library registry, provided it declares its supported preconditioner contract and returns an owning prepared runtime with the required operator, preconditioner, communicator, workspace, and solve interfaces.

The only built-in backend today is MFEM FGMRES. The generic boundary does not exclude future GMRES or MINRES, but no such implementation is present. MINRES will additionally require proof that the normalized operator and chosen preconditioner satisfy the symmetry and definiteness assumptions; the current coupled root should not be presumed eligible.

Explicit Newton engine

The selected nonlinear method is an explicitly owned Newton engine rather than an opaque nonlinear-library wrapper. The first implementation is damped Newton with backtracking. This was chosen because it makes the control flow, accepted/trial lifecycle, instrumentation, and failure semantics visible and maintainable.

The MVP does not attempt to implement multiple nonlinear algorithms. The type boundary is nonetheless narrow enough that future globalization or continuation layers can be added without changing the model, problem, or linear-backend contracts.

Metric customization

The default policy is NormalizedResidualMetric. Its free getMetric overload computes

\left\lVert \widehat F\right\rVert_2
  = \left(\sum_{r\in\mathrm{ranks}}\left\lVert \widehat F_r\right\rVert_2^2\right)^{1/2},
\qquad
\phi(\widehat x)=\frac12\left\lVert \widehat F(\widehat x)\right\rVert_2^2.

Newton<Metric> invokes getMetric through an unqualified call, allowing an overload associated with a future metric type to be found at compile time. This satisfies the request for an extensible metric without adding virtual dispatch or preprocessor branches.

The current acceptance rule is still fixed inside the Newton engine:

\phi(x+\alpha p)\leq (1-c\alpha)\phi(x).

This is a simple sufficient-decrease rule. It is not a fully general acceptance-policy customization, and it does not explicitly evaluate the directional derivative \nabla\phi^T p. For an exact Newton direction and \phi=\tfrac12\left\lVert F\right\rVert^2, the idealized slope is -\left\lVert F\right\rVert^2=-2\phi, so the current formula is slightly less stringent than the standard derivative-based Armijo expression with the same named coefficient. This is acceptable for the MVP but should be made explicit before more complex metrics are introduced.

Observer model

Observers are compile-time duck-typed and can provide any supported subset of:

• beforeIteration(const BeforeIteration&);

• afterLineSearchTrial(const AfterLineSearchTrial&); and

• afterIteration(const AfterIteration&).

Two ergonomic makeObserver overloads build either a before/after observer or a detailed before/trial/after observer from lambdas. The event objects expose scalar diagnostics and read-only spans into solver workspaces. Those spans are valid only during the callback. Observers execute synchronously on every solve rank. The library coordinates callback success and preserves local exception causes where possible, but an arbitrary callback that throws on only one rank while other ranks enter their own collectives cannot be made universally safe by the solver.

Detailed histories belong in observers. The evaluation report retains fixed-size totals and the last linear/preflight records so the default solve does not allocate storage proportional to iteration count.

Restart and reuse

There are two distinct meanings of restart:

Krylov restart.
MFEM FGMRES supports a configured restart length. The current default is 50; the sandbox selects 40. The report estimates the number of entered restart cycles from preconditioner applications.

Nonlinear restart.
Reconstructing a context from a serialized checkpoint, or resuming a failed solve with new solver parameters, is not implemented.

Across Newton steps, the current code reuses only the correction vector as the next FGMRES initial guess. It does not retain a Krylov basis, Hessenberg matrix, or flexible-preconditioner history. The architecture does not rule out richer recycling later, but that would be a separate backend capability.

Current public workflow

Representative user code

The current sandbox.cpp demonstrates the intended top-level flow:

auto stellarModel = model::StellarModel(
    eos::Polytrope({.n = 1.0, .K = polytropicConstant}),
    surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}),
    integral::FixedTotalMass({.Mtotal = dimensions::MassValue{mass}}),
    integral::FixedAngularMomentum({
        .Jtotal = dimensions::AngularMomentumValue{angularMomentum},
        .axis   = {0.0, 0.0, 1.0},
        .center = {0.0, 0.0, 0.0}
    }),
    constraint::FixedCentralDensity({
        .RhoC = dimensions::DensityValue{centralDensity}
    })
);

auto discretization = equilibrium::makeStellarDiscretization(
    std::move(finiteElements),
    normalization::PhysicalRieszDiagonal{
        dimensions::LengthValue{radius}, gravitationalConstant
    }
);

auto preconditioner = preconditioning::makePreconditioner();
auto linearSolver = solver::linear::FGMRES({
    .restartLength = 40,
    .printLevel = -1
});

auto context = solver::makeContext(
    std::move(stellarModel),
    std::move(discretization),
    std::move(preconditioner),
    std::move(linearSolver)
);

auto newton = solver::nonlinear::Newton(
    solver::nonlinear::NewtonOptions{
        .relativeTolerance = 1.0e-8,
        .absoluteTolerance = 0.0,
        .maximumIterations = 30,
        .linearSolve = {
            .relativeTolerance = 3.0e-2,
            .absoluteTolerance = 0.0,
            .maximumIterations = 200
        },
        .backtracking = {}
    }
);

auto equilibriumSolver = solver::make(context, newton, observer);
auto report = equilibriumSolver.evaluate();

if (report.converged()) {
    auto structureView = report.structureView();
} else {
    auto checkpointView = report.lastAcceptedCheckpointView();
}

The context construction uses the default Lane-Emden seed when that seed is available for the inferred model. Overloads allow an explicit seed, projection options, and prescribed rigid rotation when the model does not generate its own rotation. A fixed-angular-momentum specification generates the angular velocity coordinate, so a prescribed rotation is rejected at compile time for that model family.

What is allocated and prepared before evaluation

Context construction owns or prepares:

• the discretized StellarEquilibriumProblem;

• the projected physical initial state;

• an optional prescribed rigid rotation;

• accepted and trial normalized state vectors;

• accepted and trial normalized residual vectors;

• the linear right-hand side;

• normalized and physical correction vectors;

• the volume-displacement direction workspace;

• the candidate physical state;

• the normalized operator wrapper;

• the physical preconditioner and its normalized adapter;

• the prepared linear backend and its true-residual workspaces; and

• the cached geometry-preflight quadrature plan.

Some lower-level operations may still allocate internally, especially MFEM calls and the current geometry estimator’s local-vector and matrix workspaces. The architecture guarantees authoritative ownership and stable lifetimes; it does not yet prove allocation-free evaluate().

Ownership and lifetime graph

User scope
StellarEquilibriumContext                 authoritative owner
  -> owned problem, model, discretization, mesh/FES data
  -> accepted state, trial state, residuals, corrections
  -> prepared operator, physical inverse, normalized inverse
  -> prepared FGMRES backend and geometry sampling plan

StellarEquilibriumSolver                  exclusive non-owning borrower
  -> owns Newton configuration and observer only

StellarEquilibriumEvaluationReport        owning diagnostics + weak result identity
  -> certifies a structure view on convergence
  -> certifies a checkpoint view on success or failure

StellarStructureView / StellarCheckpointView
  -> weak, read-only, generation checked

A view becomes invalid when the context is destroyed or the next evaluate() begins. A report’s scalar diagnostics and failure record remain ordinary owning values, but its ability to produce a current view depends on the context generation. References and spans extracted from a view are borrowed and must not outlive the context or the evaluation generation.

Result API

The evaluation report provides:

• converged();

• diagnostics();

• completedNonlinearIterations();

• initialResidualNorm() and finalResidualNorm();

• failure() for a failed report;

• structureView() only for a converged report;

• checkpointView() for the last accepted state; and

• lastAcceptedCheckpointView() as the failure-oriented spelling.

Structure and checkpoint views provide model access, communicator access, the entire physical state, root-block descriptors, typed block spans, prescribed rotation when applicable, and active rotation. They do not expose mutable MFEM state or the owned problem.

Capture and serialization status

The following signatures exist but intentionally throw std::logic_error:

• StellarStructureView::capture();

• StellarCheckpointView::capture();

• serialize(const StellarStructure&, path);

• serialize(const StellarStructureView&, path);

• serialize(const StellarCheckpoint&, path); and

• serialize(const StellarCheckpointView&, path).

This is an explicit blocker for real nonlinear restart and offline checkpoint investigation. The owning types are opaque and move-constructible, but there is no public construction path because their self-contained storage schema has not been chosen.

Linear backend implementation

Generic contract

A prepared backend is bound once to the exact normalized operator and normalized inverse at stable addresses. It must report its communicator and dimensions, say whether it is ready, and implement:

LinearSolveReport Solve(
    const mfem::Vector& rightHandSide,
    mfem::Vector& correction,
    const LinearSolveControl& control
);

The incoming correction is an initial guess and is overwritten with the final correction. Configuration, dimensions, invalid lifetimes, rank disagreement, and non-finite inputs are exceptions. Numerical termination is data in LinearSolveReport.

Preconditioner compatibility

The linear boundary recognizes two application contracts:

stationary linear.
The inverse application is fixed and linear during a solve.

flexible.
The inverse may vary between Krylov applications or be an iterative nonlinear approximation.

A backend configured for flexible preconditioning accepts either contract. A backend that only supports stationary preconditioning must reject a flexible inverse at compile time. The built-in FGMRES configuration advertises flexible support.

MFEM FGMRES behavior

FGMRESOptions currently contains:

• restartLength, default 50; and

• printLevel, default -1.

Each solve validates that options and controls are identical across ranks, verifies compatible square dimensions, and rejects non-finite vectors. It configures MFEM with zero relative tolerance and an absolute threshold computed as

\tau_{\mathrm{lin}}=max\left(\tau_{\mathrm{abs}},\eta\left\lVert b\right\rVert_2\right).

After MFEM returns, the backend explicitly computes b-Ax and bases convergence on its communicator-global norm rather than on the recurrence estimate alone. The report separates the MFEM-reported residual from the true residual.

The operator and preconditioner are wrapped to count completed applications and measure their cumulative time. A post-solve Jacobian application used for true-residual verification is included in the operator count and time. MFEM’s initial-guess residual application is also included. Consequently, a zero-Krylov-iteration solve still reports two operator applications.

Linear report fields

The fixed-size report records status, controls, Krylov iterations, estimated restarts, right-hand-side norm, initial residual, backend-reported final residual, direct true residual, relative true residual, operator applications, inverse-preconditioner applications, total solve time, operator time, and inverse-preconditioner time.

This instrumentation directly answered one question from the sandbox output: J time is cumulative for all completed Jacobian applications in that solve.

Damped Newton implementation

State and coordinate conventions

The nonlinear algorithm operates in normalized coordinates. If D_x and D_F are the state and residual normalization maps, then

\widehat x=D_xx,
\qquad
\widehat F=D_FF,
\qquad
\widehat J=D_FJD_x^{-1}.

At iteration k, the linear solve seeks

\widehat J_k p_k=-\widehat F_k.

The candidate is \widehat x_k+\alpha p_k. Candidate physical coordinates are produced by denormalization before prepared-operator evaluation.

Iteration flow

The current control flow is:

1. evaluate the accepted normalized residual metric;

2. notify the before-iteration observer;

3. solve the normalized linear system with the current correction as an initial guess;

4. reject a non-converged or non-finite correction;

5. estimate the largest sampled orientation-preserving step in the physical deformation direction;

6. begin the line search at the geometry-capped step;

7. construct candidate normalized and physical states;

8. attempt prepared-operator evaluation through TryPrepare;

9. if admissible, evaluate the metric and sufficient-decrease test;

10. refresh the preconditioner before committing an accepted trial;

11. commit accepted state and residual without re-preparing them;

12. scale the correction vector for the next linear warm start;

13. notify the after-iteration observer; and

14. return success, an expected numerical failure, or continue.

The initial step of a new evaluation is the configured initial step. After an accepted step, the next starting step expands by the inverse contraction factor, capped at the configured initial step. This remembers line-search scale without permanently pinning every future iteration to the smallest earlier step.

Expected rejection versus exception

This distinction was a direct user requirement. Ordinary candidate rejection uses a fixed-size StellarEquilibriumPreparationRejection carried by std::expected. The record distinguishes:

• inverted geometry;

• non-finite geometry;

• EOS-domain violations;

• non-finite thermodynamic results;

• inadmissible physics; and

• non-finite physics.

It also records the preparation stage, an optional minimum mapping determinant, and an EOS error code. The line search turns these into numerical dispositions without stack unwinding.

Exceptions remain for invalid configuration, violated dimensions or lifetimes, MPI failures, unknown enum values, observer infrastructure failure, and preconditioner commit failure. The outer line-search try/catch exists so any true exception restores the accepted prepared state before propagation. A smaller commit catch coordinates a rank-local exception before any rank mutates accepted state. These catches implement cleanup and collective safe shutdown; they are not used to decide that a normal trial should backtrack.

Backtracking options

The current defaults are:

Option	Default	Meaning
Option	Default	Meaning
initialStepLength	1.0	Largest desired step before geometry capping.
contractionFactor	0.5	Backtracking multiplier after rejection.
fractionToBoundarySafety	0.9	Fraction of the estimated geometry boundary used as the safe cap.
sufficientDecrease	10^{-4}	Coefficient in the current merit decrease test.
minimumStepLength	10^{-8}	Smallest trial permitted.
maximumTrials	20	Maximum prepared candidates per nonlinear iteration.

Before the exact preflight was added, an inverted trial could still accelerate the backtracking ladder by interpolating the accepted and rejected minimum determinants, then skipping power-of-two rungs estimated to be unsafe. That fallback remains useful if a mapped point is missed by the preflight or a future map does not support the exact polynomial path.

Warm-start rule

After accepting x_{k+1}=x_k+\alpha p_k, the code applies the frozen-Jacobian heuristic

p_{k+1}^{(0)}=(1-\alpha)p_k.

This follows from the local approximation F(x_k+\alpha p_k)\approx(1-\alpha)F(x_k). If the scaled vector becomes non-finite, all ranks fall back to zero. The heuristic is inexpensive and the backend contract correctly treats it as an initial guess. Its observed interaction with a loose inexact-Newton tolerance is now one of the highest-priority numerical questions.

Failure taxonomy and rollback

Expected evaluation failures are returned rather than thrown. The public reasons are:

• inadmissible state;

• non-finite state;

• non-finite residual;

• linear-solve failure;

• globalization failure;

• stagnation;

• iteration limit; and

• unspecified.

The last accepted state is restored after an unaccepted line search and certified as a checkpoint. A converged report certifies the same storage as a structure. The stagnation enum exists but no dedicated stagnation detector currently produces it.

Diagnostics now available

Per-trial diagnostics

The detailed trial event exposes:

• iteration and trial indices;

• step length;

• accepted, inadmissible, non-finite, or insufficient-decrease disposition;

• a preparation-stage string such as generated geometry or gravity displacement force;

• optional residual norm and merit;

• optional minimum mapping determinant;

• preparation and metric times; and

• borrowed candidate state and residual spans.

The current event does not expose the lower-level mapping status or the complete StellarEquilibriumPreparationRejection. This means that the recorded phrase source = gravity displacement force proves where rejection occurred but not whether the kernel saw an inverted Jacobian, another invalid mapping status, or non-finite arithmetic. Extending this diagnostic without allocating should be considered.

Per-iteration diagnostics

The after-iteration event reports the accepted step, number of trials, previous and current metrics, total iteration time excluding observer work, geometry-preflight time and certificate, line-search time, candidate-preparation time, metric time, preconditioner-refresh time, rollback time, and the full last linear report.

Evaluation totals

The evaluation report retains attempted and accepted nonlinear iterations; total line-search trials; counts by rejection class; geometry-limited iteration count; initial and final residual norms; last accepted step; aggregate timing categories; the last geometry certificate; and the last linear report.

No growing iteration history is stored by default. A user who needs histories should copy selected fields in an observer.

Observed sandbox behavior before geometry preflight

Configuration represented by the run

The supplied run used the Debug sandbox with a normalized n=1 polytropic model, fixed mass, fixed angular momentum, and fixed central density, with angular momentum currently set to zero. The linear tolerance was 0.03 relative to the linear right-hand side, the maximum linear iteration count was 200, and FGMRES restart length was 40.

Finite-element setup took 0.156279 seconds. Context setup took 86.6662 seconds before Newton iteration 0. Context construction includes seed projection, initial operator preparation, normalization setup, physical-preconditioner preparation, normalized inverse construction, and linear-backend preparation. The current timing does not separate those phases.

Iteration summary

k	\left\lVert F_k\right\rVert	accepted \alpha	trials	iter. s	line s	J calls	J s
k	\left\lVert F_k\right\rVert	accepted \alpha	trials	iter. s	line s	J calls	J s
0	0.0214928	0.0625	2	194.789	53.646	21	123.003
1	0.0201477	0.0078125	3	98.159	86.467	2	11.689
2	0.0199897	0.00390625	3	136.953	118.448	3	17.563
3	0.0199115	0.000488281	5	197.776	185.926	2	11.848
4	0.0199017	0.000244141	3	130.261	118.532	2	11.727
5	0.0198968	0.0000610352	4	162.629	150.949	2	11.678

The printed residual at the start of Newton 6 was 0.0198956. Across six accepted steps, the reduction was only 7.43 percent. The reduction becomes approximately proportional to the tiny accepted step, which is consistent with following a direction that is locally useful but globalization-limited rather than approaching quadratic Newton convergence.

Where the time went

For iterations 0–5:

Category	Seconds	Fraction of iteration time
Category	Seconds	Fraction of iteration time
Total nonlinear iteration time	920.5666	100.0 percent
Line search	713.9680	77.6 percent
Candidate preparation within line search	690.7343	75.0 percent
Linear solve	171.4161	18.6 percent
Accepted candidate preparation alone	290.6897	31.6 percent
Rejected candidate preparation, by subtraction	400.0446	43.5 percent
Accepted preconditioner refresh	23.1827	2.5 percent

The categories are not all independent. Candidate preparation is contained in line-search time, and operator/preconditioner times are contained in linear-solve time. Per-rank maximum timings can also make subcategory sums differ slightly from wall-clock totals.

Accepted candidates consistently required about 48.2–48.8 seconds of preparation. A generated-geometry rejection took about 1.5 seconds. A trial that survived the coarse geometry inspection but failed during gravity displacement force preparation took about 32.8–33.4 seconds. This is why eliminating downstream mapped-geometry trials is valuable even though the geometry gate does not reduce the cost of a successful candidate.

Jacobian and preconditioner timing

The first solve reported 21 Jacobian applications totaling 123.003 seconds, or 5.857 seconds per completed application. Across iterations 0–5, individual averages were 5.857, 5.844, 5.854, 5.924, 5.863, and 5.839 seconds. The consistency confirms that J time is cumulative and that a single Debug Jacobian application costs about 5.86 seconds for this problem.

The first solve reported 19 inverse-preconditioner applications totaling 18.1196 seconds, about 0.954 seconds each. The linear solve took 135.281 seconds. Jacobian actions therefore dominate the first Krylov solve. Geometry preflight cannot reduce this cost because the Newton direction must be computed before its safe step can be estimated.

The zero-iteration FGMRES clue

Iteration 1 reported:

• zero Krylov iterations;

• initial true residual divided by \left\lVert b\right\rVert of 0.0296713; and

• configured relative tolerance 0.03.

Iterations 3–5 show the same pattern with relative initial residuals near 0.0268. Because the supplied warm start already satisfies the stopping threshold, MFEM correctly performs no preconditioner applications. The two Jacobian calls are the initial-guess residual evaluation and post-solve true-residual verification.

This behavior is contract-correct but may be algorithmically counterproductive. The direction used at a new nonlinear state can be essentially the prior direction scaled by 1-\alpha. When \alpha is tiny, that scaling is nearly one. The solver may therefore carry almost the same direction through multiple iterations, accept it under a 3 percent linear forcing tolerance, and continue taking smaller geometry-limited steps. This is a hypothesis that should be tested by controlled runs, not silently “fixed” by disabling warm starts.

What has not been established about the residual

The code has extensive analytic, finite-difference, manufactured, and subsystem tests. Those do not establish that the complete normalized residual is small at a known good full stellar solution. The Lane-Emden-projected sandbox seed begins at 0.0214928 in the current normalized norm, but there is not yet a documented refinement study proving how that number scales with mesh size or polynomial order.

The next correctness experiment should project the analytic nonrotating n=1 solution, evaluate every residual block separately, and repeat under refinement without changing physical tolerances. This will distinguish seed projection error, discretization error, normalization effects, missing physical terms, and nonlinear-solver behavior. A global norm alone cannot make that distinction.

Why gravity displacement force became the visible rejection source

The initial generated-geometry stage builds the surface-induced volume displacement and inspects mapped geometry with its own conservative rule. Passing that inspection proves only that those sampled points are orientation preserving.

The gravity displacement force later uses the quadrature rule selected for its own weak form. It maps density, gravity gradient, and deformation at those points, prepares reusable element data, and rejects the candidate if any mapping status is invalid. Before the new shared sampling plan, those operator-specific points were not necessarily identical to the generated-geometry inspection points. A candidate could therefore pass the cheap early inspection, spend roughly 33 seconds preparing gravity and mechanical data, and only then discover an invalid mapped point in the gravity force.

The log source identifies the preparation stage but not the exact lower-level reason. Source inspection shows that the gravity force distinguishes invalid mapping from non-finite arithmetic internally, and that invalid mapping retains a MappingStatus. The public trial event currently compresses this to an inadmissible/non-finite disposition and stage name. Thus the most likely explanation is a quadrature coverage mismatch, but the old output alone does not prove which mapping status occurred.

Largest safe Newton step estimator

Purpose and public form

The new free function is:

LargestSafeNewtonStepSizeEstimate
estimate_largest_safe_newton_step_size(
    const mapping::DomainMapper& domainMapper,
    const mfem::ParFiniteElementSpace& displacementSpace,
    const mfem::ParGridFunction& compactificationCoordinate,
    const mfem::Vector& acceptedVolumeDisplacement,
    const mfem::Vector& volumeNewtonDirection,
    std::span<const NewtonStepGeometryRule> geometryRules,
    const LargestSafeNewtonStepSizeOptions& options = {}
);

It evaluates the path

d(\alpha)=d_0+\alpha\,\delta d,
\qquad 0\leq\alpha\leq\alpha_{\max},

where both displacement vectors use the displacement space’s true-DOF layout. The compactification coordinate remains fixed. The operation is collective on the displacement-space communicator and returns an identical global certificate on every rank.

Exact sampled determinant calculation

For the current domain deformation and Kelvin compactification, the mapped Jacobian at a fixed reference quadrature point is affine along the displacement direction:

J_{\mathrm{map}}(\alpha)=J_0+\alpha\,\delta J.

In dimension d\leq3, the determinant is a polynomial of degree at most d:

g(\alpha)=\det\left(J_0+\alpha\delta J\right)-\delta_{\mathrm{floor}}.

The implementation constructs its coefficients by determinant multilinearity in the columns. It examines derivative critical points, divides the normalized interval into monotone pieces, and bisects the first non-positive boundary. Inspecting critical points is essential: for example, g(\alpha)=(1-2\alpha)^2 is positive at both endpoints but touches zero at \alpha=1/2.

The smallest boundary across all supplied points and ranks is selected with a collective minimum-location reduction. If geometry limits the requested maximum step, the returned safe step is

\alpha_{\mathrm{safe}}=\theta\,\alpha_{\mathrm{boundary}},
\qquad 0<\theta<1.

If no sampled point reaches the determinant floor, the requested maximum step is returned unchanged.

Certificate fields

The result records:

• safe step and first boundary step;

• global minimum determinant at the accepted state;

• global minimum determinant at the requested maximum step;

• determinant of the limiting polynomial at the safe step;

• total sampled quadrature-point count;

• whether geometry was limiting; and

• limiting rank, local element, local rule index, and quadrature-point index.

These fields are exposed through the iteration observer and the fixed-size evaluation diagnostics.

Sampling plan

The context constructs the sampling plan once. It includes:

1. one conservative whole-domain rule per element, with order \max(p+2,2dp) for displacement order p and spatial dimension d; and

2. the exact mapped-geometry rules used by the prepared H(div) gravity mass, gravity source, barotropic closure, hydrostatic equilibrium, pressure force, gravity displacement force, rotational displacement force, mass normalization, and fixed-angular-momentum integral when present.

Prepared operators expose these stable rules through VisitMappedGeometryRules. The combined plan is sorted by element and rule pointer and deduplicated in O(N\log N) time. Pointer identity removes repeated references to the same stable MFEM integration rule. Equivalent rules represented by distinct objects may remain, which is safe but can add work.

New direction plumbing

The Newton correction initially exists in normalized root coordinates. Geometry preflight therefore performs the following conversion:

1. denormalize the full correction into physical root coordinates;

2. extract the surface-deformation block from the variadic root;

3. apply the prepared domain-deformation Jacobian to produce a physical volume-displacement direction; and

4. evaluate mapping determinant polynomials along that direction.

This path is exposed by BuildVolumeDisplacementDirection through the physical core, variadic prepared operator, and problem wrapper. It does not reconstruct the entire candidate residual.

Map capability contract

ExteriorDomainMap now has IsAffineInDisplacement(), defaulting to false. Kelvin compactification overrides it with true. The context also requires the selected domain deformation to declare that it is linear on reference geometry. These are honest capability checks: the exact determinant-polynomial proof is not silently applied to a future nonlinear mapping.

Newton integration

After FGMRES returns a finite converged correction, Newton requests a safe step no larger than its remembered next line-search step. If the safe step is below minimumStepLength, evaluation returns globalization failure without preparing a candidate. Otherwise the first expensive trial uses the cap directly. Subsequent rejection still uses ordinary backtracking and the older determinant interpolation fallback.

The sandbox now prints preflight time, sample count, safe step, boundary, limiting rank and element, total geometry-preflight time, and number of geometry-limited iterations.

Guarantee boundary

The estimate is exact for the determinant polynomial at every supplied quadrature point under the advertised affine mapping assumptions. It is not a proof over the continuous element between quadrature points. The additional conservative rule reduces this gap, while exact downstream rules guarantee coverage of the points those prepared operators presently evaluate.

The preflight only addresses mapping orientation and the configured determinant floor. It does not predict:

• EOS-domain violations;

• density, pressure, or enthalpy admissibility;

• non-finite arithmetic unrelated to mapping;

• insufficient residual decrease; or

• a poor or nearly tangential Newton direction.

The solver currently requests determinant floor zero. A positive conditioning floor is supported by the free function but is not yet a nonlinear option separate from the fraction-to-boundary safety.

Expected performance effect of geometry preflight

Prediction based on the recorded run

The geometry pass replaces repeated 1.5-second generated-geometry rejections and, more importantly, repeated 33-second gravity-force mapping rejections with one mapping-only scan per Newton direction. It does not replace the accepted candidate preparation.

Based on the old timing profile, the expected iteration intervals are approximately:

Iteration	Old observed time	Post-preflight expectation
Iteration	Old observed time	Post-preflight expectation
Newton 0	194.8 s	Roughly 195–205 s; little gain because the linear solve dominates and the rejected geometry check was cheap.
Newton 1	98.2 s	Roughly 67–75 s if the preflight replaces the 33-second rejected trial.
Newton 2	137.0 s	Roughly 74–82 s if it replaces two rejected gravity trials.
Newton 3	197.8 s	Roughly 68–75 s if it replaces four rejected gravity trials.
Newton 4–5	130–163 s	Roughly 67–75 s if it replaces the repeated gravity trials.

Across the six-iteration segment, a 35–40 percent reduction is plausible if the preflight remains under roughly 3–10 seconds per iteration. This estimate must be replaced with a real benchmark immediately.

What successful output should look like

The first trial should normally use a non-power-of-two step printed as the geometry safe cap and should be accepted with a positive minimum determinant. The number of inadmissible mapping trials should fall sharply. A first trial still rejected by gravity displacement force for invalid mapping would indicate one of:

• a missing downstream quadrature rule;

• a mismatch between the volume direction used by preflight and the candidate deformation;

• a false affine-in-displacement capability claim;

• an operator evaluating mapped geometry outside its advertised rules; or

• a lower-level mapping status unrelated to determinant sign.

Such a result should be treated as a correctness bug in the preflight coverage or diagnostics, not merely as evidence to reduce the safety factor blindly.

Verification performed during this phase

Build discipline

Formatting was run before each code build. Debug targets for the main test executable and sandbox compiled successfully after the final geometry-rule deduplication and normalized-direction wiring. No physical, analytic, or finite-difference tolerance was loosened. A final whitespace/error check on the diff was clean.

Focused test evidence

The following results were recorded:

Scope	Result	What it establishes
Scope	Result	What it establishes
Two-rank safe-step estimator tests	5 cases, 32 assertions passed	First boundary, tangent root, positive floor, unconstrained step, and collective invalid-input behavior.
Standalone Newton/geometry focus	6 cases, 55 assertions passed	Geometry-related Newton and estimator paths at the then-current stage.
Final Physical Riesz geometry-cap integration test	1 case, 23 assertions passed, about 170 s	A normalized surface-only correction is denormalized, converted to volume displacement, capped before trial preparation, and accepted in one trial.
Generated rotation/angular momentum and warm-start/backtracking focus	2 cases, 231,772 assertions passed, about 417 s combined	Existing generated-rotation physics and nonlinear reuse behavior were not broken by the refactor.
Debug build of tests and sandbox	Passed	The current module graph and public sandbox compile.

The very large assertion count in the rotation/operator focus reflects element- and degree-of-freedom-level comparisons, not hundreds of thousands of independent scenarios.

What the tests cover structurally

The new and existing solver tests cover:

• generic linear-backend and preconditioner-contract detection;

• stable prepared ownership before evaluation;

• explicit/default seeds and rotation availability;

• read-only result views and opaque owning capture targets;

• failure-reason preservation;

• FGMRES restart configuration, initial guesses, and true-residual reports;

• context/view generation and solver-destruction lifetimes;

• paired observer hooks and line-search trial records;

• backtracking scale reuse and correction-vector warm starts;

• accepted final steps at the iteration limit;

• rollback after metric, observer, and preconditioner exceptions;

• recovery after non-finite or malformed backend corrections;

• end-to-end production context to MFEM FGMRES wiring; and

• the geometry cap at the normalized full-root boundary.

Why the full suite was not run

The full suite takes approximately two hours. The changes were concentrated in solver orchestration, mapped-rule exposure, and a new free geometry function, and the focused tests exercised those boundaries plus representative existing physics. Under the user’s instruction to use judgment, the full suite was not rerun solely for this handoff. It should be run before merging or checkpointing a broader follow-up that changes shared mapping, normalization, or operator semantics.

Unverified claims

The following are not established by the current tests:

• that the post-preflight real sandbox is faster;

• that every old gravity-force rejection is eliminated;

• that the physical nonrotating problem converges;

• that the residual of a known good solution decreases under refinement;

• that the current warm-start and 3 percent forcing policy are robust;

• that context setup or Jacobian actions scale acceptably;

• that a nonzero-angular-momentum sequence converges; or

• that capture, serialization, and nonlinear restart work.

Current state by subsystem

Subsystem	Status	Current meaning
Subsystem	Status	Current meaning
Variadic model/compiler	Implemented and inherited	Canonical inferred root remains the stable physics boundary.
Physical Riesz normalization	Implemented and exercised	Solver operates in normalized coordinates and denormalizes directions when physical geometry is required.
Problem object	Implemented	Owns model, discretization, compiled surface constraint, and prepared variadic operator.
User context	Implemented	Authoritative lifetime root; stable, non-copyable, non-movable; prepares all main numerical components.
Linear backend abstraction	Implemented	ADL preparation and contract-aware type checking permit future backends.
MFEM FGMRES	Implemented and focused-test verified	Restart length, warm start, true residual, counts, timings, and distributed validation exist.
Newton engine	MVP implemented	Normalized merit, backtracking, rollback, expected failures, and fixed-size diagnostics exist.
Observer API	Implemented and tested	Ergonomic two- or three-callback construction; synchronous on every rank.
Result views	Implemented and tested	Converged structure view and last-accepted checkpoint view are report-certified and generation checked.
Owning capture	Stub only	Signatures exist; calls throw.
Serialization	Stub only	Free overloads exist; calls throw.
Safe geometry estimator	Implemented in dirty tree	Exact at supplied points for affine-in-displacement maps; collective and diagnostic.
Newton geometry integration	Implemented in dirty tree and focused-test verified	Caps first trial before expensive candidate preparation.
Physical convergence	Not established	Recorded residual decrease is slow and steps collapse.
Performance quality	Not established	Context, accepted preparation, and Jacobian actions are expensive in Debug.
Continuation/restart	Future work	Architecture does not exclude it, but no persisted restart or continuation controller exists.

Known architectural limitations and issues

The context is intentionally more than a POD

The user requested a context or problem-like owner and described it as a POD. The implementation is a user-facing context but not a trivial aggregate. Its stable-address and one-active-solver invariants justify the stronger class. A future cleanup should preserve the semantic simplicity of “one user-owned lifetime root” even if internal storage is refactored.

Solver construction still moves configuration values

The undesired ownership transfer from user to solver was removed for the large problem state: the solver borrows the context. The small Newton configuration and observer are still stored by value in the solver, so move-only observers are supported. This is ordinary configuration ownership, not a transfer of the physical problem or solution.

Only one nonlinear method and one built-in linear backend exist

The interfaces are generic, but implementation breadth is deliberately narrow. Only Newton and MFEM FGMRES are built in. This matches the MVP decision. It should not be represented as demonstrated plug-and-play support for MINRES, arbitrary nonlinear methods, or basis recycling.

Metric extensibility is partial

The metric value is customizable at compile time. Acceptance, convergence threshold construction, and the concept of sufficient decrease remain hardwired. A future metric that needs a directional derivative, blockwise trust criterion, constraint filter, or multi-objective acceptance decision will need a slightly richer policy boundary.

View invalidation is strict

The next evaluation invalidates every earlier structure/checkpoint view from that context, even if an earlier report object remains alive. This is honest for mutable context-owned storage but may surprise users. Self-contained capture() is the intended escape hatch and is not implemented yet.

No complete persistence schema

A self-contained structure must decide how to store distributed fields, manifest metadata, mesh and refinement identity, finite-element order, units and normalization metadata, model specifications, generated scalar values, rotation, communicator-independent global numbering, and versioning. Saving only the raw root vector is insufficient for reliable offline investigation or restart.

Current result is not a scientific validation record

The report stores numerical solver totals, not virial error, mass error, angular-momentum error, surface residual, mesh quality distribution, or refinement comparison. Those diagnostics should be computed through separate physics observers or post-processing once the nonlinear path becomes reliable.

Distributed observer caveat

Observers run on every rank. They must not perform rank-divergent collectives or throw on one rank while peers enter unrelated collectives. The helper coordinates ordinary callback failure, but no generic library can repair arbitrary inconsistent user callback control flow.

Allocation behavior is not fully characterized

Most persistent solver workspaces are context-owned, but the geometry estimator currently creates local true-to-local vectors, element vectors, matrices, and mapping workspace objects on each call. This is acceptable under the clarified requirement but is a later optimization candidate if profiles show meaningful cost.

Known numerical and performance issues

Physical convergence is poor

The observed method is reducing the residual, but not at a useful rate. The step collapse reaches 6.1\times10^{-5} by iteration 5. Even if geometry-preflight cuts wall time, the same sequence would require far too many accepted steps unless direction quality or parameterization improves.

Warm-start and forcing interaction

The 3 percent linear tolerance permits a prior scaled correction to be accepted without new Krylov work. This may be a reasonable inexact Newton solve if the direction still produces adequate nonlinear progress, but the observed geometry-limited steps suggest otherwise. The code needs experiments that independently vary:

• cold versus warm initial correction;

• fixed tolerances such as 3\times10^{-2}, 10^{-2}, and tighter values;

• an adaptive Eisenstat-Walker-style forcing term;

• a minimum required improvement over the initial linear residual; and

• optional restart of the warm start after a severely damped nonlinear step.

Do not select one of these policies solely from intuition. Record Jacobian calls, true residual, accepted step, residual reduction, and total time.

Context construction costs 86.7 seconds

The top-level timer is too coarse to identify the source. Add subphase timings around seed projection, initial prepared-operator work, normalization construction, physical-preconditioner preparation, normalized wrapper construction, geometry-rule-plan construction, and FGMRES binding. Geometry-plan construction was added after the recorded number, so its incremental cost must also be measured.

Accepted trial preparation costs about 48 seconds

Removing rejected trials does not remove this floor. The prepared root currently rebuilds the mapped and physics-dependent caches needed by gravity, displacement forces, closure, hydrostatic balance, mass, and generated constraints. Profile each stage. Candidate failure ordering already prepares the displacement residual early after gravity so mapped-force failures occur before unrelated closure and hydrostatic rows, but successful candidates still pay for everything.

Jacobian application costs about 5.86 seconds

The first linear solve spent approximately 123 seconds in 21 Jacobian applications. Block-level timing is not yet printed. Likely candidates for measurement include repeated true/local transfers, repeated element geometry evaluation, multiple quadrature traversals over the same mapped data, temporary allocations, and Debug-only checks. Optimization should follow measurements and preserve analytic/finite-difference tolerances.

Debug timings are development evidence

The user requires Debug builds during development to avoid very long build times and preserve diagnostics. These timings are useful for identifying dominant algorithms. They are not production performance claims. A later controlled Release or profiling build may quantify absolute performance, but should not replace the Debug correctness workflow without agreement.

Geometry is only one globalization boundary

The new cap handles orientation preservation. If the next failure source becomes EOS admissibility or density/enthalpy domain constraints, an analogous cheap preflight may be useful. That should be driven by measured failure frequency and should reuse explicit free functions rather than embedding more exception-driven logic in the Newton loop.

Blockers and decision points

Immediate blockers

The following block a credible MVP result:

1. The new geometry preflight has not been benchmarked on the actual sandbox sequence.

2. The full residual has not been validated at a known nonrotating solution under refinement.

3. The zero-iteration warm-start behavior has not been compared against a cold start or adaptive forcing.

4. Context setup, successful candidate preparation, and Jacobian actions lack block-level profiles.

5. The nonlinear solve has not converged even the target nonrotating reference.

Persistence is not required to run the next experiments, but it blocks reliable long-run checkpointing and restart.

Choices that do not need to be reopened now

Unless new evidence contradicts them, retain:

• the user-owned context and borrowing solver;

• evaluation reports with structure/checkpoint views;

• free-function serialization;

• self-contained capture as the only owning capture mode;

• explicit Newton with backtracking as the reference nonlinear engine;

• MFEM FGMRES as the first built-in linear backend;

• compile-time metric selection;

• observer-based growing diagnostics; and

• expected values for ordinary candidate rejection.

Choices that should be evidence driven

The next work should gather data before selecting:

• fixed versus adaptive inexact-Newton forcing;

• when to retain or reset the correction warm start;

• whether the current surface parameterization is the primary cause of geometry limitation;

• which Jacobian block or mapped-data path to optimize first;

• whether a positive determinant conditioning floor is useful; and

• when the benefit of true Krylov recycling justifies backend complexity.

Recommended next work, in order

Step 1: run the post-preflight sandbox once

Build in Debug, run formatting before the build, and capture the full observer output. Record:

• context setup time;

• preflight sample count and time per iteration;

• requested maximum step, boundary, and safe step;

• limiting rank and element;

• first-trial disposition;

• any remaining invalid-mapping stage;

• accepted preparation time;

• FGMRES initial and true residual ratios;

• Jacobian and preconditioner calls; and

• residual reduction per accepted step.

Stop early if the first geometry-capped trial is still rejected by invalid mapped geometry. Diagnose coverage before spending hours on a full sequence.

Step 2: enrich rejection diagnostics

Expose enough of the fixed-size preparation rejection in AfterLineSearchTrial to distinguish:

• inverted determinant;

• singular/non-finite inverse mapping;

• compactification-specific invalidity;

• EOS-domain error code; and

• non-finite physics arithmetic.

Prefer enums and optional scalar data over allocated strings. Keep the existing human-readable source string for sandbox ergonomics if desired.

Step 3: validate the residual on the analytic nonrotating reference

For an n=1 Lane-Emden solution with consistent radius, mass, polytropic constant, and central density:

1. project the analytic state into each finite-element block;

2. prepare the full inferred root;

3. report normalized and physical norms by residual block;

4. verify fixed-mass and fixed-central-density rows separately;

5. compare the surface residual and minimum determinant;

6. repeat over mesh refinement and polynomial order; and

7. require convergence trends without loosening existing physical or analytic tolerances.

This is the direct answer to the unresolved question “is the residual small for a known good solution?” At present, the answer for the full root is not yet established.

Step 4: isolate warm-start and forcing behavior

Run a small controlled matrix with identical initial state and geometry:

Case	Linear initial guess	Forcing policy
Case	Linear initial guess	Forcing policy
A	zero each Newton step	current 0.03 fixed tolerance
B	scaled prior correction	current 0.03 fixed tolerance
C	zero each Newton step	tighter fixed tolerance
D	scaled prior correction	tighter fixed tolerance
E	scaled prior correction	adaptive forcing based on nonlinear progress

Compare total time to a fixed residual reduction, not time per iteration alone. A run with more Krylov work can be cheaper overall if it produces materially larger accepted Newton steps.

Step 5: add subphase profiling

Instrumentation should remain low-cost and observer-friendly. Measure:

• context construction phases;

• each prepared-operator stage for accepted and rejected candidates;

• each major Jacobian block action;

• true/local gather-scatter time;

• mapped-geometry evaluation time;

• EOS and quadrature kernel time; and

• allocation counts if a suitable low-overhead mechanism is available.

Avoid preprocessor-only instrumentation if a zero-cost empty observer or compile-time policy can provide the same release behavior.

Step 6: optimize the measured hot path

Likely interventions, to be chosen only after profiling, include:

• reusing per-element mapped geometry across operator blocks at the same prepared state;

• fusing quadrature traversals that consume identical rules and fields;

• retaining true-to-local representations through a complete Jacobian action;

• moving remaining per-application temporary storage into the context or prepared operator;

• avoiding duplicate determinant/inverse evaluations; and

• strengthening the physical structure preconditioner if fewer Jacobian applications outweigh more inverse work.

Every optimization must keep analytic and finite-difference tolerances unchanged.

Step 7: converge the nonrotating reference

Require:

• nonlinear convergence under the declared normalized norm;

• direct final residual agreement;

• blockwise residual control;

• mass and central-density constraint accuracy;

• orientation-preserving geometry with useful margin;

• reproducibility from the same seed; and

• a basic refinement trend.

Only then move to nonzero angular momentum.

Step 8: continuation in angular momentum

Use the converged nonrotating structure as the initial state for small fixed-angular-momentum increments. This requires either implementing capture/restart or adding a scoped context-to-context seed path that does not weaken ownership. Record angular velocity, moment of inertia, virial defect, mass and angular-momentum errors, surface shape, minimum determinant, nonlinear steps, and Krylov work.

Do not add a public solver.restart() API merely for this first experiment. The current architecture should first establish the checkpoint representation and continuation semantics.

Step 9: persistence

Design a versioned, communicator-independent schema, then implement:

• self-contained structure capture;

• self-contained checkpoint capture;

• free serialization overloads for owner and view;

• deserialization into an inspection value;

• compatibility checks against model/discretization metadata; and

• a controlled path to project or restore a checkpoint into a new context.

Serialization should not persist opaque MFEM pointers or assume identical MPI rank count.

File map for takeover

Committed solver architecture

File	Role
File	Role
libmeanfield/interface/solver/linear_backend.cppm	Generic backend concepts, linear controls/report, MFEM FGMRES configuration and prepared implementation.
libmeanfield/interface/solver/newton.cppm	Newton/backtracking options, metric policy, dispositions, observer events and helpers.
libmeanfield/interface/solver/stellar_context.cppm	Authoritative context ownership, construction overloads, prepared runtime, and explicit Newton loop.
libmeanfield/interface/solver/stellar_equilibrium_types.cppm	Public failure report and fixed-size evaluation diagnostics.
libmeanfield/interface/solver/stellar_structure.cppm	Report certification, views, opaque owning types, and persistence stubs.
libmeanfield/interface/solver/stellar_equilibrium.cppm	Solver module aggregation.
libmeanfield/interface/operators/stellar_equilibrium_problem.cppm	Owned discretized problem and solver-facing preparation/residual/Jacobian forwarding.
libmeanfield/interface/preconditioning/stellar_recipe.cppm	Default preconditioner recipe introduced for solver assembly.
sandbox.cpp	User API exercise and physical performance diagnostics.

Geometry-preflight working-tree files

File	Role
File	Role
libmeanfield/interface/deformation/safe_newton_step.cppm	Public rule, options, certificate, and free-function declaration.
libmeanfield/impl/deformation/safe_newton_step.cpp	Collective determinant-polynomial estimator and root search.
libmeanfield/interface/mapping/compactification/compactification.cppm	Affine-in-displacement capability contract.
libmeanfield/interface/mapping/compactification/kelvin.cppm	Kelvin capability declaration.
libmeanfield/interface/solver/stellar_context.cppm	Cached rule plan, normalized-to-physical direction conversion, cap invocation, and diagnostics.
libmeanfield/interface/solver/newton.cppm	Geometry settings and observer certificate fields.
libmeanfield/interface/solver/stellar_equilibrium_types.cppm	Evaluation-level geometry totals and last certificate.
libmeanfield/interface/mean_field.cppm	Public export of the estimator.
CMakeLists.txt	New source, interface, serial test, and MPI-test registration.

The prepared operator interfaces listed below expose their exact mapped quadrature rules:

• libmeanfield/interface/operators/prepared_hdiv_mass.cppm;

• libmeanfield/interface/operators/prepared_gravity_source.cppm;

• libmeanfield/interface/operators/prepared_barotropic_closure_operator.cppm;

• libmeanfield/interface/operators/prepared_hydrostatic_equilibrium_operator.cppm;

• libmeanfield/interface/operators/prepared_pressure_force.cppm;

• libmeanfield/interface/operators/prepared_gravity_displacement_force.cppm;

• libmeanfield/interface/operators/prepared_rotation_displacement_force.cppm;

• libmeanfield/interface/operators/prepared_mass_normalization.cppm; and

• libmeanfield/interface/operators/prepared_angular_momentum.cppm.

The barotropic, mass-normalization, and angular-momentum prepared runtimes now retain stable integration-rule pointers set by their implementations. These pointers refer to MFEM’s stable rule registry and are borrowed by the context sampling plan.

Direction forwarding files

The new physical volume-direction path touches:

• libmeanfield/interface/operators/prepared_stellar_equilibrium.cppm;

• libmeanfield/impl/operators/prepared_stellar_equilibrium.cpp;

• libmeanfield/interface/operators/prepared_variadic_stellar_equilibrium.cppm; and

• libmeanfield/interface/operators/stellar_equilibrium_problem.cppm.

Primary tests

File	Focus
File	Focus
tests/solver/stellar_equilibrium_architecture.cpp	Public compile-time and ownership contracts.
tests/solver/stellar_equilibrium_architecture_internal.cpp	Internal storage lifetime and generated-rotation report behavior.
tests/solver/stellar_equilibrium_runtime.cpp	FGMRES, Newton, observers, rollback, warm starts, geometry wiring, and production integration.
tests/deformation/safe_newton_step.cpp	Polynomial-boundary and collective estimator behavior in serial/MPI test targets.
tests/user-api/stellar_equilibrium.cpp	Simple, complete, intermediate, and advanced public API examples.

Working-tree preservation notes

1. Do not reset to HEAD; the geometry estimator and its tests are untracked.

2. Do not assume every modified line belongs to one agent. User-authored domain-mapper comments and sandbox changes coexist with the geometry work.

3. Do not delete untracked build trees, CSV experiments, notes, archives, or scripts without explicit scope confirmation.

4. Before a source build, run ./format as requested by the user.

5. Prefer Debug builds for development.

6. Do not run the approximately two-hour full suite unless the regression risk justifies it; report focused coverage honestly.

7. Never loosen physical, analytic, or finite-difference tolerances to make solver changes pass.

8. The user designated 76818f2f828ca12606a495a0b84f5d6ce4b78dc6 as the conversation checkpoint even though current HEAD is later.

Invariants future work should preserve

1. One authoritative owner. The context owns the physical problem and numerical lifetime; solvers borrow it.

2. Report-certified access. A structure view exists only after convergence; a failed solve exposes only a checkpoint view.

3. No pack-specific solver. The solver consumes the inferred normalized problem and preconditioner contracts, not named constraint combinations.

4. Linear/nonlinear separation. Linear backend choice and nonlinear parameters remain independently selectable.

5. Expected failures are data. Ordinary line-search rejection does not use exceptions.

6. Exceptions restore accepted state. Infrastructure and callback failures cannot leave a candidate silently installed.

7. Direct linear residual. Backend recurrence convergence is never the only reported correctness criterion.

8. Stable prepared addresses. Backend and preconditioner borrows cannot be invalidated by object moves.

9. Observer cost is optional. Growing history and detailed diagnostics remain outside the default fixed-size report.

10. Physical/normalized distinction. Geometry is evaluated from a correctly denormalized physical direction.

11. Capability honesty. Exact polynomial preflight applies only to maps and deformations that advertise the required affine/linear contracts.

12. Quadrature coverage. Every operator that can reject mapped geometry must advertise the points it evaluates or share a stronger common geometry plan.

13. No tolerance erosion. Optimization and globalization work cannot weaken existing physics tests.

14. Persistence is self-contained. A captured structure or checkpoint must not depend on a live solver, context, communicator, or MFEM pointer.

15. Scientific claims require independent evidence. Solver convergence, block residuals, invariants, virial error, and refinement are separate checks.

Takeover checklist

A developer resuming this work should:

1. confirm repository root, checkpoint, current HEAD, and dirty-tree status;

2. read the previous variadic handoff for the inherited model/compiler architecture;

3. read linear_backend.cppm, newton.cppm, stellar_structure.cppm, and stellar_context.cppm in that order;

4. read prepared_stellar_equilibrium.cppm/.cpp and the displacement-force rejection path;

5. read the complete safe-step interface, implementation, and tests;

6. inspect the current sandbox configuration and preserve user edits;

7. run formatting before any Debug build;

8. run the focused geometry and Newton tests after relevant changes;

9. benchmark only enough of the real sandbox to verify first-trial geometry behavior;

10. do not run a long sequence if the capped first trial still fails for invalid mapping;

11. perform the known-solution block-residual study before claiming solver error is the only problem;

12. perform the warm/cold and forcing-policy comparison before removing or expanding warm-start behavior;

13. profile by block before optimizing Jacobian or accepted candidate preparation;

14. preserve fixed-size failure records and observer-based histories; and

15. commit the geometry work separately from unrelated user notes and artifacts when the user requests a checkpoint.

Glossary

Term	Meaning in the current solver
Term	Meaning in the current solver
Accepted state	The committed normalized/physical root whose prepared operator, residual, and preconditioner are mutually current.
Candidate or trial state	A temporary line-search state that may be rejected without changing the accepted state.
Checkpoint view	Read-only, generation-checked access to the last accepted state, available after success or failure.
Context	The user-owned, stable lifetime root for the problem, state, prepared operator, preconditioner, backend, and workspaces.
Correction warm start	The previous accepted Newton correction scaled by 1-\alpha and supplied as the next FGMRES initial guess.
FGMRES restart	Starting a new Krylov cycle after the configured basis length; not nonlinear checkpoint restart.
Geometry preflight	Mapping-only calculation of the first sampled determinant boundary along a volume-displacement direction.
Inexact Newton forcing	The permitted linear residual relative to the Newton right-hand side. The sandbox currently uses 0.03.
Metric	Compile-time policy converting a normalized residual into residual norm and scalar merit.
Prepared backend	An owning linear runtime bound to stable operator, preconditioner, communicator, and internal workspaces.
Structure view	Read-only, generation-checked access certified only by a converged report.
True residual	Directly recomputed b-Ax, independent of the Krylov recurrence estimate.

Recorded raw sandbox sequence

The essential supplied output is reproduced here so later work can compare like for like. Whitespace is normalized, but numerical values are unchanged.

Finite-element setup: 0.156279 s
Solver context setup: 86.6662 s
Newton 0: |F| = 0.0214928
  trial 1: step = 1, outcome = inadmissible state,
           source = generated geometry, min(det J_map) = -10.0826,
           prepare = 1.48201 s
  trial 2: step = 0.0625, outcome = accepted, |F| = 0.0201477,
           min(det J_map) = 0.309289, prepare = 48.3029 s
  step = 0.0625, trials = 2, iteration = 194.789 s,
  line search = 53.6456 s, trial preparation = 49.7849 s,
  accepted refresh = 3.85453 s
  FGMRES: converged, iterations = 19, restarts = 0,
  initial/|b| = 1, true/|b| = 0.0267531, J calls = 21,
  P^-1 calls = 19, solve = 135.281 s, J time = 123.003 s,
  P^-1 time = 18.1196 s

Newton 1: |F| = 0.0201477
  trial 1: step = 0.125, outcome = inadmissible state,
           source = generated geometry, min(det J_map) = -0.986494,
           prepare = 1.49634 s
  trial 2: step = 0.015625, outcome = inadmissible state,
           source = gravity displacement force, prepare = 32.8684 s
  trial 3: step = 0.0078125, outcome = accepted, |F| = 0.0199897,
           min(det J_map) = 0.22833, prepare = 48.2396 s
  step = 0.0078125, trials = 3, iteration = 98.1586 s,
  line search = 86.4674 s, trial preparation = 82.6044 s,
  accepted refresh = 3.8553 s
  FGMRES: converged, iterations = 0, restarts = 0,
  initial/|b| = 0.0296713, true/|b| = 0.0296713,
  J calls = 2, P^-1 calls = 0, solve = 5.85094 s,
  J time = 11.6888 s

Newton 2: |F| = 0.0199897
  trials: 0.015625 rejected by gravity displacement force;
          0.0078125 rejected by gravity displacement force;
          0.00390625 accepted with |F| = 0.0199115
  iteration = 136.953 s, line search = 118.448 s,
  trial preparation = 114.522 s, accepted refresh = 3.91821 s
  FGMRES: converged, iterations = 1, true/|b| = 0.0266376,
  J calls = 3, P^-1 calls = 1, solve = 12.6573 s,
  J time = 17.5629 s, P^-1 time = 0.939231 s

Newton 3: |F| = 0.0199115
  trials: 0.0078125, 0.00390625, 0.00195312, and 0.000976562
          rejected by gravity displacement force;
          0.000488281 accepted with |F| = 0.0199017
  iteration = 197.776 s, line search = 185.926 s,
  trial preparation = 182.025 s, accepted refresh = 3.88939 s
  FGMRES: converged, iterations = 0, true/|b| = 0.0267994,
  J calls = 2, solve = 5.92722 s, J time = 11.8481 s

Newton 4: |F| = 0.0199017
  trials: 0.000976562 and 0.000488281 rejected by
          gravity displacement force;
          0.000244141 accepted with |F| = 0.0198968
  iteration = 130.261 s, line search = 118.532 s,
  trial preparation = 114.693 s, accepted refresh = 3.83085 s
  FGMRES: converged, iterations = 0, true/|b| = 0.0268207,
  J calls = 2, solve = 5.85722 s, J time = 11.7267 s

Newton 5: |F| = 0.0198968
  trials: 0.000488281, 0.000244141, and 0.00012207 rejected by
          gravity displacement force;
          0.0000610352 accepted with |F| = 0.0198956
  iteration = 162.629 s, line search = 150.949 s,
  trial preparation = 147.105 s, accepted refresh = 3.83446 s
  FGMRES: converged, iterations = 0, true/|b| = 0.0268315,
  J calls = 2, solve = 5.84238 s, J time = 11.6781 s

Newton 6: |F| = 0.0198956

Suggested benchmark record

For each post-preflight iteration, save one row with:

iteration, residual, linear_initial_relative, linear_true_relative,
linear_iterations, J_calls, J_seconds, P_calls, P_seconds,
preflight_seconds, preflight_samples, boundary_step, safe_step,
limiting_rank, limiting_element, trials, rejection_stages,
accepted_step, accepted_prepare_seconds, residual_after, iteration_seconds

This is sufficient to test the predicted speedup and the warm-start hypothesis without storing full vectors.

Final handoff statement

MeanField has advanced from a solver-ready inferred operator stack to an actual reference Newton-Krylov orchestration with explicit ownership, typed configuration, robust result semantics, and unusually detailed low-cost diagnostics. The user-facing design now matches the central intent: the context is the authoritative owner; the solver borrows it; evaluation produces a report; and the report grants an honest borrowed structure or checkpoint view.

The first physical run exposed two separate problems rather than one. First, the line search wasted most of its time discovering mapped-geometry limits at expensive downstream quadrature points. The new exact sampled determinant preflight is a principled response and is wired into Newton before candidate preparation. Second, the accepted Newton directions and broader operator stack remain too expensive and too weakly convergent. In particular, repeated zero-iteration FGMRES acceptance of a carried direction under the 3 percent forcing tolerance is a concrete mechanism that may help explain step collapse.

The next agent should resist both extremes: do not mistake the geometry preflight for a convergence fix, and do not discard it merely because the underlying solve is still slow. Benchmark it, verify its coverage, then use the time and residual records to separate geometry globalization, linear forcing, preconditioner quality, Jacobian cost, and discretization correctness. Establish the full residual at a known nonrotating solution before tuning toward rapid rotation. Preserve the dirty tree, the checkpoint semantics, and the user’s no-tolerance-loosening rule throughout.