formoniq: one solver interface, with subspace correction as its general form

[luiswirth]

Motivation

Every solve in the crate names its own backend. FaerLu::new(...).solve(...)
appears literally in fe.rs, hodge.rs, whitney_complex.rs, time.rs,
problems/dirac.rs and problems/elliptic.rs. There is no seam to swap at, so
the multigrid and auxiliary-space machinery of #115 cannot be reached from a
problem however good it is: solve_source<C: HilbertComplex> is built from one
mesh, and a V-cycle needs a hierarchy. The datum is missing from the signature,
which is why the third bullet of #115 never landed.

Two further symptoms of the same gap. Grade0Multigrid and GradeKHodgeHx are
each a bundle of a tower, an operator and a solve method, split on a grade
branch; that is one object written twice. And the AFW block preconditioner in
elliptic.rs builds its blocks from DirectInverse, so the mixed Hodge-Laplace
solve has a mesh-independent iteration count but a non-optimal cost per
iteration.

What is actually solved

The PDE operators are Hodge-Laplace, Hodge-Dirac and the Lie derivative. The
algebraic systems they produce differ more than the operators do:

system	where	character
M_k	fe (L2 projection), hodge::codif, whitney_complex, Leapfrog, dirac	SPD, conditioning O(1) in h
A_k = M_k + D^T M_{k+1} D	elliptic blocks, multigrid, hx	SPD, O(h^-2)
mixed Hodge-Laplace KKT, bordered by the harmonic constraint	elliptic::solve_source	symmetric indefinite
Hodge-Dirac bordered system	dirac	symmetric indefinite
IRK stage system I ⊗ M - dt (A_tab ⊗ A)	time::LinearIrk, used by heat, wave, advection	nonsymmetric, block over stages
shift-invert pencil	linalg::eigen	symmetric indefinite

Only the second row admits multigrid directly. That is not narrow, because the
others are preconditioned by it: the KKT blocks are A_{k-1} and A_k; the
Dirac operator satisfies D^2 = -Δ, so the Laplace preconditioner is the Dirac
preconditioner; and an implicit stage block is M_k + dt a_ii A_k, a shifted
Hodge operator, which is the operator Hiptmair-Xu was written for in the first
place.

Design

Subspace correction is the general interface

Multigrid, Hiptmair-Xu and the AFW block preconditioner are one theorem, Xu's
subspace correction: decompose V = Σ Π_i V_i, solve on each subspace, combine.
They differ only in the decomposition and the combinator.

	subspaces	combinator
multigrid	the coarse levels of a tower	multiplicative (VCycle)
Hiptmair-Xu	d Λ^{k-1} and the vector-nodal space	additive (AuxiliarySpace)
AFW blocks	the σ, u, p factors	additive (BlockDiagonal)
domain decomposition (absent)	overlapping patches	additive

iterative already carries both combinators, and correctly. What it lacks is a
name for the input they share:

/// A decomposition of the space into subspaces, each with its prolongation
/// Π_i: V_i -> V and an approximate inverse on it. The single input of
/// subspace correction; VCycle combines the pieces multiplicatively,
/// AuxiliarySpace additively.
pub struct Decomposition { .. }

Every FEEC-specific question then reduces to one: given this operator, which
decomposition does it admit?

The extension point

// formoniq::solve

pub trait FeecOperator {
  fn assemble(&self, complex: &impl HilbertComplex) -> CsrMatrix;

  /// The decomposition this operator admits, `None` where it admits none
  /// and a direct solve is the answer.
  fn decompose(&self, structure: &Structure) -> Option<Decomposition> { None }
}

/// What is known about the space beyond its entries: a refinement tower, the
/// de Rham complex, later a domain partition.
pub struct Structure<'a> { .. }

The default decompose is what makes this cheap to extend: implementing the
trait costs one method and yields a working direct solve; the decomposition is
opted into once it is known. Advection can land with None and be correct.

An operator family rather than an operator menu carries the concrete cases:

/// a(u,v) = <α u, v> + <β du, dv> on grade k, with α and β material weights.
struct HodgeForm<A, B> { .. }

Mass is β = 0, the HΛ(d) gram is α = β = 1, an implicit stage block is
β = dt a_ii, and Maxwell with materials and Helmholtz are the same form with
field-valued weights. Scalar coefficients would not suffice: operators.rs
already carries WeightedHodgeMassElmat, and curl (1/μ) curl - ω² ε needs
tensor-valued material data.

Its decompose is one recursion: at grade 0 the tower's levels, at grade k
the regular decomposition whose gradient subspace is the grade-(k-1) form's own
decomposition. Base case grade 0, so there is no k = 1 special case.

The solver

The choice of solver is not a property of the space, so it does not belong on
HilbertComplex, which states facts about the space and nothing else. It is a
separate object carrying the space it acts on:

pub trait Solver {
  fn complex(&self) -> &impl HilbertComplex;
  fn inverse(&self, op: &impl FeecOperator) -> Option<HodgeInverse>;
}

Direct::new(&complex)                  // entries only: a factorization
Subspace::new(&tower).with_sweeps(2)   // reads the operator's decomposition

Subspace::new(&tower) fixes the complex to the tower's finest level, so the
space and the strategy cannot disagree. Problem entry points take the solver
instead of the complex:

elliptic::solve_source(&Direct::new(&complex), source, 1)?
elliptic::solve_source(&Subspace::new(&tower), source, 1)?

solve_source's body is otherwise untouched, and its AFW blocks become
solver.inverse(&HodgeForm::hdif(k-1)) and solver.inverse(&HodgeForm::hdif(k)),
which is what makes the mixed solve optimal end to end. None keeps its present
meaning throughout: not positive definite, the Lorentzian case, falling back to a
whole-system factorization.

Krylov: the symmetry class picks the method

cg requires an SPD operator, minres a symmetric one, and neither currently
says so. Add a PositiveDefinite marker beside the existing SelfAdjoint and
require them:

fn cg    (a: impl LinearOperator + SelfAdjoint + PositiveDefinite, ..)
fn minres(a: impl LinearOperator + SelfAdjoint, ..)
fn gmres (a: impl LinearOperator, restart: usize, ..)   // new

The wrong pairing then fails to compile. Definiteness that depends on the
geometry rather than the form stays runtime and keeps the Option contract
DirectInverse::try_new already has.

GMRES is the missing method and advection is what needs it: the Lie derivative
is nonsymmetric, so neither CG nor MINRES applies. Its symmetric part is a
HodgeForm and can supply the decomposition, without pretending the operator is
symmetric.

The scalar field

iterative::InnerProductSpace is over real vectors. Time-harmonic Maxwell,
frequency-domain problems and anything Schrödinger-shaped are complex-valued.
This is not an operator that FeecOperator accommodates; it is a
generalization of iterative itself, and it touches every Krylov method, both
combinators and the direct backend.

InnerProductSpace is the right seam, and it is already abstract over the vector
type. Generalizing the scalar while the crate is being reworked is cheap;
retrofitting it afterwards is a sweep of everything. Worth doing in the same
pass, even though no problem in the crate needs it yet.

Where things live

iterative gains Decomposition, GMRES and the PositiveDefinite marker, and
stays free of FEEC. formoniq::solve is new and holds FeecOperator,
HodgeForm, Structure and the two solvers. formoniq::linalg stays the
backend bridge. Grade0Multigrid and GradeKHodgeHx::solve are deleted, leaving
multigrid as the tower and cycle and hx as a preconditioner constructor.

That leaves an open question about DirectInverse and the faer bridge. They are
purely linear-algebraic, with no grade, mesh or forms, and they sit in formoniq
only because it was the one crate running a direct solve. If iterative becomes
the crate answering "invert this operator", direct and iterative both, they
belong there, and the crate's name should change with its subject. The
eigensolver is a separate concern and can stay where the shift-invert pencil is
FEEC-shaped.

Known gaps

• Essential boundary conditions. The refinement tower's prolongation is not yet
  defined on RelativeWhitneyComplex: P must commute with the inclusion E,
  which holds for uniform refinement but is neither written nor tested. This is
  the only part of the plan that can fail, so it goes first.
• The IRK stage system is nonsymmetric even for the heat equation, because the
  tableau is. A block preconditioner over the stages built from shifted
  HodgeForm inverses driving GMRES is the general answer; a diagonally implicit
  tableau makes the stages sequential and avoids the question.
• LinearIrk factors once and solves per step for the whole integration, where a
  direct factorization genuinely wins. Direct stays a first-class choice, not a
  fallback.

Laws / validation

• The preconditioned iteration count is bounded under refinement for the mixed
  Hodge-Laplace source problem, swept over grades, against the direct-block
  baseline: the property formoniq: FEEC multigrid and auxiliary-space preconditioning #115 established for the standalone operator, now
  through the problem API.
• Direct and Subspace agree to solver tolerance on every problem: the
  preconditioner changes the path, never the fixed point.
• Subspace correction reproduces the existing V-cycle and auxiliary-space
  results once both are expressed as a Decomposition, so the refactor is
  behaviour-preserving where it is not an extension.

Follow-up to #115, which this completes; the intrinsic vector-nodal space
remains with #119.

DISCLAIMER: This is a design proposal and has not been implemented or verified.