Module linear_operator 

Layer 2 — the discretisation-agnostic LinearOperator contract and the generic Krylov drivers (cg_op, bicgstab_op, gmres_op) written against it, so a finite-element crate can reuse this solver layer without adopting the finite-volume LduMatrix layout (GitHub issue #175). Purely additive: the existing finite-volume solvers are unchanged. Discretisation-agnostic linear-algebra contract, plus Krylov solvers written against it.

What belongs in this module

The LinearOperator contract — “I am an n x n linear map you can apply to a vector, and I can hand you my diagonal” — and the OperatorPreconditioner contract, together with the three matrix-free Krylov drivers cg_op, bicgstab_op and gmres_op that need nothing more than those two contracts.

What does NOT belong here

Any storage format. This module stores no coefficients of its own: it is the seam at which a representation (LDU for the face-addressed finite-volume side, CSR for a finite-element side) meets the solvers.

Why it exists

Before this module every Krylov entry point in the crate (crate::krylov::gmres, crate::krylov::bicgstab, [crate::ldu_matrix::solvers::conjugate_gradient], GAMG, Gauss-Seidel) took a concrete &LduMatrix. LduMatrix is a good representation for a finite-volume mesh — one coefficient per internal face, owner/neighbour addressing straight off the mesh — and a poor one for a finite-element stiffness matrix, whose degree-of-freedom graph is denser than the cell graph and has no face to hang a coefficient on.

A finite-element crate therefore had exactly two bad options: reformulate structural mechanics as finite volume so it could reuse the solvers, or reimplement the solvers so it could keep its own matrix. This module is the third option (GitHub issue #175): one Krylov layer, two matrix layouts.

Additive by construction

Nothing in this module changes any existing function. The existing finite-volume solvers keep their &LduMatrix signatures and their exact behaviour; the functions here are new generic siblings. LduMatrix implements LinearOperator, so the finite-volume side may use them too, but is not obliged to.

No trait objects

Per the workspace design rules, LinearOperator and OperatorPreconditioner are used only as generic bounds (fn cg_op<A: LinearOperator, M: OperatorPreconditioner>), never as Box<dyn LinearOperator>. Each call site monomorphises to one concrete matrix type, so the matrix-vector product inlines exactly as it would have done in a hand-written solver.

Units

Everything here is dimensionless f64. A Krylov subspace mixes residuals, search directions and solution increments that share no single physical dimension, so no uom typing is applied — exactly as in crate::krylov. Apply units at the layer that assembles the matrix.

Example

use outram_foam_basic_lib::krylov::KrylovSettings;
use outram_foam_basic_lib::ldu_matrix::LduMatrix;
use outram_foam_basic_lib::linear_operator::{cg_op, JacobiOperatorPreconditioner};

// Symmetric positive-definite 3-cell chain.
let mut a = LduMatrix::new(3, vec![0, 1], vec![1, 2]);
a.diag = vec![4.0, 4.0, 4.0];
a.lower = vec![-1.0, -1.0];
a.upper = vec![-1.0, -1.0];
let b = vec![1.0, 2.0, 3.0];

let m = JacobiOperatorPreconditioner::new(&a);
let settings = KrylovSettings { tolerance: 1e-12, max_iter: 100, restart: 30 };
let (x, result) = cg_op(&a, &b, None, &m, &settings);

assert!(result.converged);
let mut ax = vec![0.0; 3];
outram_foam_basic_lib::linear_operator::LinearOperator::apply(&a, &x, &mut ax);
for i in 0..3 {
    assert!((ax[i] - b[i]).abs() < 1e-9);
}

Structs

IdentityPreconditioner

No preconditioning: M = I, so z = r.

JacobiOperatorPreconditioner

Diagonal (Jacobi) preconditioning built from any LinearOperator: z[i] = r[i] / A[i][i].

Traits

LinearOperator

A square linear map A that can be applied to a vector.

OperatorPreconditioner

An approximate inverse M^{-1} ~ A^{-1} usable by the generic Krylov drivers.

Functions

bicgstab_op

Preconditioned BiCGStab for a general (possibly non-symmetric) operator.

cg_op

Preconditioned conjugate gradients for a symmetric positive-definite operator.

gmres_op

Restarted right-preconditioned GMRES(m) for a general operator.