diff --git a/experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md b/experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md index 9138be2..fe05fe1 100644 --- a/experiments/GRAVITY_REDUCTION_IMPLEMENTATION_2026-09-17.md +++ b/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. diff --git a/experiments/plateau_numerical_strategy_diagnostics.hpp b/experiments/plateau_numerical_strategy_diagnostics.hpp index f9664a4..71f3595 100644 --- a/experiments/plateau_numerical_strategy_diagnostics.hpp +++ b/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 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::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::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 residuals{mfem::Vector(base.Size()), mfem::Vector(base.Size())}; std::array completed{}; std::array gravityNorms{nan, nan};