feat(solver/preconditioner): Implements preconditioner system

The system matrix for a problem may be poorly conditioned, preconditioners are supposed to efficiently estimate the inverse of the system matrix so that a simpler problem, with the approximate identity matrix may be solved instead. The problem is that the best approach to this varies with different fields. The preconditioner system implements infrastructure to define seperate preconditioners for different fields which are then composed together.

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

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);
+    }
+}

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

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

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;
+}

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

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

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

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

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

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

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

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

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