Resolve reduced derivative trend with a recorded four-step diagnostic grid

experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md

@@ -91,3 +91,27 @@ The focused test source `tests/solver/gravity_reduction.cpp` covers:
 These algebraic tests establish the reduction contract. They do not establish
 FEM convergence, a faster implementation, or an additional virial-qualified
 stellar equilibrium; those require the integrated runtime experiments.
+
+## Resolved FEM difference sampling
+
+The first frozen n=1.5 probe completed all gravity solves and preparations. Its
+central-difference errors at steps 0.01, 0.001, and 0.0001 were respectively
+`1.7761e-9`, `1.9661e-11`, and `7.1208e-11`. Only the first exceeded twenty times
+the recorded arithmetic scale. The observed first-interval order was about
+1.956, but the conservative checker correctly left the trend unresolved because
+the finer endpoint did not meet its resolution criterion.
+
+The diagnostic now uses the explicit, validated step array `{0.1, 0.03, 0.01,
+0.001}` and records it in the summary. The larger steps test whether truncation
+dominates sufficiently to reveal a central-difference trend; retaining 0.001
+also preserves an independent small-step accuracy check. This adds one central
+pair, hence at most four extra gravity correction solves under the existing
+two-correction completion budget. The inner tolerance, prepared completion gate,
+and checker thresholds are unchanged. There is no assumption that larger steps
+will prepare successfully or automatically pass the trend check.
+
+These remain finite-range empirical derivative checks. The recorded arithmetic
+scale is heuristic, and a gravity residual bound by itself is not a bound on its
+amplified effect in the retained equations. The larger grid is intended to
+resolve a measurable trend without claiming such an unmeasured amplification
+bound or interpreting floor-limited differences as derivative failures.

experiments/plateau_numerical_strategy_diagnostics.hpp

@@ -19,6 +19,9 @@ struct NumericalStrategyProbeOptions final {
     int restartLength{80};
     int maximumCompletionCorrections{2};
     double gravityCompletionTolerance{1.e-11};
+    // Coarser differences resolve truncation before cancellation dominates;
+    // the final small step independently checks derivative accuracy.
+    std::array<double, 4> reducedDifferenceSteps{.1, .03, .01, .001};
 };
 
 namespace numerical_strategy_detail {
@@ -55,6 +58,12 @@ void ProbeNumericalStrategies(
     if (options.restartLength <= 0 || options.maximumCompletionCorrections <= 0 ||
         !std::isfinite(options.gravityCompletionTolerance) || options.gravityCompletionTolerance <= 0.)
         throw std::invalid_argument("Numerical strategy probes require valid inner and completion controls.");
+    double previousDifferenceStep = std::numeric_limits<double>::infinity();
+    for (const double step : options.reducedDifferenceSteps) {
+        if (!std::isfinite(step) || step <= 0. || step >= previousDifferenceStep)
+            throw std::invalid_argument("Reduced finite-difference steps must be finite, positive, and strictly decreasing.");
+        previousDifferenceStep = step;
+    }
     int ranks = 0;
     MPI_Comm_size(state.Problem().GetCommunicator(), &ranks);
     if (ranks != 1 || !finiteElements.mesh || normalizedBase.Size() != state.Problem().StateSize() ||
@@ -92,6 +101,8 @@ void ProbeNumericalStrategies(
             << "\ninner_maximum_iterations," << options.inner.maximumIterations
             << "\ninner_restart_length," << options.restartLength
             << "\ngravity_completion_tolerance," << options.gravityCompletionTolerance << '\n';
+    for (std::size_t index = 0; index < options.reducedDifferenceSteps.size(); ++index)
+        summary << "reduced_difference_step_" << index << ',' << options.reducedDifferenceSteps[index] << '\n';
     const double nan = std::numeric_limits<double>::quiet_NaN();
     int preparationCount = 0;
     double restorationError = 0.;
@@ -281,7 +292,7 @@ void ProbeNumericalStrategies(
                         << "\nlifted_gravity_action_norm," << gravityResidual.Norml2()
                         << "\nlifted_retained_action_difference," << d::Difference(retainedAction, reducedAction)
                         << "\ncompleted_base_residual_norm," << completedBaseResidual.Norml2() << '\n';
-                for (const double alpha : {1.e-2, 1.e-3, 1.e-4}) {
+                for (const double alpha : options.reducedDifferenceSteps) {
                     std::array<mfem::Vector, 2> residuals{mfem::Vector(base.Size()), mfem::Vector(base.Size())};
                     std::array<bool, 2> completed{};
                     std::array<double, 2> gravityNorms{nan, nan};