Compare compact gravity applications against full FEM blocks

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
 )

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)

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