diff --git a/src/include/serif/solver/preconditioner/contract.hpp b/src/include/serif/solver/preconditioner/contract.hpp
new file mode 100644
index 00000000..77b5711d
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/contract.hpp
@@ -0,0 +1,11 @@
+#pragma once
+
+namespace serif::solver::preconditioner
+{
+    /** @brief Whether repeated applications represent one fixed linear map or a variable inner solve. */
+    enum class ApplicationContract
+    {
+        fixed_linear,
+        variable
+    };
+} // namespace serif::solver::preconditioner
diff --git a/src/include/serif/solver/preconditioner/default_preconditioners.hpp b/src/include/serif/solver/preconditioner/default_preconditioners.hpp
new file mode 100644
index 00000000..67607537
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/default_preconditioners.hpp
@@ -0,0 +1,75 @@
+#pragma once
+
+#include <concepts>
+#include <type_traits>
+
+#include "serif/discretization/blocks/concepts.hpp"
+#include "serif/discretization/blocks/fields.hpp"
+#include "serif/discretization/blocks/system_layout.hpp"
+#include "serif/solver/preconditioner/recipes/bordered.hpp"
+#include "serif/solver/preconditioner/recipes/block.hpp"
+#include "serif/solver/preconditioner/strategies/inverse/diagonal.hpp"
+#include "serif/solver/preconditioner/strategies/inverse/fixed_amg.hpp"
+#include "serif/solver/preconditioner/strategies/inverse/sparse_direct.hpp"
+
+namespace serif::solver::preconditioner {
+    namespace detail {
+        using PrimalFields = utils::dependencies::KeyList<discretization::blocks::density::mass, discretization::blocks::gravity::potential>;
+
+        /** @brief Separate the primal fields from scalar constraints and optional shape coordinates. */
+        template <typename Coordinates, typename Shape>
+        struct PrimalCoordinates;
+
+        template <typename... Coordinates, typename Shape>
+        struct PrimalCoordinates<utils::dependencies::KeyList<Coordinates...>, Shape> {
+            using Fields = PrimalFields;
+            using Scalars = utils::dependencies::CombinedKeyList<std::conditional_t<discretization::blocks::GlobalScalarTerm<Coordinates>,
+                utils::dependencies::KeyList<Coordinates>, utils::dependencies::KeyList<>>...>;
+            using Border = utils::dependencies::CombinedKeyList<std::conditional_t<
+                discretization::blocks::GlobalScalarTerm<Coordinates> || std::same_as<Coordinates, Shape>,
+                utils::dependencies::KeyList<Coordinates>, utils::dependencies::KeyList<>>...>;
+
+            static constexpr bool supported =
+                utils::dependencies::key_count_v<utils::dependencies::CombinedKeyList<Fields, utils::dependencies::KeyList<Coordinates...>>> == sizeof...(Coordinates) &&
+                ((utils::dependencies::ContainsKey<Coordinates, Fields> || discretization::blocks::GlobalScalarTerm<Coordinates> || std::same_as<Coordinates, Shape>) && ...);
+            static constexpr bool has_shape = (std::same_as<Coordinates, Shape> || ...);
+        };
+
+        template <typename Layout>
+        using PrimalValues = PrimalCoordinates<typename std::remove_cvref_t<Layout>::Values, discretization::blocks::surface_deformation::coefficients>;
+
+        template <typename Layout>
+        using PrimalRows = PrimalCoordinates<typename std::remove_cvref_t<Layout>::Residuals, discretization::blocks::surface_deformation::shape_equilibrium>;
+    }
+
+    /** @brief Require density and potential fields with paired optional shape and scalar constraints. */
+    template <typename Layout>
+    concept PrimalPreconditionerLayout = discretization::blocks::SystemLayoutType<Layout> &&
+        detail::PrimalValues<Layout>::supported && detail::PrimalRows<Layout>::supported &&
+        detail::PrimalValues<Layout>::has_shape == detail::PrimalRows<Layout>::has_shape &&
+        utils::dependencies::key_count_v<typename detail::PrimalValues<Layout>::Scalars> ==
+            utils::dependencies::key_count_v<typename detail::PrimalRows<Layout>::Scalars>;
+
+    /**
+     * @brief Couple scalar density and potential responses, then correct every shape and global constraint together.
+     *
+     * Potential responds first and its force contribution then adjusts density. The existing dense border retains
+     * surface-to-field and field-to-surface responses, including the mean radius mode; shape has no separate mass inverse.
+     * Density uses its diagonal without regularization because its pressure-force Jacobian includes drift and need not
+     * satisfy the positive-diagonal contract of the elliptic AMG strategy used for potential.
+     * @note The dense border stores one field response per shape/global coordinate and a square border matrix.
+     * This default targets bounded surface catalogs; large nodal surfaces require an explicitly chosen scalable recipe.
+     * @note This is an approximate inverse, not a guarantee that the physical constraint Schur system is nonsingular.
+     * @param potential_inverse Inverse of the potential block: three fixed AMG cycles by default, or for example
+     * strategies::inverse::SparseDirect for an exact factorization of moderate-size blocks.
+     */
+    template <PrimalPreconditionerLayout Layout, typename PotentialInverse = strategies::inverse::FixedAMG>
+    [[nodiscard]] auto make_primal_preconditioner(const Layout&, PotentialInverse potential_inverse = strategies::inverse::FixedAMG{.cycles = 3}) {
+        namespace blocks = discretization::blocks;
+        using Constraints = ConstraintCoordinates<typename detail::PrimalRows<Layout>::Border, typename detail::PrimalValues<Layout>::Border>;
+        const auto density = make_preconditioner<blocks::density::mass>(strategies::inverse::Diagonal{});
+        const auto potential = make_preconditioner<blocks::gravity::potential>(potential_inverse);
+        const auto fields = couple(density, potential, strategies::coupling::SecondThenFirst{});
+        return with_constraints<Constraints>(fields);
+    }
+}
diff --git a/src/include/serif/solver/preconditioner/failure.hpp b/src/include/serif/solver/preconditioner/failure.hpp
new file mode 100644
index 00000000..ac7e43f1
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/failure.hpp
@@ -0,0 +1,14 @@
+#pragma once
+
+#include <string>
+#include "serif/utils/error/codes.hpp"
+
+namespace serif::solver::preconditioner
+{
+    /** @brief A recoverable preparation failure, including the numerical reason and how to correct it. */
+    struct Failure final
+    {
+        utils::errors::SERiFErrorCode code{utils::errors::SERiFErrorCode::PRECONDITIONER_ERROR};
+        std::string message;
+    };
+} // namespace serif::solver::preconditioner
diff --git a/src/include/serif/solver/preconditioner/keys.hpp b/src/include/serif/solver/preconditioner/keys.hpp
new file mode 100644
index 00000000..3534e904
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/keys.hpp
@@ -0,0 +1,25 @@
+#pragma once
+
+#include "serif/operators/linearized_system_view.hpp"
+
+namespace serif::solver::preconditioner::keys
+{
+    template <typename Rows, typename Values>
+    using Block = operators::PhysicalBlock<Rows, Values>;
+
+    /** @brief The inverse maps the selected equation residuals back to their physical unknowns. */
+    template <typename Rows, typename Values>
+    struct Inverse
+    {
+    };
+
+    template <typename Rows, typename Values>
+    struct Diagonal
+    {
+    };
+
+    template <typename Rows, typename Values>
+    struct Matrix
+    {
+    };
+} // namespace serif::solver::preconditioner::keys
diff --git a/src/include/serif/solver/preconditioner/preconditioner.hpp b/src/include/serif/solver/preconditioner/preconditioner.hpp
new file mode 100644
index 00000000..a56939de
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/preconditioner.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "serif/backend/mfem/solver/preconditioner/graph.hpp"
+#include "serif/backend/mfem/solver/preconditioner/solver.hpp"
+#include "serif/solver/preconditioner/default_preconditioners.hpp"
+
+namespace serif::solver::preconditioner {
+    using backend::mfem::Solver;
+}
diff --git a/src/include/serif/solver/preconditioner/recipes/block.hpp b/src/include/serif/solver/preconditioner/recipes/block.hpp
new file mode 100644
index 00000000..ab205f3f
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/recipes/block.hpp
@@ -0,0 +1,25 @@
+#pragma once
+
+#include <utility>
+#include "serif/solver/preconditioner/keys.hpp"
+
+namespace serif::solver::preconditioner
+{
+    /** @brief Select a physical equation/unknown pair and the numerical approximation to its inverse. */
+    template <typename Row, typename Field, typename Strategy>
+    struct BlockPreconditioner final
+    {
+        using Rows = utils::dependencies::KeyList<Row>;
+        using Values = utils::dependencies::KeyList<Field>;
+        using Result = keys::Inverse<Rows, Values>;
+        Strategy inverse;
+    };
+
+    /** @brief Describe an inverse without preparing matrices or binding runtime state. */
+    template <typename Field, typename Row = Field, typename Strategy>
+    [[nodiscard]] auto make_preconditioner(Strategy inverse)
+    {
+        return BlockPreconditioner<Row, Field, Strategy>{.inverse = std::move(inverse)};
+    }
+
+} // namespace serif::solver::preconditioner
diff --git a/src/include/serif/solver/preconditioner/recipes/bordered.hpp b/src/include/serif/solver/preconditioner/recipes/bordered.hpp
new file mode 100644
index 00000000..ca94acfc
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/recipes/bordered.hpp
@@ -0,0 +1,43 @@
+#pragma once
+
+#include <utility>
+#include "serif/solver/preconditioner/recipes/coupled.hpp"
+#include "serif/solver/preconditioner/strategies/inverse/dense.hpp"
+
+namespace serif::solver::preconditioner
+{
+    /** @brief Select arbitrary constraint rows and unknowns, including coupled generated scalars. */
+    template <utils::dependencies::IsKeyList RowsT, utils::dependencies::IsKeyList ValuesT>
+    struct ConstraintCoordinates final
+    {
+        using Rows = RowsT;
+        using Values = ValuesT;
+    };
+
+    /** @brief Correct a physical field inverse through its dense constraint Schur system. */
+    template <typename Fields, typename Constraints>
+        requires DisjointSubsystems<Fields, Constraints>
+    struct BorderedPreconditioner final
+    {
+        using Rows = utils::dependencies::CombinedKeyList<typename Fields::Rows, typename Constraints::Rows>;
+        using Values = utils::dependencies::CombinedKeyList<typename Fields::Values, typename Constraints::Values>;
+        using Result = keys::Inverse<Rows, Values>;
+        Fields fields;
+        strategies::inverse::DenseDirect constraint_inverse;
+    };
+
+    /** @brief Attach constraints without a registry of built-in mass, rotation, or central-density providers. */
+    template <typename Constraints, typename Fields>
+        requires DisjointSubsystems<Fields, Constraints>
+    [[nodiscard]] auto with_constraints(Fields fields, strategies::inverse::DenseDirect inverse = {})
+    {
+        if constexpr (utils::dependencies::key_count_v<typename Constraints::Rows> == 0 && utils::dependencies::key_count_v<typename Constraints::Values> == 0)
+        {
+            return fields;
+        }
+        else
+        {
+            return BorderedPreconditioner<Fields, Constraints>{.fields = std::move(fields), .constraint_inverse = inverse};
+        }
+    }
+} // namespace serif::solver::preconditioner
diff --git a/src/include/serif/solver/preconditioner/recipes/coupled.hpp b/src/include/serif/solver/preconditioner/recipes/coupled.hpp
new file mode 100644
index 00000000..b7bacc89
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/recipes/coupled.hpp
@@ -0,0 +1,37 @@
+#pragma once
+
+#include <utility>
+#include "serif/solver/preconditioner/keys.hpp"
+#include "serif/solver/preconditioner/strategies/coupling/policies.hpp"
+
+namespace serif::solver::preconditioner
+{
+    /** @brief Require disjoint ownership on both sides rather than silently deduplicating physical coordinates. */
+    template <typename First, typename Second>
+    concept DisjointSubsystems =
+        utils::dependencies::key_count_v<utils::dependencies::CombinedKeyList<typename First::Rows, typename Second::Rows>> ==
+            utils::dependencies::key_count_v<typename First::Rows> + utils::dependencies::key_count_v<typename Second::Rows> &&
+        utils::dependencies::key_count_v<utils::dependencies::CombinedKeyList<typename First::Values, typename Second::Values>> ==
+            utils::dependencies::key_count_v<typename First::Values> + utils::dependencies::key_count_v<typename Second::Values>;
+
+    /** @brief Compose physical subsystem descriptions without creating an additional preparation owner. */
+    template <typename First, typename Second, typename Policy>
+        requires DisjointSubsystems<First, Second>
+    struct CoupledPreconditioner final
+    {
+        using Rows = utils::dependencies::CombinedKeyList<typename First::Rows, typename Second::Rows>;
+        using Values = utils::dependencies::CombinedKeyList<typename First::Values, typename Second::Values>;
+        using Result = keys::Inverse<Rows, Values>;
+        First first;
+        Second second;
+        Policy coupling;
+    };
+
+    /** @brief Choose the physical subsystem ordering explicitly; the same algebra serves material/surface and material/gravity. */
+    template <typename First, typename Second, typename Policy>
+        requires DisjointSubsystems<First, Second>
+    [[nodiscard]] auto couple(First first, Second second, Policy coupling)
+    {
+        return CoupledPreconditioner<First, Second, Policy>{.first = std::move(first), .second = std::move(second), .coupling = coupling};
+    }
+} // namespace serif::solver::preconditioner
diff --git a/src/include/serif/solver/preconditioner/strategies/coupling/order.hpp b/src/include/serif/solver/preconditioner/strategies/coupling/order.hpp
new file mode 100644
index 00000000..c14396af
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/strategies/coupling/order.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+namespace serif::solver::preconditioner::strategies::coupling
+{
+    /** @brief Choose how the two groups of physical unknowns respond to each other's corrections. */
+    enum class CouplingOrder
+    {
+        independent,
+        second_then_first,
+        lower_diagonal_upper
+    };
+} // namespace serif::solver::preconditioner::strategies::coupling
diff --git a/src/include/serif/solver/preconditioner/strategies/coupling/policies.hpp b/src/include/serif/solver/preconditioner/strategies/coupling/policies.hpp
new file mode 100644
index 00000000..a7787662
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/strategies/coupling/policies.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+#include "serif/solver/preconditioner/strategies/coupling/order.hpp"
+
+namespace serif::solver::preconditioner::strategies::coupling
+{
+    /** @brief Correct the first subsystem using the second subsystem's update. */
+    struct SecondThenFirst final
+    {
+        static constexpr CouplingOrder order = CouplingOrder::second_then_first;
+    };
+} // namespace serif::solver::preconditioner::strategies::coupling
diff --git a/src/include/serif/solver/preconditioner/strategies/inverse/dense.hpp b/src/include/serif/solver/preconditioner/strategies/inverse/dense.hpp
new file mode 100644
index 00000000..9158e527
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/strategies/inverse/dense.hpp
@@ -0,0 +1,13 @@
+#pragma once
+
+#include "serif/solver/preconditioner/contract.hpp"
+
+namespace serif::solver::preconditioner::strategies::inverse
+{
+    /** @brief Directly solve a small dense system; optionally reject an unreliable infinity-norm condition estimate. */
+    struct DenseDirect final
+    {
+        static constexpr ApplicationContract application_contract{ApplicationContract::fixed_linear};
+        double minimum_reciprocal_condition{0.0};
+    };
+} // namespace serif::solver::preconditioner::strategies::inverse
diff --git a/src/include/serif/solver/preconditioner/strategies/inverse/diagonal.hpp b/src/include/serif/solver/preconditioner/strategies/inverse/diagonal.hpp
new file mode 100644
index 00000000..0e8b2298
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/strategies/inverse/diagonal.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+#include "serif/solver/preconditioner/contract.hpp"
+
+namespace serif::solver::preconditioner::strategies::inverse
+{
+    /** @brief Invert each diagonal entry independently. */
+    struct Diagonal final
+    {
+        static constexpr ApplicationContract application_contract{ApplicationContract::fixed_linear};
+    };
+} // namespace serif::solver::preconditioner::strategies::inverse
diff --git a/src/include/serif/solver/preconditioner/strategies/inverse/fixed_amg.hpp b/src/include/serif/solver/preconditioner/strategies/inverse/fixed_amg.hpp
new file mode 100644
index 00000000..8e89dedb
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/strategies/inverse/fixed_amg.hpp
@@ -0,0 +1,13 @@
+#pragma once
+
+#include "serif/solver/preconditioner/contract.hpp"
+
+namespace serif::solver::preconditioner::strategies::inverse
+{
+    /** @brief Apply a fixed number of BoomerAMG cycles from zero, preserving a fixed linear inverse. */
+    struct FixedAMG final
+    {
+        static constexpr ApplicationContract application_contract{ApplicationContract::fixed_linear};
+        int cycles{1};
+    };
+} // namespace serif::solver::preconditioner::strategies::inverse
diff --git a/src/include/serif/solver/preconditioner/strategies/inverse/sparse_direct.hpp b/src/include/serif/solver/preconditioner/strategies/inverse/sparse_direct.hpp
new file mode 100644
index 00000000..06764446
--- /dev/null
+++ b/src/include/serif/solver/preconditioner/strategies/inverse/sparse_direct.hpp
@@ -0,0 +1,21 @@
+#pragma once
+
+#include "serif/solver/preconditioner/contract.hpp"
+
+namespace serif::solver::preconditioner::strategies::inverse
+{
+    /**
+     * @brief Exact inverse of the assembled block through a sparse LU factorization, prepared once per state.
+     * @note The backend is chosen at build time with the sparse_direct Meson option (Apple Accelerate on macOS 15.5+,
+     * SuiteSparse UMFPACK elsewhere, or none); without one, preparation fails rather than falling back to another
+     * inverse. backend::mfem::sparse_direct_available() reports it at run time. Suited to blocks of tens of
+     * thousands of rows, where one factorization and its triangular solves are much cheaper than repeated
+     * multigrid cycles. The symbolic analysis is reused across Newton steps while the sparsity pattern is
+     * unchanged (bitwise the same factorization), and batched applications (ArrayMult) solve their right-hand
+     * sides together, which agrees with single solves to rounding.
+     */
+    struct SparseDirect final
+    {
+        static constexpr ApplicationContract application_contract{ApplicationContract::fixed_linear};
+    };
+} // namespace serif::solver::preconditioner::strategies::inverse