Linear Solvers¶
Once a model has assembled and constrained its system \(A\,\mathbf{x} = \mathbf{b}\), a linear solver computes the solution \(\mathbf{x}\). This is the innermost and often most expensive part of a simulation. mufem offers two families, direct and iterative, selected per model on its Solver page (for example the Time-Harmonic Magnetic Solver).
Direct solvers¶
A direct solver computes \(\mathbf{x}\) by factorizing the system matrix (a sparse \(LU\) or Cholesky decomposition) and then performing forward and backward substitution. mufem uses MUMPS, a parallel sparse direct solver, for both real- and complex-valued systems.
Direct solvers are robust. They solve any non-singular system in a single pass, with no tolerances to tune and no convergence to monitor. Their cost is memory and time. The factorization fills in many of the zero entries of a sparse matrix, so storage and work grow rapidly with problem size, especially in three dimensions. So direct solvers are the natural choice for small to medium problems, or a dependable reference when an iterative solver struggles to converge.
Iterative solvers¶
An iterative solver never forms a factorization. From an initial guess it builds a sequence of better approximations. Each step needs only a few matrix-vector products with \(A\). mufem uses Krylov subspace methods: the conjugate gradient method (CG) for symmetric positive-definite systems, and GMRES / flexible GMRES (FGMRES) for the general case.
Because they touch the matrix only through matrix-vector products, iterative solvers keep the sparsity of \(A\) and have a much smaller memory footprint. They scale to large meshes where a direct factorization would be infeasible. The trade-off is that convergence is not guaranteed in a fixed number of steps. It depends strongly on the conditioning of the system, which is where preconditioning enters.
The iteration is governed by a stopping criterion and an iteration budget:
| Control | Meaning |
|---|---|
| Relative tolerance | Stop once the residual norm has dropped by this factor relative to the initial residual. |
| Absolute tolerance | Stop once the residual norm falls below this floor, regardless of the initial value. |
| Maximum iterations | Cap on the number of iterations, so a poorly conditioned solve terminates. |
The exact defaults and setters are documented on each model's Solver page.
Preconditioning¶
The conditioning of the system sets the convergence rate of a Krylov solver. A preconditioner \(M\) approximates \(A\) but is cheap to invert. It transforms the system into an equivalent one that iterates far faster:
A good preconditioner makes \(M^{-1} A\) close to the identity, so the solver converges in few iterations. The art is to balance that quality against the cost of applying \(M^{-1}\) every step. mufem provides several, matched to the structure of the problem:
| Preconditioner | Typical use |
|---|---|
| Jacobi / diagonal | Cheapest option; rescales by the matrix diagonal. |
| Algebraic multigrid (BoomerAMG) | Strong general-purpose preconditioner for scalar elliptic problems (for example thermal conduction). |
| Auxiliary-space Maxwell (AMS) | Tailored to \(H(\mathrm{curl})\) edge-element systems that arise in electromagnetics. |
| Block-diagonal | For coupled multi-field systems: preconditions each physical field block separately. |
The multigrid and auxiliary-space preconditioners let the iterative solvers keep a nearly mesh-independent iteration count as you refine the problem. This is the key to scaling to large three-dimensional models.