From f77c7d2b292882bce362206d5154225d046d21c0 Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <thomas@boudreauxmail.com>
Date: Thu, 17 Sep 2026 10:19:17 -0400
Subject: [PATCH] Compare compact gravity applications against full FEM blocks

---
 CMakeLists.txt                                |  1 +
 experiments/plateau_diagnostic.cpp            | 10 ++++++-
 tests/solver/stellar_numerical_strategies.cpp | 30 +++++++++++++++++++
 3 files changed, 40 insertions(+), 1 deletion(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 50969f1..776bc21 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -373,6 +373,7 @@ add_executable(tests
         tests/solver/stellar_equilibrium_runtime.cpp
         tests/solver/enthalpy_admissibility.cpp
         tests/solver/gravity_subsystem.cpp
+        tests/solver/gravity_compact_applications.cpp
         tests/solver/gravity_reduction.cpp
         tests/solver/stellar_numerical_strategies.cpp
 )
diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp
index 23df102..6b395b2 100644
--- a/experiments/plateau_diagnostic.cpp
+++ b/experiments/plateau_diagnostic.cpp
@@ -21,6 +21,7 @@ import mean_field;
 #include "plateau_pressure_obstruction_diagnostics.hpp"
 #include "plateau_shape_curvature_diagnostics.hpp"
 #include "plateau_numerical_strategy_diagnostics.hpp"
+#include "plateau_gravity_application_diagnostics.hpp"
 #include "plateau_positivity_strategy_diagnostics.hpp"
 
 using namespace mean_field;
@@ -250,8 +251,15 @@ void FocusedProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Ve
         for (const auto &[key, value] : experiment::plateau_diagnostics::MeasureEnthalpyStep(s, f, o.n, &physicalBase))
             bounds << key << ',' << value << '\n';
         bounds.flush();
-        if (o.strategyProbe)
+        if (o.strategyProbe) {
+            s.physicalInverse->Refresh();
+            using Form = typename std::remove_cvref_t<decltype(s.Problem())>::FormType;
+            solver::PreparedGravitySubsystem gravity(*s.normalizedOperator, *s.normalizedInverse,
+                solver::MakeGravityBlockPartition<Form>(s.Problem().GetManifest().layout()),
+                s.Problem().GetCommunicator());
+            plateau_diagnostics::ProbeGravityApplications(gravity, o.output);
             experiment::plateau_diagnostics::ProbeNumericalStrategies(s, f, base, direction, o.output);
+        }
         if (o.positivityProbe)
             experiment::plateau_diagnostics::ProbePolynomialEnthalpy(s, f, base, direction, o.output);
         if (o.surfaceTrace)
diff --git a/tests/solver/stellar_numerical_strategies.cpp b/tests/solver/stellar_numerical_strategies.cpp
index ef9f8ee..3c43735 100644
--- a/tests/solver/stellar_numerical_strategies.cpp
+++ b/tests/solver/stellar_numerical_strategies.cpp
@@ -110,6 +110,34 @@ namespace stellar_numerical_strategy_test {
         CHECK(relativeDifference(before.residual, after.residual) < 1.e-13);
     }
 
+    template <typename Context>
+    void checkGravityApplications(Context &context) {
+        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &state, const auto &) {
+            using Form = typename std::remove_cvref_t<decltype(state.Problem())>::FormType;
+            const auto partition = solver::MakeGravityBlockPartition<Form>(state.Problem().GetManifest().layout());
+            solver::PreparedGravitySubsystem compact(*state.normalizedOperator, *state.normalizedInverse,
+                                                      partition, state.Problem().GetCommunicator());
+            solver::PreparedGravitySubsystem reference(*state.normalizedOperator, *state.normalizedInverse,
+                partition, state.Problem().GetCommunicator(), solver::linear::FGMRES{},
+                solver::RootRestrictedGravityApplications{});
+            REQUIRE(compact.UsesCompactOperator());
+            REQUIRE(compact.UsesCompactPreconditioner());
+            mfem::Vector direction(partition.GravityValueSize()), actual(direction.Size()), expected(direction.Size());
+            // Dense directions expose ordering/scaling mistakes which a few
+            // isolated basis vectors may miss in a distributed field block.
+            for (int trial = 0; trial < 3; ++trial) {
+                for (int i = 0; i < direction.Size(); ++i)
+                    direction(i) = std::sin(.37*(i+1)*(trial+1)) + .4*std::cos(.11*(i+1)*(trial+2));
+                compact.GetOperator().Mult(direction, actual);
+                reference.GetOperator().Mult(direction, expected);
+                CHECK(relativeDifference(expected, actual) < 2.e-12);
+                compact.GetPreconditioner().Mult(direction, actual);
+                reference.GetPreconditioner().Mult(direction, expected);
+                CHECK(relativeDifference(expected, actual) < 2.e-12);
+            }
+        });
+    }
+
     [[nodiscard]] nonlinear::NewtonOptions oneAttempt() {
         // A relative threshold of exactly ||b|| need not stop MFEM before
         // its first iteration. Use a separated absolute threshold so this
@@ -151,6 +179,8 @@ TEST_CASE("Production gravity strategies retain optional model borders and resto
         CAPTURE(CentralConstraint);
         auto context = makeContext<CentralConstraint>();
         const auto original = snapshot(context);
+        checkGravityApplications(context);
+        checkPreserved(original, snapshot(context));
         REQUIRE(original.gravitySize > 0);
         REQUIRE(original.retainedSize > 0);
         // The fixed mass and angular momentum are always retained. Central