Linear solvers

A nonlinear solver solves one linear system per Newton iteration. Which linear solver it uses is its linsolver option: the name of a *linearSolver block, or a type of linear solver with its defaults.

*solver, solver=NIST, name=theSolver
linsolver=pardiso

A *linearSolver block defines a linear solver with options, as *material defines a material: the keyword line holds only its identity, name and type; the data lines are the options of that type, one key=value per line (an option the type does not have is an error):

*linearSolver, name=fast, type=blockamg
cornerNodeInterpolation=false
verbosity=info

*solver, solver=NIST, name=theSolver
linsolver=fast

A type name in linsolver= means that type with its defaults, so a block must not be named like a type. linsolver unset means pardiso if it was built, else superlu. From Python, the same linear solver is getLinSolverByName("blockamg", {"rigidBodyModesPerAggregate": False}) (edelweissfe.config.linsolve).

The types of linear solver
Type	Direct	What it is
`superlu`	✓	SciPy’s sparse LU (SuperLU), always available
`umfpack`	✓	UMFPACK through SciPy (SuperLU if `scikit-umfpack` is missing)
`pardiso`	✓	Intel MKL PARDISO (`edelweissfe.linsolve.pardiso.pardiso`), the default
`panuapardiso`	✓	Panua PARDISO (`edelweissfe.linsolve.panuapardiso.panuapardiso`)
`klu`	✓	SuiteSparse’s KLU (`edelweissfe.linsolve.klu.klu`)
`petsclu`	✓	PETSc’s LU, through `petsc4py`
`mumps`	✓	MUMPS, through its Python bindings
`blockamg`	✗	the preset `blockAMG()`: LGMRES around a block Gauss–Seidel sweep over the fields, with one AMG hierarchy (or two-level cycle) per field
`custom`	—	a composed solver defined by you, in a Python script (see Composing a linear solver)

Direct or iterative

Direct solvers factorize the matrix and back-substitute: superlu and umfpack (SciPy, always available), pardiso (Intel MKL, the production default), panuapardiso, klu, petsclu and mumps. Pick one of these unless the factorization has become the run’s bottleneck or exceeds available memory. A direct solver keeps its backend for the whole step: PARDISO with reuseSymbolicFactorization=true keeps its reordering and symbolic factorization while only the values of the matrix change. A singular matrix gives NaN, which the nonlinear solver turns into a cutback.

Iterative solvers never factorize the full matrix, trading exactness for O(n) memory: past roughly a million DOFs the fill-in of a direct factorization exceeds memory, whereas algebraic multigrid stays linear. blockamg is the field-split solver for large coupled multi-field systems (e.g. displacement + gradient-enhanced damage): one AMG hierarchy per physical field, combined by a block Gauss–Seidel sweep to precondition an outer LGMRES. The block structure – which DOFs belong to which field, and each field’s dimension – is discovered automatically from the model. It requires the optional amgcl extension.

*linearSolver, name=fieldSplit, type=blockamg
sweeps=1
symmetric=false

*solver, solver=NIST, name=theSolver
linsolver=fieldSplit

A custom solver script composes any other tree: a Schur complement that condenses the dofs on which AMG alone stalls (bars tied to their host elements by bond springs, a softening zone), other smoothers, other coarse spaces. How to write one, and the tree of blockamg, are on Composing a linear solver; why the field-split method works, and each of its steps, on Theory: the field-split AMG preconditioner.

One adapter for every type

Every type is a composed solver (Composing a linear solver) behind one adapter, ComposedLinearSolver: blockamg is the tree of the blockAMG() preset, and a direct solver is a tree of one leaf, type=pardiso being DirectSolver("pardiso"). What the adapter does therefore works the same for every type, with the same data lines (LinearSolverSchema, described in The adapter):

The data lines of every type
Data line	Default	What it sets
`verbosity`	`warning`	`silent`, `warning`, `info` or `debug` (Verbosity levels)
`warningOuterIterations`	100	a Krylov solve with more outer iterations is warned about, and is a slow solve
`residualGrowthFactor`	4	a right-hand side this many times larger than the previous one marks a new increment
`dumps`, `dumpDirectory`, `maxDumps`, `dumpContextSolves`	`none`, `linsolveDumps`, 10, 0	which solves are written to disk, as a sequence of authentic systems for an offline replay (Dumps of the solves)

The options of each type

``pardiso``:
Option	Type / Default	Description
`verbosity`	str, default=’warning’	silent, warning (a warning for a slow or unconverged solve), info (one line per solve, the structure of the solver, and the ‘linear solve’ column of the nonlinear solver) or debug (also why a preconditioner was rebuilt or reused).
`warningOuterIterations`	int, default=100	A solve with more outer Krylov iterations than this is warned about, and is a slow solve.
`residualGrowthFactor`	float, default=4.0	A right-hand side this many times larger than the previous one marks a new increment.
`dumps`	str, default=’none’	Which solves are written to disk for an offline replay: none, slow (above warningOuterIterations, or stopped at the iteration cap), every, or the numbers of the solves separated by spaces, e.g. ‘12 13 14’.
`dumpDirectory`	str, default=’linsolveDumps’	Where the dumps go.
`maxDumps`	int, default=10	The most systems dumped to dumpDirectory in the run (a disk guard).
`dumpContextSolves`	int, default=0	How many of the solves preceding a dumped one are dumped with it.
`reuseSymbolicFactorization`	bool, default=False	Keep PARDISO’s reordering and symbolic factorization while the sparsity pattern stays the same. Off by default: it has been observed to give silently wrong results for some coupled problems.

The other direct solvers take only the common data lines above.

``blockamg``:
Option	Type / Default	Description
`verbosity`	str, default=’warning’	silent, warning (a warning for a slow or unconverged solve), info (one line per solve, the structure of the solver, and the ‘linear solve’ column of the nonlinear solver) or debug (also why a preconditioner was rebuilt or reused).
`warningOuterIterations`	int, default=100	A solve with more outer Krylov iterations than this is warned about, and is a slow solve.
`residualGrowthFactor`	float, default=4.0	A right-hand side this many times larger than the previous one marks a new increment.
`dumps`	str, default=’none’	Which solves are written to disk for an offline replay: none, slow (above warningOuterIterations, or stopped at the iteration cap), every, or the numbers of the solves separated by spaces, e.g. ‘12 13 14’.
`dumpDirectory`	str, default=’linsolveDumps’	Where the dumps go.
`maxDumps`	int, default=10	The most systems dumped to dumpDirectory in the run (a disk guard).
`dumpContextSolves`	int, default=0	How many of the solves preceding a dumped one are dumped with it.
`rigidBodyModesPerAggregate`	bool, default=True	A (linear) 3D displacement field gets the two-level cycle on its rigid-body modes (else one AMG hierarchy, on the translations).
`cornerNodeInterpolation`	bool, default=True	A quadratic 3D displacement field gets the two-level cycle on its corner nodes (else the case of a linear one).
`tolerance`	float, default=None	A fixed relative tolerance per solve; unset: the Eisenstat-Walker forcing tolerance.
`maxIterations`	int, default=100	The outer iteration cap.
`reusePreconditioner`	bool, default=False	Keep the preconditioner across solves until a new increment or a stale preconditioner asks for a rebuild (false: a new preconditioner for every solve).
`sweeps`	int, default=1	The block Gauss-Seidel sweeps over the fields.
`symmetric`	bool, default=False	A symmetric block Gauss-Seidel sweep (forward, then backward).
`krylovOnScaledSystem`	bool, default=False	Run LGMRES on the diagonally scaled matrix, with true-residual restarts.
`trueResidualRestarts`	int, default=2	The true-residual restarts (only with krylovOnScaledSystem).
`residualGapCompensation`	bool, default=False	Tighten the first pass by the learnt gap between the scaled and the true residual (only with krylovOnScaledSystem).
`dropTolerance`	float, default=0.0	Build every field’s preconditioner on its block without the entries below this (relative).
`lumpDroppedEntries`	bool, default=False	Add the dropped entries to the diagonal (keeps the row sums).

``custom``:
Option	Type / Default	Description
`verbosity`	str, default=’warning’	silent, warning (a warning for a slow or unconverged solve), info (one line per solve, the structure of the solver, and the ‘linear solve’ column of the nonlinear solver) or debug (also why a preconditioner was rebuilt or reused).
`warningOuterIterations`	int, default=100	A solve with more outer Krylov iterations than this is warned about, and is a slow solve.
`residualGrowthFactor`	float, default=4.0	A right-hand side this many times larger than the previous one marks a new increment.
`dumps`	str, default=’none’	Which solves are written to disk for an offline replay: none, slow (above warningOuterIterations, or stopped at the iteration cap), every, or the numbers of the solves separated by spaces, e.g. ‘12 13 14’.
`dumpDirectory`	str, default=’linsolveDumps’	Where the dumps go.
`maxDumps`	int, default=10	The most systems dumped to dumpDirectory in the run (a disk guard).
`dumpContextSolves`	int, default=0	How many of the solves preceding a dumped one are dumped with it.
`script`	str (required)	The solver script: a Python file that binds the composed solver, a tree of components, to its global linearSolver (a path relative to the working directory).
`hotReload`	bool, default=False	Build the solver again whenever the script changes on disk, before a solve.

Linear solvers

Direct or iterative

One adapter for every type

The options of each type

The pages on linear solvers