Scale prepared gravity consistency by uncancelled field contributions

libmeanfield/interface/solver/newton_numerics.cppm

@@ -1,7 +1,10 @@
 module;
 
 #include <concepts>
+#include <cmath>
 #include <cstdint>
+#include <limits>
+#include <span>
 #include <type_traits>
 #include <utility>
 
@@ -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<double>::quiet_NaN()};
+        double targetThreshold{std::numeric_limits<double>::quiet_NaN()};
+    };
+
+    [[nodiscard]] inline GravityPreparedResidualAssessment AssessGravityPreparedResidual(
+        const std::span<const double> 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<double>(256.0L * std::numeric_limits<double>::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;
+    }
+}