Composing a linear solver
Every linear solver of EdelweissFE is a composed solver behind one adapter (Linear solvers): type=blockamg is
the blockAMG() preset below, type=pardiso the one leaf DirectSolver("pardiso"), and type=custom the
tree of a solver script of your own.
A composed linear solver is a tree of small objects, built with plain constructors in a Python script next to the
input file. Each object is one textbook component – a Krylov solver, a diagonal scaling, a correction on one field, a
multiplicative combination of corrections, an AMG hierarchy, a Schur complement – and its __call__ is the formula
of the textbook. There is one way to build each method, and it is the fast one – with one documented exception, the
two-level cycle (see The two-level cycle: two ways, on purpose).
This page teaches how to build a tree. The arguments and data lines of every component and of the adapter are on Composed linear solvers: reference, the method of the field-split preconditioner on Theory: the field-split AMG preconditioner, and the Python API on Linear solvers: Python API.
A solver script
# solver.py -- the linear solver of a displacement / nonlocal damage model
from edelweissfe.linsolve.composition import (
AMG, LGMRES, DiagonalScaling, Multiplicative, OnField, RigidBodyModesPerAggregate, TwoLevel,
)
# the displacement field: a two-level cycle on the rigid-body modes, AMG below the first coarse level
displacement = OnField("displacement", TwoLevel(coarseSpace=RigidBodyModesPerAggregate()))
# the nonlocal damage field: one AMG hierarchy
damage = OnField("nonlocal damage", AMG())
# both fields, one block Gauss-Seidel sweep: displacement first, then the nonlocal damage
blockGaussSeidel = Multiplicative([displacement, damage])
# accelerated by LGMRES
linearSolver = LGMRES(DiagonalScaling(blockGaussSeidel))
and in the input file:
*linearSolver, name=mySolver, type=custom
script=solver.py
verbosity=info
*solver, solver=NIST, name=theSolver
linsolver=mySolver
LGMRES’s defaults are the production ones: an Eisenstat–Walker tolerance per solve (EisenstatWalkerTolerance), a
RebuildPolicy that builds the preconditioner anew for every solve (RebuildPolicy(reusePreconditioner=True), in
an input file the blockamg data line reusePreconditioner=true, keeps it across Newton iterations,
Rebuilding the preconditioner for every solve), and the detection of a new increment
through the adapter’s residualGrowthFactor; for a plain solve to a fixed tolerance, write
LGMRES(..., tolerance=FixedTolerance(1e-8)).
script is a path relative to the working directory, like the output files of the run (only *include is
relative to the including input file): script=solver.py finds a script next to the input file when the simulation
is started in the directory of the input file, as run_tests_edelweissfe does.
The script defines the global linearSolver and holds only the math: the tree. What the adapter does around it –
the verbosity, the dumps, the hot reload – are data lines of the block (a ComposedLinearSolver bound in the script
is refused; see The adapter). The journal shows the structure of the solver and a hash of the script in a frame at the
first solve (with verbosity=info), above the header of the iteration table.
How to read the nesting
The script defines the components bottom-up, with named variables: first the basic method (the block Gauss–Seidel sweep over the fields), then what is wrapped around it, and the Krylov solver last – the textbook order “build a basic iteration, then accelerate it”.
At run time the nesting means delegation: the outer component calls its children, many times, before and after its
own steps. LGMRES calls its preconditioner once per iteration; the diagonal scaling multiplies by \(D\) before
and after the sweep; a Schur complement solves with \(A_{bb}\) before and after its child. These are
sandwiches, not a pipeline – which is why there is no | syntax.
Where a component sits in the tree says what it acts on. SchurComplement(..., DiagonalScaling(Multiplicative(...)))
scales only what the sweep sees, while the outer Krylov solve still converges the true residual of \(A\).
Components
Component |
\(M^{-1}\) |
|---|---|
|
the outer Krylov solve \(A x = b\), right-preconditioned (AMGCL); owns when the preconditioner is rebuilt
( |
|
flexible GMRES (AMGCL): \(x = Z_m y_m\) from the preconditioned vectors themselves, so the preconditioner may change between iterations – e.g. an inner Krylov solve |
|
the tolerance per solve: Eisenstat–Walker, or fixed |
|
\(D \, \hat M^{-1} D\), \(\hat M\) built on \(\hat A = D A D\), \(D = \operatorname{diag}(\lvert a_{ii} \rvert^{-1/2})\) |
|
a correction on one field, \(R^T M^{-1} R\), \(M\) built on the field’s block \(A_{ff}\) |
|
successive corrections \(z \leftarrow z + M_k^{-1}(r - A z)\), or summed ones \(\sum_k M_k^{-1} r\);
over |
|
one AMGCL hierarchy on a field, one cycle per application. |
|
a two-level cycle with a coarse space of its own (a |
|
block LU: \(z_k = M_{\tilde S}^{-1}(r_k - A_{kb} A_{bb}^{-1} r_b)\),
\(z_b = A_{bb}^{-1}(r_b - A_{bk} z_k)\); |
|
\(A^{-1}\) by a registered direct solver package ( |
|
\(P A_c^{-1} P^T\) with \(A_c = P^T A P\), and \(\omega \operatorname{diag}(A)^{-1}\): the pieces
of a two-level cycle spelled out, |
The combinators show that methods which look unrelated share one structure. Block Gauss–Seidel over the fields is
Multiplicative([OnField("displacement", M_u), OnField("nonlocal damage", M_d)]), and a two-level cycle is
Multiplicative([smoother, CoarseCorrection(P, coarseSolve), smoother]). They are also fast: a combination asks each
correction which field it acts on, and forms the residual of a correction on field \(i\) only on that field’s
rows, \(r_i - \sum_{j: z_j \neq 0} A_{ij} z_j\), from the field blocks of the matrix (formed once, in one
threaded pass, for all corrections).
Wrong compositions fail as early as possible: at construction when the tree alone shows it (a child that is not a
component, an empty combination, LGMRES around a Krylov solve, a Schur complement with nothing to condense), at the
first set-up otherwise, with a message that names the cause (a field the system does not have, a one-field component
such as AMG on several fields, rigid-body modes on a scalar field or without node coordinates).
The two-level cycle: two ways, on purpose
A two-level cycle can be written in two ways, and both are kept deliberately:
TwoLevel(coarseSpace=RigidBodyModesPerAggregate()) # the fast way
Multiplicative([Jacobi(), CoarseCorrection(P, AMG()), Jacobi()]) # the spelled-out way
TwoLevel is a leaf, like AMG: pre-smoothing (Chebyshev), the coarse correction
\(P A_c^{-1} P^T (r - A z)\) with \(A_c = P^T A P\), and post-smoothing run as one cycle in C++. Use it for
production runs; blockAMG() (type=blockamg) uses it for 3D displacement fields.
The spelled-out way builds the same cycle from general components, each of which can be replaced – another smoother,
a coarse space of your own (any \(P\)), a direct coarse solve. Use it to learn how a two-level cycle works and to
try a new coarse space or smoother. It is correct, but slower than TwoLevel, at the same iterations
(The two-level cycle written out). Not because it is Python: its products all run threaded in C++ –
Multiplicative forms its residuals in the threaded pass of the field split, and CoarseCorrection forms
\(P^T r\), \(P z\) and \(P^T A P\) with AMGCL’s threaded products – but because of the smoother.
TwoLevel keeps one compact copy of the field block (single precision, 3x3 node blocks) for its Chebyshev smoother,
which reads half the bytes per product; a smoother component works on the general double-precision matrix.
The testfile deck testfiles/edelweiss-only/TwoLevelSpelledOut runs the spelled-out cycle in a simulation, on
quadratic elements with the P1 interpolant as the coarse space. Its prolongation callable returns \(P\) together
with the context of the coarse dofs, context.ofCoarseNodes(field, isCorner) – the displacement field on the corner
nodes, with their coordinates – so its coarse solve is the rigid-body coarse level,
DiagonalScaling(TwoLevel(coarseSpace=RigidBodyModesPerAggregate())), the coarse solve the blockAMG() preset
gives a 3D displacement field but for its smoother: this deck keeps the field’s default tree (degree 16), the preset
takes the weaker smoother of the corner-node two-level cycle (degree 8, Default field trees). (A
callable whose
coarse dofs are no nodes returns context.ofCoarseDofs(nCoarse), one field of scalar dofs; a fixed matrix \(P\)
gets that context.) \(P\) acts on the matrix the correction is handed: the deck has no DiagonalScaling above
the cycle, so its variables are the displacements, and the plain weights \(\tfrac12, \tfrac12\) are right. Under a
DiagonalScaling the variables are \(D^{-1} u\), and the same interpolation of displacements is
\(D^{-1} P_{geo} D_c\) – the scaled weights of CornerNodeInterpolation (p-multigrid: TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve))), with \(D\)
from context.scaling and \(D_c\) from the coarse context’s scaling. Unlike TwoLevel, the deck’s \(P\)
keeps the Dirichlet rows, so a Dirichlet corner dof is no identity row of \(A_c\); the coarse solve treats it as
free, which is still a valid preconditioner. Its smoother is Jacobi(omega=0.4): damped Jacobi
damps the fast error modes only for \(\omega < 2 / \lambda_{max}(\operatorname{diag}(A)^{-1} A)\), and on
quadratic elements \(\lambda_{max}\) is larger than for the Laplacian the default \(\omega = 2/3\) is made
for (about 3.7 there), so the default amplifies them and the cycle stalls.
There is no native Chebyshev component yet, so the spelled-out way cannot reproduce TwoLevel exactly; a
Chebyshev component is issue #237 (https://github.com/Edelweiss-Numerics/EdelweissFE/issues/237). The exception
goes away once the remaining two pieces exist natively: the compact field-block copy kept by the system matrix, and a
Chebyshev smoother component on it (CoarseCorrection already multiplies with \(P\) threaded). Then TwoLevel
is replaced by its spelled-out form.
Which form to use: what each costs
Every component composes with every other, and many methods can be written in more than one way. The ways do not cost the same. The rule: a written-out form is for learning and for trying something new; for production runs, take the native component where one exists. What each form costs, as far as it was measured:
Written as |
Costs |
For production runs |
|---|---|---|
a two-level cycle spelled out, |
more wall time than |
|
\(P\) built by a callable of your own, in NumPy or in a Python loop over the nodes (the deck
|
one thread, at every set-up: on a large model a serial pass over its nodes or its matrix (e.g. a search for
the Dirichlet rows in NumPy, The two-level cycle written out) is a noticeable share of a set-up of
|
|
|
a weaker smoother than Chebyshev, and only one that damps at all for \(\omega < 2 / \lambda_{max}\) (0.4, not the textbook 2/3, on quadratic elements); not measured against Chebyshev on large systems |
the Chebyshev smoother inside |
|
double-precision scalar entries, more bytes per product: the corner-node two-level cycle was slower per application than on the default single-precision 3x3 blocks (The two-level cycle written out) |
the defaults of The AMGCL backend of a field; the scalar backend only where its near null-space is needed (all rigid-body modes in one AMG hierarchy) |
a component of your own in NumPy or SciPy (Writing your own component) |
NumPy and SciPy run on one thread: SciPy’s sparse product too, whatever |
fine for trying an idea; keep its products with \(A\) few, or make them native once it pays |
a preset copied into a solver script and changed |
nothing: the components are the same as the preset’s, and the composition of |
– |
The combinators themselves cost nothing extra: Multiplicative and Additive form their residuals with the
threaded products of the field split, and CoarseCorrection its products with \(P\) and \(P^T A P\)
threaded too. What makes a written-out form slower is what it is written out of.
Presets
Presets are plain functions that return a tree. blockAMG() is the blockamg type, its arguments the data lines
of that type (*linearSolver, type=blockamg), and in a script linearSolver = blockAMG(...). It builds:
LGMRES(
DiagonalScaling(
Multiplicative([OnField(f, M_f) for every field f]) # M_f: TwoLevel(RigidBodyModesPerAggregate())
) # for a 3D displacement field, else AMG()
)
The corrections are made at the first set-up, for the fields of the model; sweeps and symmetric order them
with gaussSeidelOrder. blockAMG(krylovOnScaledSystem=True) puts the scaling around LGMRES instead
(A Krylov solve on the scaled system), trueResidualRestarts and residualGapCompensation go to LGMRES, and dropTolerance and
lumpDroppedEntries to every field (A preconditioner on the filtered matrix).
Each field gets the preconditioner chosen for it by measurement (blockAMGFieldPreconditioner,
Why the defaults are what they are):
a quadratic 3D displacement field (one with midside nodes) the two-level cycle on its corner nodes,
TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve))(p-multigrid: TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve))), whose coarse solve, under a diagonal scaling of its own, is that of the next case, with the smoother of the corner-node two-level cycle: three levels, the quadratic field, its linear part on the corner nodes, and the rigid-body modes of aggregates of those; the corner and midside nodes come from the shape each element declares (ensightType:hexa20,quad8andbar3have midside nodes;point,bar2,tria3,quad4,tetra4,penta6andhexa8corners only). An element of any other shape (e.g.tetra10) is an error, never skipped: switch the corner-node interpolation off for such a model (cornerNodeInterpolation=False). An element may have more nodes than its shape – a bond element has the nodes of its bar, then those of the host element it couples to –; each node after its shape’s must be a node of another element’s shape, else that is an error too;a linear 3D displacement field the rigid-body coarse level,
TwoLevel(coarseSpace=RigidBodyModesPerAggregate()): a first coarse level built from all rigid-body modes, rotations included, which a bending-dominated structure needs (Rigid-body coarse level: TwoLevel(coarseSpace=RigidBodyModesPerAggregate()));any other field one AMG hierarchy,
AMG().
Each two-level cycle can be switched off; the rigid-body coarse level where the solves converge in a few outer iterations without it – there it saves little, for a larger setup:
*linearSolver, name=translationsOnly, type=blockamg
rigidBodyModesPerAggregate=false
cornerNodeInterpolation=false
A 2D displacement field, or any other vector field, gets the coarse level in a custom solver script,
OnField("displacement", TwoLevel(coarseSpace=RigidBodyModesPerAggregate())).
Everything else – another tree for a field, the rebuild policy, the parameters of the forcing tolerance or of LGMRES,
the options of the rigid-body coarse level, the p-multigrid variant, a Schur complement – is a custom solver
script: copy the preset (edelweissfe.linsolve.composition.presets, a plain function) into the script and change
what you need.
A Schur complement
Some dofs make the field-split AMG stall: the dofs of embedded bars with bond-slip, for instance, tied to the
nodes of their host elements by stiff bond springs. A SchurComplement condenses such dofs exactly and builds its
child on the approximate Schur complement of all other dofs (the method: Schur complement condensation: SchurComplement (opt-in)). The
condensed dofs are chosen with the markers of adaptive mesh refinement: the dofs of the interior nodes of the
marked elements – the nodes all of whose elements are marked. For embedded bars, mark the bar and the bond elements:
# solver.py
from edelweissfe.config.markerlibrary import getMarkerClass
from edelweissfe.linsolve.composition import (
LGMRES, DiagonalScaling, Multiplicative, OnField, SchurComplement, gaussSeidelOrder,
)
from edelweissfe.linsolve.composition.presets import blockAMGFieldPreconditioner
bars = getMarkerClass("elementSet").fromOptions({"elSet": "bars"})
bonds = getMarkerClass("elementSet").fromOptions({"elSet": "bondElements"})
def fieldCorrections(context):
corrections = [OnField(f.name, blockAMGFieldPreconditioner(f, context)) for f in context.fields]
return gaussSeidelOrder(corrections, sweeps=2)
blockGaussSeidel = Multiplicative(fieldCorrections)
linearSolver = LGMRES(SchurComplement(condensedDofMarkers=[bars, bonds], schurComplementSolve=DiagonalScaling(blockGaussSeidel)))
*linearSolver, name=barSchur, type=custom
script=solver.py
*solver, solver=NISTParallel, name=theSolver
linsolver=barSchur
The corrections here are a function of the context: they are made at the first set-up, one per field of the remaining
dofs, with the preconditioner the preset would give that field. A fieldOutput marker condenses, e.g., a softening
zone instead (testfile GradientEnhancedDamageSchurComplement). The other arguments of SchurComplement are
listed in The arguments of SchurComplement.
Writing your own component
A component of your own subclasses Preconditioner – in the solver
script itself, if you like – and calls super().__init__() in its constructor. It then composes with every other
component. The matrix it is given, A, is a SystemMatrix,
not a SciPy matrix: A.matrix() gives the SciPy CSR matrix (the system matrix also forms the field blocks and the
scaled matrix for the components that need them).
Note
A component of your own in NumPy or SciPy runs on one thread (SciPy’s sparse products too): fine for trying an idea, but slower than the native components on a large model (Which form to use: what each costs).
Four hooks make a component:
_setUp(A, context)Build what the application needs for the matrix
A, e.g. the inverse of its diagonal.context(aSetUpContext) holds the fields of the sub-problem and their node coordinates._forget()Drop what
_setUpbuilt: set each of those attributes to its value before the first set-up (self.inverseDiagonal = None). The constructor of the base class calls it, so this is where they are declared. The publicsetUpcallsforget()(this hook and those of the children) before_setUp: the old factorization or hierarchy is freed before the new one is built, never both at once, which on a large model is gigabytes. What is reused on purpose (a symbolic factorization, an estimate carried from one build to the next) is not forgotten._apply(r)Return \(z = M^{-1} r\).
_update(A, change)Take the matrix of the next Newton iteration. By default it builds anew (it calls
self.setUp); override it to reuse what you built, rebuilding (again byself.setUp) only ifchange.requiresRebuild(a new increment, a new sparsity pattern, a stale preconditioner, as theMatrixChangeand the Krylov solver above decide) or if nothing is built – afterforget(), which a Krylov solver calls on its preconditioner before it rebuilds it, and the adapter on the whole tree when the model changes.
Then describe() returns a Description – the name, the formula
(unicode and LaTeX), the parameters and the steps of one application – which the views show. The public
setUp, update and __call__ of the base class record each call and its time in the component’s statistics
and call the hook; a component calls its children through those public methods.
A complete component, run as a deck of the test suite (testfiles/edelweiss-only/OwnComponentSolverScript):
# A component of our own, composed with those of the library: a symmetric Gauss-Seidel smoother before and after one
# AMG cycle, accelerated by LGMRES.
from scipy.sparse import tril, triu
from scipy.sparse.linalg import spsolve_triangular
from edelweissfe.linsolve.composition import (
AMG,
LGMRES,
Description,
FixedTolerance,
Multiplicative,
Preconditioner,
)
class SymmetricGaussSeidel(Preconditioner):
"""One forward and one backward Gauss-Seidel sweep: z = (D + U)⁻¹ D (D + L)⁻¹ r, with A = L + D + U split into its
strictly lower triangle, its diagonal and its strictly upper triangle."""
def _forget(self):
# what _setUp builds, before it is built (also called by the constructor of Preconditioner)
self.lower = self.upper = self.diagonal = None
def _setUp(self, A, context):
# A is a SystemMatrix, as every component gets it; A.matrix() is the SciPy matrix
matrix = A.matrix()
self.lower = tril(matrix, format="csr") # D + L
self.upper = triu(matrix, format="csr") # D + U
self.diagonal = matrix.diagonal() # D
def _update(self, A, change):
# only the values changed: the triangles of the last matrix still make a good smoother for the new one
if change.requiresRebuild or self.lower is None:
self.setUp(A, change.context)
def _apply(self, r):
y = spsolve_triangular(self.lower, r, lower=True) # forward sweep: (D + L) y = r
return spsolve_triangular(self.upper, self.diagonal * y, lower=False) # backward sweep: (D + U) z = D y
def describe(self):
return Description(
name="SymmetricGaussSeidel", formula="(D + U)⁻¹ D (D + L)⁻¹", latex="(D + U)^{-1} D (D + L)^{-1}"
)
# smooth, one AMG cycle, smooth again: one smoother, applied twice
gaussSeidel = SymmetricGaussSeidel()
cycle = Multiplicative([gaussSeidel, AMG(), gaussSeidel])
linearSolver = LGMRES(cycle, tolerance=FixedTolerance(1e-8))
The rest of the protocol of Preconditioner has defaults that fit a
leaf; override a method only in the case named:
children(): a component with children returns its(role, child)pairs, in the order it uses them. The tree is walked through it (the model, the end of a solve, the questions below, the views), and the role is the placeholder of the child in the component’s formula.fieldIndexand_applyRestricted(r_f): a component that acts on one field returns its index in the context; a combination then forms the residual of that field alone and callsapplyRestrictedwith it.withRebuildReasons(change): a component that adds reasons for a rebuild of its own (a Krylov solver’sRebuildPolicy), or keeps something built for the matrix above such a child (the \(D\) of aDiagonalScalingaround a Krylov solve), so that it rebuilds together with the child.setModel(model, dofManager): a component that derives something from the model, e.g. which dofs a marker selects.finishSolve(): a component that holds something only one solve needs (e.g. the field blocks of its matrix) releases it here, so that no copy of the last matrix is kept while the nonlinear solver assembles the next one.isFixedLinearOperator: false for a component whose \(M^{-1} r\) is not one fixed linear operator between two set-ups (a Krylov solve inside);LGMRESrefuses such a preconditioner,FGMREStakes it.containsKrylovSolver: true for a Krylov solver (the solve is recorded in the “linear solve” column).containsDirectSolver: true for a direct solver (a sparse LU factorization); the nonlinear solver then prunes the explicitly stored zeros out of the multi-point-condensed system matrix (pruneCondensedMatrixZeros), which such a factorization wants, unless the tree reuses its preconditioner.reusesPreconditioner: true for a Krylov solver that keeps its preconditioner across solves (RebuildPolicy(reusePreconditioner=True)); the nonlinear solver then does not prune, which would change the sparsity pattern in almost every iteration.builtMatrices(),solveRecord(),condensedDofs(): the matrices a component formed (forviews.spy), what its last Krylov solve did, and the dofs a Schur complement condensed (for the dumps).
Views
Every component describes itself (describe()); the views in edelweissfe.linsolve.composition.views draw
the descriptions along the tree of children(). For the blockAMG() preset on a displacement / nonlocal damage
system, after two solves:
print(solver) – the structure:
LGMRES(tolerance=EisenstatWalkerTolerance(minimum=0.0001, maximum=0.001), maxIterations=100, restart=50, augmentationVectors=3)
└─ DiagonalScaling(D=diag(1/√|aᵢᵢ|))
└─ Multiplicative
├─ M1: OnField(field=displacement)
│ └─ TwoLevel(coarse space=rigid-body modes, prolongation=smoothed, coarse solve=AMG V-cycle)
└─ M2: OnField(field=nonlocal damage)
└─ AMG(amgclTree=default)
views.formula(solver) – the composed approximate inverse (latex=True for LaTeX):
LGMRES(A, D · (Rᵀ AMG(nonlocal damage)⁻¹ R ∘ Rᵀ TwoLevel(displacement; rigid-body modes)⁻¹ R) · D)⁻¹
views.runStatistics(solver) – what each component recorded (inclusive times):
LGMRES(...) updated 2× 0.07s · applied 2× 0.05s
└─ DiagonalScaling(...) updated 2× 0.07s · applied 18× 0.04s
└─ Multiplicative set up 1× 0.06s · updated 2× 0.06s · applied 18× 0.04s
├─ M1: OnField(...) set up 1× 0.06s · updated 1× 0.00s · applied 18× 0.03s
│ └─ TwoLevel(...) set up 1× 0.06s · ... · applied 18× 0.03s · rows 9555 · aggregates 198 · coarse dofs 1188
└─ M2: OnField(...) set up 1× 0.00s · updated 1× 0.00s · applied 18× 0.00s
└─ AMG(...) set up 1× 0.00s · ... · applied 18× 0.00s · rows 3185
The adapter only ever updates the root; a combination whose corrections are made per field sets them up at its first update (a set-up within that update), and reuses them at the next.
views.dataFlow(solver) – one application as a graphviz graph (DOT, or SVG with the graphviz package): each
component a box around its steps, a step that delegates leading into its child’s box. views.spy(solver, A) draws
the sparsity patterns of the matrices the components were built on, views.writeHtml(solver, "solver.html") writes
all views into one page, and in Jupyter a solver shows itself as that page.