Theory: the field-split AMG preconditioner
This page explains a field-split algebraic-multigrid preconditioner for the coupled multi-field equation systems EdelweissFE’s implicit solvers produce – displacement coupled to a scalar nonlocal (gradient-enhanced) damage field, typically with penalty contact, tie constraints and adaptive refinement on top – and the outer Krylov solve it preconditions.
The method is a tree of components (Composing a linear solver), LGMRES(DiagonalScaling(Multiplicative([OnField(f,
M_f), ...]))), and every step below is one of them, named in the heading of its section:
DiagonalScaling (the equilibration),
OnField (the field split),
AMG (one hierarchy per field),
Multiplicative (the block Gauss–Seidel sweep),
LGMRES with its
EisenstatWalkerTolerance and
RebuildPolicy (the outer solve, its tolerance, the reuse of the
hierarchies), TwoLevel (the p-multigrid variant and the rigid-body
coarse level) and SchurComplement (the condensation). The
near null-spaces are built in nullspace. The AMG hierarchies and the outer solve need the
optional amgcl extension (edelweissfe.linsolve.amgcl.amgcl).
The preset blockAMG() builds this tree with its defaults; it is the
linear solver type blockamg (linsolver=blockamg, or a *linearSolver, type=blockamg block, see
Linear solvers). Each setting below is named where it lives: an argument of a component in a custom solver
script, and for the settings the preset exposes, a data line of type=blockamg
(Presets). The measurements behind the defaults are collected on their own page,
Field-split AMG: measurements behind the defaults, so that this page reads as the method alone; the arguments and defaults of the
components are on Composed linear solvers: reference.
Why a field-split preconditioner
On large coupled fracture models a direct factorization eventually stops being an option: its fill-in grows superlinearly, and past roughly \(10^6\) DOFs it exceeds the available memory regardless of how much time one is willing to spend. Algebraic multigrid (AMG) has \(O(n)\) memory and is the route to those sizes.
Applied monolithically to the coupled system, however, AMG does not work: a single hierarchy cannot represent the disparate physics and scales of the fields at once. Measured directly on a real coupled system, a monolithic hierarchy’s residual reduction plateaus around 0.2 instead of converging, and it does so regardless of how the near-null-space is chosen – the missing ingredient is the field-block structure itself, not a better single hierarchy.
The remedy, following the block-preconditioning strategy of Alkmim et al. (IJNME 2026), is to precondition blockwise: build one AMG hierarchy per physical field, combine them with a block Gauss–Seidel sweep, and use that sweep as the preconditioner of an outer Krylov solve over the full coupled system. Each field’s operator is individually AMG-friendly – elasticity for a displacement field, a screened-Poisson (Helmholtz) operator for a gradient-enhanced damage field (see Adaptive refinement: error estimation and marking for that operator) – even though their monolithic coupling is not.
The field-split preconditioner buys reach, not necessarily speed.
On the reference models AMGCL’s smoothed aggregation converges, but not tightly, on the
non-symmetric (contact-condensed) displacement block, so the outer Krylov solve needs
many iterations – with the former defaults typically \(O(100)\); with the current
defaults (see Why the defaults are what they are), mostly 5–50. The hardest damage-onset
solves stop at its cap of 100 iterations
(maxIterations), return their iterate, and the hierarchies are rebuilt for
the next solve. The cap is a deliberate budget: a solve that needs more is taken as a sign that the
hierarchies no longer fit the Jacobian, not given more iterations. That is an acceptable price where the
alternative is not fitting in memory at all; at problem sizes a direct
factorization still handles comfortably, a direct solver (e.g. pardiso) can be competitive or
faster. Which one wins depends on problem size, conditioning, and
how severely damage and contact nonlinearity degrade the per-field hierarchies on a given increment.
Where the block structure comes from
Nothing about the block layout is configured by hand. A sparse matrix does not carry the information
of which DOF belongs to which physical field, so it is pushed in from the model side:
setModel() hands every registered linear solver the
live FEModel and
DofManager, and the base class derives from them the
ordered list of FieldBlock records – field name, half-open DOF
range [start, stop), and nodal dimension (3 for a 3D displacement field, 1 for a scalar damage
field). The nonlinear solver calls it whenever the equation system is (re)built: on the first solve,
and again after every adaptive-refinement or connectivity change, so the layout tracks the mesh.
Two consequences worth knowing:
The fields must tile the DOF vector contiguously in DOF order. They do, by construction, for every model the DOF manager builds; a gap or overlap raises an error rather than being silently worked around. DOFs past the last node field (scalar variables) are folded into one trailing scalar block named
scalar variables. In EdelweissFE these are at present the Lagrange multipliers of constraints (rigid bodies, equal-value constraints), whose rows have a zero diagonal: such a system is refused (see the equilibration step below), so the block only serves scalar variables with a diagonal of their own.Beyond the field ranges, the components read node coordinates from the model on their own, lazily, for the near-null-space construction described below (and, for the opt-in p-multigrid variant, the corner/midside element topology). Nothing has to be pushed in ahead of it.
setFieldStructure() is a lower-level escape hatch for a
caller that knows the field blocks but has no model to hand over – e.g. an offline script replaying
a captured system through the solver. On that path node coordinates are unavailable unless given
(The adapter), and an AMG hierarchy degrades gracefully to a translations-only near-null-space.
The solve pipeline
The linear solver sees only the condensed, Dirichlet-eliminated system: the nonlinear solver applies the multi-point-constraint transform \(\hat K = T^\mathsf{T} K T + I_s\) (hanging nodes, ties; the identity \(I_s\) fills the structurally zero slave rows and columns, and the slave DOFs are recovered from their masters after the solve, so the transform keeps a symmetric \(K\) symmetric and adds no slave-to-master coupling that the equilibration below could blow up) and then the symmetric Dirichlet elimination before calling the linear solver: each constrained row becomes the identity row, and each constrained column is zeroed after its coupling to the prescribed value has been moved to the right-hand side – so a symmetric operator stays symmetric, as smoothed aggregation and the Chebyshev smoother assume. Both steps are size-preserving and leave the DOF manager’s ordering intact, so the field blocks still describe the matrix that actually arrives.
The MPC condensation is skipped entirely while no multi-point constraint has a record – e.g. an
hAdaptivity before its first refinement, whose hanging-node constraint is registered up front.
This matters for hierarchy reuse: the condensed matrix would be pruned of exact zeros
(pruneCondensedMatrixZeros, a nonlinear-solver option), which changes its sparsity pattern from
one Newton iteration to the next, and a changed pattern forces a hierarchy rebuild (step (4) below).
Left unset, that option is automatic: it prunes for a linear solver that factorizes the matrix (a direct solver wants
the stored zeros gone) and does not reuse a preconditioner, and leaves the matrix alone for the others
(needsPrunedSystemMatrix): a tree of Krylov solvers and AMG gets
the same iterations with and without the stored zeros, so pruning would only cost it a copy of the matrix and a
pass over it in every Newton iteration (0.8 s of a 1 s call on a model with 376 million stored nonzeros), and a Krylov
solve that reuses its preconditioner (RebuildPolicy(reusePreconditioner=True)) must not have its pattern changed
anew in every iteration. An explicit
pruneCondensedMatrixZeros = true or false in the nonlinear solver’s own configuration overrides the choice.
The uncondensed matrix keeps the assembler’s fixed pattern.
Beyond that matrix and right-hand side, the linear solver is told nothing about the Newton loop – like
every registered linear solver it is called as (A, b) -> x. Both of its stateful mechanisms
(adaptive tolerance and hierarchy reuse) therefore reconstruct the loop state from
\(\lVert b \rVert\) and from the matrix’s own bookkeeping across calls.
The symbols of this page, their variables in the worked example (A worked example), and the components that hold them:
Symbol |
Meaning |
Worked example |
Component |
|---|---|---|---|
\(A, b, x\) |
the system as the nonlinear solver hands it over, and its solution |
|
the adapter’s call |
\(i = 1 \dots n\) |
the fields (DOF ranges, in DOF order) |
|
|
\(D\) |
the diagonal scaling, entries \(\lvert a_{ii} \rvert^{-1/2}\) |
|
|
\(\hat A = D A D,\ \hat b = D b\) |
the scaled (equilibrated) system, \(x = D z\) |
|
formed for the components below a |
\(\hat A_{ij}\) |
its field blocks |
|
the field blocks of the |
\(M_i^{-1}\) |
one AMG cycle on \(\hat A_{ii}\) |
|
|
\(M^{-1}\) |
the block Gauss–Seidel sweep, \(v = M^{-1} r\) |
|
|
\(\hat A M^{-1} y = \hat b\) |
the outer Krylov solve, right-preconditioned |
|
|
\(\eta\) |
the requested relative tolerance of one solve |
|
|
Newton iteration k, in the nonlinear solver
assembled K (VIJ -> CSR) and residual R
|
| MPC condensation K <- T' K T + I_s both size preserving:
| Dirichlet elimination row, column i -> e_i DOF ordering unchanged
|
| setModel(model, dofManager) -- on a (re)build only
| -> field blocks, nodal dimensions, node coordinates
v
(A, b) ---> the composed linear solver
|
(1) Newton state, from ||b|| and the solver's own bookkeeping
||b|| > residualGrowthFactor * ||b||_prev ?
yes -> a new increment: eta = its maximum,
refresh (with rebuildOnNewIncrement)
no -> eta from Eisenstat-Walker forcing
refresh the hierarchies? first solve | block layout | size |
nnz differs from the rebuild's | new increment | setModel |
hot reload | previous solve stale or at the iteration cap
|
(2) equilibration Ahat = D A D, bhat = D b, D = diag(1/sqrt|A_ii|)
(no refresh: the D of the rebuild solve)
|
(3) field split diagonal blocks Ahat_ii + couplings Ahat_ij
|
(4) hierarchies refresh -> build one AMG hierarchy per field
(smoothed aggregation + Chebyshev,
per-field near-null-space)
otherwise -> reuse the standing ones
|
(5a) preconditioner M = block Gauss-Seidel sweep over the fields,
one AMG cycle per diagonal block
|
(5b) outer Krylov solve, preconditioned by M
krylovOnScaledSystem = true: on Ahat z = bhat, x = D z
krylovOnScaledSystem = false: on A x = b, with D M^-1 D as the preconditioner
|
(6) measure the TRUE residual ||A x - b|| / ||b||
> eta and true-residual restarts left (scaled system only)?
yes -> warm restart at a tighter tolerance, back to (5b)
no -> done
|
(7) record ||b||, eta and the outer iteration count for the next solve
(after a rebuild also nnz and the staleness reference count)
v
x ------> nonlinear solver
|
| MPC slave recovery x_s <- x_s + W x (master-interpolated part)
v
the Newton correction
In matrix form
With the fields numbered in DOF order, \(\hat A\) splits into blocks \(\hat A_{ij}\). The preconditioner of step (5a), \(v = M^{-1} r\), is one forward block Gauss–Seidel sweep with the per-field AMG cycles \(M_i^{-1}\) in place of exact block solves,
applied by forward substitution, \(v_i = M_i^{-1} \bigl( r_i - \sum_{j < i} \hat A_{ij} v_j \bigr)\).
With exact \(M_i = \hat A_{ii}\) and a block lower triangular \(\hat A\), it is \(\hat A^{-1}\).
The outer Krylov method (LGMRES, AMGCL’s implementation) solves the right-preconditioned system for its own
unknown \(y\),
whose residual \(\hat b - \hat A M^{-1} y = \hat b - \hat A z = D (b - A x)\) is the scaled one –
hence the true-residual check of step (6). With
krylovOnScaledSystem = false (the default, LGMRES(DiagonalScaling(...)) rather than
DiagonalScaling(LGMRES(...))) the outer method solves the unscaled
system instead, with the same \(M\) between two scalings,
and its own (right-preconditioned) residual \(b - A x\) is the true one. Steps (1)–(7) are the components of the tree, one section each below.
Diagonal equilibration: DiagonalScaling
A coupled multi-field system has a large diagonal dynamic range: elastic stiffnesses, a Helmholtz-like damage operator and penalty contact stiffnesses sit in the same matrix, and the Dirichlet-eliminated rows carry a unit diagonal among them. AMG’s strength-of-connection test is scale-sensitive, so every solve begins with a symmetric diagonal (Jacobi) equilibration
(in the code, the diagonal of \(D\) is the vector scaling), a two-sided diagonal scaling:
\(\hat A z = \hat b\) with \(x = D z\) is the same system in other variables, whatever
\(D\) is. That is what makes it possible to reuse an older \(D\) together with a reused
hierarchy (see below): a stale scaling or hierarchy costs outer iterations, and the accuracy of a solve
only if it drives that solve to the iteration cap.
The field split: OnField
Each OnField(name, M) acts on the dofs of one field, \(R^T M^{-1} R\). The matrix the tree is set up on
(SystemMatrix) splits
\(\hat A\) along the field ranges into diagonal blocks \(\hat A_{ii}\) (one per
field, the operators the AMG hierarchies are built for) and the off-diagonal coupling blocks
\(\hat A_{ij}\), which are needed on every solve regardless of whether the hierarchies are
rebuilt. All of these blocks, and the full-system operator used for the outer matrix-vector products
and the true-residual check, are handed to AMGCL’s OpenMP-threaded sparse matrix-vector kernels:
SciPy’s own sparse matvec is single-threaded C code irrespective of OMP_NUM_THREADS, and on a
solver whose inner loop is nothing but sparse matrix-vector products that is a substantial share of
the wall-clock. They are rebuilt on every call, because the values of \(A\) – and with them those
of \(D A D\) – change on every Newton iteration; a stable sparsity pattern would only save the
counting pass.
Per-field AMG hierarchies: AMG
Each diagonal block gets its own AMGCL hierarchy, applied on every visit of the sweep to its field
(once per outer Krylov iteration with symmetric = false; with symmetric = true the reverse sweep
revisits all fields but the last), built anew for every solve by default or, with RebuildPolicy(reusePreconditioner=True),
reused across Newton iterations while it still fits the Jacobian (see the reuse rules below). The defaults differ by nodal dimension (defaultFieldTree() in
the code): a 3D vector field (e.g. displacement) is stored in \(3\times3\) nodal blocks in
float, coarsened by smoothed aggregation with a strong-connection threshold of 0.01 and
estimate_spectral_radius (see Near-null-space and the Chebyshev smoother) from the translations of its equilibrated block (see
below), and smoothed by a degree-16 Chebyshev polynomial
(one post-smoothing sweep, no pre-smoothing) – and in the blockAMG() preset its hierarchy starts with the
rigid-body coarse level, a first coarse level of all its rigid-body modes, rotations included (see
below), or, for quadratic elements, with the corner-node level
(p-multigrid: TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve)));
a 2D vector field gets the same coarsening and
smoothing but in \(2\times2\) blocks, in double (this wrapper has no float 2x2 backend); a vector
field of any other dimension gets the same recipe on the plain scalar backend – the only backend that
takes the full rigid-body near-null-space (see below), which keeps that path available for a field
shape the block backend was never tuned for; a scalar field takes the plain scalar tree (smoothed
aggregation, Chebyshev at AMGCL’s own defaults). Any field can be given another tree in a custom solver
script, OnField(name, AMG(amgclTree=...)); the default trees are listed in Default field trees.
The tree is AMGCL’s parameter tree, plus one key of EdelweissFE’s own, spectralRadiusSteps, which its Chebyshev
smoother and its smoothed aggregation read (Near-null-space and the Chebyshev smoother); a Chebyshev relax rejects keys it does
not know. Three keyword arguments of AMG (the first two also of TwoLevel) select
the backend the hierarchy is built on; left unset, backendOf() chooses
them from the field:
backendPrecision"double"or"float"; by default"float"on \(3\times3\) blocks and"double"otherwise. A single-precision backend should in principle roughly halve the hierarchy’s memory traffic, but measured across a full set of captured reference systems it inflated outer iteration counts and netted only a small gain in aggregate – plausibly because the Chebyshev smoother’s own spectral-radius estimate loses accuracy infloat32, degrading smoother quality enough to eat the bandwidth saving. That holds for the scalar backend, which builds the whole hierarchy in single precision. WithbackendBlockSize = 3,"float"instead sets the hierarchy up in double – coarsening, Galerkin products, every spectral-radius estimate, the coarsest direct solve – and stores only the level matrices the cycle reads in float, with double vectors and double accumulation.backendBlockSizeDefaults to the field’s own nodal dimension, for a block-valued backend operating on \(B \times B\) nodal blocks:
3for a 3D vector field,2for a 2D one,1(scalar) otherwise.3for a displacement field was measured faster per outer iteration (less index traffic) and validated inside the block Gauss–Seidel sweep on captured and live runs. Note that a block hierarchy has the identity block as its near-null-space, unlessnearNullspacesays otherwise: constants per component in the equilibrated variables. Those are \(D\)-weighted translations, \(D\) times a constant, not the translations \(D^{-1} t\) of \(\hat A\) (see below); the two agree only where the diagonal is uniform, which it is not on quadratic elements (corner and midside nodes) or across a damage band. A block hierarchy also has no rotational rigid-body modes, which the scalar backend carries (described next);nearNullspace(on a block backend the translations, by default) gives it the translations, and the rigid-body coarse level (TwoLevel(coarseSpace=RigidBodyModesPerAggregate()), the default for the 3D displacement field) a first coarse level of all rigid-body modes (see below).nearNullspaceWhat smoothed aggregation starts from:
"rigidBodyModes"(the default),"translations"orNone. The default takes the field’s rigid-body modes as far as the backend holds them. On a block backend (backendBlockSize2 or 3), which takes the translations at most, that is the translations \(D^{-1} t\) of \(\hat A_{ii}\), instead of AMGCL’s identity blocks – see the near-null-space of a block backend below; the aggregation, the prolongation smoothing and the smoother stay AMGCL’s. On the scalar backend it is all rigid-body modes of a displacement field (translations and rotations, from the node coordinates; described next) and the translations of any other vector field, soAMG(backendBlockSize=1)needs nothing else."translations"takes the translations on either backend;Nonegives AMGCL no near-null-space (its identity blocks, or its constant). A near-null-space on a block backend needs smoothed aggregation, which is checked: a tree with another coarsening needsnearNullspace=None.
Near-null-space and the Chebyshev smoother
Smoothed aggregation and the Chebyshev smoother have to be understood together, because they play deliberately complementary roles and only work if they agree on what “smooth error” means.
The smoother damps error components whose eigenvalues lie inside an estimated window \([\,\mathrm{lower},\ \mathrm{higher}\,]\cdot\rho\) of the operator’s spectral radius \(\rho\). With the vector-field default
lower = 0.01, eigenvalues below 1% of the spectral radius are by construction not targeted by the smoother at all.Everything below that window – including the operator’s near-null-space, the directions in which it has almost no stiffness – is left to the coarse-grid correction, which removes it only to the extent that smoothed aggregation’s prolongation operator can represent it. Smoothed aggregation builds that prolongation from a supplied near-null-space basis.
The Chebyshev smoother of degree \(k\) is \(k\) steps of the Chebyshev iteration on \(\hat A_{ii} u = r_i\) for the interval \([\ell, h] = [\mathrm{lower}, \mathrm{higher}]\cdot\rho\), with centre \(d = (h + \ell)/2\) and half-width \(c = (h - \ell)/2\):
with \(\alpha_1 = 2d / (2d^2 - c^2)\) and \(\alpha_k = (d - \alpha_{k-1} c^2 / 4)^{-1}\)
(AMGCL’s relaxation::chebyshev). It damps every eigencomponent in \([\ell, h]\) by the
Chebyshev polynomial’s bound, and each step costs one matrix-vector product with the block: degree
16 with npre = 0, npost = 1 (the vector-field default) costs 16 products per cycle on the
finest level, degree 5 with npre = npost = 1 (the former vector-field default) 10.
For 3D linear elasticity the physical near-null-space is the six-dimensional space of rigid-body
motions: three translations and three infinitesimal rotations, each of which costs zero elastic
strain energy. Of the whole, constrained operator they are not eigenvectors – the Dirichlet rows hold
them – but they are locally near-null: restricted to an aggregate away from the constrained boundary,
a rigid-body motion is almost free of energy, and that local property is what smoothed aggregation
needs. rigidBodyNullspace() builds that full basis from the
field’s node coordinates, about the field’s own coordinate centroid, node-major to match the DOF
vector’s layout:
with \((x,y,z)\) the coordinates relative to the centroid. A 2D field gets two translations plus
the single in-plane rotation \((-y,x)\); a field of any other nodal dimension has no rotational
rigid-body mode and takes the translations alone
(translationNullspace()). A scalar field, e.g. nonlocal
damage, has no rigid-body rotation either and takes AMGCL’s default constant.
Because the hierarchy is built for the equilibrated block, the basis has to be transformed consistently: the near-null-space of \(D A D\) is \(D^{-1}\) times that of \(A\), which is why both functions divide the raw modes by the field’s own slice of the scaling vector.
Supplying only the translations – the common simplification, since rotations need node coordinates rather than just the DOF-block structure – leaves the coarse grid unable to represent rotational rigid-body error exactly, and by the construction above the Chebyshev smoother is not targeting that error class either, so it has nowhere efficient to go.
The full basis was measured to need fewer outer iterations than translations alone, and to be robust to the thread count (The full rigid-body near-null-space).
The spectral radius is the smoother’s second, independent subtlety. The Chebyshev polynomial damps the eigenvalues in
[lower, higher] times an estimate of \(\rho\); an estimate too far below \(\rho\) leaves the largest
eigenvalues outside that interval, where the polynomial amplifies them. Every spectral radius of EdelweissFE – of
the Chebyshev smoother of every AMG level, of the fine level of TwoLevel, and of \(D^{-1} A\) for the weight of
the prolongation smoothing – is estimated one way (edelweissfe/linsolve/amgcl/amgcl-spectral-radius.hpp):
Arnoldi steps, not the power iteration. The power iteration converges with the ratio of the two largest eigenvalues, close to one for finite element operators, so a short power iteration can stop far below \(\rho\). The largest Ritz value \(\theta\) of the Arnoldi process from the same number of products converges much faster: a few steps reach what takes the power iteration many. Our operators are not symmetric (contact, softening, non-associated plasticity, the condensation of multi-point constraints), so the symmetric Lanczos process has no theory for them; Arnoldi holds for any matrix, and its orthogonalization costs \(k^2/2\) vector operations for \(k\) steps, little against the \(k\) sparse products. On the captured systems, markedly nonsymmetric ones among them, the two give the same estimates (The spectral radius estimate of the Chebyshev smoother).
A start vector seeded per row. AMGCL seeds its start vector per OpenMP thread, so its estimate – and with it the outer iterations, on a hard system by far – depends on the thread count (The spectral radius estimate of the Chebyshev smoother). The start vector is computed from the row index alone, zero on the Dirichlet rows of a two-level cycle, and the inner products sum over a fixed partition of the rows in order: the estimate is the same, bit for bit, at any thread count.
The Ritz value of largest modulus plus its residual norm, \(|\theta| + r\) with \(r = |h_{k+1,k} y_k|\) for the unit eigenvector \(y\) of the small \(k \times k\) Arnoldi matrix, which Eigen’s dense eigensolver computes. The modulus, because the spectral radius is the largest modulus of the eigenvalues: a dominant eigenvalue that is negative or complex, as an indefinite or strongly nonsymmetric level (contact, softening) may have, is not missed, as the largest real part would miss it. For a symmetric matrix there is an eigenvalue within \(r\) of \(\theta\), so \(|\theta| + r \geq \rho\) once \(\theta\) is the Ritz value of the dominant eigenvalue. For a nonsymmetric matrix \((\theta, x)\) is an exact eigenpair of a matrix within \(r\) of \(A\), which bounds the error only together with the conditioning of the eigenvectors: \(r\) is an indicator there, not a bound, and the smoother keeps its own margin,
higher1.1 – the default of every Chebyshev smoother in EdelweissFE, a tree that leaves it out included. AMGCL’s own default is 1.0, which sits on its Gershgorin bound, never below \(\rho\); the Arnoldi estimate can be, slightly, and an eigenvalue above the interval is amplified, not damped.A fixed number of steps,
spectralRadiusSteps, 10 by default: on the captured systems the estimate then lies slightly above \(\rho\), well within the smoother’s margin (The spectral radius estimate of the Chebyshev smoother).Its memory. Arnoldi keeps its \(k\) basis vectors while it runs, \(k\) doubles per DOF (80 bytes per DOF for the 10 default steps); the power iteration needed two. Against the operator itself (about 80 stored entries per DOF of a 3D field) that is small, and it is freed before the next level of a hierarchy is built.
The key spectralRadiusSteps is EdelweissFE’s own, not AMGCL’s: it sits in a tree’s relax for the smoother and
in its coarsening for the prolongation smoothing (with estimate_spectral_radius). AMGCL’s power_iters, which
it replaces, is refused, also power_iters 0, AMGCL’s Gershgorin bound: there is one estimate. The Chebyshev
relaxation and the smoothed aggregation that use it are EdelweissFE’s own (amgcl-chebyshev.hpp), step for step
AMGCL’s otherwise, with one Chebyshev polynomial for every Chebyshev smoother; every other relaxation and coarsening is
AMGCL’s.
The near-null-space of a block backend
A block backend (backendBlockSize 2 or 3, the default for 2D and 3D vector fields) stores a field’s
block in \(B \times B\) nodal blocks. Without a near-null-space, its tentative prolongation gives every
node \(i\) of an aggregate \(a\) the identity block, \(P^{\mathrm{tent}}_{ia} = I_B\), so the coarse
space holds the constants per component of the variables the hierarchy is built on. For the equilibrated block \(\hat A_{ii} = D A_{ii} D\) these
are \(z = D^{-1} x\), and a constant \(z\) is the displacement \(x = D\,c\): a
\(D\)-weighted translation. The translations of \(\hat A_{ii}\) are \(D^{-1} t\) instead, the
same transformation the scalar backend applies to its rigid-body modes (see above).
nearNullspace (the default "rigidBodyModes" and "translations") builds the hierarchy on a block backend
from those translations
(amgcl-block-translations.hpp of the AMGCL wrapper edelweissfe.linsolve.amgcl.amgcl). Per node
they form the diagonal block \(W_i = \operatorname{diag}(w_{i1}, \dots, w_{iB})\),
\(w_{ik} = 1/s_{ik}\) with \(s\) the equilibration
(translationNullspace()). Smoothed aggregation orthonormalizes its
near-null-space on every aggregate by a QR factorization of the stacked blocks. Column \(k\) of that
stack is nonzero only in the rows of component \(k\), so the columns are orthogonal and the
factorization is a normalization per component:
So \((P^{\mathrm{tent}})^\mathsf{T} P^{\mathrm{tent}} = I\) on every aggregate, and
\(P^{\mathrm{tent}} R = W\) reproduces the translations exactly; \(R\) is the near-null-space of
the coarse level, which repeats the construction. Every block stays a diagonal \(B \times B\) block,
so the hierarchy keeps AMGCL’s block value type. The rotations do not fit that way (six modes on
\(3 \times 3\) blocks would need \(6 \times 6\) coarse blocks): this hierarchy takes the translations only
and refuses any other near-null-space, and the rotations stay with the scalar backend – or with the rigid-body coarse
level, the default for the 3D displacement field, which builds the first coarse level of a block field itself (its
prolongation in \(B \times M\) blocks, the Galerkin product in \(M \times M\) blocks) and hands the unblocked
coarse operator to a scalar AMGCL hierarchy with aggr.block_size \(M\) (see
Rigid-body coarse level: TwoLevel(coarseSpace=RigidBodyModesPerAggregate()), which also says why AMGCL’s as_scalar coarsening, the route of AMGCL’s
runtime coarsening for a near-null-space on a block backend, does not fit).
Everything else is AMGCL’s smoothed aggregation, unchanged: the aggregation on \(\hat A_{ii}\), the
prolongation smoothing and the Galerkin product. Unit weights \(w = 1\) give AMGCL’s own hierarchy back,
up to a scaling of the coarse variables.
A hierarchy built on \(A_{ii}\) itself, whose constant blocks are its translations \(t\), gives the
same coarse space: with the same aggregates, it and the hierarchy of \(\hat A_{ii}\) on \(D^{-1} t\)
are one hierarchy in other variables (a similarity transform),
\(M_{\hat A}^{-1} = D^{-1} M_A^{-1} D^{-1}\), provided the smoother is invariant under the scaling, as a
Jacobi-scaled Chebyshev smoother ("relax": {"scale": true}) is and an unscaled one is not. Only the
aggregates can differ: for blocks, AMGCL’s strength test compares traces of block products,
\(\operatorname{tr}(\varepsilon^2 A_{ii} A_{jj}) < \operatorname{tr}(A_{ij} A_{ij})\), which a symmetric
diagonal scaling leaves unchanged only where \(D\) is uniform within each node.
Measured offline and live, the translations need fewer outer iterations on most models; they help less once softening makes the operator indefinite (The translations of a block backend).
The block Gauss–Seidel sweep: Multiplicative
The preconditioner \(v = M^{-1} r\) applied at each outer Krylov iteration is a block Gauss–Seidel sweep over the fields. Starting from \(v_j = 0\) for all fields, each field \(i\) in turn forms its local residual against the current parts of \(v\) of all fields and corrects its own part by one cycle of its own AMG hierarchy:
where \(M_i^{-1}\) denotes one AMG cycle on \(\hat A_{ii}\). On a field’s first visit
\(v_i = 0\), and this is the plain \(v_i = M_i^{-1}(r_i - \sum_{j \neq i} \hat A_{ij} v_j)\);
the correction form matters on revisits, which an inexact \(M_i^{-1}\) would otherwise throw
away. With symmetric = true (default false) each forward sweep is followed by a reverse-order sweep
of the same correction form, \(i = n-1, \dots, 1\) – starting one field before the one just
updated. With exact field solves that field has nothing left to correct; with AMG cycles it has, and
skipping it saves one cycle of that field. With exact \(M_i = \hat A_{ii}\) this is the classic block
symmetric Gauss–Seidel preconditioner,
with inexact \(M_i\), a reverse sweep in correction form over all fields would give \((B + U)^{-1} (2B - \hat D) (B + L)^{-1}\), \(\hat D = \operatorname{blockdiag}(\hat A_{ii})\), the classic operator for \(B = \hat D\); skipping the last field changes its last block row.
sweeps repeats the (possibly reversed) sweep. Both are arguments of the blockAMG() preset (data lines of
type=blockamg); in a tree, Multiplicative(gaussSeidelOrder([OnField(...), ...], sweeps, symmetric))
(gaussSeidelOrder()).
Neither is needed for the outer solve: right-preconditioned GMRES/LGMRES only needs \(M^{-1}\)
to be a fixed linear operator, not a symmetric one – the measured configurations use
sweeps = 1, symmetric = false.
What that one forward sweep does depends on how the fields couple. On the gradient-enhanced damage model of the reference runs the coupling is one-sided: \(K_{ud}\) (displacement rows, damage columns) is a noticeable fraction of \(K_{uu}\), while \(K_{du}\) is many orders of magnitude smaller than \(K_{dd}\). With the displacement field first in DOF order, the only coupling product the forward sweep forms is \(\hat A_{du} v_u\), which is numerically zero there: on that model the default sweep acts as block Jacobi, and still pays that product on every outer iteration. Sweeping the damage field first would instead make it the exact block back substitution of the (then block upper triangular) system, up to the inexact field solves. Measured, the order made no difference to the outer iterations on the reference models (The order of the block Gauss–Seidel sweep).
This is deliberately only an approximate inverse of the coupled operator – which is exactly why the true-residual bookkeeping below exists.
The outer Krylov solve: LGMRES
The outer solve, LGMRES, is AMGCL’s own native
amgcl::solver::lgmres, driven through
PyAMGCLKrylovSolver; the block sweep is passed straight in
as the preconditioner callback. (FGMRES is the alternative for a
variable preconditioner, one with a Krylov solve inside, which LGMRES refuses.) A former fallback on
scipy.sparse.linalg.gmres() was removed once the native path was gated live against
it at several thread counts and found to produce the same Newton trajectory. Orchestrating GMRES from
Python – the Arnoldi steps, orthogonalization and restart bookkeeping – was measured as a
substantial share of this solver’s wall-clock on a reference model when run that way, and, being serial Python,
anti-scaled past one NUMA node on multi-socket hardware, which the native implementation does not.
LGMRES(..., maxIterations) (the data line maxIterations) bounds all Arnoldi steps of a solve, across all
its internal restarts.
LGMRES augments each restart with a few vectors carried over from the previous restart cycle of the same
solve. LGMRES(..., augmentationVectors) sets how many, and a restart cycle is restart + augmentationVectors
iterations, every one of them a preconditioner application: restart Krylov steps plus one step per augmentation
vector (amgcl/solver/lgmres.hpp). The first cycle of a solve has no augmentation vectors yet and is plain
GMRES(restart + augmentationVectors). Every solve starts from an empty augmentation subspace. Keeping it alive
across solves (“recycling”) sounds attractive – a hard solve’s converged subspace as a warm start for
the next – but measured on captured systems it contributed nothing (bit-identical iteration counts), and
on a live run it was actively harmful: a solve that had already struggled left a poorly conditioned
subspace which compounded into a following solve needing far more iterations than a fresh one, and in
one run an outright NaN. It was therefore removed.
Adaptive outer tolerance: EisenstatWalkerTolerance
Most Newton iterates do not need their linear solve solved tightly, and loosening the requested
tolerance for those cuts outer iterations substantially. By default LGMRES therefore chooses
its own outer tolerance per solve by Eisenstat–Walker forcing (“choice 2”; EisenstatWalkerTolerance(minimum, maximum,
gamma, alpha), EisenstatWalkerTolerance, the default tolerance of
LGMRES), from the ratio of successive residual norms:
clamped to \([\eta_\mathrm{min}, \eta_\mathrm{max}]\) (minimum, maximum). The classic safeguard against
over-loosening is applied: a large tightening step is only trusted if the previous tolerance was
already small, i.e. \(\eta_k\) is raised back to \(\gamma \eta_{k-1}^{\alpha}\) whenever that
quantity exceeds 0.1. At the classic parameters that needs \(\eta_{k-1} > 0.26\), so the safeguard
never acts with the default \(\eta_\mathrm{max}\) (well below that). With no meaningful history to compare against – the first solve, a residual of
exactly zero, or the solve right after a detected new increment, where the ratio across the
jump says nothing about Newton convergence – the tolerance falls back to \(\eta_\mathrm{max}\).
A solve is treated as the first of a new increment when its residual norm exceeds
residualGrowthFactor times the previous solve’s. That single test drives both the tolerance reset
and a hierarchy refresh. It is a guess, and it has a blind spot: a cutback after a diverged attempt
usually starts from a residual smaller than the last, large one of that attempt, so the retry is not
recognized – its first solve, typically the hardest, gets a tolerance from the ratio to the diverged
residual (clamped to \(\eta_\mathrm{min}\)) and, unless the sparsity pattern changed, the
hierarchies of the discarded attempt. Only a cutback after an attempt whose residual had decreased is
detected. (A notification from the nonlinear solver at the start of every increment would replace the
guess; that is planned.)
A fixed relative tolerance – the data line tolerance, LGMRES(..., tolerance=FixedTolerance(eta)) – switches
the scheme off entirely. Adaptive forcing is safe to use as the default only because the requested tolerance is
enforced on the true residual, which is the subject of the next section: without that enforcement,
a loosened request combined with an imperfect preconditioner can quietly let the true residual run
looser still, which does change the Newton path.
True-residual enforcement: LGMRES on the scaled system
On the scaled system (krylovOnScaledSystem = true; not the default), the outer Krylov method’s own
stopping test is not on the true residual but on the residual of the scaled (equilibrated) system it
iterates on, \(\lVert \hat b - \hat A z \rVert / \lVert \hat b
\rVert\): AMGCL’s LGMRES is right-preconditioned and tests it directly. The
scaled residual can differ substantially from the true relative residual (stiff rows are down-weighted
by \(D\)), so “converged” can mean a true relative residual well above the requested \(\eta\);
measured on a reference model, one solve’s achieved true residual was found far looser than what had
been asked for.
Every solve on the scaled system therefore measures the true relative residual \(\lVert A x - b \rVert / \lVert b
\rVert\) explicitly (computed as \(D^{-1}(\hat A z - \hat b)\), which is exact and rides on the
threaded operator already built, rather than paying a second unscaled matrix-vector product). If it
still exceeds \(\eta\), the solve is restarted on the true residual: a warm restart from the current iterate
at a tighter requested tolerance, up to trueResidualRestarts times (the data line, LGMRES(...,
trueResidualRestarts)). By
default the tolerance
is tightened geometrically by a fixed factor of 0.01 per round – deliberately expressed purely in
the units of the requested relative tolerance, which the Krylov solver interprets consistently by
construction, rather than being derived from the solver’s internal residual bookkeeping. Setting
trueResidualRestarts = 0 disables enforcement: the outer solve’s own test decides.
On the unscaled system (krylovOnScaledSystem = false) none of this is needed: the residual the
outer solve converges is the true one, so a pass that misses \(\eta\) has stopped at its iteration
cap, and a true-residual restart would only extend the cap (maxIterations is the knob for more iterations). What
it trades is the norm itself. \(\lVert b - A x \rVert\) sums the residuals of all
fields in their own units, while the scaled residual \(D (b - A x)\) weighs every row by its own
diagonal. In a system whose fields differ widely in residual scale – a displacement field in newtons
next to a damage field with entries of order one, say – a field with small residual entries can thus
be left less resolved at the same \(\eta\) on the unscaled system than on the scaled one. On the
reference models the unscaled system was the faster of the two, with no more Newton iterations (see
Why the defaults are what they are), but that is a measurement on those models, not a guarantee.
Gap-compensated tolerance: ResidualGapCompensation
Because the system is equilibrated, the outer Krylov method minimises the scaled residual \(\lVert \hat b - \hat A z \rVert\) while acceptance is judged on the true relative residual. The ratio between the two – the gap – has been measured to vary widely, from close to one on some models to considerably larger on a contact model with multi-point constraints (stiff rows are down-weighted by the equilibration), and it is why almost every solve needs a true-residual restart on the scaled system (the former default): on a reference gradient-enhanced damage model, nearly every solve ran one, i.e. effectively every linear solve was performed twice. A true-residual restart is disproportionately expensive, because restarting discards most of the accumulated Krylov subspace.
With gap compensation (the data line residualGapCompensation=true; in a tree
LGMRES(..., residualGapCompensation=ResidualGapCompensation(safetyFactor, maxFactor))) the gap is measured on each pass (as the
achieved true
relative residual over the achieved scaled relative residual of the same solution – not over the
tolerance the pass was given, which a pass stopped at the iteration cap does not reach; that version
inflated the gap by the unreached factor and made every following first pass ask for that much more
accuracy), smoothed across solves by an exponential moving
average, floored at 1 and capped at maxFactor, and used to pre-compensate the first pass: it is asked for
safetyFactor\(\cdot\,\eta\,/\,\text{gap}\) so that it lands on target in one pass. Any
true-residual restart that is still needed then tightens by the measured gap rather than by the fixed factor
(never looser, and always at least a factor of two per round so the loop still terminates).
Gap compensation is opt-in, and need not be the default: the default already runs the outer solve on the unscaled
system (krylovOnScaledSystem = false), where there is no gap to compensate at all. Why it is not the default on
the scaled system either – the uncompensated solve over-solves, and Newton relies on it – is measured in
Gap compensation.
Reuse across Newton iterations: RebuildPolicy
Building the per-field hierarchies is a large fraction of a solve, and it is avoidable when the
Jacobian has barely moved since the last one. LGMRES can therefore keep the hierarchies below it (and the
equilibration they were built for) standing across solves, RebuildPolicy(reusePreconditioner=True) (the data line
reusePreconditioner=true of type=blockamg); when they are rebuilt is the outer solve’s
RebuildPolicy,
LGMRES(..., rebuild=RebuildPolicy(stalenessFactor, rebuildOnNewIncrement, reusePreconditioner)).
By default it does not (reusePreconditioner=False): every solve builds its own preconditioner. A reused hierarchy
was built for an earlier Jacobian, and the outer iterations it costs grow faster than its build saves: rebuilding for
every solve needed far fewer outer iterations per solve and was faster (Rebuilding the preconditioner for every solve).
The outer Krylov solve always operates on the current matrix (equilibrated with the scaling the standing hierarchies were built for, which is the same system whatever the scaling); only the preconditioner may be stale. A stale preconditioner costs outer iterations and, if it drives a solve to the iteration cap, the accuracy of that one Newton step; the next solve is then rebuilt (see below).
The hierarchies are rebuilt when
reusePreconditioneris off: then every solve builds its own, and none of the triggers below matters,there is nothing to reuse yet (the first solve),
the field-block layout or the system size changed (e.g. an adaptive-refinement event),
the number of stored nonzeros changed – a free \(O(1)\) proxy for “the sparsity pattern changed”. It is not exact (two different patterns could coincidentally share an
nnz) but it errs safely, since an unnecessary refresh costs time and never correctness. A pattern change is not a graceful degradation: measured on one transition, a hierarchy built for the previous pattern needed far more outer iterations than a fresh one, an outright wall-clock regression,a residual jump marks a new increment (see
residualGrowthFactorand its blind spot above) andrebuildOnNewIncrementis set (the default), orthe previous solve’s outer iteration count exceeded both
stalenessFactortimes and 15 more than the count of the solve the hierarchies were built for – a growing count being the signal that the reused hierarchies are drifting away from the current Jacobian. The reference is that first solve, not the previous one, so a slow drift across several solves is caught; the absolute margin keeps a factor alone from firing on noise at small counts, orthe previous solve hit the outer iteration cap. A missed tolerance alone does not trigger a rebuild: on the scaled system, with true-residual restarts off, it typically comes from the gap between the scaled residual the outer solve converges and the true one (stiff rows down-weighted by the equilibration, e.g. on a model with multi-point constraints), which a fresh hierarchy cannot close – measured slower on such a model when it did trigger rebuilds; on the unscaled system (the default) a miss without a cap hit happens only at a rounding tie, since the outer solve then converges the true residual itself. A cap hit asks for a rebuild even when the capped solve had just rebuilt its hierarchies: not rebuilding again was measured neutral on a live run (the same Newton path, at a comparable cost), so the simpler rule stays.
On a model with contact or tie constraints the active constraint set can change between Newton
iterations, and then the sparsity pattern – and with it the nonzero count – changes and the
hierarchies are rebuilt. How often that happens depends on the model: the multi-point condensation
would otherwise drop exact zeros (pruneCondensedMatrixZeros, a nonlinear-solver option), which
changes the count on almost every Newton iteration even when the constraint set does not; left unset,
that option’s “auto” default leaves pruning off for every tree without a direct solver, and so whenever the
hierarchies are reused (needsPrunedSystemMatrix), so reuse pays on
such models without any extra configuration (on a contact model with many multi-point constraint
slaves, a configuration relying on it ran noticeably faster than with pruning on).
The native LGMRES instance is deliberately not torn down on a hierarchy refresh – only on an
actual system-size change, which invalidates AMGCL’s preallocated scratch vectors. (A hot reload of a custom
solver script builds a new LGMRES, and with it a new instance.)
p-multigrid: TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve))
For a field discretized with quadratic elements, a two-level cycle whose coarse level is the field’s linear part
replaces the field’s single-level hierarchy, OnField(name, TwoLevel(coarseSpace=CornerNodeInterpolation())): the
blockAMG() preset gives it to a quadratic 3D displacement field (cornerNodeInterpolation=false switches it
off), a custom solver script to any quadratic field (CornerNodeTwoLevelCycle):
npre Chebyshev sweeps on the field’s own scaled block \(\hat A = D A D\), restriction of the residual through
\(P^\mathsf{T}\), the coarse solve on the Galerkin-projected coarse operator \(A_c = P^\mathsf{T} \hat A P\),
prolongation through \(P\), and npost more Chebyshev sweeps. The coarse space is purely topological – the
mesh’s corner nodes – so no re-discretization is involved. The coarse solve is a component; its default,
DiagonalScaling(AMG()), works on any field. For a 3D displacement field the blockAMG() preset passes the faster
rigid-body coarse level, three levels built from components:
# each level its own tree, here each a copy of the corner-node two-level cycle's (Chebyshev, degree 8)
coarse = DiagonalScaling(
TwoLevel(coarseSpace=RigidBodyModesPerAggregate(), amgclTree=defaultCornerNodeTwoLevelTree())
)
displacement = OnField(
"displacement",
TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve=coarse), amgclTree=defaultCornerNodeTwoLevelTree()),
)
The corner-node two-level cycle smooths with a weaker polynomial. On a 3D displacement field given no tree, its fine
level, the quadratic field, takes a Chebyshev polynomial of degree 8 on the interval from 0.03 to 1.1 times the spectral
radius (defaultCornerNodeTwoLevelTree()), not the degree 16 on 0.01 to 1.1
of the default field trees. The coarse solve has a tree of its own; the blockAMG() preset gives the rigid-body
coarse level of the corner-node field the same. A sweep of degree 8 costs half the matrix products of one of degree 16;
the corner-node two-level cycle then needs more outer iterations, but fewer than twice as many, so it is faster.
AMG() and the rigid-body coarse level alone gain nothing from degree 8 – their iterations grow about as much as
each one gets cheaper – and keep degree 16 (The smoother of the corner-node two-level cycle).
The prolongation interpolates the scaled variables. The P1 interpolant \(P_{geo}\) keeps a corner node’s displacement and gives a midside node the mean of its two edge-endpoint corners. But the cycle works on \(\hat A = D A D\), whose variables are \(\hat u = D^{-1} u\), not \(u\). Interpolating \(\hat u\) with the weights \(\tfrac12, \tfrac12\) is interpolating something else than the displacement: a scaled translation \(D^{-1} t\) – the smoothest error there is, which the coarse level exists to remove – is reproduced at a midside node only where its diagonal equals its corners’; on captured systems the plain weights missed it at most free midside dofs, and by far where the diagonals along an edge differ much (The corner-node two-level cycle with scaled weights and a composed coarse solve). The interpolation of \(u\), written in the scaled variables, is
for a midside dof \(m\) on the edge \((a, b)\), \(s\) the diagonal of \(D\) and \(D_c\) its restriction to the corners. It reproduces \(D^{-1} t\) to rounding, and the coarse operator \(A_c = P^\mathsf{T} \hat A P = D_c P_{geo}^\mathsf{T} A P_{geo} D_c\) is the P1 operator of the physical problem, in the corners’ scaled variables. There is one way to do it, so there is no switch back to the plain weights; without a scaling above the cycle (\(D = I\)), \(P = P_{geo}\). The scaled weights need far fewer outer iterations than the plain ones (The corner-node two-level cycle with scaled weights and a composed coarse solve).
The coarse problem is a field of its own. \(A_c\) is handed to the component coarseSolve with the context of
the corner-node field (ofCoarseNodes()): the same field name
and dimension on the corner nodes, their coordinates, \(D_c\) as its scaling, and an identity row for each Dirichlet
corner dof (a zero column of \(P\)). A midside node with a free dof whose two edge corners are both constrained in
that component (a symmetry plane and a tie, issue #223) would get an empty row of \(P\), no coarse correction at
all: it joins the coarse nodes, as a node the midside corners cannot classify becomes a corner. This is decided at the
build, where the Dirichlet rows are known. So every component works on it unchanged: RigidBodyModesPerAggregate
builds the rigid-body modes of the corners, AMG(backendBlockSize=1) takes them as its near null-space,
DirectSolver() solves exactly. It is called back from the C++ cycle once per application.
Why the coarse solve scales the coarse operator once more. Where a midside diagonal is tiny against its corners’,
the weights \(s_a / s_m\) are large, and \(A_c\) inherits entries of their square: a smoother on \(A_c\)
would work with a spectral radius far above that of a unit-diagonal operator. DiagonalScaling around the coarse
solve rescales \(A_c\) to a unit diagonal, \(D_2 A_c D_2\) with
\(D_2 = |\operatorname{diag} A_c|^{-1/2}\); on a model with such diagonals, it saves outer iterations and time. Of
the coarse solves measured on a 3D displacement field,
DiagonalScaling(TwoLevel(coarseSpace=RigidBodyModesPerAggregate())) was the fastest, ahead of an AMG hierarchy
with the rigid-body modes of the corners, DiagonalScaling(AMG(backendBlockSize=1))
(The corner-node two-level cycle with scaled weights and a composed coarse solve).
Its preconditioner is rebuilt exactly when the single-level hierarchies are (see the reuse rules above): while only
the values change, the cycle and its coarse solve are both kept. Enable it per field with OnField(name,
TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve), amgclTree=...)); the topology map is computed lazily from the model and
recomputed whenever the model changes (setModel, e.g. after adaptive refinement). There is no fallback to the
single-level default: on linear elements every node is a corner, so the coarse operator is the field’s block itself and
the cycle gains nothing, and an element topology the map does not know raises. A node that two elements classify
inconsistently (seen at non-conforming refinement interfaces) is made a corner; the log says once per model how many.
The whole cycle, and every stage of its build, runs in C++ and is OpenMP-threaded. The Dirichlet
DOFs of the field’s block – rows whose only nonzero is the diagonal, which includes the identity
slave rows of the multi-point-constraint condensation – are decoupled exactly rather than removed:
the fine right-hand side and iterate are kept zero on them, which a Dirichlet row
preserves under every residual and every smoothing step, and \(P\) has no rows for them and zero columns
for the Dirichlet corner dofs, which are identity rows of \(A_c\). This is the free-submatrix cycle up to floating-point summation order, without extracting the
free submatrix, and it allows TwoLevel(..., backendBlockSize=3): the fine operator is stored in node-wise
\(3\times 3\) blocks even where a node is only partially constrained, which roughly halves the
memory traffic of the fine smoother and residual – the dominant cost of the apply.
backendPrecision="float" (the default on \(3\times 3\) blocks) halves the bytes per stored entry once more.
Only the stored entries
are rounded: vectors and every sum stay double (a float matrix times a double vector), and the
spectral radius is estimated on the double operator before it is rounded – unlike
AMG(backendPrecision="float") on the scalar backend, whose whole setup, the estimate included, runs in single
precision.
The fine Chebyshev smoother is the tree’s relax (degree, lower, higher and the spectral radius estimate,
Near-null-space and the Chebyshev smoother), applied npre times before and npost times after the coarse-grid correction – the
same spelling as for the rigid-body coarse level. It is the Chebyshev polynomial of every Chebyshev smoother, with the
bounds of its tree (on a 3D displacement field by default lower 0.03, higher 1.1 times the spectral radius;
the defaults of relax are 1/30 and 1.1, AMGCL’s 1/30 and 1 but the margin, Near-null-space and the Chebyshev smoother); its
estimate runs on the free rows (the start vector is zero on the Dirichlet rows) and is built anew with the cycle.
The same cycle on the free submatrix, written out in Python: the free submatrix, the scaled \(P\), an exact coarse solve and the Chebyshev polynomial. A test checks the C++ cycle against it, to rounding, for several numbers of sweeps and degrees:
def _geometricProlongation(isCorner, edgeEndpoints, nDim):
"""P_geo: a corner node keeps its value, a midside node takes the mean of its two edge-endpoint corners."""
nNodes = len(isCorner)
cornerIndex = np.cumsum(isCorner) - 1
rows, cols, vals = [], [], []
for node in range(nNodes):
corners = [node] if isCorner[node] else edgeEndpoints[node]
for corner in corners:
rows, cols, vals = rows + [node], cols + [cornerIndex[corner]], vals + [1.0 / len(corners)]
P_node = sp.csr_matrix((vals, (rows, cols)), shape=(nNodes, int(np.count_nonzero(isCorner))))
return sp.kron(P_node, sp.identity(nDim), format="csr")
def _cornerDofs(isCorner, nDim):
cornerNodes = np.nonzero(isCorner)[0]
return (cornerNodes[:, None] * nDim + np.arange(nDim)).ravel()
def _referenceCycle(As, scaling, isCorner, edgeEndpoints, spectralRadius, nuPre, nuPost, degree, r):
"""A pure-Python cycle on the extracted free submatrix, with P = D^-1 P_geo D_c and an exact coarse solve."""
n, nNodes = As.shape[0], len(isCorner)
nDim = n // nNodes
isDirichlet = np.diff(As.indptr) == 1
free = np.nonzero(~isDirichlet)[0]
A_free = As[free, :][:, free].tocsr()
cornerDofs = _cornerDofs(isCorner, nDim)
P = sp.diags(1.0 / scaling) @ _geometricProlongation(isCorner, edgeEndpoints, nDim) @ sp.diags(scaling[cornerDofs])
freeCoarse = np.nonzero(~isDirichlet[cornerDofs])[0]
P_free = P.tocsr()[free, :][:, freeCoarse].tocsr()
coarseSolve = spla.factorized((P_free.T @ A_free @ P_free).tocsc())
# amgcl::relaxation::chebyshev with its float lower/higher parameters
lo = spectralRadius * float(np.float32(0.01))
hi = spectralRadius * float(np.float32(1.1))
d, c = 0.5 * (hi + lo), 0.5 * (hi - lo)
def smooth(x, b):
p = np.zeros_like(x)
alpha = 0.0
for k in range(degree):
residual = b - A_free @ x
if k == 0:
alpha, beta = 1.0 / d, 0.0
elif k == 1:
alpha = 2 * d * (1.0 / (2 * d * d - c * c))
beta = alpha * d - 1.0
else:
alpha = 1.0 / (d - 0.25 * alpha * c * c)
beta = alpha * d - 1.0
p = alpha * residual + beta * p
x = x + p
return x
rFree = r[free]
xFree = np.zeros_like(rFree)
for _ in range(nuPre):
xFree = smooth(xFree, rFree)
xFree = xFree + P_free @ coarseSolve(P_free.T @ (rFree - A_free @ xFree))
for _ in range(nuPost):
xFree = smooth(xFree, rFree)
x = r.copy()
x[free] = xFree
return x
Rigid-body coarse level: TwoLevel(coarseSpace=RigidBodyModesPerAggregate())
By default, the hierarchy of the 3D displacement field starts with a coarse level built from all of its rigid-body
modes – the translations and the rotations – in place of AMGCL’s hierarchy on the translations
(TwoLevel(coarseSpace=RigidBodyModesPerAggregate()), rigidbodymodestwolevelcycle); a 2D
displacement field gets it in a custom solver script. It keeps what makes the block backend fast, the fine level
in node-wise blocks stored in float, and adds what that backend cannot take, the rotations.
Why the translations are not enough
The coarse space of a block hierarchy holds the translations of every aggregate (see the near-null-space of a block backend): displacements that are constant on each aggregate. Smoothed aggregation represents a smooth error well only where it is close to such a combination on every aggregate, and the Chebyshev smoother leaves smooth error alone by construction: its window starts at 1% of the spectral radius (see the near-null-space and the Chebyshev smoother above). In a bending-dominated structure – a slender beam, a tunnel lining, a plate – the smooth error is bending: every aggregate rotates almost as a rigid body, its displacement varying linearly across it. Neither part of the cycle reduces such an error, and one application of the hierarchy on the translations can even amplify it. Measured on the displacement blocks \(A_{uu}\) of captured systems, \(\lVert v - M A_{uu} v \rVert / \lVert v \rVert\) after one application \(M\) of the field’s preconditioner, for the rigid rotations \(v\) about the field’s centroid and for a smooth quadratic field (\(u_x = y^2 + z^2\) about the centroid), both zero on the Dirichlet rows:
error |
tunnel segment, elastic |
tunnel segment, damaged |
damage cylinder |
edge breakout |
|---|---|---|---|---|
rotations about x, y, z |
1.12 2.06 0.88 -> 0.15 0.27 0.47 |
1.21 0.97 0.81 -> 0.12 0.10 0.18 |
0.28 0.28 0.80 -> 0.12 0.11 0.12 |
0.19 0.37 0.50 -> 0.03 0.03 0.02 |
quadratic field |
1.53 -> 0.62 |
1.34 -> 0.46 |
1.50 -> 1.06 |
0.30 -> 0.17 |
A value above 1 means that one application amplifies that error. On the tunnel segment, a shell-like lining of quadratic hexahedra, the hierarchy on the translations amplifies the rotations and the quadratic field, and the outer solve pays with iterations: 67 of them, against 22 with the rotations in the coarse space (see The rigid-body coarse level). Not every error fares better: the rigid translations, which both coarse spaces hold, are reduced less on the elastic tunnel system (0.11 0.04 0.20 -> 0.70 0.18 0.41). The outer iteration counts of The rigid-body coarse level are the net effect.
How the coarse level is built
On the field’s equilibrated block \(\hat A = D A D\), with its near-null-space \(B\) – the field’s rigid-body
modes divided by \(D\) (the near-null-space of \(D A D\) is \(D^{-1}\) times that of \(A\), as for the
translations), which Python builds with rigidBodyNullspace(), the one definition
of the modes, and hands to C++ – in C++, amgcl-rigid-body-coarse-level.hpp, one function per step:
Aggregates. AMGCL’s plain aggregation of the nodes: the aggregates of AMGCL’s own first level on the block backend, with its strength test on the \(B \times B\) blocks (see above) and the tree’s
coarsening.aggr.eps_strong. A node without a strong connection – a fully constrained one, for instance – is in no aggregate and left to the smoother.Tentative prolongation. Per aggregate \(a\), the rows of \(B\) on its dofs, orthonormalized – what AMGCL’s own tentative prolongation does with a near-null-space:
\[V_a = B \big|_a = \begin{bmatrix} b_1 & \cdots & b_M \end{bmatrix}_a = Q_a R_a, \qquad P^{\mathrm{tent}} \big|_a = Q_a ,\]with \(b_1, \dots, b_B\) the translations and \(b_{B+1}, \dots, b_M\) the rotations about the field’s centroid. The rows of Dirichlet dofs and multi-point-constraint slaves are zero. The QR factorization (modified Gram–Schmidt with a second sweep) makes the columns of \(P^{\mathrm{tent}}\) orthonormal, and \(P^{\mathrm{tent}} R_a\) reproduces the modes on the aggregate exactly. A rotation about the field’s centroid \(g\) carries a translation of the order of \(|c_a - g|\) across an aggregate of extent \(h_a\) around its centroid \(c_a\); projecting the translations out loses about \(\log_{10}(|c_a - g| / h_a)\) digits, 2 to 4 on a mesh of elements a hundred to ten thousand times smaller than the structure, far from the threshold of the dependent modes below, so the modes need no recentring on every aggregate. The Dirichlet rows are zero because the cycle solves them exactly, apart from the coarse-grid correction (as the corner-node two-level cycle does), and a rigid-body mode is near-null only on the free dofs: so \(P^\mathsf{T} \hat A P = P_F^\mathsf{T} \hat A_{FF} P_F\).
Dependent modes. A mode the Gram–Schmidt sweeps reduce to at most \(\sqrt{\varepsilon_{\mathrm{mach}}}\) (1.5e-8) of its norm lies in the span of the modes before it: the rotations of an aggregate of a single node, the rotation about the line through a collinear aggregate, a mode whose components are all constrained (an aggregate on a symmetry plane). Its normalized remainder would be rounding noise, not a direction of the aggregate, so it is dropped: its column of \(P^{\mathrm{tent}}\) is zero, and its coarse dof gets a diagonal entry and nothing else. The threshold is not tuned: an exactly dependent mode leaves a remainder of about 1e-16 of its norm, an independent one of the order of \(h_a / |c_a - g|\) (see above), or of the aggregate’s aspect ratio.
Smoothed prolongation. \(P = (I - \omega\, \mathrm{diag}_B(\hat A_F)^{-1} \hat A_F) P^{\mathrm{tent}}\), AMGCL’s own, shared with the translations of the block backend: \(\hat A_F\) holds the strong connections, the weak ones lumped onto the diagonal blocks, \(\mathrm{diag}_B\) takes its node-wise \(B \times B\) diagonal blocks, and \(\omega\) follows the tree’s
coarsening(relax,estimate_spectral_radius, and the spectral radius estimate, Near-null-space and the Chebyshev smoother). One damped block-Jacobi step widens every column by a layer of nodes and smooths away its jump at the aggregate’s boundary. With the option"prolongation": "tentative", \(P = P^{\mathrm{tent}}\).Galerkin coarse operator. \(A_c = P^\mathsf{T} \hat A P\), with \(M = B (B + 1) / 2\) coarse dofs per aggregate, 6 in 3D and 3 in 2D, formed block by block, in \(M \times M\) blocks, and handed to AMGCL unblocked to a scalar matrix: a dropped mode’s coarse dof gets the mean of the kept diagonal entries as its diagonal entry, and the zeros of its row and column stay stored – AMGCL’s aggregation below takes the largest entry of every \(M \times M\) block, which they do not change.
The hierarchy below. On aggregate \(a\), the field’s rigid-body modes are \(V_a = P^{\mathrm{tent}} R_a\): their coefficients in the coarse basis, \(R_a\), are the near-null-space of an AMGCL hierarchy on \(A_c\), as in AMGCL’s own recursion. The field’s tree configures it as AMGCL would configure its next level: the strong-connection threshold halved, the \(M\) dofs of an aggregate one point (
aggr.block_size), in a scalar hierarchy.
What runs outside AMGCL is steps 2 to 5: the tentative prolongation (in \(B \times M\) blocks), its smoothing
(AMGCL’s own algorithm, on those blocks), the Galerkin product (in \(M \times M\) blocks) and its unblocking; the
aggregates of step 1 and the hierarchy below are AMGCL’s. AMGCL’s as_scalar coarsening, which would run a scalar
smoothed aggregation with all modes on the block backend, is not used, for three reasons: it re-blocks the prolongation
into \(B \times B\) blocks, which the \(M = 3\) modes of a 2D field do not fit; its Householder QR keeps an
orthonormal completion of a dependent mode instead of dropping it; and its modes are not zero on the Dirichlet rows.
And AMGCL’s sparse products take one block type, where the Galerkin product chains \(B \times B\),
\(B \times M\) and \(M \times B\) blocks: so the blocks have kernels of their own.
One application is the two-level cycle the corner-node interpolation uses too: restriction \(P^\mathsf{T} r\), one V-cycle of the
hierarchy below, prolongation, and the fine Chebyshev sweeps of the tree (npre before, npost after the
coarse-grid correction: 0 and 1, degree 16, in the default trees), with the Dirichlet rows decoupled exactly.
\(P\) and \(P^\mathsf{T}\) are stored like the fine level (in float with the 3D default tree), every sum in
double. With the translations alone in place of all modes, the construction is AMGCL’s first level on the block
backend – the same aggregates, tentative prolongation, smoothing and Galerkin product, \(M = B\) – except on the
Dirichlet rows, where AMGCL’s translations are not zero. A test of the suite checks that, on a system without Dirichlet rows,
the two preconditioners then agree to rounding: what the rigid-body coarse level adds is the rotations.
When it pays, and when to switch it off
The coarse level pays where the hierarchy on the translations needs many outer iterations although the operator is
elastic or mildly damaged: typically a bending-dominated structure. It is the default for the 3D displacement field, of
any element type: the rigid-body modes are the near null-space of elasticity on the displacement field (the field of
type vector in phenomena). A rotation field – of beams, or of the reference
point of a rigid body – or any other vector field has another near null-space (its translations), and
RigidBodyModesPerAggregate rejects it (carriesRigidBodyModes() decides, by the field’s
phenomenon). A 2D displacement field gets the coarse level in a custom script, where every field’s preconditioner
is written out, OnField(name, TwoLevel(coarseSpace=RigidBodyModesPerAggregate(), amgclTree=...)), and nothing is added to it. A 3D
displacement field without node coordinates (see below) keeps the hierarchy on the translations in the preset, logged
once at info. The field’s tree, the default one or one given to TwoLevel(amgclTree=...), configures the smoother,
the aggregation and the hierarchy below.
Switch it off on a model whose solves converge in a few outer iterations with the hierarchy on the translations: the
coarse level saves little there, and its setup is larger wherever its coarse operator is – up to 1.2 times that of
the hierarchy on the translations on the larger models of The rigid-body coarse level. On the edge breakout, where the translations need 22
and 23 outer iterations, the coarse level saves 2 and 1, at a setup of 0.94 and 0.99 times theirs and a solve time,
setup included, of 0.92 and 0.95 – a wash. The data line rigidBodyModesPerAggregate=false of type=blockamg
switches it off; its options are those of a custom script (see its arguments).
The modes are built from the node coordinates of the reference configuration, as the model holds them, which
the coarse level reads from the model: after adaptive refinement those of the new model, under a SchurComplement those
of the remaining nodes. Rigid-body modes about the reference configuration suit small rotations; under large rotations
they represent the rigid motions of the deformed structure less exactly, which costs convergence speed only, not
accuracy. Without a model – a script driving the solver itself – the coordinates must be given as
nodeCoordinates (see The adapter): without them, TwoLevel(coarseSpace=RigidBodyModesPerAggregate())
raises an error rather than switch off the coarse level it asks for, while the preset keeps the hierarchy on the
translations.
What it does not fix:
Softening. The rigid-body modes are the near-null-space of the elastic operator. Once damage makes the displacement block indefinite, other directions of low energy appear that no coarse space of rigid bodies holds: on the damaged tunnel system (The rigid-body coarse level), the iterations fall from 61 to 36, against 67 to 22 on the elastic one.
Scalar fields, and the scalar backend. A scalar field has no rotations, and the scalar backend (
AMG(backendBlockSize=1)) takes the rigid-body modes as its near-null-space itself (AMG(nearNullspace="rigidBodyModes"), the default);TwoLevel(coarseSpace=RigidBodyModesPerAggregate())on either is rejected.
Schur complement condensation: SchurComplement (opt-in)
Aggregation AMG fails on some dofs whatever the field split: the dofs of embedded bars with bond-slip are the prominent case. A bar node is tied by bond springs to the nodes of the host element it lies in, so a low-energy error has the bar follow the host, \(u_{bar} \approx \sum_i N_i u_i\) – an interpolation of nodes at other positions, which an aggregate moving as a rigid body cannot represent. The bars are also stiff 1D chains inside a softer 3D continuum, a strong anisotropy. The field-split preconditioner alone stalls on such systems. (With perfect bond, the same relation is a multi-point constraint, eliminated before the solve, and the field-split preconditioner without condensation is the right choice.)
SchurComplement condenses these dofs exactly – we call them \(b\), the condensed (“bar”) dofs – and leaves
AMG the remaining (kept) dofs \(k\):
The Schur complement \(S\) is the operator of the remaining dofs once \(x_b\) is eliminated – for bars, the concrete stiffened by the bars as seen through their bond. It is again an operator of a continuum, which is what the per-field hierarchies handle well. The preconditioner of the outer Krylov solve is the block LU factorization,
which is nothing but one symmetric block Gauss–Seidel sweep over the two blocks \([b, k]\): the first and the last step solve the \(b\) block exactly (a direct factorization of \(A_{bb}\)), the middle step is the usual field-split sweep \(M_{\tilde S}^{-1}\), built on an approximation \(\tilde S\) of \(S\) instead of on \(A\). With \(\tilde S = S\) and exact field solves, the preconditioner is \(A^{-1}\); the outer Krylov solve corrects for both approximations.
This is the inverse of \(A\) written as the product of both block factors (blockFactors = "both"), with two
solves with \(A_{bb}\). By default (blockFactors = "upper") only the block upper triangular factor
\(U = \begin{bmatrix} A_{bb} & A_{bk} \\ 0 & S \end{bmatrix}\) of
\(A = \begin{bmatrix} I & 0 \\ A_{kb} A_{bb}^{-1} & I \end{bmatrix} U\) is inverted, with one solve:
The remaining dofs are solved first, and the condensed ones follow them – for bars, the bars follow the concrete, which is exactly the relation AMG cannot represent.
Why one factor is enough: its last step solves the rows of \(b\) exactly, so after the first outer iteration the residual in these rows is zero, and from then on the outer Krylov solve iterates on the rows of \(k\) alone, with the operator \(S \tilde S^{-1}\) – the same as with the full block LU factorization, which saves one iteration more. With \(\tilde S = S\) and exact field solves, GMRES converges in two iterations instead of one: \(A U^{-1}\) is the unit block lower triangular factor \(L\) of \(A = L U\), and \((L - I)^2 = 0\). With the approximations, the outer iterations were found to hardly grow, and every one saves a solve with \(A_{bb}\). The block lower triangular factor (bars first, concrete held fixed) was tried too: its last step leaves a residual in the rows of \(b\), coupled to \(k\) by the stiff bond, and it needed considerably more outer iterations.
The approximate Schur complement. \(S\) is dense wherever \(A_{bb}^{-1}\) is – along a whole bar. But the entries of the inverse of a banded, well-conditioned matrix decay exponentially away from the diagonal (Demko, Moss and Smith 1984); physically, a bar on its bond is a beam on an elastic foundation, and a load on it is felt only a few host elements away. \(\tilde S\) therefore keeps only the significant entries:
a banded component of \(A_{bb}\) – a chain: a bar, ordered along the bar by the reverse Cuthill–McKee ordering – is cut into overlapping windows, cores of
windowNodesnodes with halos ofhaloNodesnodes on either side; for the dofs of a core, \(A_{bb}^{-1}\) is replaced by the inverse of the window’s block, and only the core rows are kept (restricted additive Schwarz, Cai and Sarkis 1999). Windows count nodes, as the decay length is a length along the bar, whatever the dofs per node (a truss, a beam). A short component (a short bar) is a single window, which is exact;a component that is not banded (a zone of a 2D or 3D mesh, e.g. a softening zone) gets its exact correction from a factorization of its own block, solved for the coupled remaining dofs in batches of right-hand sides – windows along a chain would cut such a zone into pieces far smaller than the distance over which its inverse decays;
entries \(|\tilde S_{ij}| < \tau \sqrt{|\tilde S_{ii} \tilde S_{jj}|}\) are dropped (
dropTolerance).
Choosing the condensed dofs. The markers of adaptive mesh refinement select elements
(edelweissfe.adaptivity.marking; any type with the options of its >>marker line), and \(b\) holds
all dofs of the interior nodes of the marked elements: the nodes that have at least one element, all of them
marked. One rule serves every use:
bar node: bar, bond -> all marked -> condensed
host node: bond, concrete, concrete -> not all -> kept
embedded bars: mark the bar elements and the bond elements (
elementSetmarkers for the bar and the bond element sets). A bar node belongs only to them; a host node also belongs to continuum elements, so it stays in \(k\);a softening zone: a
fieldOutputmarker on aperElementfieldOutput (damage, a hardening variable) with a threshold, and ahaloso that the nodes of the zone itself become interior. All dofs of a node are condensed, e.g. displacement and nonlocal damage. A fieldOutput has no result before the first converged increment, so nothing is condensed until then.
The markers are evaluated again whenever what they read can change (initialOnly is rejected): a fieldOutput
marker when its fieldOutput is updated, after a converged increment; a marker of sets or surfaces when the model
changes; a marker that cannot tell (e.g. recoveryError) on every solve. A changed selection rebuilds
\(\tilde S\) and the hierarchies. Should the markers select every
dof, nothing is condensed (with a warning): there would be nothing left for the field split, and a direct solver is
the better choice; more than half of the dofs is warned about too. Only the dofs of nodes can be condensed, and
always all fields of a node together: dofs that belong to no node (e.g. Lagrange multipliers) cannot be selected.
Robustness. If \(A_{bb}\) cannot be factorized (a zone or a bond that lost all stiffness), or \(\tilde S\) cannot be formed (a singular window block of an \(A_{bb}\) that factorizes as a whole), the solve does without the condensation, for all dofs, with a warning: the condensed dofs come from the last converged increment, so a cutback would meet the same block again. The condensation then keeps the state of its last successful update. An approximate Schur complement with non-positive diagonal entries (the condensed zone does not hold all of the softening, and AMG on the rest may struggle) is warned about as well.
Within a solve. Every solve cuts \(A_{bb}\), \(A_{bk}\) and \(A_{kb}\) out of the current matrix
(at positions found once per sparsity pattern) and factorizes \(A_{bb}\); the operator of its outer solve,
\(A\) itself, is kept per sparsity pattern too, its values refreshed in place. \(\tilde S\) is computed when
the hierarchies are rebuilt – under the same reuse rules as without condensation – and the hierarchies are
built on it: the field blocks are those of the remaining dofs (a field whose dofs are all condensed, such as the
rotations of embedded beams, has none), and the node coordinates of the near null-space and of the rigid-body
coarse level and the midside corners of p-multigrid are restricted to the remaining nodes (a midside node whose edge endpoint is condensed becomes a
corner). Between rebuilds, the sweep keeps the couplings of the standing \(\tilde S\), while \(A_{bb}\),
\(A_{bk}\) and \(A_{kb}\) are those of the current matrix: unlike without condensation, where the couplings
between the fields are taken from the current matrix on every solve. The outer Krylov solve runs on \(A\)
itself: the Schur complement sits under the outer Krylov solve, LGMRES(SchurComplement(...)).
Two sweeps. The sweep in residual-correction form makes sweeps = 2 exactly two steps of the stationary
iteration \(z \leftarrow z + M^{-1} (r - \tilde S z)\) with the one-sweep preconditioner \(M^{-1}\),
which on bar systems was found a good balance between AMG cycles and outer iterations; sweeps = 1 is
cheaper per iteration and needs more of them.
The direct solver. condensedBlockSolve=DirectSolver(package) takes any registered direct solver package; superlu
and pardiso factorize once for several right-hand sides, the others (klu, mumps, petsclu,
umfpack, panuapardiso) factorize again for every one, which makes them a poor choice here. Every outer
iteration solves with \(A_{bb}\), so the cost of a back substitution counts most:
for the block-tridiagonal chains of bars, superlu (the default) is the cheaper one (with blockFactors =
"both", twice per iteration). A large compact zone (a
softening zone of a 3D mesh) shifts the weight to the exact correction of \(\tilde S\), which solves for many
right-hand sides at once – there a threaded factorization (pardiso) does better. For a softening zone, mind
that condensation pays only where AMG struggles with the zone: the field-split preconditioner alone handles a typical
gradient-enhanced damage band well, and the condensation then costs more per iteration and per rebuild of
\(\tilde S\) (rebuilt whenever the zone grows) than it saves in iterations.
The condensation is a custom solver script, LGMRES(SchurComplement(condensedDofMarkers=[...], schurComplementSolve=...)): a
complete one is A Schur complement, the arguments of
SchurComplement are listed in
The arguments of SchurComplement.
A worked example
The whole method, written out with NumPy, SciPy and one AMGCL hierarchy per field on a small
displacement / damage system – equilibration, field split, per-field AMG, one block Gauss–Seidel sweep as
the preconditioner, and the outer Krylov solve – and the same solve as a composed solver. It is a test
of the test suite; displacementDamageSystem, ExactFieldSolve and hierarchyLevels come from the
suite’s tests/_linsolvesystems.py:
import numpy as np
import scipy.sparse as sp
from _linsolvesystems import ExactFieldSolve, displacementDamageSystem, hierarchyLevels
from scipy.sparse.linalg import LinearOperator, gmres
from edelweissfe.linsolve.composition import (
AMG,
LGMRES,
ComposedLinearSolver,
DiagonalScaling,
FixedTolerance,
Multiplicative,
OnField,
)
#: One AMGCL hierarchy per field: smoothed aggregation, Chebyshev smoothing.
FIELD_AMG = {"coarsening": {"type": "smoothed_aggregation"}, "relax": {"type": "chebyshev"}}
def test_the_algorithm_step_by_step():
# A coupled system with the field structure of a gradient-enhanced damage model: a "displacement" field (3 DOFs
# per node) and a nonlocal damage field (1 DOF per node) on a 15^3 node grid, coupled both ways, the displacement
# rows scaled over two orders of magnitude. Not its physics: the displacement block is a vector Laplacian (three
# copies of a grid Laplacian) standing in for elasticity, so it shows the method, not why elasticity needs the
# rigid-body modes.
A, b, fieldBlocks = displacementDamageSystem(15, coupling="symmetric", rowScaling=True)
fields = [slice(block.start, block.stop) for block in fieldBlocks]
# The steps are numbered as in "The solve pipeline" of the theory page. Steps (1) and (7) carry state from one
# Newton iteration to the next and have no part in a single solve.
# (2) Equilibrate: solve (D A D) z = D b with D = diag(scaling), scaling = 1 / sqrt|a_ii|, then x = D z.
scaling = 1.0 / np.sqrt(np.abs(A.diagonal()))
D = sp.diags(scaling)
A_scaled = (D @ A @ D).tocsr()
b_scaled = scaling * b
# (3) Split by field into the blocks A_ij, and (4) build one AMG hierarchy per diagonal block A_ii, with AMGCL
# (an optional extension of EdelweissFE, hence imported here).
A_blocks = [[A_scaled[row, column] for column in fields] for row in fields]
from edelweissfe.linsolve.amgcl.amgcl import PyAMGCLSolver
fieldAMG = []
for i in range(len(fields)):
amg = PyAMGCLSolver({"precond": FIELD_AMG})
amg.build(A_blocks[i][i])
assert hierarchyLevels(amg) > 1 # a real hierarchy, not a direct solve of a small block
fieldAMG.append(amg)
# (5a) The preconditioner v = M^-1 r: one forward block Gauss-Seidel sweep over the fields, each field's part of
# v from one AMG cycle on its own residual, v_i = M_i^-1 (r_i - sum_{j < i} A_ij v_j). Without the coupling
# sum it is block Jacobi.
def blockGaussSeidel(r, fieldSolves=fieldAMG, coupled=True):
v = np.zeros_like(r)
for i, field in enumerate(fields):
residual = r[field].copy()
if coupled:
for j in range(i):
residual -= A_blocks[i][j] @ v[fields[j]]
v[field] = fieldSolves[i].applyPreconditioner(residual)
return v
# With exact field solves, the sweep is the block lower triangular solve M^-1 = (B + L)^-1, B the
# diagonal blocks and L the blocks below them: what the AMG cycles approximate. Check: (B + L) v = r.
exact = [ExactFieldSolve(A_blocks[i][i]) for i in range(len(fields))]
B_plus_L = sp.bmat([[A_blocks[i][j] if j <= i else None for j in range(len(fields))] for i in range(len(fields))])
r = np.random.default_rng(1).standard_normal(A.shape[0])
assert np.linalg.norm(B_plus_L @ blockGaussSeidel(r, exact) - r) <= 1e-10 * np.linalg.norm(r)
# (5b) The outer Krylov solve on the coupled system, right-preconditioned by the sweep: GMRES solves
# (A_scaled M^-1) y = b_scaled for y, and z = M^-1 y. GMRES converges the residual of the operator it is
# given, b_scaled - A_scaled M^-1 y = b_scaled - A_scaled z, the scaled residual D (b - A x).
A_scaled_Minv = LinearOperator(A.shape, matvec=lambda y: A_scaled @ blockGaussSeidel(y))
y, info = gmres(A_scaled_Minv, b_scaled, rtol=1e-10, atol=0.0, restart=200, maxiter=1)
z = blockGaussSeidel(y)
x = scaling * z
# (6) The true residual b - A x = D^-1 (b_scaled - A_scaled z) is then at most max(D) / min(D) times larger,
# relative to b: here 100 (sqrt of the 1e4 spread of the diagonal), hence 1e-8 after a 1e-10 solve.
assert info == 0 and np.linalg.norm(b - A @ x) <= 1e-8 * np.linalg.norm(b)
# (5b) on the unscaled system, with the same preconditioner between two scalings (krylovOnScaledSystem = false,
# the default): GMRES solves (A D M^-1 D) y = b for y, and x = D M^-1 D y. The residual it
# converges, b - A D M^-1 D y = b - A x, is then the true one, and step (6) has nothing left to check.
A_DMinvD = LinearOperator(A.shape, matvec=lambda y: A @ (scaling * blockGaussSeidel(scaling * y)))
y, info = gmres(A_DMinvD, b, rtol=1e-10, atol=0.0, restart=200, maxiter=1)
x_unscaled = scaling * blockGaussSeidel(scaling * y)
assert info == 0 and np.linalg.norm(b - A @ x_unscaled) <= 1.0001e-10 * np.linalg.norm(b)
assert np.linalg.norm(x_unscaled - x) <= 1e-6 * np.linalg.norm(x)
# What the preconditioner buys: GMRES iterations to the same 1e-10, without a preconditioner, with the sweep
# (AMG cycles), and with exact field solves, block Jacobi against block Gauss-Seidel. Measured: far fewer
# iterations with the sweep than without a preconditioner at all; with AMG-cycle field solves, block Jacobi
# needs about as many iterations as block Gauss-Seidel (the coupling is weak, and the AMG cycles' own error
# dominates); with exact field solves, the coupling term of the sweep becomes visible and Gauss-Seidel needs
# fewer iterations than Jacobi.
def iterations(preconditioner):
history = []
A_scaled_Minv = LinearOperator(A.shape, matvec=lambda y: A_scaled @ preconditioner(y))
gmres(
A_scaled_Minv,
b_scaled,
rtol=1e-10,
atol=0.0,
restart=200,
maxiter=1,
callback=history.append,
callback_type="pr_norm",
)
return len(history)
unpreconditioned = iterations(lambda y: y)
assert unpreconditioned > 5 * iterations(blockGaussSeidel)
assert iterations(lambda y: blockGaussSeidel(y, exact)) < iterations(lambda y: blockGaussSeidel(y, exact, False))
# The same solve as a composed solver: the same method, with the field split in threaded C++ and AMGCL's
# right-preconditioned LGMRES as the outer solver (5b), on the unscaled system, LGMRES(DiagonalScaling(...)) as
# the blockamg preset builds it, or on the scaled one, DiagonalScaling(LGMRES(...)). The preset adds the reuse of
# the hierarchies across Newton iterations, an adaptive (Eisenstat--Walker) tolerance, and a near null-space for
# the displacement field: its translations here, given only the field layout, and all rigid-body modes given the
# model (setModel).
def blockGaussSeidelTree():
return Multiplicative([OnField(field.name, AMG(amgclTree=dict(FIELD_AMG))) for field in fieldBlocks])
for tree in (
LGMRES(DiagonalScaling(blockGaussSeidelTree()), tolerance=FixedTolerance(1e-10)),
DiagonalScaling(LGMRES(blockGaussSeidelTree(), tolerance=FixedTolerance(1e-10))),
):
solver = ComposedLinearSolver(tree, verbosity="silent")
solver.setFieldStructure(fieldBlocks)
assert np.linalg.norm(solver(A, b) - x) <= 1e-6 * np.linalg.norm(x)
The written-out half is for reading: its products run in SciPy, on one thread. For a real model use the composed solver (Which form to use: what each costs).
Practical guidance
Reach for the field-split preconditioner (
type=blockamg) when a direct factorization no longer fits in memory, or when the factorization dominates the run. Below that, measure againstpardisorather than assuming.Run the first characterisation at
verbosity=infoand read the job’s performance table. Outer iteration counts up to the low hundreds were the expected regime at the former defaults; with the current ones, mostly 5–50, with the hardest damage-onset solves reaching the cap; counts triggering the warning threshold, or cap hits, repeatedly point at the preconditioner, not the tolerance.Set
OMP_NUM_THREADSfor the AMGCL kernels as usual (see Parallelization). Results are reproducible at a fixed thread count up to last-bit summation order; across thread counts the hierarchies differ slightly (see Near-null-space and the Chebyshev smoother), and lowering the power-iteration setting of the vector-field default makes that difference large.Reach for the knobs in roughly this order: the tree of a field whose hierarchy is clearly the weak one (a
customscript), thensweeps/symmetricfor the sweep itself, then the tolerance.The corner-node two-level cycle is on for a quadratic 3D displacement field by default, the rigid-body coarse level for a linear one: give a 2D displacement field either in a
customscript, and switch the rigid-body coarse level off (rigidBodyModesPerAggregate=false) on a model whose solves converge in a few outer iterations without it, see p-multigrid: TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve)) and Rigid-body coarse level: TwoLevel(coarseSpace=RigidBodyModesPerAggregate()).When a run degrades reproducibly, capture it with
dumps=slowplus a fewdumpContextSolvesrather than trying to reproduce it by rerunning the simulation (Dumps of the solves).
Further reading
S. C. Eisenstat, H. F. Walker, Choosing the forcing terms in an inexact Newton method, SIAM J. Sci. Comput. 17 (1996) 16–32 – the adaptive outer tolerance, “choice 2”.
P. Vaněk, J. Mandel, M. Brezina, Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems, Computing 56 (1996) 179–196 – the coarsening, and why it needs a near-null-space basis.
A. H. Baker, E. R. Jessup, T. Manteuffel, A technique for accelerating the convergence of restarted GMRES, SIAM J. Matrix Anal. Appl. 26 (2005) 962–984 – LGMRES and its recycled vectors.
D. Demidov, AMGCL: an efficient, flexible, and extensible algebraic multigrid implementation, Lobachevskii J. Math. 40 (2019) 535–546 – the library providing the hierarchies, the relaxation kernels and the native outer solver.
N. Alkmim et al., Int. J. Numer. Methods Eng. (2026) – the field-split block-preconditioning strategy for coupled multi-field fracture models that this solver follows.
S. Demko, W. F. Moss, P. W. Smith, Decay rates for inverses of band matrices, Math. Comp. 43 (1984) 491–499 – why the Schur complement of condensed bars is local.
X.-C. Cai, M. Sarkis, A restricted additive Schwarz preconditioner for general sparse linear systems, SIAM J. Sci. Comput. 21 (1999) 792–797 – the windows of the approximate Schur complement.
H. C. Elman, D. J. Silvester, A. J. Wathen, Finite Elements and Fast Iterative Solvers, Oxford University Press (2014) – block (Schur complement) preconditioners.