Add frozen surface and scalar Schur preconditioner experiment

CMakeLists.txt

@@ -389,6 +389,7 @@ add_executable(tests
         tests/solver/stellar_equilibrium_runtime.cpp
         tests/solver/enthalpy_admissibility.cpp
         tests/solver/plateau_coupled_coarse.cpp
+        tests/solver/plateau_surface_border_schur.cpp
         tests/solver/gravity_subsystem.cpp
         tests/solver/gravity_compact_applications.cpp
         tests/solver/gravity_reduction.cpp
@@ -417,6 +418,10 @@ add_executable(continuation_contract_tests EXCLUDE_FROM_ALL
         tests/solver/stellar_numerical_strategies.cpp
         tests/solver/continuation_snapshot.cpp)
 target_link_libraries(continuation_contract_tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost)
+add_executable(surface_border_schur_contract_tests EXCLUDE_FROM_ALL
+        tests/test_main.cpp
+        tests/solver/plateau_surface_border_schur.cpp)
+target_link_libraries(surface_border_schur_contract_tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost)
 
 # A deliberately opt-in API workbench. It is compiled explicitly during API
 # verification but temporary user experiments do not break the default build.

experiments/plateau_preconditioner_benchmark.hpp

@@ -14,6 +14,7 @@
 #include <memory>
 
 #include "plateau_coupled_coarse.hpp"
+#include "plateau_surface_border_schur.hpp"
 #include <stdexcept>
 #include <string>
 #include <type_traits>
@@ -31,7 +32,9 @@ enum class PreconditionerBenchmarkVariant {
     gravity_schur,
     stellar_ldu_gravity_schur,
     coupled_shape_coarse,
-    gravity_schur_coarse
+    gravity_schur_coarse,
+    surface_border_schur,
+    gravity_schur_surface_border_schur
 };
 
 inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const std::string &name) {
@@ -44,6 +47,9 @@ inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const
     if (name == "stellar-ldu-gravity-schur") return PreconditionerBenchmarkVariant::stellar_ldu_gravity_schur;
     if (name == "coupled-shape-coarse") return PreconditionerBenchmarkVariant::coupled_shape_coarse;
     if (name == "gravity-schur-coarse") return PreconditionerBenchmarkVariant::gravity_schur_coarse;
+    if (name == "surface-border-schur") return PreconditionerBenchmarkVariant::surface_border_schur;
+    if (name == "gravity-schur-surface-border-schur")
+        return PreconditionerBenchmarkVariant::gravity_schur_surface_border_schur;
     throw std::invalid_argument("Unknown preconditioner benchmark variant: " + name);
 }
 
@@ -271,6 +277,12 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
                "true_check_jacobian_applications,krylov_other_seconds\n";
     auto coarseReport = detail::Open(output / "preconditioner_coarse.csv");
     coarseReport << "role,variant,stage,metric,value\n";
+    auto schurReport = detail::Open(output / "preconditioner_surface_border_schur.csv");
+    schurReport << "role,variant,stage,metric,value\n";
+    auto schurCoordinates = detail::Open(output / "preconditioner_surface_border_coordinates.csv");
+    schurCoordinates << "role,variant,coarse_index,residual_index,correction_index\n";
+    auto schurMatrices = detail::Open(output / "preconditioner_surface_border_matrices.csv");
+    schurMatrices << "role,variant,row,column,uneliminated,schur\n";
     auto work = detail::Open(output / "preconditioner_work.csv");
     work << "role,variant,stage,metric,count\n";
     auto actions = detail::Open(output / "preconditioner_actions.csv");
@@ -285,7 +297,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
     history << "role,variant,iteration,reported_norm,final\n";
 
     const auto run = [&]<pc::StellarStructureFactorizationPolicy Policy, bool MaterialLDU = false,
-                         bool RefinedGravity = false, bool UseCoarse = false>(
+                         bool RefinedGravity = false, bool UseCoarse = false, bool UseSurfaceBorderSchur = false>(
                          const std::string &role, const std::string &variant) {
         const auto candidateStart = detail::Clock::now();
         auto physical = detail::Prepare<Policy, MaterialLDU, RefinedGravity>(state.Problem(), options);
@@ -319,8 +331,38 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         const double bindingSeconds = detail::Seconds(sectionStart);
         using NormalizedOperation = std::remove_reference_t<decltype(*state.normalizedOperator)>;
         using Coarse = PreparedCoupledShapeCoarse<NormalizedOperation, decltype(normalized)>;
+        using UnborderedStructure = NormalizedUnborderedStructureInverse<NormalizedOperation, decltype(physical)>;
+        using SurfaceBorderSchur = PreparedSurfaceBorderSchur<NormalizedOperation, UnborderedStructure>;
         std::unique_ptr<Coarse> coarse;
+        std::unique_ptr<UnborderedStructure> unborderedStructure;
+        std::unique_ptr<SurfaceBorderSchur> surfaceBorderSchur;
+        static_assert(!(UseCoarse && UseSurfaceBorderSchur));
         mfem::Solver *applicationInverse = &normalized;
+        const auto recordSurfaceBorderSchur = [&](const char *stage) {
+            if (!surfaceBorderSchur) return;
+            const auto &stats = surfaceBorderSchur->GetStatistics();
+            const auto metric = [&](const char *name, auto value) {
+                schurReport << role << ',' << variant << ',' << stage << ',' << name << ',' << value << '\n';
+            };
+            metric("status", static_cast<int>(stats.status));
+            metric("coarse_size", stats.coarseSize);
+            metric("reciprocal_condition", stats.reciprocalCondition);
+            metric("minimum_scaled_pivot", stats.minimumScaledPivot);
+            metric("setup_jacobian_applications", stats.setupOperatorApplications);
+            metric("setup_inverse_applications", stats.setupInverseApplications);
+            metric("setup_seconds", stats.lastSetupSeconds);
+            metric("setup_jacobian_seconds", stats.setupOperatorSeconds);
+            metric("setup_inverse_seconds", stats.setupInverseSeconds);
+            metric("persistent_vector_bytes", stats.persistentVectorBytes);
+            metric("applications", stats.applications);
+            metric("internal_jacobian_applications", stats.operatorApplications);
+            metric("internal_inverse_applications", stats.inverseApplications);
+            metric("internal_jacobian_seconds", stats.operatorSeconds);
+            metric("internal_inverse_seconds", stats.inverseSeconds);
+            metric("application_seconds", stats.applicationSeconds);
+            metric("failed_applications", stats.failedApplications);
+            schurReport.flush();
+        };
         const auto recordCoarse = [&](const char *stage) {
             if (!coarse) return;
             const auto &stats = coarse->GetStatistics();
@@ -398,6 +440,34 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
                 throw std::runtime_error("The actual FEM coarse inverse failed its exact-subspace check.");
         }
 
+        if constexpr (UseSurfaceBorderSchur) {
+            sectionStart = detail::Clock::now();
+            unborderedStructure = std::make_unique<UnborderedStructure>(*state.normalizedOperator, physical);
+            surfaceBorderSchur = std::make_unique<SurfaceBorderSchur>(
+                *state.normalizedOperator, *unborderedStructure, MakeSurfaceBorderSchurCoordinates(physical),
+                state.Problem().GetCommunicator());
+            const double extraConstructionSeconds = detail::Seconds(sectionStart);
+            constructionSeconds += extraConstructionSeconds;
+            const auto &coordinates = surfaceBorderSchur->GetCoordinates();
+            const auto &schur = surfaceBorderSchur->GetSchurMatrix();
+            const auto &uneliminated = surfaceBorderSchur->GetUneliminatedCoarseMatrix();
+            for (int row = 0; row < schur.Height(); ++row) {
+                schurCoordinates << role << ',' << variant << ',' << row << ','
+                                 << coordinates.residualIndices[row] << ','
+                                 << coordinates.correctionIndices[row] << '\n';
+                for (int column = 0; column < schur.Width(); ++column)
+                    schurMatrices << role << ',' << variant << ',' << row << ',' << column << ','
+                                  << uneliminated(row, column) << ',' << schur(row, column) << '\n';
+            }
+            schurCoordinates.flush();
+            schurMatrices.flush();
+            recordSurfaceBorderSchur("constructed");
+            recordWork("after_surface_border_setup");
+            if (!surfaceBorderSchur->IsCurrent())
+                throw std::runtime_error("The measured surface/scalar Schur matrix is singular or ill-conditioned; no fallback was applied.");
+            applicationInverse = surfaceBorderSchur.get();
+        }
+
         mfem::Vector inverseAction(rhs.Size()), inverseJacobianAction(rhs.Size());
         sectionStart = detail::Clock::now();
         applicationInverse->Mult(rhs, inverseAction);
@@ -412,11 +482,23 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         if (!detail::Finite(inverseJacobianAction) || !std::isfinite(actionNorm))
             throw std::runtime_error("The diagnostic inverse Jacobian image or its norm is nonfinite.");
         inverseJacobianAction -= rhs;
+        if (surfaceBorderSchur) {
+            double squared = 0.0;
+            for (int row : surfaceBorderSchur->GetCoordinates().residualIndices)
+                squared += inverseJacobianAction(row) * inverseJacobianAction(row);
+            const double coarseRelative = std::sqrt(squared) / rhsNorm;
+            schurReport << role << ',' << variant
+                        << ",validation,coarse_equation_defect_relative_to_full_rhs," << coarseRelative << '\n';
+            schurReport.flush();
+            if (!std::isfinite(coarseRelative) || coarseRelative > 1.e-8)
+                throw std::runtime_error("The Schur inverse failed the authoritative-J coarse-equation check.");
+        }
         mfem::Vector referenceDifference(inverseAction);
         referenceDifference -= referenceCorrection;
         const double referenceRelative = detail::Relative(referenceDifference.Norml2(),
             std::max(inverseAction.Norml2(), referenceCorrection.Norml2()));
         const bool requireEquivalent = std::same_as<Policy, pc::IndependentStellarSubsystems> && !MaterialLDU && !UseCoarse &&
+            !UseSurfaceBorderSchur &&
             (!RefinedGravity || options.schurCorrectionSweeps == 0);
         actions << role << ',' << variant << ',' << inverseAction.Norml2() << ',' << actionNorm << ','
                 << inverseJacobianAction.Norml2() << ',' << inverseJacobianAction.Norml2() / rhsNorm << ','
@@ -447,6 +529,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         const double krylovSetupSeconds = detail::Seconds(sectionStart);
         recordWork("before_krylov");
         recordCoarse("before_krylov");
+        recordSurfaceBorderSchur("before_krylov");
         countedJacobian.ResetStatistics();
         countedInverse.ResetStatistics();
         sectionStart = detail::Clock::now();
@@ -454,6 +537,7 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         const double krylovSeconds = detail::Seconds(sectionStart);
         recordWork("after_krylov");
         recordCoarse("after_krylov");
+        recordSurfaceBorderSchur("after_krylov");
         const auto jacobianStats = countedJacobian.GetStatistics();
         const auto inverseStats = countedInverse.GetStatistics();
         const auto lifecycle = countedInverse.GetLifecycleStatistics();
@@ -537,6 +621,10 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
         run.template operator()<pc::IndependentStellarSubsystems, false, false, true>("candidate", options.variant); break;
     case PreconditionerBenchmarkVariant::gravity_schur_coarse:
         run.template operator()<pc::IndependentStellarSubsystems, false, true, true>("candidate", options.variant); break;
+    case PreconditionerBenchmarkVariant::surface_border_schur:
+        run.template operator()<pc::IndependentStellarSubsystems, false, false, false, true>("candidate", options.variant); break;
+    case PreconditionerBenchmarkVariant::gravity_schur_surface_border_schur:
+        run.template operator()<pc::IndependentStellarSubsystems, false, true, false, true>("candidate", options.variant); break;
     }
     auto scope = detail::Open(output / "preconditioner_scope.txt");
     scope << "Both solves use fresh physical inverses, separately completed scalar borders, the same physical J,\n"
@@ -549,6 +637,12 @@ void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
           << "unchanged_state_refresh_seconds measures an immediate Refresh; it is not a changed-state refresh benchmark.\n"
           << "Each candidate receives one diagnostic inverse/Jv application before the timed cold Krylov solve.\n"
           << "Application timers/counts in the summary include only Krylov work; the independent true check is separate.\n"
+          << "Surface/scalar Schur setup is included in construction; its internal J applications are separately counted.\n"
+          << "For that inverse, outer Krylov J counts exclude internal Schur J work; inverse seconds include it.\n"
+          << "Schur status codes are 0=unprepared, 1=ready, 2=singular, 3=ill-conditioned; rejected setup has no fallback.\n"
+          << "The Schur partition uses every native surface parameter and generated scalar block, including a central border if present.\n"
+          << "Its field inverse restricts the unbordered structure inverse, not an exact field elimination.\n"
+          << "Existing scalar-border construction is retained in setup cost but its inverse is not applied by the new wrapper.\n"
           << "setup_plus_solve_seconds includes construction, normalization binding, Krylov setup, and Krylov solve.\n"
           << "setup_solve_and_true_check_seconds adds the separately measured original-J true-residual gate.\n"
           << "measured_candidate_seconds also includes diagnostic refresh, coarse exactness checks, and inverse/Jv comparison.\n"

experiments/plateau_surface_border_schur.hpp

@@ -0,0 +1,432 @@
+#pragma once
+
+// Include after import mean_field. The experiment changes only the inverse
+// supplied to Krylov; every column below comes from the authoritative Jacobian.
+#include <algorithm>
+#include <cmath>
+#include <concepts>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <span>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include <mfem.hpp>
+#include <mpi.h>
+
+namespace experiment::plateau_diagnostics {
+    struct SurfaceBorderSchurCoordinates final {
+        // Equation and unknown ordering need not coincide in a generated form.
+        std::vector<int> residualIndices;
+        std::vector<int> correctionIndices;
+    };
+
+    template <class PhysicalInverse>
+    [[nodiscard]] SurfaceBorderSchurCoordinates MakeSurfaceBorderSchurCoordinates(const PhysicalInverse &inverse) {
+        using namespace mean_field;
+        using namespace utils::blocks;
+        using Root = std::remove_cvref_t<decltype(inverse.GetBlock())>;
+        using Structure = typename Root::StructureComponent;
+        constexpr auto shapeColumn = type_index_v<surface_deformation::parameters::value,
+                                                  typename Root::CorrectionBlocks>;
+        constexpr auto shapeRow = type_index_v<surface_deformation::shape_equilibrium::residual,
+                                               typename Root::ResidualBlocks>;
+        static_assert(Root::borderValueArity == Root::borderResidualArity);
+        static_assert(shapeColumn < Structure::CorrectionBlocks::size);
+        static_assert(shapeRow < Structure::ResidualBlocks::size);
+        const auto &map = inverse.GetCoordinateMap();
+        const auto &columns = map.GetCorrectionRanges();
+        const auto &rows = map.GetResidualRanges();
+        SurfaceBorderSchurCoordinates result;
+        const auto append = [](const auto &range, std::vector<int> &indices) {
+            for (int i = 0; i < range.size; ++i)
+                indices.push_back(range.equilibriumOffset + i);
+        };
+        append(columns[shapeColumn], result.correctionIndices);
+        append(rows[shapeRow], result.residualIndices);
+        // Generated scalar blocks follow the complete structure descriptor.
+        // This includes the optional central phase condition without naming a
+        // particular model or interpreting a runtime stable-ID string.
+        for (std::size_t i = Structure::CorrectionBlocks::size; i < columns.size(); ++i)
+            append(columns[i], result.correctionIndices);
+        for (std::size_t i = Structure::ResidualBlocks::size; i < rows.size(); ++i)
+            append(rows[i], result.residualIndices);
+        return result;
+    }
+
+    template <class Candidate>
+    concept SurfaceBorderSchurOperation = std::derived_from<Candidate, mfem::Operator> &&
+        requires(const Candidate &operation) {
+            { operation.IsPrepared() } -> std::same_as<bool>;
+            { operation.GetProblem().GetPreparationGeneration() } -> std::convertible_to<std::uint64_t>;
+        };
+
+    template <class Candidate>
+    concept SurfaceBorderSchurInverse = std::derived_from<Candidate, mfem::Solver> &&
+        mean_field::solver::LinearPreconditionerApplicationContractAvailable<Candidate> &&
+        (mean_field::solver::linearPreconditionerApplicationContract<Candidate> ==
+         mean_field::preconditioning::ApplicationContract::stationary_linear) &&
+        requires(const Candidate &inverse) {
+            { inverse.IsCurrent() } -> std::same_as<bool>;
+        };
+
+    // Reusing the complete scalar-bordered inverse here would impose zero
+    // scalar-constraint RHS before elimination. Its field restriction can be
+    // singular even when A_FF is invertible (a zero scalar diagonal suffices).
+    // Borrow only the original structure inverse, preserving its typed packing
+    // and the exact physical/normalized coordinate maps.
+    template <SurfaceBorderSchurOperation Operation, SurfaceBorderSchurInverse PhysicalInverse>
+    class NormalizedUnborderedStructureInverse final : public mfem::Solver {
+        using Root = std::remove_cvref_t<decltype(std::declval<const PhysicalInverse &>().GetBlock())>;
+        using Backend = typename Root::StructureComponent::BackendType;
+    public:
+        static constexpr auto applicationContract = mean_field::preconditioning::backend::applicationContract<Backend>;
+        static_assert(applicationContract == mean_field::preconditioning::ApplicationContract::stationary_linear);
+
+        NormalizedUnborderedStructureInverse(const Operation &operation, const PhysicalInverse &inverse)
+            : mfem::Solver(operation.Width(), operation.Height()), m_operation(&operation), m_inverse(&inverse),
+              m_generation(operation.GetProblem().GetPreparationGeneration()), m_physicalRhs(Width()),
+              m_physicalAction(Height()), m_groupedRhs(Width()), m_groupedAction(Height()) {
+            if (&operation.GetProblem() != &inverse.GetProblem() || inverse.Height() != Height() ||
+                inverse.Width() != Width() || !IsCurrent())
+                throw std::invalid_argument("The structure adapter requires the current inverse of the same problem.");
+        }
+        void SetOperator(const mfem::Operator &operation) override {
+            if (&operation != m_operation || !IsCurrent())
+                throw std::invalid_argument("The structure adapter cannot bind another or stale Jacobian.");
+        }
+        [[nodiscard]] bool IsCurrent() const {
+            return m_operation->IsPrepared() && m_inverse->IsCurrent() && !m_inverse->iterative_mode &&
+                   !m_inverse->GetGroupedPreconditioner().GetStructurePreconditioner().iterative_mode &&
+                   m_generation == m_operation->GetProblem().GetPreparationGeneration();
+        }
+        [[nodiscard]] const Operation &GetNormalizedJacobian() const noexcept { return *m_operation; }
+        void Mult(const mfem::Vector &rhs, mfem::Vector &action) const override {
+            if (!IsCurrent() || rhs.Size() != Width() || action.Size() != Height())
+                throw std::logic_error("The structure adapter requires current objects and compatible vectors.");
+            const auto &structure = m_inverse->GetGroupedPreconditioner().GetStructurePreconditioner();
+            const auto &map = m_inverse->GetCoordinateMap();
+            m_operation->DenormalizeResidual(rhs, m_physicalRhs);
+            map.PackResidual(m_physicalRhs, m_groupedRhs);
+            m_groupedAction = 0.0;
+            const mfem::Vector fieldRhs(m_groupedRhs, 0, structure.Width());
+            mfem::Vector fieldAction(m_groupedAction, 0, structure.Height());
+            structure.Mult(fieldRhs, fieldAction);
+            fieldAction.SyncAliasMemory(m_groupedAction);
+            map.UnpackCorrection(m_groupedAction, m_physicalAction);
+            m_operation->NormalizeState(m_physicalAction, action);
+        }
+    private:
+        const Operation *m_operation;
+        const PhysicalInverse *m_inverse;
+        std::uint64_t m_generation;
+        mutable mfem::Vector m_physicalRhs, m_physicalAction, m_groupedRhs, m_groupedAction;
+    };
+
+    enum class SurfaceBorderSchurStatus { unprepared, ready, singular, ill_conditioned };
+
+    struct SurfaceBorderSchurOptions final {
+        double minimumReciprocalCondition{1.e-14};
+    };
+
+    struct SurfaceBorderSchurStatistics final {
+        SurfaceBorderSchurStatus status{SurfaceBorderSchurStatus::unprepared};
+        int coarseSize{0};
+        std::uint64_t setups{0}, setupOperatorApplications{0}, setupInverseApplications{0};
+        double reciprocalCondition{0.0};
+        double minimumScaledPivot{0.0};
+        double lastSetupSeconds{0.0}, setupOperatorSeconds{0.0}, setupInverseSeconds{0.0};
+        std::uint64_t applications{0}, operatorApplications{0}, inverseApplications{0}, failedApplications{0};
+        double applicationSeconds{0.0}, operatorSeconds{0.0}, inverseSeconds{0.0};
+        std::size_t persistentVectorBytes{0};
+    };
+
+    /* With C the complete surface/scalar block and F every remaining field,
+     * B = R_F P0 E_F approximates A_FF^{-1}. Store
+     *
+     *   Z = E_C - E_F B A_FC,       S = R_C A Z.
+     *
+     * Application gives y = E_F B r_F and x = y + Z S^{-1}(r_C-R_C A y).
+     * The coarse equations are satisfied exactly (up to arithmetic), even
+     * when B is approximate. The complete inverse is exact only when B is.
+     * Restricting the existing inverse defines B; its internal border response
+     * is not claimed to be an exact elimination of the frozen field block.
+     * No symmetry, positive definiteness, or reduced physical equation is
+     * assumed. A singular/unsafe S is reported and cannot be applied.
+     */
+    template <SurfaceBorderSchurOperation Operation, SurfaceBorderSchurInverse Inverse>
+    class PreparedSurfaceBorderSchur final : public mfem::Solver {
+    public:
+        static constexpr auto applicationContract = mean_field::preconditioning::ApplicationContract::stationary_linear;
+        static constexpr std::size_t maximumCoarseSize = 128;
+
+        PreparedSurfaceBorderSchur(const Operation &operation, const Inverse &inverse,
+                                   SurfaceBorderSchurCoordinates coordinates, MPI_Comm communicator,
+                                   SurfaceBorderSchurOptions options = {})
+            : mfem::Solver(operation.Width(), operation.Height()), m_operation(&operation), m_inverse(&inverse),
+              m_coordinates(std::move(coordinates)), m_communicator(communicator), m_options(options),
+              m_rhs(operation.Height()), m_fieldAction(operation.Width()), m_image(operation.Height()) {
+            Validate();
+            Refresh();
+        }
+
+        PreparedSurfaceBorderSchur(const PreparedSurfaceBorderSchur &) = delete;
+        PreparedSurfaceBorderSchur &operator=(const PreparedSurfaceBorderSchur &) = delete;
+        PreparedSurfaceBorderSchur(PreparedSurfaceBorderSchur &&) = delete;
+        PreparedSurfaceBorderSchur &operator=(PreparedSurfaceBorderSchur &&) = delete;
+
+        void SetOperator(const mfem::Operator &operation) override {
+            if (&operation != m_operation)
+                throw std::invalid_argument("The surface-border Schur inverse cannot bind a different Jacobian.");
+            RequireCurrent();
+        }
+
+        [[nodiscard]] bool IsCurrent() const {
+            return m_statistics.status == SurfaceBorderSchurStatus::ready && DependenciesCurrent();
+        }
+        [[nodiscard]] const SurfaceBorderSchurStatistics &GetStatistics() const noexcept { return m_statistics; }
+        [[nodiscard]] const SurfaceBorderSchurCoordinates &GetCoordinates() const noexcept { return m_coordinates; }
+        [[nodiscard]] const mfem::DenseMatrix &GetSchurMatrix() const noexcept { return m_schur; }
+        [[nodiscard]] const mfem::DenseMatrix &GetUneliminatedCoarseMatrix() const noexcept { return m_coarse; }
+        [[nodiscard]] std::span<const mfem::Vector> GetCompletedColumns() const noexcept { return m_columns; }
+
+        void Refresh() {
+            const double start = MPI_Wtime();
+            Validate();
+            m_statistics.status = SurfaceBorderSchurStatus::unprepared;
+            m_generation = m_operation->GetProblem().GetPreparationGeneration();
+            const int count = static_cast<int>(m_coordinates.correctionIndices.size());
+            m_statistics.coarseSize = count;
+            m_columns.clear();
+            m_columns.reserve(count);
+            m_schur.SetSize(count);
+            m_coarse.SetSize(count);
+            m_solveWork.SetSize(count);
+            mfem::Vector basis(Height()), image(Width()), response(Height());
+            for (int column = 0; column < count; ++column) {
+                basis = 0.0;
+                basis(m_coordinates.correctionIndices[column]) = 1.0;
+                ApplyOperation(basis, image, true);
+                for (int row = 0; row < count; ++row)
+                    m_coarse(row, column) = image(m_coordinates.residualIndices[row]);
+                ApplyFields(image, response, true);
+                response *= -1.0;
+                response += basis;
+                m_columns.emplace_back(response);
+                ApplyOperation(response, image, true);
+                for (int row = 0; row < count; ++row)
+                    m_schur(row, column) = image(m_coordinates.residualIndices[row]);
+            }
+            Factor();
+            if (!DependenciesCurrent())
+                throw std::logic_error("The surface-border Schur dependencies changed during preparation.");
+            m_statistics.lastSetupSeconds = MPI_Wtime() - start;
+            m_statistics.persistentVectorBytes = sizeof(mfem::real_t) *
+                (static_cast<std::size_t>(count) * Height() + 3 * static_cast<std::size_t>(Height()) + count);
+            ++m_statistics.setups;
+        }
+
+        void Mult(const mfem::Vector &rightHandSide, mfem::Vector &action) const override {
+            const double start = MPI_Wtime();
+            ++m_statistics.applications;
+            try {
+                RequireCurrent();
+                if (rightHandSide.Size() != Width() || action.Size() != Height())
+                    throw std::invalid_argument("The surface-border Schur inverse requires preallocated root vectors.");
+                RequireFinite(rightHandSide);
+                ApplyFields(rightHandSide, m_fieldAction, false);
+                ApplyOperation(m_fieldAction, m_image, false);
+                for (int row = 0; row < m_solveWork.Size(); ++row)
+                    m_solveWork(row) = rightHandSide(m_coordinates.residualIndices[row]) -
+                                       m_image(m_coordinates.residualIndices[row]);
+                SolveFactored(m_solveWork);
+                // No write to the caller's output occurs before all RHS uses,
+                // so in-place application has the same algebra as separate buffers.
+                action = m_fieldAction;
+                for (int column = 0; column < m_solveWork.Size(); ++column)
+                    action.Add(m_solveWork(column), m_columns[column]);
+                RequireFinite(action);
+                RequireCurrent();
+                m_statistics.applicationSeconds += MPI_Wtime() - start;
+            } catch (...) {
+                ++m_statistics.failedApplications;
+                m_statistics.applicationSeconds += MPI_Wtime() - start;
+                throw;
+            }
+        }
+
+    private:
+        void Validate() const {
+            int ranks = 0;
+            if (m_communicator == MPI_COMM_NULL || MPI_Comm_size(m_communicator, &ranks) != MPI_SUCCESS || ranks != 1)
+                throw std::invalid_argument("The surface-border Schur experiment currently requires one MPI rank.");
+            if (!m_operation->IsPrepared() || !m_inverse->IsCurrent() || m_inverse->iterative_mode ||
+                m_operation->Height() != m_operation->Width() || m_inverse->Height() != Height() ||
+                m_inverse->Width() != Width())
+                throw std::invalid_argument("The surface-border Schur experiment requires current compatible objects.");
+            if constexpr (requires { m_inverse->GetNormalizedJacobian(); }) {
+                if (&m_inverse->GetNormalizedJacobian() != m_operation)
+                    throw std::invalid_argument("The inverse belongs to a different normalized Jacobian.");
+            }
+            if (!std::isfinite(m_options.minimumReciprocalCondition) ||
+                m_options.minimumReciprocalCondition <= 0.0 || m_options.minimumReciprocalCondition >= 1.0)
+                throw std::invalid_argument("The Schur reciprocal-condition threshold must lie strictly between zero and one.");
+            if (m_coordinates.residualIndices.empty() ||
+                m_coordinates.residualIndices.size() != m_coordinates.correctionIndices.size() ||
+                m_coordinates.residualIndices.size() >= static_cast<std::size_t>(Height()))
+                throw std::invalid_argument("The Schur partition requires nonempty fields and square coarse coordinates.");
+            if (m_coordinates.residualIndices.size() > maximumCoarseSize)
+                throw std::invalid_argument("The bounded surface-border experiment supports at most 128 coarse coordinates.");
+            const auto verify = [this](const auto &indices) {
+                std::vector<bool> used(Height(), false);
+                for (int index : indices) {
+                    if (index < 0 || index >= Height() || used[index])
+                        throw std::invalid_argument("Schur coordinate indices must be distinct and in range.");
+                    used[index] = true;
+                }
+            };
+            verify(m_coordinates.residualIndices);
+            verify(m_coordinates.correctionIndices);
+        }
+
+        [[nodiscard]] bool DependenciesCurrent() const {
+            return m_operation->IsPrepared() && m_inverse->IsCurrent() && !m_inverse->iterative_mode &&
+                   m_generation == m_operation->GetProblem().GetPreparationGeneration();
+        }
+        void RequireCurrent() const {
+            if (!IsCurrent())
+                throw std::logic_error("The surface-border Schur inverse is stale, singular, or ill-conditioned.");
+        }
+        static void RequireFinite(const mfem::Vector &vector) {
+            for (int i = 0; i < vector.Size(); ++i)
+                if (!std::isfinite(vector(i)))
+                    throw std::runtime_error("The surface-border Schur experiment encountered a nonfinite vector.");
+        }
+        void ApplyOperation(const mfem::Vector &input, mfem::Vector &output, bool setup) const {
+            const double start = MPI_Wtime();
+            m_operation->Mult(input, output);
+            const double seconds = MPI_Wtime() - start;
+            if (setup) {
+                ++m_statistics.setupOperatorApplications;
+                m_statistics.setupOperatorSeconds += seconds;
+            } else {
+                ++m_statistics.operatorApplications;
+                m_statistics.operatorSeconds += seconds;
+            }
+            RequireFinite(output);
+        }
+        void ApplyFields(const mfem::Vector &input, mfem::Vector &output, bool setup) const {
+            m_rhs = input;
+            for (int index : m_coordinates.residualIndices)
+                m_rhs(index) = 0.0;
+            const double start = MPI_Wtime();
+            m_inverse->Mult(m_rhs, output);
+            const double seconds = MPI_Wtime() - start;
+            if (setup) {
+                ++m_statistics.setupInverseApplications;
+                m_statistics.setupInverseSeconds += seconds;
+            } else {
+                ++m_statistics.inverseApplications;
+                m_statistics.inverseSeconds += seconds;
+            }
+            RequireFinite(output);
+            for (int index : m_coordinates.correctionIndices)
+                output(index) = 0.0;
+        }
+
+        void Factor() {
+            const int count = m_schur.Height();
+            m_lu = m_schur;
+            m_pivots.resize(count);
+            m_scale = 0.0;
+            for (int i = 0; i < count * count; ++i) {
+                if (!std::isfinite(m_lu.Data()[i]))
+                    throw std::runtime_error("The coarse Schur matrix is nonfinite.");
+                m_scale = std::max(m_scale, std::abs(m_lu.Data()[i]));
+            }
+            m_statistics.reciprocalCondition = 0.0;
+            m_statistics.minimumScaledPivot = 0.0;
+            m_statistics.status = SurfaceBorderSchurStatus::singular;
+            if (m_scale == 0.0)
+                return;
+            for (int i = 0; i < count * count; ++i)
+                m_lu.Data()[i] /= m_scale;
+            double norm = 0.0;
+            for (int row = 0; row < count; ++row) {
+                double sum = 0.0;
+                for (int column = 0; column < count; ++column)
+                    sum += std::abs(m_lu(row, column));
+                norm = std::max(norm, sum);
+            }
+            m_statistics.minimumScaledPivot = std::numeric_limits<double>::infinity();
+            for (int column = 0; column < count; ++column) {
+                int pivot = column;
+                for (int row = column + 1; row < count; ++row)
+                    if (std::abs(m_lu(row, column)) > std::abs(m_lu(pivot, column)))
+                        pivot = row;
+                m_pivots[column] = pivot;
+                const double magnitude = std::abs(m_lu(pivot, column));
+                m_statistics.minimumScaledPivot = std::min(m_statistics.minimumScaledPivot, magnitude);
+                if (!std::isfinite(magnitude) || magnitude == 0.0)
+                    return;
+                if (pivot != column)
+                    for (int j = 0; j < count; ++j)
+                        std::swap(m_lu(column, j), m_lu(pivot, j));
+                for (int row = column + 1; row < count; ++row) {
+                    m_lu(row, column) /= m_lu(column, column);
+                    for (int j = column + 1; j < count; ++j)
+                        m_lu(row, j) -= m_lu(row, column) * m_lu(column, j);
+                }
+            }
+            // A small explicit inverse supplies an infinity-norm condition
+            // estimate without regularizing a physically important coarse mode.
+            mfem::Vector column(count), sums(count);
+            sums = 0.0;
+            for (int j = 0; j < count; ++j) {
+                column = 0.0;
+                column(j) = m_scale;
+                SolveFactored(column);
+                for (int i = 0; i < count; ++i)
+                    sums(i) += std::abs(column(i));
+            }
+            const double inverseNorm = sums.Max();
+            m_statistics.reciprocalCondition = 1.0 / (norm * inverseNorm);
+            m_statistics.status = std::isfinite(m_statistics.reciprocalCondition) &&
+                                  m_statistics.reciprocalCondition >= m_options.minimumReciprocalCondition
+                ? SurfaceBorderSchurStatus::ready : SurfaceBorderSchurStatus::ill_conditioned;
+        }
+
+        void SolveFactored(mfem::Vector &values) const {
+            const int count = values.Size();
+            values /= m_scale;
+            for (int i = 0; i < count; ++i) {
+                if (m_pivots[i] != i)
+                    std::swap(values(i), values(m_pivots[i]));
+            }
+            for (int i = 0; i < count; ++i)
+                for (int j = 0; j < i; ++j)
+                    values(i) -= m_lu(i, j) * values(j);
+            for (int i = count - 1; i >= 0; --i) {
+                for (int j = i + 1; j < count; ++j)
+                    values(i) -= m_lu(i, j) * values(j);
+                values(i) /= m_lu(i, i);
+            }
+        }
+
+        const Operation *m_operation;
+        const Inverse *m_inverse;
+        SurfaceBorderSchurCoordinates m_coordinates;
+        MPI_Comm m_communicator;
+        SurfaceBorderSchurOptions m_options;
+        std::uint64_t m_generation{0};
+        mfem::DenseMatrix m_schur, m_coarse, m_lu;
+        std::vector<mfem::Vector> m_columns;
+        std::vector<int> m_pivots;
+        double m_scale{0.0};
+        mutable mfem::Vector m_rhs, m_fieldAction, m_image, m_solveWork;
+        mutable SurfaceBorderSchurStatistics m_statistics;
+    };
+} // namespace experiment::plateau_diagnostics

tests/solver/plateau_surface_border_schur.cpp

@@ -0,0 +1,336 @@
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <initializer_list>
+#include <limits>
+#include <type_traits>
+#include <vector>
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <mfem.hpp>
+#include <mpi.h>
+
+import mean_field;
+
+#include "../../experiments/plateau_surface_border_schur.hpp"
+
+namespace surface_border_schur_test {
+    namespace probe = experiment::plateau_diagnostics;
+    using namespace mean_field;
+
+    mfem::DenseMatrix Matrix(int size, std::initializer_list<double> entries) {
+        mfem::DenseMatrix result(size);
+        auto entry = entries.begin();
+        for (int row = 0; row < size; ++row)
+            for (int column = 0; column < size; ++column)
+                result(row, column) = *entry++;
+        return result;
+    }
+
+    mfem::Vector Vector(std::initializer_list<double> values) {
+        mfem::Vector result(values.size());
+        int i = 0;
+        for (double value : values)
+            result(i++) = value;
+        return result;
+    }
+
+    double Difference(const mfem::Vector &left, const mfem::Vector &right) {
+        mfem::Vector difference(left);
+        difference -= right;
+        return difference.Norml2();
+    }
+
+    class Operation final : public mfem::Operator {
+    public:
+        explicit Operation(mfem::DenseMatrix matrix)
+            : mfem::Operator(matrix.Height()), matrix(std::move(matrix)) {}
+        void Mult(const mfem::Vector &input, mfem::Vector &output) const override {
+            ++applications;
+            matrix.Mult(input, output);
+        }
+        bool IsPrepared() const { return prepared; }
+        const Operation &GetProblem() const { return *this; }
+        std::uint64_t GetPreparationGeneration() const { return generation; }
+        mfem::DenseMatrix matrix;
+        bool prepared{true};
+        std::uint64_t generation{1};
+        mutable std::uint64_t applications{0};
+    };
+
+    class Inverse : public mfem::Solver {
+    public:
+        static constexpr auto applicationContract = preconditioning::ApplicationContract::stationary_linear;
+        Inverse(const Operation &operation, mfem::DenseMatrix matrix)
+            : mfem::Solver(matrix.Height()), operation(operation), matrix(std::move(matrix)),
+              generation(operation.generation) {}
+        void SetOperator(const mfem::Operator &) override {}
+        void Mult(const mfem::Vector &input, mfem::Vector &output) const override {
+            ++applications;
+            matrix.Mult(input, output);
+            if (nonfinite)
+                output(0) = std::numeric_limits<double>::infinity();
+        }
+        bool IsCurrent() const { return generation == operation.generation; }
+        const Operation &GetNormalizedJacobian() const { return operation; }
+        const Operation &operation;
+        mfem::DenseMatrix matrix;
+        std::uint64_t generation;
+        bool nonfinite{false};
+        mutable std::uint64_t applications{0};
+    };
+
+    class FlexibleInverse : public Inverse {
+    public:
+        using Inverse::Inverse;
+        static constexpr auto applicationContract = preconditioning::ApplicationContract::flexible;
+    };
+    static_assert(probe::SurfaceBorderSchurOperation<Operation>);
+    static_assert(probe::SurfaceBorderSchurInverse<Inverse>);
+    static_assert(!probe::SurfaceBorderSchurInverse<FlexibleInverse>);
+    static_assert(!probe::SurfaceBorderSchurOperation<mfem::Operator>);
+
+    using Schur = probe::PreparedSurfaceBorderSchur<Operation, Inverse>;
+    static_assert(probe::SurfaceBorderSchurInverse<Schur>);
+    static_assert(!std::is_copy_constructible_v<Schur>);
+
+    // Fields occupy different equation/unknown positions. A contiguous or
+    // same-offset assumption would silently factor the wrong block here.
+    probe::SurfaceBorderSchurCoordinates Coordinates() { return {{1, 3}, {0, 2}}; }
+
+    mfem::DenseMatrix CoupledMatrix() {
+        return Matrix(4, {0.3, 4, -1, 1,
+                          0, 2,  2, 0.4,
+                        0.5, 2, 0.7, 3,
+                          3, -0.2, 1, 0.9});
+    }
+
+    mfem::DenseMatrix FieldInverse(bool exact) {
+        mfem::DenseMatrix result(4);
+        result = 0.0;
+        result(1, 0) = exact ? 0.3 : 0.2;
+        result(1, 2) = exact ? -0.1 : 0.0;
+        result(3, 0) = exact ? -0.2 : 0.0;
+        result(3, 2) = exact ? 0.4 : 0.2;
+        // These coarse contributions must be excluded, not accidentally used
+        // to complete the field response twice.
+        result(0, 0) = 11;
+        result(2, 2) = -17;
+        result(1, 1) = 23;
+        result(3, 3) = -31;
+        return result;
+    }
+
+    namespace blocks = utils::blocks;
+    struct StructureDescriptor {
+        using CorrectionBlocks = blocks::type_list<blocks::density::mass::value,
+            blocks::surface_deformation::parameters::value, blocks::enthalpy::specific::value>;
+        using ResidualBlocks = blocks::type_list<blocks::enthalpy::specific::residual,
+            blocks::surface_deformation::shape_equilibrium::residual, blocks::density::mass::residual>;
+    };
+    struct RootDescriptor {
+        using StructureComponent = StructureDescriptor;
+        using CorrectionBlocks = blocks::type_list<blocks::density::mass::value,
+            blocks::surface_deformation::parameters::value, blocks::enthalpy::specific::value,
+            blocks::fixed_total_mass::mass_normalization::value,
+            blocks::fixed_central_density::central_value::value>;
+        using ResidualBlocks = blocks::type_list<blocks::enthalpy::specific::residual,
+            blocks::surface_deformation::shape_equilibrium::residual, blocks::density::mass::residual,
+            blocks::fixed_total_mass::mass_normalization::residual,
+            blocks::fixed_central_density::central_value::residual>;
+        static constexpr std::size_t borderValueArity = 2, borderResidualArity = 2;
+    };
+    struct CoordinateMap {
+        using Range = preconditioning::EquilibriumCoordinateRange;
+        std::array<Range, 5> columns{{{5, 0, 2}, {1, 2, 3}, {0, 5, 1}, {4, 6, 1}, {7, 7, 1}}};
+        std::array<Range, 5> rows{{{6, 0, 1}, {2, 1, 3}, {0, 4, 2}, {5, 6, 1}, {7, 7, 1}}};
+        const auto &GetCorrectionRanges() const { return columns; }
+        const auto &GetResidualRanges() const { return rows; }
+    };
+    struct CoordinateBinding {
+        RootDescriptor root;
+        CoordinateMap coordinates;
+        const auto &GetBlock() const { return root; }
+        const auto &GetCoordinateMap() const { return coordinates; }
+    };
+}
+
+TEST_CASE("Surface Border Schur Coordinates Include All Native Modes And Generated Central Border",
+          "[solver][surface_border_schur][serial]") {
+    using namespace surface_border_schur_test;
+    const auto coordinates = probe::MakeSurfaceBorderSchurCoordinates(CoordinateBinding{});
+    CHECK(coordinates.correctionIndices == std::vector<int>{1, 2, 3, 4, 7});
+    CHECK(coordinates.residualIndices == std::vector<int>{2, 3, 4, 5, 7});
+}
+
+TEST_CASE("Surface Border Schur With Exact Field Inverse Solves Full Nonsymmetric System",
+          "[solver][surface_border_schur][serial]") {
+    using namespace surface_border_schur_test;
+    Operation operation(CoupledMatrix());
+    Inverse inverse(operation, FieldInverse(true));
+    Schur schur(operation, inverse, Coordinates(), MPI_COMM_WORLD);
+    REQUIRE(schur.IsCurrent());
+    CHECK(schur.GetStatistics().setupOperatorApplications == 4);
+    CHECK(schur.GetStatistics().setupInverseApplications == 2);
+    const auto rhs = Vector({1, -2, 3, 4});
+    mfem::Vector action(4), expected(4), defect(4);
+    mfem::DenseMatrixInverse direct(operation.matrix);
+    direct.Mult(rhs, expected);
+    schur.Mult(rhs, action);
+    CHECK(Difference(action, expected) < 1.e-13);
+    operation.Mult(action, defect);
+    CHECK(Difference(defect, rhs) < 1.e-13);
+    CHECK(schur.GetStatistics().operatorApplications == 1);
+    CHECK(schur.GetStatistics().inverseApplications == 1);
+    for (int column = 0; column < 2; ++column) {
+        const auto &completed = schur.GetCompletedColumns()[column];
+        operation.Mult(completed, defect);
+        CHECK(std::abs(defect(0)) < 1.e-13);
+        CHECK(std::abs(defect(2)) < 1.e-13);
+        CHECK(completed(Coordinates().correctionIndices[column]) == 1.0);
+    }
+}
+
+TEST_CASE("Surface Border Schur With Approximate Fields Enforces Every Coarse Equation",
+          "[solver][surface_border_schur][serial]") {
+    using namespace surface_border_schur_test;
+    Operation operation(CoupledMatrix());
+    Inverse inverse(operation, FieldInverse(false));
+    Schur schur(operation, inverse, Coordinates(), MPI_COMM_WORLD);
+    REQUIRE(schur.IsCurrent());
+    const auto rhs = Vector({1, -2, 3, 4});
+    mfem::Vector action(4), image(4);
+    schur.Mult(rhs, action);
+    operation.Mult(action, image);
+    for (int index : Coordinates().residualIndices)
+        CHECK(std::abs(image(index) - rhs(index)) < 1.e-13);
+    CHECK(Difference(image, rhs) > 1.e-3);
+    // Verify S = D - A_CF B A_FC independently in this explicit packing.
+    const auto &s = schur.GetSchurMatrix();
+    CHECK(s(0, 0) == Catch::Approx(-0.16));
+    CHECK(s(0, 1) == Catch::Approx(2.344));
+    CHECK(s(1, 0) == Catch::Approx(2.922));
+    CHECK(s(1, 1) == Catch::Approx(0.834));
+}
+
+TEST_CASE("Surface Border Schur Is Stationary Linear And Supports Aliasing",
+          "[solver][surface_border_schur][serial]") {
+    using namespace surface_border_schur_test;
+    Operation operation(CoupledMatrix());
+    Inverse inverse(operation, FieldInverse(false));
+    Schur schur(operation, inverse, Coordinates(), MPI_COMM_WORLD);
+    const auto first = Vector({1, 2, -1, 3}), second = Vector({-2, 4, 0.5, 1});
+    mfem::Vector a(4), b(4), combined(first), expected(4), observed(4);
+    schur.Mult(first, a);
+    schur.Mult(second, b);
+    combined *= 2.3;
+    combined.Add(-0.7, second);
+    expected = a;
+    expected *= 2.3;
+    expected.Add(-0.7, b);
+    schur.Mult(combined, observed);
+    CHECK(Difference(observed, expected) < 1.e-13);
+    schur.Mult(combined, combined);
+    CHECK(Difference(combined, expected) < 1.e-13);
+    schur.Mult(first, observed);
+    CHECK(Difference(observed, a) == 0.0);
+}
+
+TEST_CASE("Surface Border Schur Reports Singular And Ill Conditioned Coarse Blocks Without Regularization",
+          "[solver][surface_border_schur][serial]") {
+    using namespace surface_border_schur_test;
+    Operation operation(Matrix(3, {1, 0, 0, 0, 0, 0, 0, 0, 1}));
+    Inverse inverse(operation, Matrix(3, {1, 0, 0, 0, 0, 0, 0, 0, 0}));
+    Schur singular(operation, inverse, {{1, 2}, {1, 2}}, MPI_COMM_WORLD);
+    CHECK(singular.GetStatistics().status == probe::SurfaceBorderSchurStatus::singular);
+    CHECK_FALSE(singular.IsCurrent());
+    auto rhs = Vector({1, 2, 3});
+    mfem::Vector action(3);
+    CHECK_THROWS_AS(singular.Mult(rhs, action), std::logic_error);
+    operation.matrix(1, 1) = 1.e-16;
+    ++operation.generation;
+    inverse.generation = operation.generation;
+    singular.Refresh();
+    CHECK(singular.GetStatistics().status == probe::SurfaceBorderSchurStatus::ill_conditioned);
+    CHECK(singular.GetStatistics().reciprocalCondition == Catch::Approx(1.e-16));
+    CHECK_FALSE(singular.IsCurrent());
+    CHECK_THROWS_AS(singular.SetOperator(operation), std::logic_error);
+    CHECK(singular.GetSchurMatrix()(0, 0) == 1.e-16);
+}
+
+TEST_CASE("Surface Border Schur Handles Pivoting And Uniform Matrix Scale",
+          "[solver][surface_border_schur][serial]") {
+    using namespace surface_border_schur_test;
+    for (double scale : {1.e-100, 1.0, 1.e100}) {
+        auto matrix = Matrix(3, {2, 0, 0, 0, 0, 2, 0, 3, 1});
+        matrix *= scale;
+        Operation operation(std::move(matrix));
+        auto base = Matrix(3, {0.5, 0, 0, 0, 0, 0, 0, 0, 0});
+        base *= 1.0 / scale;
+        Inverse inverse(operation, std::move(base));
+        Schur schur(operation, inverse, {{1, 2}, {1, 2}}, MPI_COMM_WORLD);
+        REQUIRE(schur.IsCurrent());
+        auto rhs = Vector({1, 2, 4});
+        rhs *= scale;
+        mfem::Vector action(3);
+        schur.Mult(rhs, action);
+        CHECK(Difference(action, Vector({0.5, 1, 1})) < 1.e-13);
+        CHECK(schur.GetStatistics().reciprocalCondition == Catch::Approx(0.5));
+    }
+}
+
+TEST_CASE("Surface Border Schur Requires Unbordered Fields For A Zero Scalar Diagonal",
+          "[solver][surface_border_schur][serial]") {
+    using namespace surface_border_schur_test;
+    Operation operation(Matrix(2, {1, 1, 1, 0}));
+    Inverse alreadyBordered(operation, Matrix(2, {0, 1, 1, -1}));
+    Schur wrongRestriction(operation, alreadyBordered, {{1}, {1}}, MPI_COMM_WORLD);
+    CHECK(wrongRestriction.GetStatistics().status == probe::SurfaceBorderSchurStatus::singular);
+    Inverse unbordered(operation, Matrix(2, {1, 0, 0, 0}));
+    Schur schur(operation, unbordered, {{1}, {1}}, MPI_COMM_WORLD);
+    REQUIRE(schur.IsCurrent());
+    auto rhs = Vector({3, 2});
+    mfem::Vector action(2);
+    schur.Mult(rhs, action);
+    CHECK(Difference(action, Vector({2, 1})) < 1.e-13);
+}
+
+TEST_CASE("Surface Border Schur Rejects Invalid Partitions And Stale Or Nonfinite Dependencies",
+          "[solver][surface_border_schur][serial]") {
+    using namespace surface_border_schur_test;
+    Operation operation(CoupledMatrix());
+    Inverse inverse(operation, FieldInverse(true));
+    const auto reject = [&](probe::SurfaceBorderSchurCoordinates coordinates) {
+        CHECK_THROWS_AS(Schur(operation, inverse, std::move(coordinates), MPI_COMM_WORLD), std::invalid_argument);
+    };
+    reject({{}, {}});
+    reject({{1}, {0, 2}});
+    reject({{1, 1}, {0, 2}});
+    reject({{1, 3}, {0, 0}});
+    reject({{1, 4}, {0, 2}});
+    reject({{0, 1, 2, 3}, {0, 1, 2, 3}});
+    CHECK_THROWS_AS(Schur(operation, inverse, Coordinates(), MPI_COMM_NULL), std::invalid_argument);
+    Schur schur(operation, inverse, Coordinates(), MPI_COMM_WORLD);
+    auto rhs = Vector({1, 2, 3, 4});
+    mfem::Vector action(4);
+    Operation other(CoupledMatrix());
+    CHECK_THROWS_AS(schur.SetOperator(other), std::invalid_argument);
+    CHECK_THROWS_AS(Schur(other, inverse, Coordinates(), MPI_COMM_WORLD), std::invalid_argument);
+    ++operation.generation;
+    CHECK_FALSE(schur.IsCurrent());
+    CHECK_THROWS_AS(schur.Mult(rhs, action), std::logic_error);
+    inverse.generation = operation.generation;
+    CHECK_FALSE(schur.IsCurrent());
+    schur.Refresh();
+    REQUIRE(schur.IsCurrent());
+    inverse.iterative_mode = true;
+    CHECK_FALSE(schur.IsCurrent());
+    CHECK_THROWS_AS(schur.Mult(rhs, action), std::logic_error);
+    inverse.iterative_mode = false;
+    inverse.nonfinite = true;
+    CHECK_THROWS_AS(schur.Mult(rhs, action), std::runtime_error);
+    inverse.nonfinite = false;
+    rhs(0) = std::numeric_limits<double>::quiet_NaN();
+    CHECK_THROWS_AS(schur.Mult(rhs, action), std::runtime_error);
+}