diff --git a/experiments/NUMERICAL_STRATEGIES_2026-09-17.md b/experiments/NUMERICAL_STRATEGIES_2026-09-17.md new file mode 100644 index 0000000..a2d1e29 --- /dev/null +++ b/experiments/NUMERICAL_STRATEGIES_2026-09-17.md @@ -0,0 +1,85 @@ +# Optional Newton strategies: implementation and validation + +The experiments start from `7a07cae`. Each hypothesis has an isolated, build-configurable branch; the integration branch combines the implementations behind explicit typed policies. The default physical model, normalization, preconditioner prescription, and ordinary Newton path remain available. Central density is still an optional model specification. + +## Selection and ownership + +`Newton(options)` retains its existing meaning. An expert numerical configuration can supply a third constructor argument: + +```cpp +namespace nl = mean_field::solver::nonlinear; +auto numerics = nl::NewtonNumerics{ + nl::PolynomialEnthalpy{}, + nl::GravityCurvatureCorrection{} +}; +auto configuration = nl::Newton(options, nl::NormalizedResidualMetric{}, numerics); +auto solve = mean_field::solver::make(context, configuration); +auto report = solve.evaluate(); +``` + +The enthalpy choices are `SampledEnthalpy` and `PolynomialEnthalpy`. Gravity choices are `MonolithicGravity`, `GravityCurvatureCorrection`, and `EliminatedGravity`. Their combinations are concrete C++ types. Unregistered policy types are rejected by concepts. Block identities come from the compiled equation form; mesh sizes, FE eligibility, arithmetic values, and distributed preparation generations are checked at runtime. + +The context owns each prepared workspace. Numerical calls are statically dispatched; only destruction of the selected context-owned storage is erased. Workspaces are constructed when a solve selects them and are destroyed before their borrowed operator/backend dependencies. Reduced gravity recovers the outer backend prescription through `GetConfiguration()`; custom backends need that capability only when selecting elimination. Polynomial direction repair currently requires the normalized residual metric, because its derivative supplies the line-search slope. + +## 1. Polynomial enthalpy feasibility + +The implementation transforms the actual nodal polynomial to a tensor Bernstein representation using outward-rounded interval arithmetic. Adaptive subdivision either certifies nonnegativity, proves negativity through a negative polynomial value, or reports unresolved bounds. A negative Bernstein coefficient alone is never classified as a negative polynomial. Exact zero traces are retained algebraically. + +The production adapter works in physical enthalpy coordinates; a general diagonal normalization need not preserve the polynomial encoded by normalized coefficients. It supports conforming scalar H1 hexahedral elements with the supported closed nodal tensor bases and runtime polynomial order. Unsupported spaces fail explicitly. Unsafe fast-math compilation is rejected. Gradual underflow and immutable FE spaces remain prerequisites. + +An invalid or uncertifiable initial polynomial is rejected before a residual-based convergence declaration. For a valid base, a certified initial enthalpy field provides an immutable anchor. The repair family is + +\[ +p_h(\beta)=(1-\beta)p_h+\beta(h_{\rm anchor}-h),\qquad 0\le\beta\le1. +\] + +Shared true DOFs and prescribed zero traces are preserved. Bernstein inequalities provide a sufficient lower bound on the repair fraction, and the full linear residual model selects a fraction within that feasible interval. The actual represented candidate is certified again after normalization and state addition. A non-descent repair fails explicitly. Backtracking uses the repaired direction's merit derivative and still evaluates the complete nonlinear residual. + +This is a conservative feasible-direction method, not the solution of a constrained Newton optimization problem. It may be too restrictive or have no descent direction. It does not prove positivity of the independently represented density field. Unchanged feasible directions require one extra Jacobian action; a model-selected enthalpy repair requires three. Combining it with elimination also relifts gravity and recomputes the actual merit derivative. + +## 2. Gravity curvature correction + +Write the normalized state as retained variables `u` and gravity variables `z`, with gravity equations `G(u,z)`. The gravity partition includes both gradient and potential equations and retains their boundary rows. Given a base `x`, Newton direction `p`, and trial `y=x+alpha*p`, solve + +\[ +G_z(y)\,\delta z=-\{G(y)-G(x)-\alpha G'(x)p\}. +\] + +The right side removes only the nonlinear remainder. It vanishes on an affine gravity path even when the base gravity residual is nonzero; it is second order in the trial length when the operator is smooth. The base residual and Jacobian action are captured before trial preparation. The candidate inverse is refreshed before solving. Only the candidate's gravity coordinates change. + +A fresh complete preparation verifies the intended gravity target and evaluates the full merit. A failed inner solve, failed preparation, or unacceptable corrected merit does not commit a state. Trial observers must publish the corrected physical and normalized vectors together. The tangent is distinct from the secant of this curved trial path. + +## 3. Gravity elimination + +The nonlinear trial path completes gravity at fixed retained variables. The initial state is also explicitly completed before defining the nonlinear convergence reference. Retained coordinates include density, enthalpy, shape, and all model specification borders; optional central-density borders remain optional. + +For the partitioned Jacobian, the reduced action is + +\[ +S=F_u-F_zG_z^{-1}G_u. +\] + +Right-hand-side condensation includes any remaining gravity defect, and lifting reconstructs the corresponding full correction. The retained projection of the supplied full inverse preconditions the reduced system. Every lifted outer correction must satisfy the original full normalized Jacobian residual tolerance. Candidate gravity completion is separately checked after fresh nonlinear preparation. + +Tolerance-stopped inner Krylov solves approximate a linear Schur action. Flexible outer preconditioning does not make this action exact. Tight inner tolerances, directional derivative checks, and the full lifted residual gate control this limitation. Inner failure propagates explicitly. Operator actions never prepare or mutate the physical state. + +The generic reference implementation restricts the full operator and supplied inverse. That preserves normalization, boundary conditions, model borders, and custom preconditioner recipes, but can do unnecessary work in non-gravity blocks. A specialized restricted inverse is a possible later performance improvement once equivalence is tested. + +## Interpreting reports + +The complete nonlinear residual remains the merit and convergence authority. Gravity residual reduction alone does not establish a stellar equilibrium. The existing physical checker separately requires solver convergence and all eight energy/force, enthalpy/density, and quadrature-order virial measurements to pass `1e-5`; it also reports local EOS, force, sign, and conservation diagnostics. + +For elimination, a successful linear backend report describes the reduced system; the lifted full-system true residual is a separate diagnostic. A propagated inner failure describes the gravity subsystem. Nested gravity solve times overlap the outer reduced solve time, and backend solve time excludes some condensation/lifting work. Use complete iteration and campaign wall times for performance comparisons, not a sum of these overlapping measurements. + +## Reproducibility and evidence + +Branches: + +- `research/polynomial-positive-newton-20260917` +- `research/gravity-curvature-correction-20260917` +- `research/gravity-reduced-newton-20260917` +- `research/numerical-strategies-20260917` (combined production integration) + +The campaign runner builds a committed clean source revision, copies an immutable executable, and records source, mesh, executable, input-state, and output hashes. Serial execution keeps timing comparisons interpretable. Each new run uses a fresh output directory. Frozen probes preserve accepted buffers and verify restoration; they never certify a final physical model. + +Validation results and measured limits are recorded below after execution. diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 66b40f6..23df102 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -374,7 +374,7 @@ void FrozenProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vec template int RunContext(Context &context, const Options &o, const Scales &scales) { std::vector rows, values; - mfem::Vector residualFactors; + mfem::Vector residualFactors, stateFactors; std::vector centralRows; int centralOffset = -1; auto metadata = File(o.output / "metadata.txt"); @@ -441,6 +441,7 @@ template int RunContext(Context &context, const Options &o, cons rows.assign(r.begin(), r.end()); values.assign(v.begin(), v.end()); residualFactors = s.normalizedOperator->GetNormalization().ResidualFactors(); + stateFactors = s.normalizedOperator->GetNormalization().StateFactors(); for (const auto &b : values) if (b.stableId == "fixed_central_density.border") centralOffset = b.offset; @@ -462,6 +463,8 @@ template int RunContext(Context &context, const Options &o, cons blocklog << "iteration,stage,block,normalized_norm,physical_norm,max_abs_normalized\n"; auto triallog = File(o.output / "trials.csv"); triallog << "iteration,trial,alpha,disposition,source,F,min_jacobian\n"; + auto stateConsistency = File(o.output / "trial_state_consistency.csv"); + stateConsistency << "iteration,trial,normalized_physical_roundtrip_error,relative_error\n"; mfem::Vector captureBase; auto captureDirectory = [&](int iteration) { return o.output / ("capture_" + std::to_string(iteration)); @@ -503,6 +506,20 @@ template int RunContext(Context &context, const Options &o, cons std::cout << "Newton " << e.iteration << " F=" << e.residualNorm << std::endl; }, [&](const solver::nonlinear::AfterLineSearchTrial &e) { + if (e.candidatePhysicalState.size() != e.candidateNormalizedState.size() || + e.candidateNormalizedState.size() != static_cast(stateFactors.Size())) + throw std::runtime_error("Trial observer exposed incompatible state representations"); + long double stateError = 0.0L, stateNorm = 0.0L; + for (int i = 0; i < stateFactors.Size(); ++i) { + const long double difference = e.candidatePhysicalState[i] * stateFactors(i) - e.candidateNormalizedState[i]; + stateError += difference*difference; + stateNorm += static_cast(e.candidateNormalizedState[i])*e.candidateNormalizedState[i]; + } + const double relativeStateError = std::sqrt(stateError)/std::max(1.0L,std::sqrt(stateNorm)); + stateConsistency << e.iteration << ',' << e.trial << ',' << std::sqrt(stateError) << ',' << relativeStateError << '\n'; + stateConsistency.flush(); + if (std::isfinite(relativeStateError) && relativeStateError > 64.0*std::numeric_limits::epsilon()) + throw std::runtime_error("Trial observer physical state disagrees with its normalized state"); triallog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << int(e.disposition) << ',' << e.rejectionSource << ',' << (e.metric ? e.metric->residualNorm : NAN) << ',' << e.minimumJacobianDeterminant.value_or(NAN) << '\n'; @@ -520,7 +537,8 @@ template int RunContext(Context &context, const Options &o, cons direction /= e.stepLength; reconstructed.Add(e.stepLength, direction); reconstructed -= candidate; - SaveVector(directory / "normalized_direction.txt", direction); + SaveVector(directory / "normalized_effective_trial_direction.txt", direction); + SaveVector(directory / "normalized_direction.txt", copySpan(e.normalizedDirection)); auto info = File(directory / "capture.txt"); info << "iteration=" << e.iteration << "\nfirst_trial_alpha=" << e.stepLength << "\ndisposition=" << int(e.disposition) << "\nsource=" << e.rejectionSource @@ -529,7 +547,9 @@ template int RunContext(Context &context, const Options &o, cons << "\ndirection_roundoff_scale_estimate=" << 8. * std::numeric_limits::epsilon() * (captureBase.Norml2() + candidate.Norml2()) / e.stepLength - << "\ndirection_source=first_trial_difference_divided_by_alpha\n" + << "\ndirection_source=exact_observer_tangent\n" + << "effective_trial_direction_source=candidate_difference_divided_by_alpha\n" + << "roundoff_scale_applies_to_effective_trial_direction_only=1\n" << "roundoff_scale_is_estimate_not_certificate=1\n"; } }, @@ -573,7 +593,7 @@ template int RunContext(Context &context, const Options &o, cons const nl::GravitySolveOptions gravity{ .control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14, .maximumIterations = o.gravityIterations}, - .backend = {.restartLength = o.gravityRestart}}; + .backend = {.restartLength = o.gravityRestart, .printLevel = -1}}; auto run = [&](auto numerics) { auto solve = solver::make(context, nl::Newton(options, nl::NormalizedResidualMetric{}, numerics), observer); return solve.evaluate(); @@ -597,6 +617,8 @@ template int RunContext(Context &context, const Options &o, cons metric("gravity_solves", strategyDiagnostics.gravitySolves); metric("gravity_iterations", strategyDiagnostics.gravityIterations); metric("gravity_failures", strategyDiagnostics.gravityFailures); + metric("gravity_prepared_failures", strategyDiagnostics.gravityPreparedFailures); + metric("linear_report_space", static_cast(strategyDiagnostics.lastLinearReportSpace)); metric("gravity_solve_seconds", strategyDiagnostics.gravitySolveSeconds); metric("last_gravity_true_residual", strategyDiagnostics.lastGravityTrueResidual); metric("last_gravity_correction_norm", strategyDiagnostics.lastGravityCorrectionNorm); diff --git a/libmeanfield/interface/solver/newton.cppm b/libmeanfield/interface/solver/newton.cppm index 97eaf7a..2a13485 100644 --- a/libmeanfield/interface/solver/newton.cppm +++ b/libmeanfield/interface/solver/newton.cppm @@ -265,6 +265,9 @@ export namespace mean_field::solver::nonlinear { std::span candidatePhysicalState{}; std::span candidateNormalizedState{}; std::span candidateNormalizedResidual{}; + // A curved trial path cannot recover its Newton tangent by dividing + // the candidate displacement by alpha. Observers borrow it directly. + std::span normalizedDirection{}; }; struct AfterIteration final { diff --git a/libmeanfield/interface/solver/newton_numerics.cppm b/libmeanfield/interface/solver/newton_numerics.cppm index fd724fd..7b3d59f 100644 --- a/libmeanfield/interface/solver/newton_numerics.cppm +++ b/libmeanfield/interface/solver/newton_numerics.cppm @@ -1,6 +1,7 @@ module; #include +#include #include #include @@ -22,7 +23,7 @@ export namespace mean_field::solver::nonlinear { struct GravitySolveOptions final { LinearSolveControl control{.relativeTolerance = 1.0e-10, .absoluteTolerance = 1.0e-14, .maximumIterations = 300}; - linear::FGMRESOptions backend{.restartLength = 80}; + linear::FGMRESOptions backend{.restartLength = 80, .printLevel = -1}; void Validate() const { control.Validate(); backend.Validate(); } }; @@ -64,14 +65,18 @@ export namespace mean_field::solver::nonlinear { template concept NewtonNumericalConfiguration = IsNewtonNumerics>::value; + enum class LinearReportSpace : std::uint8_t { full_system, reduced_system, gravity_subsystem }; + struct NumericalStrategyDiagnostics final { int polynomialChecks{0}; int directionRepairs{0}; double lastRepairFraction{0.0}; double lastMeritDirectionalDerivative{0.0}; - int gravitySolves{0}; - int gravityIterations{0}; - int gravityFailures{0}; + std::uint64_t gravitySolves{0}; + std::uint64_t gravityIterations{0}; + std::uint64_t gravityFailures{0}; + std::uint64_t gravityPreparedFailures{0}; + LinearReportSpace lastLinearReportSpace{LinearReportSpace::full_system}; double gravitySolveSeconds{0.0}; double lastGravityTrueResidual{0.0}; double lastGravityCorrectionNorm{0.0}; diff --git a/libmeanfield/interface/solver/stellar_context.cppm b/libmeanfield/interface/solver/stellar_context.cppm index 9947b8b..92b3689 100644 --- a/libmeanfield/interface/solver/stellar_context.cppm +++ b/libmeanfield/interface/solver/stellar_context.cppm @@ -1307,9 +1307,14 @@ export namespace mean_field::solver { } if (!strategy.CaptureDirection(state, diagnostics.numericalStrategies)) { - NotifyAfter(iteration, nonlinear::IterationDisposition::globalization_failure, false, 0.0, 0, + const bool linearFailure = strategy.DirectionSolveFailed(); + NotifyAfter(iteration, linearFailure ? nonlinear::IterationDisposition::linear_solve_failure : + nonlinear::IterationDisposition::globalization_failure, false, 0.0, 0, diagnostics.initialResidualNorm, acceptedMetric, acceptedMetric, linearReport, timings, state); - return Failure(state, std::move(diagnostics), StellarEquilibriumFailureReason::globalization_failure, + return Failure(state, std::move(diagnostics), linearFailure ? + StellarEquilibriumFailureReason::linear_solve_failure : + StellarEquilibriumFailureReason::globalization_failure, + linearFailure ? "Lifting the repaired direction failed its checked gravity solve." : "No certified enthalpy repair with a decreasing linear merit was available."); } @@ -1780,7 +1785,8 @@ export namespace mean_field::solver { .candidateNormalizedState = detail::ReadOnlySpan(state.trialNormalizedState), .candidateNormalizedResidual = residualAvailable ? detail::ReadOnlySpan(state.trialNormalizedResidual) - : std::span{} + : std::span{}, + .normalizedDirection = detail::ReadOnlySpan(state.normalizedCorrection) }; const Clock::time_point observerStart = Clock::now(); nonlinear::detail::InvokeAfterLineSearchTrial(m_observer, event); diff --git a/libmeanfield/interface/solver/stellar_newton_strategy.cppm b/libmeanfield/interface/solver/stellar_newton_strategy.cppm index 448527b..38f0ea1 100644 --- a/libmeanfield/interface/solver/stellar_newton_strategy.cppm +++ b/libmeanfield/interface/solver/stellar_newton_strategy.cppm @@ -118,6 +118,8 @@ namespace mean_field::solver::detail { report.gravitySolves = subsystem.GetSolveCount(); report.gravityIterations = subsystem.GetIterationCount(); report.gravityFailures = subsystem.GetFailureCount(); + report.gravityPreparedFailures = preparedFailures; + report.lastLinearReportSpace = lastLinearReportSpace; report.gravitySolveSeconds = subsystem.GetSolveSeconds(); if (const auto &last = subsystem.GetLastSolveReport(); last.has_value()) { report.lastGravityTrueResidual = last->trueResidualNorm; @@ -174,6 +176,7 @@ namespace mean_field::solver::detail { std::isfinite(targetDefect) && targetDefect <= tolerance; lastPreparedDefect = targetDefect; } + if (!preparation || !trialSolveSucceeded) { ++preparedFailures; } return preparation; } @@ -186,6 +189,8 @@ namespace mean_field::solver::detail { double lastCorrectionNorm{0.0}; double lastPreparedDefect{0.0}; double liftedRelativeDefect{0.0}; + std::uint64_t preparedFailures{0}; + nonlinear::LinearReportSpace lastLinearReportSpace{nonlinear::LinearReportSpace::full_system}; }; template @@ -259,13 +264,18 @@ namespace mean_field::solver::detail { LinearSolveReport Solve(State &state, const LinearSolveControl &control) { this->subsystem.Refresh(); try { + this->lastLinearReportSpace = nonlinear::LinearReportSpace::gravity_subsystem; auto inner = m_operator.TryCondenseRightHandSide(state.linearRightHandSide, m_rhs); if (!inner.Converged()) { return inner; } this->subsystem.GetPartition().GatherRetainedValues(state.normalizedCorrection, m_correction); + this->lastLinearReportSpace = nonlinear::LinearReportSpace::reduced_system; auto report = m_backend.Solve(m_rhs, m_correction, control); if (!report.Converged()) { return report; } inner = m_operator.TryLiftCorrection(m_correction, state.linearRightHandSide, state.normalizedCorrection); - if (!inner.Converged()) { return inner; } + if (!inner.Converged()) { + this->lastLinearReportSpace = nonlinear::LinearReportSpace::gravity_subsystem; + return inner; + } state.normalizedOperator->Mult(state.normalizedCorrection, this->fullAction); this->fullAction -= state.linearRightHandSide; const auto communicator = this->subsystem.GetCommunicator(); @@ -283,6 +293,7 @@ namespace mean_field::solver::detail { } return report; } catch (const GravityReductionFailure &failure) { + this->lastLinearReportSpace = nonlinear::LinearReportSpace::gravity_subsystem; return failure.GetReport(); } } @@ -368,9 +379,14 @@ namespace mean_field::solver::detail { return report; } bool CaptureDirection(State &state, nonlinear::NumericalStrategyDiagnostics &diagnostics) { + m_directionSolveFailed = false; if (!m_enthalpy.Repair(state, diagnostics)) { return false; } if constexpr (!std::same_as || altersDirection) { - if (!m_gravity.CaptureDirection(state)) { m_gravity.CopyDiagnostics(diagnostics); return false; } + if (!m_gravity.CaptureDirection(state)) { + m_directionSolveFailed = true; + m_gravity.CopyDiagnostics(diagnostics); + return false; + } } if constexpr (altersDirection && std::same_as) { mfem::Vector action(state.Problem().EquationSize()); @@ -389,7 +405,9 @@ namespace mean_field::solver::detail { return preparation; } bool TrialSolveSucceeded() const noexcept { return m_gravity.TrialSolveSucceeded(); } + bool DirectionSolveFailed() const noexcept { return m_directionSolveFailed; } private: + bool m_directionSolveFailed{false}; PreparedEnthalpyStrategy m_enthalpy; PreparedGravityStrategy m_gravity; };