diff --git a/experiments/GRAVITY_PREPARED_ROUNDOFF_2026-09-17.md b/experiments/GRAVITY_PREPARED_ROUNDOFF_2026-09-17.md new file mode 100644 index 0000000..ed615e0 --- /dev/null +++ b/experiments/GRAVITY_PREPARED_ROUNDOFF_2026-09-17.md @@ -0,0 +1,45 @@ +# Gravity prepared-residual arithmetic experiment + +This hypothesis branch changes only the arithmetic assessment used after a +gravity correction has passed its inner true-residual gate. It does not change +that inner tolerance, the full-system lifted Newton gate, nonlinear merit +acceptance, physical equations, or virial validation. + +For fixed retained state, write the normalized gravity residual as +`G = A_q q + A_phi phi + b`. Its norm can be tiny while the separate terms are +large. Estimating arithmetic from `||G||` and the combined correction action +therefore misses cancellation and the floating-point addition of corrected +gravity coefficients. + +Before changing the candidate, the workspace now measures the normalized +gradient and potential actions separately, for both the old state and its +correction. It infers the fixed source as `b = G_old - A_q q - A_phi phi`. +Selection uses the compiled field block types, and applications reuse the +prepared gravity operator. Scratch vectors are allocated with the workspace. + +The assessment uses `256 epsilon` times the sum of those five contribution +norms and the norms of the old, desired, and freshly prepared gravity residuals. +This is an explicit arithmetic estimate, **not a rigorous forward-error bound**: +it does not bound cancellation inside each individual field action. The factor +does not create an absolute residual floor, so a common scaling of the gravity +equations scales the assessment equally. + +Two separate requirements remain: + +1. The discrepancy between the affine prediction and freshly prepared gravity + residual must fit the arithmetic estimate. +2. The prepared residual's distance from its target must not exceed the + **measured** inner true residual plus that estimate. + +Diagnostics retain the raw affine discrepancy, raw target defect, arithmetic +estimate, and resulting target threshold separately. Nonfinite or negative +norm inputs fail the assessment. A failed preparation leaves these measurements +unavailable rather than reusing the preceding trial's values. Diagnostics also +report whether compact gravity operator and inverse applications were selected. + +Three focused tests, appended to `tests/solver/gravity_reduction.cpp`, cover +large cancelling contributions, rejection of defects beyond either gate, +nonfinite measurements, and scaling over twenty-four orders of magnitude. +The real stellar rejection must still be reproduced and reassessed; passing a +synthetic arithmetic test alone establishes neither the hypothesis nor stellar +convergence. diff --git a/libmeanfield/interface/solver/newton_numerics.cppm b/libmeanfield/interface/solver/newton_numerics.cppm index 7b3d59f..42e37fc 100644 --- a/libmeanfield/interface/solver/newton_numerics.cppm +++ b/libmeanfield/interface/solver/newton_numerics.cppm @@ -1,7 +1,10 @@ module; #include +#include #include +#include +#include #include #include @@ -76,13 +79,53 @@ export namespace mean_field::solver::nonlinear { std::uint64_t gravityIterations{0}; std::uint64_t gravityFailures{0}; std::uint64_t gravityPreparedFailures{0}; + bool compactGravityOperator{false}; + bool compactGravityInverse{false}; LinearReportSpace lastLinearReportSpace{LinearReportSpace::full_system}; double gravitySolveSeconds{0.0}; double lastGravityTrueResidual{0.0}; double lastGravityCorrectionNorm{0.0}; double lastGravityPreparedDefect{0.0}; + double lastGravityAffineDisagreement{0.0}; + double lastGravityArithmeticAllowance{0.0}; + double lastGravityPreparedTargetThreshold{0.0}; double initialGravityResidual{0.0}; double completedInitialGravityResidual{0.0}; double lastLiftedRelativeTrueResidual{0.0}; }; } + +export namespace mean_field::solver::detail { + struct GravityPreparedResidualAssessment final { + bool admissible{false}; + double arithmeticAllowance{std::numeric_limits::quiet_NaN()}; + double targetThreshold{std::numeric_limits::quiet_NaN()}; + }; + + [[nodiscard]] inline GravityPreparedResidualAssessment AssessGravityPreparedResidual( + const std::span contributionNorms, + const double affineDisagreement, + const double targetDefect, + const double measuredInnerResidual + ) { + GravityPreparedResidualAssessment result; + if (!std::isfinite(affineDisagreement) || affineDisagreement < 0.0 || + !std::isfinite(targetDefect) || targetDefect < 0.0 || + !std::isfinite(measuredInnerResidual) || measuredInnerResidual < 0.0) + return result; + long double scale = 0.0; + for (const double norm : contributionNorms) { + if (!std::isfinite(norm) || norm < 0.0) return result; + scale += norm; + } + // This arithmetic estimate includes large terms that cancel in G. + // It is not a rigorous error bound or a replacement for either the + // independently checked inner residual or the full nonlinear merit. + result.arithmeticAllowance = static_cast(256.0L * std::numeric_limits::epsilon() * scale); + result.targetThreshold = measuredInnerResidual + result.arithmeticAllowance; + result.admissible = std::isfinite(result.arithmeticAllowance) && std::isfinite(result.targetThreshold) && + affineDisagreement <= result.arithmeticAllowance && + targetDefect <= result.targetThreshold; + return result; + } +} diff --git a/libmeanfield/interface/solver/stellar_newton_strategy.cppm b/libmeanfield/interface/solver/stellar_newton_strategy.cppm index 38f0ea1..b77c8bc 100644 --- a/libmeanfield/interface/solver/stellar_newton_strategy.cppm +++ b/libmeanfield/interface/solver/stellar_newton_strategy.cppm @@ -108,6 +108,10 @@ namespace mean_field::solver::detail { expectedResidual(subsystem.GetPartition().GravityResidualSize()), gravityAction(subsystem.GetPartition().GravityResidualSize()), desiredResidual(subsystem.GetPartition().GravityResidualSize()), + sourceResidual(subsystem.GetPartition().GravityResidualSize()), + scaleGravityValues(subsystem.GetPartition().GravityValueSize()), + scaleGravityAction(subsystem.GetPartition().GravityResidualSize()), + scaleRootValues(state.Problem().StateSize()), fullCorrection(state.Problem().StateSize()), fullAction(state.Problem().EquationSize()) { m_options.Validate(); } @@ -119,6 +123,8 @@ namespace mean_field::solver::detail { report.gravityIterations = subsystem.GetIterationCount(); report.gravityFailures = subsystem.GetFailureCount(); report.gravityPreparedFailures = preparedFailures; + report.compactGravityOperator = subsystem.UsesCompactOperator(); + report.compactGravityInverse = subsystem.UsesCompactPreconditioner(); report.lastLinearReportSpace = lastLinearReportSpace; report.gravitySolveSeconds = subsystem.GetSolveSeconds(); if (const auto &last = subsystem.GetLastSolveReport(); last.has_value()) { @@ -126,6 +132,9 @@ namespace mean_field::solver::detail { } report.lastGravityCorrectionNorm = lastCorrectionNorm; report.lastGravityPreparedDefect = lastPreparedDefect; + report.lastGravityAffineDisagreement = lastAffineDisagreement; + report.lastGravityArithmeticAllowance = lastArithmeticAllowance; + report.lastGravityPreparedTargetThreshold = lastPreparedTargetThreshold; report.lastLiftedRelativeTrueResidual = liftedRelativeDefect; } @@ -147,11 +156,37 @@ namespace mean_field::solver::detail { } auto ApplyCorrection(State &state) { + lastPreparedDefect = lastAffineDisagreement = lastArithmeticAllowance = + lastPreparedTargetThreshold = std::numeric_limits::quiet_NaN(); const auto &partition = subsystem.GetPartition(); partition.InjectGravityValues(gravityCorrection, fullCorrection); - state.normalizedOperator->Mult(fullCorrection, fullAction); + subsystem.GetOperator().Mult(gravityCorrection, gravityAction); partition.GatherGravityResiduals(state.trialNormalizedResidual, expectedResidual); - partition.GatherGravityResiduals(fullAction, gravityAction); + sourceResidual = expectedResidual; + // Split the two gravity fields before taking norms: their large + // cancellation is precisely what the residual alone cannot expose. + using Form = typename Problem::FormType; + using namespace utils::blocks; + constexpr value_block> gradient{}; + constexpr value_block> potential{}; + const auto &layout = state.Problem().GetManifest().layout(); + std::array contributionNorms{}; + const auto measure = [&](const mfem::Vector &values, auto block, const bool source) { + scaleRootValues = 0.0; + for (int index = layout.offset(block); index < layout.offset(block) + layout.size(block); ++index) + scaleRootValues(index) = values(index); + partition.GatherGravityValues(scaleRootValues, scaleGravityValues); + subsystem.GetOperator().Mult(scaleGravityValues, scaleGravityAction); + if (source) sourceResidual -= scaleGravityAction; + return StrategyNorm(scaleGravityAction, subsystem.GetCommunicator()); + }; + contributionNorms[0] = measure(state.trialNormalizedState, gradient, true); + contributionNorms[1] = measure(state.trialNormalizedState, potential, true); + contributionNorms[2] = StrategyNorm(sourceResidual, subsystem.GetCommunicator()); + contributionNorms[3] = measure(fullCorrection, gradient, false); + contributionNorms[4] = measure(fullCorrection, potential, false); + contributionNorms[5] = StrategyNorm(expectedResidual, subsystem.GetCommunicator()); + contributionNorms[6] = StrategyNorm(desiredResidual, subsystem.GetCommunicator()); expectedResidual += gravityAction; lastCorrectionNorm = StrategyNorm(gravityCorrection, subsystem.GetCommunicator()); partition.AddGravityValues(gravityCorrection, state.trialNormalizedState); @@ -159,22 +194,21 @@ namespace mean_field::solver::detail { auto preparation = state.PrepareTrial(); if (preparation) { partition.GatherGravityResiduals(state.trialNormalizedResidual, gravityResidual); - const double scale = StrategyNorm(gravityResidual, subsystem.GetCommunicator()) + - StrategyNorm(expectedResidual, subsystem.GetCommunicator()) + - StrategyNorm(desiredResidual, subsystem.GetCommunicator()) + - StrategyNorm(gravityAction, subsystem.GetCommunicator()); + contributionNorms[7] = StrategyNorm(gravityResidual, subsystem.GetCommunicator()); expectedResidual -= gravityResidual; const double disagreement = StrategyNorm(expectedResidual, subsystem.GetCommunicator()); gravityResidual -= desiredResidual; const double targetDefect = StrategyNorm(gravityResidual, subsystem.GetCommunicator()); - const double tolerance = m_options.control.ConvergenceThreshold(lastSolve.rightHandSideNorm) + - 256.0 * std::numeric_limits::epsilon() * scale; // G is affine in gravity for fixed retained variables. Verify // that property against a fresh preparation, including its // boundary rows, before trusting a completed candidate. - trialSolveSucceeded = std::isfinite(disagreement) && disagreement <= tolerance && - std::isfinite(targetDefect) && targetDefect <= tolerance; + const auto assessment = AssessGravityPreparedResidual( + contributionNorms, disagreement, targetDefect, lastSolve.trueResidualNorm); + trialSolveSucceeded = assessment.admissible; lastPreparedDefect = targetDefect; + lastAffineDisagreement = disagreement; + lastArithmeticAllowance = assessment.arithmeticAllowance; + lastPreparedTargetThreshold = assessment.targetThreshold; } if (!preparation || !trialSolveSucceeded) { ++preparedFailures; } return preparation; @@ -183,11 +217,15 @@ namespace mean_field::solver::detail { nonlinear::GravitySolveOptions m_options; Subsystem subsystem; mfem::Vector gravityResidual, gravityCorrection, expectedResidual, gravityAction, desiredResidual; + mfem::Vector sourceResidual, scaleGravityValues, scaleGravityAction, scaleRootValues; mfem::Vector fullCorrection, fullAction; LinearSolveReport lastSolve{}; bool trialSolveSucceeded{true}; double lastCorrectionNorm{0.0}; double lastPreparedDefect{0.0}; + double lastAffineDisagreement{0.0}; + double lastArithmeticAllowance{0.0}; + double lastPreparedTargetThreshold{0.0}; double liftedRelativeDefect{0.0}; std::uint64_t preparedFailures{0}; nonlinear::LinearReportSpace lastLinearReportSpace{nonlinear::LinearReportSpace::full_system}; diff --git a/tests/solver/gravity_reduction.cpp b/tests/solver/gravity_reduction.cpp index 485ea16..08f044b 100644 --- a/tests/solver/gravity_reduction.cpp +++ b/tests/solver/gravity_reduction.cpp @@ -421,3 +421,59 @@ TEST_CASE("Reduced gravity action is repeatable, linear, and alias safe with an operation.Mult(first, secondAction); CHECK(relativeError(firstAction, secondAction) == 0.0); } + +TEST_CASE("Prepared gravity assessment retains the scale of cancelling field contributions", + "[solver][gravity_prepared_residual]") { + using mean_field::solver::detail::AssessGravityPreparedResidual; + const std::array contributions{1.e6, 1.e6, 1.e-8, 2.0, 3.0, 1.e-8, 1.e-8, 1.e-8}; + const double affineDisagreement = 1.e-8; + const double measuredInnerResidual = 5.e-15; + const auto result = AssessGravityPreparedResidual(contributions, affineDisagreement, + affineDisagreement + measuredInnerResidual, measuredInnerResidual); + REQUIRE(result.admissible); + CHECK(result.arithmeticAllowance > affineDisagreement); + CHECK(result.targetThreshold == result.arithmeticAllowance + measuredInnerResidual); + + const std::array residualOnly{1.e-8, 1.e-8, 1.e-8}; + CHECK_FALSE(AssessGravityPreparedResidual(residualOnly, affineDisagreement, + affineDisagreement + measuredInnerResidual, measuredInnerResidual).admissible); + CHECK_FALSE(AssessGravityPreparedResidual(contributions, 2.0 * result.arithmeticAllowance, + measuredInnerResidual, measuredInnerResidual).admissible); + CHECK_FALSE(AssessGravityPreparedResidual(contributions, 0.0, + 2.0 * result.targetThreshold, measuredInnerResidual).admissible); +} + +TEST_CASE("Prepared gravity assessment scales with equation units without an absolute floor", + "[solver][gravity_prepared_residual]") { + using mean_field::solver::detail::AssessGravityPreparedResidual; + const std::array contributions{1.e8, 1.e8, 5.0, 3.0, 2.0}; + const auto reference = AssessGravityPreparedResidual(contributions, 1.e-7, 1.5e-7, 1.e-8); + REQUIRE(reference.admissible); + for (const double factor : {1.e-12, 1.e12}) { + auto scaled = contributions; + for (double &value : scaled) value *= factor; + const auto result = AssessGravityPreparedResidual(scaled, factor * 1.e-7, factor * 1.5e-7, factor * 1.e-8); + CHECK(result.admissible == reference.admissible); + CHECK(result.arithmeticAllowance / factor == Catch::Approx(reference.arithmeticAllowance).epsilon(1.e-14)); + CHECK(result.targetThreshold / factor == Catch::Approx(reference.targetThreshold).epsilon(1.e-14)); + } + const std::array zero{0.0}; + const auto exactZero = AssessGravityPreparedResidual(zero, 0.0, 0.0, 0.0); + CHECK(exactZero.admissible); + CHECK(exactZero.arithmeticAllowance == 0.0); + CHECK_FALSE(AssessGravityPreparedResidual(zero, 1.e-30, 0.0, 0.0).admissible); +} + +TEST_CASE("Prepared gravity assessment rejects nonfinite or invalid measurements", + "[solver][gravity_prepared_residual]") { + using mean_field::solver::detail::AssessGravityPreparedResidual; + const std::array contributions{1.0, 2.0}; + for (const double invalid : {-1.0, std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN()}) { + CHECK_FALSE(AssessGravityPreparedResidual(contributions, invalid, 0.0, 0.0).admissible); + CHECK_FALSE(AssessGravityPreparedResidual(contributions, 0.0, invalid, 0.0).admissible); + CHECK_FALSE(AssessGravityPreparedResidual(contributions, 0.0, 0.0, invalid).admissible); + const std::array invalidContributions{1.0, invalid}; + CHECK_FALSE(AssessGravityPreparedResidual(invalidContributions, 0.0, 0.0, 0.0).admissible); + } +}