Linear solvers: Python API

The classes and functions of the linear solvers. How to choose a linear solver is on Linear solvers, how to compose one on Composing a linear solver, the arguments and data lines at a glance on Composed linear solvers: reference.

The types and the configuration

The types of linear solver: what linsolver=<type> (and, with its data lines, *linearSolver, type=<type>) names.

Every type is a composed solver behind the one adapter, ComposedLinearSolver: a type is the schema of its data lines and the tree they define (LinearSolverType.tree()). A direct solver is a tree of one leaf, type=pardiso is DirectSolver("pardiso"). So what the adapter does – the verbosity, the dumps, the “linear solve” column – works the same for every type, and the data lines of LinearSolverSchema are those of every type.

class edelweissfe.linsolve.types.BlockAMG[source]

blockamg: the preset blockAMG(), block Gauss–Seidel over the fields of the diagonally scaled matrix, inside LGMRES.

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

schema: alias of BlockAMGSchema

classmethod tree(configuration)[source]

The composed solver of these data lines.

Parameters:: configuration (BlockAMGSchema)
Return type:: Preconditioner

class edelweissfe.linsolve.types.BlockAMGSchema(verbosity='warning', warningOuterIterations=100, residualGrowthFactor=4.0, dumps='none', dumpDirectory='linsolveDumps', maxDumps=10, dumpContextSolves=0, rigidBodyModesPerAggregate=True, cornerNodeInterpolation=True, tolerance=None, maxIterations=100, reusePreconditioner=False, sweeps=1, symmetric=False, krylovOnScaledSystem=False, trueResidualRestarts=2, residualGapCompensation=False, dropTolerance=0.0, lumpDroppedEntries=False)[source]

The data lines of blockamg: the arguments of the preset blockAMG().

Attributes:

tolerance

Parameters:

verbosity (str)
warningOuterIterations (int)
residualGrowthFactor (float)
dumps (str)
dumpDirectory (str)
maxDumps (int)
dumpContextSolves (int)
rigidBodyModesPerAggregate (bool)
cornerNodeInterpolation (bool)
tolerance (float)
maxIterations (int)
reusePreconditioner (bool)
sweeps (int)
symmetric (bool)
krylovOnScaledSystem (bool)
trueResidualRestarts (int)
residualGapCompensation (bool)
dropTolerance (float)
lumpDroppedEntries (bool)

class edelweissfe.linsolve.types.Custom[source]

custom: the composed solver of a solver script, defined by you (see loadSolverScript()).

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

schema: alias of CustomSchema

classmethod create(configuration)[source]

The linear solver of these data lines: the tree behind the adapter.

Parameters:: configuration (CustomSchema)
Return type:: ComposedLinearSolver

class edelweissfe.linsolve.types.CustomSchema(verbosity='warning', warningOuterIterations=100, residualGrowthFactor=4.0, dumps='none', dumpDirectory='linsolveDumps', maxDumps=10, dumpContextSolves=0, script=None, hotReload=False)[source]

The data lines of custom.

Attributes:

script

Parameters:

verbosity (str)
warningOuterIterations (int)
residualGrowthFactor (float)
dumps (str)
dumpDirectory (str)
maxDumps (int)
dumpContextSolves (int)
script (str)
hotReload (bool)

class edelweissfe.linsolve.types.DirectSolverType[source]

A direct solver: the tree of one DirectSolver leaf on the package package, factorized anew on every solve.

Attributes:

package

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

classmethod tree(configuration)[source]

The composed solver of these data lines.

Return type:: Preconditioner

package: str = None: The registered name of the direct solver package (registry category directsolver).

class edelweissfe.linsolve.types.KLU[source]

klu: SuiteSparse’s KLU, if built.

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

package: str = 'klu': The registered name of the direct solver package (registry category directsolver).

class edelweissfe.linsolve.types.LinearSolverSchema(verbosity='warning', warningOuterIterations=100, residualGrowthFactor=4.0, dumps='none', dumpDirectory='linsolveDumps', maxDumps=10, dumpContextSolves=0)[source]

The data lines of every type: how its solves are reported and dumped.

Parameters:

verbosity (str)
warningOuterIterations (int)
residualGrowthFactor (float)
dumps (str)
dumpDirectory (str)
maxDumps (int)
dumpContextSolves (int)

class edelweissfe.linsolve.types.LinearSolverType[source]

A type of linear solver: the schema of its data lines, and the composed solver they define.

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

schema: alias of LinearSolverSchema

classmethod create(configuration)[source]

The linear solver of these data lines: the tree behind the adapter.

Return type:: ComposedLinearSolver

classmethod fromOptions(options)[source]

The linear solver of these options (the data lines, by name; unknown names raise).

Parameters:: options (dict)
Return type:: ComposedLinearSolver

classmethod tree(configuration)[source]

The composed solver of these data lines.

Return type:: Preconditioner

static adapterSettings(configuration)[source]

The settings of the adapter (ComposedLinearSolver) in the data lines every type shares (LinearSolverSchema), as its keyword arguments.

Return type:: dict

class edelweissfe.linsolve.types.Mumps[source]

mumps: MUMPS, if its Python bindings are installed.

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

package: str = 'mumps': The registered name of the direct solver package (registry category directsolver).

class edelweissfe.linsolve.types.PanuaPardiso[source]

panuapardiso: Panua PARDISO, if built.

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

package: str = 'panuapardiso': The registered name of the direct solver package (registry category directsolver).

class edelweissfe.linsolve.types.Pardiso[source]

pardiso: Intel MKL PARDISO (the default linear solver, if built).

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

schema: alias of PardisoSchema

classmethod tree(configuration)[source]

The composed solver of these data lines.

Parameters:: configuration (PardisoSchema)
Return type:: Preconditioner

package: str = 'pardiso': The registered name of the direct solver package (registry category directsolver).

class edelweissfe.linsolve.types.PardisoSchema(verbosity='warning', warningOuterIterations=100, residualGrowthFactor=4.0, dumps='none', dumpDirectory='linsolveDumps', maxDumps=10, dumpContextSolves=0, reuseSymbolicFactorization=False)[source]

The data lines of pardiso.

Parameters:

verbosity (str)
warningOuterIterations (int)
residualGrowthFactor (float)
dumps (str)
dumpDirectory (str)
maxDumps (int)
dumpContextSolves (int)
reuseSymbolicFactorization (bool)

class edelweissfe.linsolve.types.PetscLU[source]

petsclu: PETSc’s LU, if petsc4py is installed.

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

package: str = 'petsclu': The registered name of the direct solver package (registry category directsolver).

class edelweissfe.linsolve.types.SuperLU[source]

superlu: SciPy’s own sparse LU, always available.

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

package: str = 'superlu': The registered name of the direct solver package (registry category directsolver).

class edelweissfe.linsolve.types.UMFPACK[source]

umfpack: UMFPACK through SciPy (SuperLU if scikit-umfpack is missing).

Methods

`adapterSettings`(configuration)	The settings of the adapter (`ComposedLinearSolver`) in the data lines every type shares (`LinearSolverSchema`), as its keyword arguments.
`create`(configuration)	The linear solver of these data lines: the tree behind the adapter.
`fromOptions`(options)	The linear solver of these options (the data lines, by name; unknown names raise).
`schema`
`tree`(configuration)	The composed solver of these data lines.

package: str = 'umfpack': The registered name of the direct solver package (registry category directsolver).

edelweissfe.linsolve.types.solveDumpsOf(configuration)[source]

The dumps the data lines ask for, or None.

Parameters:: configuration (LinearSolverSchema)
Return type:: SolveDumps | None

Created on Sat Feb 10 10:27:25 2018

@author: Matthias Neuner

edelweissfe.config.linsolve.createLinearSolver(linsolver, definitions)[source]

The linear solver a nonlinear solver’s linsolver= names: a *linearSolver block, else a type with its defaults, else (empty) the default linear solver.

Parameters:

linsolver (str) – The name of a block of definitions, or of a type; "" for the default.
definitions (dict) – The *linearSolver blocks by name (linearSolverDefinitionsFromInputFile()).

Return type:

LinearSolver

edelweissfe.config.linsolve.getDefaultLinSolver()[source]

Get the linear solver to use when a step declares no linsolver option.

Returns:: pardiso if its optional extension was built, else SciPy’s own sparse LU (superlu), with their defaults.
Return type:: LinearSolver

edelweissfe.config.linsolve.getLinSolverByName(linsolverName, options)[source]

Return a linear solver of the type linsolverName, configured with options.

Parameters:

linsolverName (str) – The type of the linear solver, case-insensitively (e.g. "pardiso", "blockamg", "custom"; registry category linsolver, see edelweissfe.linsolve.types).
options (dict) – The data lines of that type by name, e.g. {"reuseSymbolicFactorization": True} or {"script": "solver.py"}; {} for its defaults.

Returns:

A callable (A, b) -> x solving A x = b.

Return type:

LinearSolver

Raises:

edelweissfe.config.registry.RegistryLookupError – If no linear solver is registered under linsolverName.
ValueError – If options names an option the type does not have, or a value it cannot take.
ImportError – If the requested solver’s optional backend is not available in this installation (see getDefaultLinSolver(), which relies on it).

edelweissfe.config.linsolve.linearSolverDefinitionsFromInputFile(inputfile)[source]

The *linearSolver blocks of an input file, by name.

Raises:: ValueError – For a block named like a type (linsolver=<type> means that type with its defaults), a name used twice, a data line that is not key=value, or an option the type does not have.
Parameters:: inputfile (dict)
Return type:: CaseInsensitiveDict

The common base class every registered linsolver inherits, so the nonlinear solver can treat all of them uniformly instead of special-casing per capability.

Every linsolve/<name>/__init__.py’s createSolver(opts) factory returns an instance of a LinearSolver subclass, callable as (A, b) -> x. LinearSolver.setJournal() and LinearSolver.setModel() are part of that base contract with safe defaults, so a caller can call either on any registered solver unconditionally, whether or not that particular solver actually uses the information – a solver that needs more than the default overrides it; everything else inherits a default that does nothing harmful.

setModel() hands a solver the live FEModel and DofManager it is solving for. Everything a solver may want is derivable from those two objects (the field layout from the DofManager, node coordinates and element topology from the FEModel), so a solver that wants more than the default field-structure bookkeeping (e.g. the composed solver reading node coordinates for a rigid-body near null-space, or building a P1 topology map) keeps the references and computes what it needs on its own schedule – the driver need not know in advance what a solver might want. setFieldStructure() remains as a lower-level way in for callers that know the field-block layout directly but have no full model to hand over (e.g. an offline replay of a dumped system).

class edelweissfe.linsolve.base.DirectSolverBackend[source]

A sparse direct solver library (SuperLU, PARDISO, KLU, …): factorizes a matrix once, then solves for as many right-hand sides as needed by back substitution.

Not a linear solver of its own: the composed leaf DirectSolver applies it, \(z = A^{-1} r\), and the type=pardiso (etc.) linear solvers are that one leaf. Backends are registered in the directsolver registry category, and take their options as keyword arguments.

The contract every backend keeps: factorize() factorizes (once per set-up of its leaf), and reports a matrix it cannot factorize there; solveFactorized() only applies the factors, as often as asked – inside a Krylov solve, once per iteration – and solves with the factorization of the last factorize(), also when the assembly has since overwritten that matrix in place. A backend that keeps its symbolic factorization (the reordering) while only the values of the matrix change does so inside factorize(): one backend instance lives as long as its DirectSolver leaf, which releases the numeric factorization (releaseFactorization()) before it factorizes anew.

factorize(A)[source]

Factorize A (CSR) for the following solveFactorized() calls. A matrix that cannot be factorized (singular) raises numpy.linalg.LinAlgError.

Return type:: None

releaseFactorization()[source]

Free the numeric factorization of the last factorize() (nothing to free before the first), so that it is gone before the next factorize() computes another: never two at once. A symbolic factorization the backend reuses stays.

Return type:: None

solveFactorized(b)[source]: Solve A x = b with the factorization of the last factorize(), for one right-hand side (a vector) or several (the columns of a matrix), returned in the shape of b.

class edelweissfe.linsolve.base.LinearSolver[source]

Common base for every linsolver registry entry. Callable as (A, b) -> x.

Subclasses implement __call__(). setJournal(), setModel() and setFieldStructure() have safe defaults here so the nonlinear solver can call any of them on any solver without asking first which ones care.

setFieldStructure(fields)[source]

Receive the ordered field blocks of the DOF vector directly (in DOF order), bypassing setModel().

An escape hatch for a caller that knows the field-block layout but has no full FEModel/DofManager to hand over – e.g. an offline probe script driving a solver directly on a captured (A, b) system. A live run goes through setModel() instead, whose default implementation computes exactly this and stores it the same way; a solver that only ever needs the field-block structure (not the model/dofManager references themselves) does not need to care which path supplied it.

Only field-split solvers (e.g. a composed solver with OnField corrections) read it; every other solver simply ignores it.

Parameters:: fields (list[FieldBlock])
Return type:: None

setJournal(journal)[source]

Receive the shared Journal instance.

Default just stores it on self._journal for solvers that want to log through it; solvers with no logging needs simply never read the attribute.

Return type:: None

setModel(model, dofManager)[source]

Receive the live model and DOF manager this solver is being asked to solve for.

Called by the nonlinear solver whenever the equation system is (re)built – i.e. on the first solve and again after any AMR/connectivity change. The default implementation derives and stores the per-field DOF-block structure (name, DOF range, nodal dimension) any field-split solver needs – equivalent to calling setFieldStructure() with those blocks – and keeps model/dofManager themselves on self._model/self._dofManager for a solver that wants to go further (e.g. node coordinates for a rigid-body near null-space, or a P1 topology map) without the driver having to know about it.

A solver overriding this should normally call super().setModel(model, dofManager) first to get the field-structure bookkeeping for free, then do whatever else it needs.

Parameters:

model – The live FEModel.
dofManager – The live DofManager describing the current equation system’s layout.

Return type:

None

identification = 'LinearSolver': Display name for this solver’s own log lines; overridden by solvers with their own name.

lastSolveSummary: LinearSolveSummary | None = None: The most recent call’s diagnostics as a LinearSolveSummary, or None; read by the nonlinear solver right after each call.

needsPrunedSystemMatrix = False: Whether this solver wants the explicitly stored zeros compacted out of the multi-point-condensed system matrix before it solves. A direct factorization does (its reordering is sensitive to the extra structural entries); a Krylov solver with an AMG preconditioner does not (a stored zero is no strong connection: the same iterations, and the compaction costs a copy of the matrix and a pass over it in every Newton iteration); one that reuses its preconditioner must not have it (the pruned pattern changes on almost every iteration). Read by the “auto” default of pruneCondensedMatrixZeros; an explicit true/false in the input file overrides it.

reportsSolveSummary = False: Whether this solver reports per-solve diagnostics (iterations, residual) for an extra “linear solve” column in the nonlinear solver’s convergence table. If True, lastSolveSummary must be populated after every __call__().

The components

The building block of a composed linear solver: Preconditioner, an approximate inverse \(M^{-1}\) of a matrix \(A\), and what is handed to it – the SetUpContext of its (sub-)problem and the MatrixChange since its last set-up.

A composed solver is a tree of preconditioners. Each one solves its part of the problem and hands a smaller problem to its children: a Krylov solver calls its preconditioner once per iteration, a Schur complement solves the condensed dofs exactly and hands the remaining ones to its child, a diagonal scaling multiplies by \(D\) before and after its child. A Krylov solver is itself a preconditioner – an approximate inverse of \(A\) – so solvers nest like any other component.

A component implements what makes it what it is, in protected hooks:

Preconditioner._setUp() builds it for a matrix (and its children for their sub-problems);
Preconditioner._forget() drops what it built, and so declares what that is (the constructor calls it);
Preconditioner._update() takes a changed matrix and decides what to keep (by default: set up again);
Preconditioner._apply() applies it, \(z = M^{-1} r\) – the formula of the textbook;
Preconditioner.children() names its children, if it has any (the tree is walked through it);
Preconditioner.describe() says what it is, for the views (views).

It is used through the public methods of the base class, setUp(), update() and __call__() (a template method each): they record the call and its time in the component’s Statistics and call the hook. So a component’s _apply holds nothing but its formula, and a hook that calls another hook of its own component (a Krylov solver’s _setUp is its _update for a first matrix) counts as one call.

A component builds only through setUp(), also when an update decides to rebuild (it calls self.setUp). setUp() first forgets what was built before (forget(): the _forget() of the component and of its children), so that the old and the new hierarchies, factorizations and matrices, gigabytes for a large model, are never held at once. An update that does not rebuild keeps what was built, and what is reused on purpose is never forgotten: the selection of a Schur complement’s condensed dofs, a symbolic factorization, the statistics. A component with nothing built builds at its next update, whatever the update says.

class edelweissfe.linsolve.composition.base.Description(name, formula, latex='', parameters=<factory>, steps=())[source]

What a component is, independent of how it is shown: the views (views) draw it.

Parameters:

name (str)
formula (str)
latex (str)
parameters (dict)
steps (tuple)

name

The component’s name, e.g. "SchurComplement".

Type:: str

formula

Its approximate inverse in unicode, with {role} for the formula of each of its children(), e.g. "D · {inner} · D".

Type:: str

latex

The same in LaTeX, with {role} placeholders.

Type:: str

parameters

The settings that make it what it is, e.g. {"restart": 50}.

Type:: dict

steps

The steps of one application, in order.

Type:: tuple

class edelweissfe.linsolve.composition.base.MatrixChange(context, rebuildReason=None, newIncrement=False, residualNorm=0.0, previousResidualNorm=None)[source]

What changed since a component’s last set-up, and whether what it built must be built again.

The facts about the matrix (a new pattern, a resized system) are found by the composed solver (ComposedLinearSolver); a Krylov solver adds its own reasons (an increment, a stale preconditioner, see RebuildPolicy); a Schur complement adds its own when its condensed dofs change. Each component then decides in Preconditioner.update() what to keep.

context

The context of the changed matrix.

Type:: edelweissfe.linsolve.composition.base.SetUpContext

rebuildReason

Why what was built must be built again ("first", "model", "amr", "pattern", "increment", "stale", "always", "condensed dofs"), or None: only the values changed, and what was built for the previous matrix is reused.

Type:: str | None

newIncrement

Whether the right-hand side jumped, which marks a new increment of the nonlinear solver.

Type:: bool

residualNorm

\(\lVert b \rVert\) of this solve.

Type:: float

previousResidualNorm

\(\lVert b \rVert\) of the previous solve, None before the first one (or after a change of the model).

Type:: float | None

Attributes:

previousResidualNorm
rebuildReason
requiresRebuild: Whether what was built must be built again.

Parameters:

context (SetUpContext)
rebuildReason (str | None)
newIncrement (bool)
residualNorm (float)
previousResidualNorm (float | None)

Methods

`withContext`(context)	This change, for a child with its own context.
`withRebuild`(reason)	This change, asking for a rebuild for `reason` -- unless it asks for one already, whose reason stays.

withContext(context)[source]

This change, for a child with its own context.

Parameters:: context (SetUpContext)
Return type:: MatrixChange

withRebuild(reason)[source]

This change, asking for a rebuild for reason – unless it asks for one already, whose reason stays.

Parameters:: reason (str)
Return type:: MatrixChange

property requiresRebuild: bool: Whether what was built must be built again.

class edelweissfe.linsolve.composition.base.Preconditioner[source]

An approximate inverse \(M^{-1}\) of a matrix \(A\): built by setUp(), applied as \(z = M^{-1} r\).

See the module documentation. A subclass implements the hooks _setUp(), _apply(), _forget() if it builds anything, and, to reuse what it built, _update(), and calls super().__init__() in its constructor (which creates its statistics and declares what it builds, by _forget()).

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

applyRestricted(r_f)[source]

\(M^{-1} r_f\) on the dofs of fieldIndex alone: _applyRestricted(), recorded as an application.

Parameters:: r_f (ndarray)
Return type:: ndarray

builtMatrices()[source]

(name, scipy matrix) pairs of the matrices this component formed and keeps (for the matrix view); empty for most.

Return type:: tuple

children()[source]

The (role, child) pairs of this component, in the order the children are used (empty for a leaf). A component with children overrides this: the tree is walked through it (setModel(), finishSolve(), isFixedLinearOperator, …, and the views), and the role is the placeholder of the child in the component’s formula (Description.formula).

Return type:: tuple

condensedDofs()[source]

The dofs a Schur complement in this component condensed in the last solve, or None (for a degradation dump, which records them for an offline replay without the model).

Return type:: ndarray | None

describe()[source]

What this component is, for the views.

Return type:: Description

finishSolve()[source]

The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children. What the next solve may reuse (a preconditioner) is kept. Without this, a solver would hold a copy of the last matrix between solves – while the nonlinear solver assembles the next one.

Return type:: None

forget()[source]

Drop what this component and its children built (_forget()), so that nothing is held twice while it is built again. What they reuse on purpose stays (see the module documentation).

Return type:: None

setModel(model, dofManager)[source]

Take a new model (the first, or one after a change of the mesh), and hand it to the children. A component that derives something from the model (e.g. which dofs a marker selects) overrides this.

Return type:: None

setUp(A, context)[source]

Build this component for the matrix A (a SystemMatrix), and its children for their sub-problems: forget() what was built before, then _setUp(), recorded.

Parameters:: context (SetUpContext)
Return type:: None

setUpNotes()[source]

Short facts about what this component and its children built in the last set-up or update, one string each, for the “linear solve” column of the nonlinear solver’s table at verbosity "debug": at most 24 characters, the width of that column. Most components have none.

Return type:: tuple

solveRecord()[source]: What the last Krylov solve in this component did, a SolveRecord: a Krylov solver’s own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else None.

update(A, change)[source]

Take the changed matrix A: _update(), recorded.

Parameters:: change (MatrixChange)
Return type:: None

withRebuildReasons(change)[source]

change, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver’s RebuildPolicy: a new increment, a stale preconditioner). A component that keeps something built for the matrix above such a child (the \(D\) of a DiagonalScaling around a Krylov solve) asks for them first, so that it rebuilds together with the child.

Parameters:: change (MatrixChange)
Return type:: MatrixChange

property containsDirectSolver: bool: Whether this component is, or has among its children, a direct solver (a sparse LU factorization). Such a factorization wants the explicitly stored zeros pruned out of the system matrix (needsPrunedSystemMatrix).

property containsKrylovSolver: bool: Whether this component is, or has among its children, a Krylov solver. A Krylov solver records each solve (the “linear solve” column of the nonlinear solver); a direct solve does not.

property fieldIndex: int | None: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs. A combination (combinators) then forms the residual of that field alone and calls applyRestricted().

property isFixedLinearOperator: bool

true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds). A Krylov method that relies on a fixed preconditioner (LGMRES, CG) refuses one that is not.

Type:: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups

property reusesPreconditioner: bool: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse). A new sparsity pattern forces a rebuild that such a solver would have saved, so the nonlinear solver then does not prune the system matrix (needsPrunedSystemMatrix), which changes its pattern on almost every iteration.

statistics: What this component recorded while it ran, see Statistics.

class edelweissfe.linsolve.composition.base.SetUpContext(fields, scaling=None, coordinatesOf=None, midsideCornersOf=None, model=None, journal=None, log=None)[source]

What a component may need besides its matrix: the fields of its (sub-)problem and where their nodes are.

fields

The field blocks of the dofs of this (sub-)problem, contiguous and in dof order.

Type:: tuple

scaling

The diagonal scaling between the physical dofs and the variables of this (sub-)problem: the matrix a component gets is \(D A D\) of the physical one, None for \(D = I\). A near null-space, physical by nature (the rigid-body modes of a displacement field), is \(D^{-1}\) times its physical modes in these variables.

Type:: numpy.ndarray | None

coordinatesOf

fieldName -> (nNodes, dimension) node coordinates of the field, in its dof order, or None if unknown (an offline replay without a model).

Type:: object

midsideCornersOf

fieldName -> (isCorner, edgeEndpoints), the corner/midside topology of a quadratic field, or None.

Type:: object

model

The model, if there is one (markers read it), else None.

Type:: object

journal

The journal, if there is one, else None.

Type:: object

log

(level, message) -> None: the composed solver’s log, which shows a message if its verbosity reaches level ("warning", "info", "debug"), or None (then message() prints warnings only).

Type:: object

Attributes:

coordinatesOf
journal
log
midsideCornersOf
model
scaling
size: The number of dofs of this (sub-)problem.

Parameters:

fields (tuple)
scaling (ndarray | None)
coordinatesOf (object)
midsideCornersOf (object)
model (object)
journal (object)
log (object)

Methods

`coordinates`(fieldName)	The node coordinates of a field, or `None`.
`message`(level, text)	Log `text` at `level` through the composed solver (see `log`).
`midsideCorners`(fieldName)	The midside corners of a field (which nodes are corners, and the corners of each midside node), or `None`.
`ofCoarseDofs`(nCoarse)	The context of a coarse problem whose dofs are no nodes (the coefficients of the modes of aggregates, say): one field `"coarse"` of `nCoarse` scalar dofs, without coordinates, midside corners or scaling.
`ofCoarseNodes`(field, isCoarseNode)	The context of the coarse problem of a corner-node interpolation of this problem's one field `field`: the same field on its coarse nodes alone (`isCoarseNode`, per node of the field: the corner nodes, and the midside nodes the interpolation keeps as nodes of their own), its dofs numbered node by node.
`onlyField`(owner)	The one field of this (sub-)problem, for a component that is built on one field's matrix; raises, naming `owner`, if there are several.
`restrictedTo`(fields, dofs[, nodesOf])	The context of a subset `dofs` of this problem's dofs, whose fields are `fields` (as the Schur complement leaves them).
`restrictedToField`(block)	The context of one field block of this problem: its dofs renumbered from zero.
`scaledBy`(D)	The context of \(D \hat A D\), for a child that gets this problem scaled once more by `D`.
`scalingOf`(dofs)	The scaling of some dofs (ones without a scaling).

coordinates(fieldName)[source]

The node coordinates of a field, or None.

Parameters:: fieldName (str)
Return type:: ndarray | None

message(level, text)[source]

Log text at level through the composed solver (see log).

Parameters:

level (str)
text (str)

Return type:

None

midsideCorners(fieldName)[source]

The midside corners of a field (which nodes are corners, and the corners of each midside node), or None.

Parameters:: fieldName (str)

ofCoarseDofs(nCoarse)[source]

The context of a coarse problem whose dofs are no nodes (the coefficients of the modes of aggregates, say): one field "coarse" of nCoarse scalar dofs, without coordinates, midside corners or scaling. The journal and the log stay.

Parameters:: nCoarse (int)
Return type:: SetUpContext

ofCoarseNodes(field, isCoarseNode)[source]

The context of the coarse problem of a corner-node interpolation of this problem’s one field field: the same field on its coarse nodes alone (isCoarseNode, per node of the field: the corner nodes, and the midside nodes the interpolation keeps as nodes of their own), its dofs numbered node by node. The coarse nodes keep their coordinates and their scaling \(D_c\); the coarse field is linear, so it has no midside corners, and it is no field of the model.

Parameters:

field (FieldBlock)
isCoarseNode (ndarray)

Return type:

SetUpContext

onlyField(owner)[source]

The one field of this (sub-)problem, for a component that is built on one field’s matrix; raises, naming owner, if there are several.

Parameters:: owner (str)
Return type:: FieldBlock

restrictedTo(fields, dofs, nodesOf=None)[source]

The context of a subset dofs of this problem’s dofs, whose fields are fields (as the Schur complement leaves them). nodesOf: fieldName -> boolean mask of the nodes of the field that remain, to restrict the coordinates and midside corners.

Parameters:

fields (tuple)
dofs (ndarray)

Return type:

SetUpContext

restrictedToField(block)[source]

The context of one field block of this problem: its dofs renumbered from zero.

Parameters:: block (FieldBlock)
Return type:: SetUpContext

scaledBy(D)[source]

The context of \(D \hat A D\), for a child that gets this problem scaled once more by D.

Parameters:: D (ndarray)
Return type:: SetUpContext

scalingOf(dofs)[source]

The scaling of some dofs (ones without a scaling).

Parameters:: dofs (slice | ndarray)
Return type:: ndarray

property size: int: The number of dofs of this (sub-)problem.

class edelweissfe.linsolve.composition.base.Statistics(setUps=0, updates=0, applications=0, setUpSeconds=0.0, updateSeconds=0.0, applySeconds=0.0, sizes=<factory>)[source]

What a component recorded while it ran: filled at run time, shown by the run-statistics view.

The times are inclusive: a component’s apply time contains its children’s.

Parameters:

setUps (int)
updates (int)
applications (int)
setUpSeconds (float)
updateSeconds (float)
applySeconds (float)
sizes (dict)

setUps

How often it was built from scratch (Preconditioner.setUp()).

Type:: int

updates

How often it took a changed matrix (Preconditioner.update()), whether it rebuilt or not.

Type:: int

applications

How often it was applied (Preconditioner.__call__()).

Type:: int

setUpSeconds, updateSeconds, applySeconds: The wall time spent in each.

sizes

What the component wants to report about what it built, e.g. {"rows": 841260, "levels": 4}.

Type:: dict

Methods

record(kind, seconds)

Count one call of kind ("setUp", "update" or "apply") that took seconds.

record(kind, seconds)[source]

Count one call of kind ("setUp", "update" or "apply") that took seconds.

Parameters:

kind (str)
seconds (float)

Return type:

None

class edelweissfe.linsolve.composition.base.Step(operation, child=None)[source]

One step of a component’s application (Preconditioner._apply()), for the data-flow view.

operation

What the step computes, in the symbols of the component’s formula, e.g. "y_b = A_bb⁻¹ r_b".

Type:: str

child

The role of the child the step delegates to (a role of Preconditioner.children()), or None for a step the component does itself.

Type:: str | None

Attributes:

child

Parameters:

operation (str)
child (str | None)

edelweissfe.linsolve.composition.base.requireEachComponentOnce(root)[source]

Raise if one component (one instance) sits in two places of the tree: it is set up for one matrix, so the second place would overwrite what the first one built, e.g. AMG() bound to a name and put under two OnField. A combination that applies one correction several times (a symmetric sweep) lists it once among its children, so that is no second place.

Parameters:: root (Preconditioner)
Return type:: None

edelweissfe.linsolve.composition.base.requirePreconditioner(owner, role, child)[source]

Raise at construction if child, the role of owner, is not a component (e.g. a class instead of an instance, or a bare SciPy solver).

Parameters:

owner (str)
role (str)

Return type:

None

The outer Krylov solve: LGMRES or FGMRES, its tolerance per solve (EisenstatWalkerTolerance, FixedTolerance) and when it rebuilds its preconditioner (RebuildPolicy).

A Krylov method accelerates a basic iteration \(x \leftarrow x + M^{-1}(b - A x)\): from the preconditioned residuals it builds a subspace and takes the best solution in it. LGMRES is restarted GMRES, augmented at each restart with a few error approximations of the previous cycles. AMGCL’s implementation is used, right-preconditioned, so the residual it converges is the true one, \(\lVert b - A x \rVert / \lVert b \rVert\).

Inside Newton’s method, a linear solve need not be accurate while Newton is still far from converged: the forcing tolerance of Eisenstat and Walker asks for as much accuracy as the last Newton step earned. And the preconditioner need not be rebuilt for every Newton iteration: while the matrix changes little, the old one does nearly as well, and a rebuild costs more than the few extra iterations.

class edelweissfe.linsolve.composition.krylov.EisenstatWalkerTolerance(minimum=0.0001, maximum=0.001, gamma=0.9, alpha=1.618033988749895)[source]

The Eisenstat–Walker forcing tolerance per solve, see eisenstatWalkerTolerance().

Parameters:

minimum (float) – The bounds of \(\eta\).
maximum (float) – The bounds of \(\eta\).
gamma (float) – The parameters of “choice 2”.
alpha (float) – The parameters of “choice 2”.

Methods

forSolve(change)

The tolerance of a solve, from its right-hand side norm and the previous one.

forSolve(change)[source]

The tolerance of a solve, from its right-hand side norm and the previous one.

Parameters:: change (MatrixChange)
Return type:: float

class edelweissfe.linsolve.composition.krylov.FGMRES(preconditioner, tolerance=None, maxIterations=100, restart=30, rebuild=None, trueResidualRestarts=2, residualGapCompensation=None)[source]

Flexible GMRES (Saad 1993), AMGCL’s implementation: GMRES that keeps each preconditioned vector \(z_k = M_k^{-1} v_k\) and takes the solution from their span, \(x = Z_m y_m\), so the preconditioner may change from one iteration to the next – e.g. an inner Krylov solve to a loose tolerance. It stores twice the vectors of GMRES (\(V\) and \(Z\)).

Parameters:

preconditioner (Preconditioner) – See KrylovSolver.
tolerance – See KrylovSolver.
maxIterations (int) – See KrylovSolver.
rebuild (RebuildPolicy) – See KrylovSolver.
trueResidualRestarts (int) – See KrylovSolver.
residualGapCompensation (ResidualGapCompensation) – See KrylovSolver.
restart (int) – The iterations of a cycle (the vectors kept before a restart).

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

method = 'fgmres': AMGCL’s name of the method.

class edelweissfe.linsolve.composition.krylov.FixedTolerance(eta)[source]

The same relative tolerance \(\eta\) for every solve.

Parameters:: eta (float) – \(\lVert b - A x \rVert \leq \eta \lVert b \rVert\).

Methods

forSolve(change)

The tolerance of a solve.

forSolve(change)[source]

The tolerance of a solve.

Parameters:: change (MatrixChange)
Return type:: float

class edelweissfe.linsolve.composition.krylov.KrylovSolver(preconditioner, tolerance=None, maxIterations=100, rebuild=None, trueResidualRestarts=2, residualGapCompensation=None)[source]

An outer Krylov solve \(A x = b\), right-preconditioned by preconditioner, with AMGCL’s implementation of the method: what LGMRES and FGMRES share. As a component it is an approximate inverse of \(A\) like any other, solved to the tolerance of the solve – but not a fixed linear operator (isFixedLinearOperator).

Parameters:

preconditioner (Preconditioner) – \(M^{-1}\).
tolerance – The tolerance per solve: EisenstatWalkerTolerance (default) or FixedTolerance. May be replaced between solves (an offline replay pins it to the live one).
maxIterations (int) – The iteration cap of a solve (across restarts). A solve that reaches it returns its iterate.
rebuild – When the preconditioner is rebuilt, see RebuildPolicy.
trueResidualRestarts (int) – On a diagonally scaled matrix (DiagonalScaling(LGMRES(...))) the method converges the scaled residual, while the solve is accepted on the true one, which can be larger (stiff rows are weighted down by \(D\)): a pass that misses \(\eta\) on the true residual is continued from its iterate, at most this many times, with a tighter tolerance (0.01 times the previous, or see ResidualGapCompensation). On an unscaled matrix the residual the method converges is the true one, and a pass that misses \(\eta\) has hit the iteration cap, which a warm restart would only extend: no true-residual restarts there.
residualGapCompensation (ResidualGapCompensation) – A ResidualGapCompensation, or None (the default): only on a diagonally scaled matrix.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
method
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

children()[source]

The (role, child) pairs of this component, in the order the children are used (empty for a leaf). A component with children overrides this: the tree is walked through it (setModel(), finishSolve(), isFixedLinearOperator, …, and the views), and the role is the placeholder of the child in the component’s formula (Description.formula).

Return type:: tuple

describe()[source]

What this component is, for the views.

Return type:: Description

finishSolve()[source]

The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children. What the next solve may reuse (a preconditioner) is kept. Without this, a solver would hold a copy of the last matrix between solves – while the nonlinear solver assembles the next one.

Return type:: None

solveRecord()[source]

What the last Krylov solve in this component did, a SolveRecord: a Krylov solver’s own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else None.

Return type:: SolveRecord

withRebuildReasons(change)[source]

change, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver’s RebuildPolicy: a new increment, a stale preconditioner). A component that keeps something built for the matrix above such a child (the \(D\) of a DiagonalScaling around a Krylov solve) asks for them first, so that it rebuilds together with the child.

Parameters:: change (MatrixChange)
Return type:: MatrixChange

property containsKrylovSolver: bool: Whether this component is, or has among its children, a Krylov solver. A Krylov solver records each solve (the “linear solve” column of the nonlinear solver); a direct solve does not.

property isFixedLinearOperator: bool

true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds). A Krylov method that relies on a fixed preconditioner (LGMRES, CG) refuses one that is not.

Type:: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups

lastIterations

its outer iterations, its true relative residual, whether it stopped at the cap, and whether the preconditioner was rebuilt for it (and why).

Type:: Of the last solve

method = None: AMGCL’s name of the method.

property reusesPreconditioner: bool: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse). A new sparsity pattern forces a rebuild that such a solver would have saved, so the nonlinear solver then does not prune the system matrix (needsPrunedSystemMatrix), which changes its pattern on almost every iteration.

class edelweissfe.linsolve.composition.krylov.LGMRES(preconditioner, tolerance=None, maxIterations=100, restart=50, augmentationVectors=3, rebuild=None, trueResidualRestarts=2, residualGapCompensation=None)[source]

LGMRES (Baker, Jessup and Manteuffel 2005): GMRES restarted every restart iterations, each cycle augmented with augmentationVectors error approximations of the previous cycles, AMGCL’s implementation. Right-preconditioned: \(x = M^{-1} y\) with \(y\) from the Krylov space of \(A M^{-1}\), which is a Krylov space only for a fixed linear \(M^{-1}\) – so LGMRES refuses a variable preconditioner (a Krylov solve inside), for which FGMRES is made.

Parameters:

preconditioner (Preconditioner) – See KrylovSolver.
tolerance – See KrylovSolver.
maxIterations (int) – See KrylovSolver.
rebuild (RebuildPolicy) – See KrylovSolver.
trueResidualRestarts (int) – See KrylovSolver.
residualGapCompensation (ResidualGapCompensation) – See KrylovSolver.
restart (int) – The iterations of a cycle, \(m\) (the Krylov vectors kept before a restart).
augmentationVectors (int) – The error approximations of the previous cycles each cycle is augmented with, \(k\).

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

method = 'lgmres': AMGCL’s name of the method.

class edelweissfe.linsolve.composition.krylov.RebuildPolicy(stalenessFactor=3.0, rebuildOnNewIncrement=True, reusePreconditioner=False)[source]

When a Krylov solve rebuilds its preconditioner, beyond what the matrix itself demands (the first solve, a new model, a resized system, a new sparsity pattern).

Parameters:

stalenessFactor (float)
rebuildOnNewIncrement (bool)
reusePreconditioner (bool)

stalenessFactor

A solve that needs more outer iterations than this factor times those of the solve the preconditioner was built for (see outerIterationsDrifted()) has it rebuilt for the next solve; so has a solve that hit the iteration cap.

Type:: float

rebuildOnNewIncrement

Whether a new increment of the nonlinear solver rebuilds it.

Type:: bool

reusePreconditioner

Whether it is reused across solves at all. False (the default) builds it anew for every solve: a preconditioner reused from an earlier Jacobian costs more outer iterations than its build saves (Rebuilding the preconditioner for every solve). True keeps it until one of the other reasons asks for a rebuild.

Type:: bool

Methods

decide(change, stale)

change, with a rebuild added if this policy asks for one (in this order of reasons).

decide(change, stale)[source]

change, with a rebuild added if this policy asks for one (in this order of reasons).

Parameters:

change (MatrixChange)
stale (bool)

Return type:

MatrixChange

class edelweissfe.linsolve.composition.krylov.ResidualGapCompensation(safetyFactor=0.3, maxFactor=1000.0)[source]

For a Krylov solve on a diagonally scaled matrix \(\hat A = D A D\) (DiagonalScaling(LGMRES(...))): ask the first pass for \(s \, \eta / g\) instead of \(\eta\), with \(g\) the gap between the true and the scaled relative residual of the same iterate, learned from the previous solves, so that the pass lands on the true-residual target without a true-residual restart; such a restart then tightens by the measured gap instead of the fixed factor 0.01.

The Krylov method converges the scaled residual \(\lVert D (b - A x) \rVert\), while the solve is accepted on the true one, \(\lVert b - A x \rVert / \lVert b \rVert\); their ratio was measured from close to one (gradient-enhanced damage) to far larger (contact with multi-point constraints). Measured to save true-residual restarts and iterations per increment, but the uncompensated scaled solve over-solves, and Newton relies on that: per unit of simulated time the saving was largely given back. Hence it is opt-in.

Parameters:

safetyFactor (float) – \(s\): smaller is more accurate (closer to the uncompensated Newton path).
maxFactor (float) – The upper bound of \(g\): one non-converged solve with an outsized gap would otherwise compound, through the smoothing, into an ever tighter, unmeetable tolerance.

Methods

`firstPassTolerance`(eta)	The tolerance of the first pass: \(s \, \eta / g\).
`learn`(trueResidual, scaledResidual)	Take the gap of a solve's first pass into \(g\).
`restartTolerance`(eta, previousTolerance, ...)	The tolerance of a true-residual restart: aimed at \(s \, \eta / g\) of the current gap, and at least half the previous one, so that the restarts end.

firstPassTolerance(eta)[source]

The tolerance of the first pass: \(s \, \eta / g\).

Parameters:: eta (float)
Return type:: float

learn(trueResidual, scaledResidual)[source]

Take the gap of a solve’s first pass into \(g\).

Parameters:

trueResidual (float)
scaledResidual (float)

Return type:

None

restartTolerance(eta, previousTolerance, trueResidual, scaledResidual)[source]

The tolerance of a true-residual restart: aimed at \(s \, \eta / g\) of the current gap, and at least half the previous one, so that the restarts end.

Parameters:

eta (float)
previousTolerance (float)
trueResidual (float)

gap: \(g\), smoothed over the solves (an exponential moving average with weight 1/2).

class edelweissfe.linsolve.composition.krylov.SolveRecord(iterations, residual, eta, capHit, trueResidualRestarts, rebuildReason, scaledSystem)[source]

What a Krylov solve did, for the journal, the convergence table and a degradation dump.

iterations

Its outer iterations, true-residual restarts included.

Type:: int

residual

Its true relative residual \(\lVert b - A x \rVert / \lVert b \rVert\), in the physical variables of its (sub-)problem.

Type:: float

eta

Its tolerance.

Type:: float

capHit

Whether a pass stopped at the iteration cap short of its tolerance.

Type:: bool

trueResidualRestarts

How many warm-restarted passes followed the first (only on a diagonally scaled matrix).

Type:: int

rebuildReason

Why the preconditioner was rebuilt for it, or None (reused).

Type:: str | None

scaledSystem

Whether it ran on a diagonally scaled matrix (its own residual then is not the true one).

Type:: bool

Attributes:

converged: Whether the true residual met the tolerance (written so that NaN is a miss).

Parameters:

iterations (int)
residual (float)
eta (float)
capHit (bool)
trueResidualRestarts (int)
rebuildReason (str | None)
scaledSystem (bool)

property converged: bool: Whether the true residual met the tolerance (written so that NaN is a miss).

edelweissfe.linsolve.composition.krylov.eisenstatWalkerTolerance(residualNorm, previousResidualNorm, previousEta, newIncrement, minimum, maximum, gamma, alpha)[source]

The Eisenstat–Walker “choice 2” forcing tolerance, clamped and safeguarded.

\(\eta_k = \gamma (\lVert b_k \rVert / \lVert b_{k-1} \rVert)^\alpha\), with the classic safeguard against over-solving (a large tightening step is only trusted if the previous tolerance was already small), clamped to [minimum, maximum]. The safeguard only acts while \(\gamma \eta_{k-1}^\alpha\) exceeds 0.1, i.e. for a previous tolerance above about 0.26 at the classic parameters: never with a maximum below that. Falls back to maximum with no history to compare against (the first solve, or the one right after a new-increment jump – the ratio across that jump does not reflect Newton convergence and is not meaningful). A previous norm of exactly zero (typically an already-converged, e.g. linear-elastic, prior iterate) makes the ratio undefined rather than just uninformative, and falls back the same way.

Parameters:

residualNorm (float)
previousResidualNorm (float | None)
previousEta (float)
newIncrement (bool)
minimum (float)
maximum (float)
gamma (float)
alpha (float)

Return type:

float

edelweissfe.linsolve.composition.krylov.outerIterationsDrifted(iterations, referenceIterations, stalenessFactor)[source]

Whether a solve with a reused preconditioner needed so many more outer iterations than the solve it was built for (referenceIterations) that it should be rebuilt: more than both stalenessFactor times the reference and the reference plus STALENESS_MIN_EXCESS. The reference is the solve the preconditioner was built for, not the previous one: compared with its predecessor only, a slow drift across several solves never trips.

Parameters:

iterations (int)
referenceIterations (int)
stalenessFactor (float)

Return type:

bool

edelweissfe.linsolve.composition.krylov.STALENESS_MIN_EXCESS = 15: The least excess over the reference outer iteration count that makes a reused preconditioner stale, on top of the staleness factor: at small counts a factor alone fires on noise.

Symmetric diagonal scaling: DiagonalScaling.

\[M^{-1} = D \, \hat M^{-1} D, \qquad \hat M \text{ built on } \hat A = D A D, \qquad D = \operatorname{diag}(\lvert a_{ii} \rvert^{-1/2})\]

Every diagonal entry of \(\hat A\) is \(\pm 1\). The fields of a coupled system have very different units – the rows of a displacement field carry stiffnesses, those of a nonlocal damage field entries of order one – and AMG and a block Gauss–Seidel sweep compare matrix entries with each other (the strength-of-connection test \(\lvert a_{ij} \rvert \geq \theta \sqrt{\lvert a_{ii} a_{jj} \rvert}\), the couplings between fields, the Chebyshev eigenvalue bounds). Scaled to a unit diagonal, those comparisons measure the matrix, not the choice of units.

Where the scaling sits in the tree says what it acts on: SchurComplement(..., DiagonalScaling(Multiplicative([...]))) scales only what the block Gauss–Seidel sweep sees, while the outer Krylov solve still converges the true residual of \(A\).

class edelweissfe.linsolve.composition.diagonalscaling.DiagonalScaling(inner)[source]

\(M^{-1} = D \hat M^{-1} D\), with \(\hat M\) (the child inner) built on \(\hat A = D A D\), \(D = \operatorname{diag}(\lvert a_{ii} \rvert^{-1/2})\). See the module documentation.

When only the values of \(A\) change (no rebuild asked for), \(D\) of the last build is kept: any \(D\) gives the same system in other variables, and the child’s reused preconditioner was built for this one.

Parameters:

inner (Preconditioner) – The preconditioner \(\hat M^{-1}\) of the scaled matrix.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

children()[source]

The (role, child) pairs of this component, in the order the children are used (empty for a leaf). A component with children overrides this: the tree is walked through it (setModel(), finishSolve(), isFixedLinearOperator, …, and the views), and the role is the placeholder of the child in the component’s formula (Description.formula).

Return type:: tuple

describe()[source]

What this component is, for the views.

Return type:: Description

edelweissfe.linsolve.composition.diagonalscaling.diagonalScalingOf(absoluteDiagonal)[source]

\(D = \operatorname{diag}(\lvert a_{ii} \rvert^{-1/2})\) of a matrix, from its absoluteDiagonal \(\lvert a_{ii} \rvert\), for DiagonalScaling.

A zero diagonal entry cannot be scaled: \(1/\sqrt{0}\) is infinite, which would poison \(\hat A\) with non-finite values before the preconditioner build, and the failure would surface much later as a NaN correction (i.e. as a diverged solution, which reads like a physics problem rather than an unsupported system). Flooring the divisor is not a fix either – a 1e-300 floor merely turns inf into ~1e150 and blows the scaled operator up to ~1e149, which is finite and equally useless. A system with a structurally zero diagonal is a saddle-point system (e.g. the multiplier block of a Lagrangian constraint, which never writes its own diagonal), and smoothed-aggregation AMG on a zero block is not something these preconditioners do: it raises a ValueError, and names the offending dofs.

A diagonal entry that is NaN or infinite is no better: it makes \(D\), and with it a near null-space \(D^{-1} t\), non-finite, which AMGCL’s hierarchy build rejects with a RuntimeError – which nothing in the nonlinear solvers catches, whereas a NaN correction is a diverged solution and a cutback. Such a system is declined like a non-finite right-hand side: a FloatingPointError, which the composed solver turns into a declined solve (checked after the zero diagonal, which is the user’s error to see first).

Parameters:: absoluteDiagonal (ndarray) – \(\lvert a_{ii} \rvert\).
Returns:: The diagonal of \(D\).
Return type:: np.ndarray

Algebraic multigrid on one field: AMG, one AMGCL hierarchy, applied as one cycle per application.

Smoothed aggregation builds coarse levels from a near null-space: the error components a smoother such as Chebyshev does not reduce, because the operator has (almost) no stiffness in their direction. For a displacement field those are the rigid-body modes – 3 translations and 3 rotations in 3D; for a scalar field the constant. In the variables of a scaled block \(\hat A_{ii} = D_i A_{ii} D_i\) they are \(D_i^{-1}\) times the physical modes (see nullspace).

AMGCL stores a 3D displacement field’s hierarchy in node-wise 3x3 blocks by default – faster per cycle, but a block hierarchy takes the 3 translations as its near null-space at most. The rotations then need a coarse level of their own, see TwoLevel.

class edelweissfe.linsolve.composition.amg.AMG(amgclTree=None, backendBlockSize=None, backendPrecision=None, nearNullspace='rigidBodyModes', dropTolerance=0.0, lumpDroppedEntries=False)[source]

One AMGCL hierarchy on the matrix of a field, applied as one cycle: \(z = \mathrm{AMG}(A)^{-1} r\).

Reused while only the values of the matrix change: the hierarchy is built again only if a rebuild is asked for (see MatrixChange).

Parameters:

amgclTree (dict) – The AMGCL parameter tree, with EdelweissFE’s own key spectralRadiusSteps (see _VECTOR_FIELD_TREE); None: the default tree of the field’s dimension (defaultFieldTree()).
backendBlockSize (int) – The AMGCL backend: scalar entries (1) or node-wise blocks (2, 3), stored in "double" or "float"; None: the field’s default, see backendOf().
backendPrecision (str) – The AMGCL backend: scalar entries (1) or node-wise blocks (2, 3), stored in "double" or "float"; None: the field’s default, see backendOf().
nearNullspace (str | None) – What smoothed aggregation starts from: "rigidBodyModes" (the default; the rotations of a displacement field on the scalar backend need node coordinates, a block backend takes the translations), "translations", or None (AMGCL’s own); see backendOf().
dropTolerance (float) – Build the hierarchy on the filtered matrix, see droppedSmallEntries(); 0 (the default) builds it on the matrix itself.
lumpDroppedEntries (bool) – Build the hierarchy on the filtered matrix, see droppedSmallEntries(); 0 (the default) builds it on the matrix itself.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

describe()[source]

What this component is, for the views.

Return type:: Description

edelweissfe.linsolve.composition.amg.backendOf(owner, field, amgclTree, backendBlockSize=None, backendPrecision=None, nearNullspace='rigidBodyModes')[source]

The AMGCL backend a one-field component builds on for the field (AMG, TwoLevel): what it was given, the field’s default for what it was not – the one place of those defaults.

backendBlockSize: the field’s dimension for a 2D or 3D field (node-wise blocks: less index and value traffic per product than scalar entries), else 1;
backendPrecision: "float" on 3x3 blocks, which keeps the setup in double and stores only the matrices the cycle reads in float (measured faster at the same iterations), else "double" – on the scalar backend float would build the whole hierarchy in single precision, its spectral-radius estimates included;
nearNullspace, what smoothed aggregation starts from, in the variables of the scaled block (divided by \(D\), see nullspace): "rigidBodyModes" (the default) the field’s rigid-body modes as far as the backend and the field hold them – translations and rotations of a displacement field on the scalar backend (carriesRigidBodyModes()), the translations on a block backend (which takes the translations at most) and for any other vector field; "translations" the translations \(D^{-1} t\) (translationNullspace()); None none given, AMGCL’s own: constant blocks on a block backend, which are \(D\)-weighted translations, the constant on the scalar backend. A scalar field always takes AMGCL’s constant.

Raises, naming owner and the field, for blocks that do not match the field, and for a near null-space on a block backend without smoothed aggregation.

Returns:

(backendBlockSize, backendPrecision, nearNullspace), the near null-space as built: "rigidBodyModes", "translations" or None.

Return type:

tuple

Parameters:

owner (str)
amgclTree (dict)

edelweissfe.linsolve.composition.amg.buildFieldHierarchy(block, amgclTree, backend, scaling, coords, diagonalBlock)[source]

One AMGCL hierarchy for a field’s scaled diagonal block \(\hat A_{ii}\) (diagonalBlock, a callable that returns it, so that it is cut out only now), with the field’s near null-space.

The near null-space is the one backendOf() resolved: the rigid-body modes, with coords (else the translations), the translations, or none given (AMGCL’s constant blocks, or its constant).

Parameters:

block – The field’s FieldBlock.
amgclTree (dict) – The field’s AMGCL tree, with EdelweissFE’s own key spectralRadiusSteps (see _VECTOR_FIELD_TREE).
backend (tuple) – (backendBlockSize, backendPrecision, nearNullspace), see backendOf().
scaling – The field’s scaling \(D_i\).
coords – The field’s node coordinates for its rigid-body rotations, or None.
diagonalBlock – () -> scipy CSR, the scaled diagonal block.

Returns:

The hierarchy, applied by its applyPreconditioner.

Return type:

PyAMGCLSolver

edelweissfe.linsolve.composition.amg.checkBackend(owner, backendBlockSize, backendPrecision, nearNullspace='rigidBodyModes')[source]

Raise, naming owner, for a backend AMGCL does not have, at construction (None is the default of the field, see backendOf()): backendBlockSize 1 (scalar entries), 2 or 3 (node-wise blocks), backendPrecision "double" or "float" (2x2 blocks only in double), nearNullspace one of NEAR_NULLSPACES.

Parameters:: owner (str)
Return type:: None

edelweissfe.linsolve.composition.amg.defaultCornerNodeTwoLevelTree()[source]

The default AMGCL tree of the corner-node two-level cycle on a 3D displacement field, TwoLevel(coarseSpace=CornerNodeInterpolation()), and the tree the blockAMG() preset gives the rigid-body coarse level of its coarse solve (any other field keeps defaultFieldTree()): the default tree of a vector field (_VECTOR_FIELD_TREE) with the smoother of the corner-node two-level cycle (_CORNER_NODE_TWO_LEVEL_SMOOTHER). The corner-node two-level cycle reads only the smoother of its tree (relax, npre, npost), the rigid-body coarse level the whole tree. A copy, as defaultFieldTree().

Return type:: dict

edelweissfe.linsolve.composition.amg.defaultFieldTree(dimension)[source]

The default AMGCL preconditioner tree of a field, from its nodal dimension: _VECTOR_FIELD_TREE for a vector field, a Chebyshev-smoothed smoothed aggregation with AMGCL’s defaults for a scalar field (e.g. nonlocal damage), but the smoother’s margin higher 1.1 above the estimated spectral radius, the default of EdelweissFE’s Chebyshev smoother (see amgcl-chebyshev.hpp). A copy: the trees are nested, and a caller may adapt the one it gets. The backend the tree is built on (block size, precision) is chosen apart from it, see backendOf().

Parameters:: dimension (int)
Return type:: dict

edelweissfe.linsolve.composition.amg.droppedSmallEntries(A, dropTolerance, lumpDroppedEntries=False)[source]

The filtered matrix \(\tilde A\) of A that a preconditioner is built on: the off-diagonal entries that fail the relative drop test (isSignificant()) at \(\tau\) = dropTolerance dropped. The diagonal is always kept, so no row is emptied. Shared by AMG and TwoLevel.

Why: some operators (e.g. of RKPM) store many entries per row whose median off-diagonal is orders of magnitude below the diagonal – numerically far sparser than structurally. Smoothed aggregation then weighs every one of them per row, and the Galerkin products carry them to every coarse level. Only the preconditioner is built on \(\tilde A\); the Krylov solve still applies \(A\), so the converged solution is that of \(A\), and only the number of iterations to reach it changes. The criterion is scale-invariant and symmetric in \((i, j)\): a symmetric matrix stays symmetric, as smoothed aggregation and the Chebyshev smoother assume.

Parameters:

A (csr_matrix) – The matrix.
dropTolerance (float) – \(\tau\).
lumpDroppedEntries (bool) – Add the dropped entries of each row onto its diagonal instead of discarding them: the row sums are kept, \(\tilde A \mathbf 1 = A \mathbf 1\), so the constant near null-space vector survives (the “filtered matrix” of the AMG literature). For elasticity, whose off-diagonals are mostly negative, this weakens the diagonal dominance; whether that helps depends on the operator.

Return type:

csr_matrix

edelweissfe.linsolve.composition.amg.fieldTreeOf(owner, amgclTree, field, defaultTree)[source]

The AMGCL tree a one-field component builds on: amgclTree, or the component’s default for the field, defaultTree(field) (built only then); raises, naming owner and the field, if amgclTree is not a tree.

Parameters:

owner (str)
amgclTree (dict | None)

Return type:

dict

edelweissfe.linsolve.composition.amg.isSignificant(values, rows, cols, sqrtAbsDiagonal, tau)[source]

Whether the entries values at the positions (rows, cols) pass the relative drop test of AMG’s strength of connection,

\[\lvert a_{ij} \rvert \geq \tau \sqrt{\lvert a_{ii} a_{jj} \rvert},\]

with sqrtAbsDiagonal holding \(\sqrt{\lvert a_{ii} \rvert}\). The one definition of the test: the filtered matrix of droppedSmallEntries() and the drop and thinning of the approximate Schur complement (schurcomplement) use it. rows and cols broadcast against values: the row and the column of every stored entry of a sparse matrix, or rows[:, None] and cols[None, :] for a dense block.

Returns:

The keep mask, of the shape of values.

Return type:

np.ndarray

Parameters:

sqrtAbsDiagonal (ndarray)
tau (float)

edelweissfe.linsolve.composition.amg.preconditionerMatrix(A, dropTolerance, lumpDroppedEntries, statistics)[source]

The matrix a one-field preconditioner (AMG, TwoLevel) is built on: \(A\) itself, or with dropTolerance > 0 its filtered copy (droppedSmallEntries()), with the stored entries before and after recorded in the component’s statistics.

Parameters:

dropTolerance (float)
lumpDroppedEntries (bool)

Return type:

csr_matrix

edelweissfe.linsolve.composition.amg.NEAR_NULLSPACES = ('rigidBodyModes', 'translations', None): The near null-spaces AMG builds its hierarchy from, see backendOf().

A two-level cycle with a coarse space of its own: TwoLevel.

One cycle, for \(z \approx A^{-1} r\) with a prolongation \(P\) from the coarse space:

\[\begin{split}z &\leftarrow S^{\nu_1}(0, r) && \text{pre-smoothing (Chebyshev)} \\ z &\leftarrow z + P \, A_c^{-1} P^T (r - A z) && \text{coarse correction}, \quad A_c = P^T A P \\ z &\leftarrow S^{\nu_2}(z, r) && \text{post-smoothing}\end{split}\]

with \(A_c^{-1}\) approximated. The coarse space decides what the cycle can correct:

RigidBodyModesPerAggregate: all rigid-body modes of each aggregate of nodes (3 translations and 3 rotations in 3D), where an AMG hierarchy on node-wise 3x3 blocks takes the translations alone; \(A_c^{-1}\) is one V-cycle of an AMG hierarchy;
CornerNodeInterpolation: the linear interpolant of a quadratic field from its corner nodes (P1); \(A_c^{-1}\) is a component of its own on the corner-node field, by default DiagonalScaling(AMG()). For a 3D displacement field, DiagonalScaling(TwoLevel(coarseSpace=RigidBodyModesPerAggregate())) makes three levels: quadratic, linear, rigid-body modes of aggregates. On such a field the fine level smooths with the smoother of the corner-node two-level cycle, a Chebyshev polynomial of degree 8 instead of 16 (defaultCornerNodeTwoLevelTree()), and the blockamg preset gives the coarse level the same.

Why “two-level”: the component is the fine level and one coarse correction, and its coarse problem is solved inexactly. The coarse space need not be a coarser mesh. With RigidBodyModesPerAggregate it is spanned by aggregates of nodes, and the whole cycle is smoothed aggregation with all rigid-body modes on every level, the method Trilinos ML and MueLu use for elasticity (Vaněk, Mandel and Brezina, Computing 56, 1996). The first coarse level is built here, on the node-wise blocks of the fine level, with the modes scaled by the diagonal scaling and zero on the Dirichlet rows. The levels below are built by AMGCL, which takes the coefficients of the modes in the coarse space as its near null-space. prolongation="tentative" makes the first level plain aggregation. With CornerNodeInterpolation the coarse level is a coarser grid, the corner nodes (p-multigrid): a field of its own, which its coarse component solves like any field (cornernodetwolevelcycle).

The cycle runs in C++ (rigidbodymodestwolevelcycle, cornernodetwolevelcycle; the coarse component of CornerNodeInterpolation is called back once per cycle). Its steps, written out in Python, are Multiplicative([smoother, CoarseCorrection(P, coarseSolve), smoother]) (combinators). Both ways are kept on purpose: this one is faster per application, the spelled-out one shows the steps and takes any smoother or coarse space. Why, and what would end the exception: the “two-level cycle: two ways” section of doc/source/documentation/composedlinsolvers.rst.

class edelweissfe.linsolve.composition.twolevel.CoarseSpace[source]

What spans the coarse level of a TwoLevel cycle: it checks that it can be built for a field, and builds the cycle on it. The cycle’s own parameters (the fine smoother of the field’s tree, the backend) come from TwoLevel, as cycle, the keywords of cycleParameters().

Attributes:

name

Methods

`build`(field, matrix, context, cycle)	The two-level cycle on this coarse space for a field's matrix (SciPy CSR), built: a `TwoLevelCycle`.
`check`(field, context, cycle)	Raise, naming the field, if the cycle cannot be built on this coarse space for the field: before anything is built.
`children`()	The `(role, component)` pairs of the coarse space, children of the `TwoLevel` that uses it: its coarse solve, if that is a component (`CornerNodeInterpolation`); none if the cycle solves on its own.
`coarseSolveStep`()	The coarse correction of the cycle, for the data-flow view.
`defaultTree`(field)	The AMGCL tree of a `TwoLevel` cycle on this coarse space that is given none: the default tree of the field's dimension (`defaultFieldTree()`).
`describe`()	The settings of the coarse space, for the views.
`sizesOf`(cycle)	What a built cycle reports about its coarse level (for the run statistics).

build(field, matrix, context, cycle)[source]

The two-level cycle on this coarse space for a field’s matrix (SciPy CSR), built: a TwoLevelCycle.

Parameters:

field – The field.
matrix – Its (scaled) matrix.
context (SetUpContext) – Its context: the node coordinates, the midside corners, the scaling.
cycle (dict) – The cycle’s parameters, see cycleParameters().

check(field, context, cycle)[source]

Raise, naming the field, if the cycle cannot be built on this coarse space for the field: before anything is built.

Parameters:

context (SetUpContext)
cycle (dict)

Return type:

None

children()[source]

The (role, component) pairs of the coarse space, children of the TwoLevel that uses it: its coarse solve, if that is a component (CornerNodeInterpolation); none if the cycle solves on its own.

Return type:: tuple

coarseSolveStep()[source]

The coarse correction of the cycle, for the data-flow view.

Return type:: Step

defaultTree(field)[source]

The AMGCL tree of a TwoLevel cycle on this coarse space that is given none: the default tree of the field’s dimension (defaultFieldTree()).

Return type:: dict

describe()[source]

The settings of the coarse space, for the views.

Return type:: dict

sizesOf(cycle)[source]

What a built cycle reports about its coarse level (for the run statistics).

Return type:: dict

name = None: The name of the coarse space in the views.

class edelweissfe.linsolve.composition.twolevel.CornerNodeInterpolation(coarseSolve=None)[source]

The coarse space of the linear interpolant of a quadratic field from its corner nodes (P1): TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve)). Needs the field’s corner/midside topology.

The prolongation interpolates the scaled variables of the field, \(P = D^{-1} P_{geo} D_c\), so that the coarse operator \(A_c = P^T \hat A P\) is the P1 operator of the corner-node field, scaled by \(D_c\) (see cornernodetwolevelcycle). The coarse problem is set up as that field: the same field on the coarse nodes – the corner nodes, and a midside node with a free dof between two constrained ones of its edge’s corners (issue #223), which no interpolation from them could correct – their coordinates, \(D_c\) as its scaling, an identity row for each Dirichlet coarse dof (ofCoarseNodes()).

The coarse solve belongs to this coarse space for its whole life: each build sets it up anew on that build’s \(A_c\). So one CornerNodeInterpolation serves one TwoLevel; under two, its coarse solve would sit in two places of the tree, which the composed solver refuses.

Each level has its own tree. A TwoLevel on this coarse space that is given no tree smooths a 3D displacement field with a weaker Chebyshev polynomial than an AMG hierarchy (defaultTree(), degree 8 instead of 16); any other field keeps the default tree of its dimension. The coarse solve has the tree of its own component; the blockamg preset gives it the same smoother (blockAMGFieldPreconditioner()).

Parameters:: coarseSolve (Preconditioner) – The coarse solve \(A_c^{-1}\), a component applied once per cycle. None: DiagonalScaling(AMG()), which works on any field. For a 3D displacement field DiagonalScaling(TwoLevel(coarseSpace=RigidBodyModesPerAggregate())) is faster (The corner-node two-level cycle with scaled weights and a composed coarse solve).

Methods

`build`(field, matrix, context, cycle)	Build the cycle in C++, which returns \(A_c\), set up the coarse solve on it, and hand the coarse solve to the cycle.
`check`(field, context, cycle)	It needs the field's midside corners.
`children`()	The `(role, component)` pairs of the coarse space, children of the `TwoLevel` that uses it: its coarse solve, if that is a component (`CornerNodeInterpolation`); none if the cycle solves on its own.
`coarseSolveStep`()	The coarse correction of the cycle, for the data-flow view.
`defaultTree`(field)	For a 3D displacement field, the smoother of the corner-node two-level cycle (`defaultCornerNodeTwoLevelTree()`), where it was measured; any other field keeps the default tree of its dimension, as on any coarse space.
`describe`()	The settings of the coarse space, for the views.
`sizesOf`(cycle)	What a built cycle reports about its coarse level (for the run statistics).

build(field, matrix, context, cycle)[source]

Build the cycle in C++, which returns \(A_c\), set up the coarse solve on it, and hand the coarse solve to the cycle.

Parameters:

context (SetUpContext)
cycle (dict)

check(field, context, cycle)[source]

It needs the field’s midside corners.

Parameters:

context (SetUpContext)
cycle (dict)

Return type:

None

children()[source]

The (role, component) pairs of the coarse space, children of the TwoLevel that uses it: its coarse solve, if that is a component (CornerNodeInterpolation); none if the cycle solves on its own.

Return type:: tuple

coarseSolveStep()[source]

The coarse correction of the cycle, for the data-flow view.

Return type:: Step

defaultTree(field)[source]

For a 3D displacement field, the smoother of the corner-node two-level cycle (defaultCornerNodeTwoLevelTree()), where it was measured; any other field keeps the default tree of its dimension, as on any coarse space.

Return type:: dict

describe()[source]

The settings of the coarse space, for the views.

Return type:: dict

sizesOf(cycle)[source]

What a built cycle reports about its coarse level (for the run statistics).

Return type:: dict

name = 'corner-node interpolation': The name of the coarse space in the views.

class edelweissfe.linsolve.composition.twolevel.RigidBodyModesPerAggregate(prolongation='smoothed', coarseHierarchyPrecision='double')[source]

The coarse space of all rigid-body modes per aggregate of nodes: TwoLevel(coarseSpace=RigidBodyModesPerAggregate()).

Parameters:

prolongation (str) – "smoothed" (\(P = (I - \omega D_F^{-1} A_F) P_{\mathrm{tent}}\)) or "tentative" (\(P_{\mathrm{tent}}\), the modes orthonormalized per aggregate).
coarseHierarchyPrecision (str) – The value type of the AMG hierarchy that solves the coarse level: "double" or "float".

Methods

`build`(field, matrix, context, cycle)	The two-level cycle on this coarse space for a field's matrix (SciPy CSR), built: a `TwoLevelCycle`.
`check`(field, context, cycle)	Rigid-body modes are those of a 2D or 3D displacement field with node coordinates (`carriesRigidBodyModes()`), on a block backend (the scalar backend takes them as the near null-space of its hierarchy itself, `AMG()`).
`children`()	The `(role, component)` pairs of the coarse space, children of the `TwoLevel` that uses it: its coarse solve, if that is a component (`CornerNodeInterpolation`); none if the cycle solves on its own.
`coarseSolveStep`()	The coarse correction of the cycle, for the data-flow view.
`defaultTree`(field)	The AMGCL tree of a `TwoLevel` cycle on this coarse space that is given none: the default tree of the field's dimension (`defaultFieldTree()`).
`describe`()	The settings of the coarse space, for the views.
`sizesOf`(cycle)	What a built cycle reports about its coarse level (for the run statistics).

build(field, matrix, context, cycle)[source]

The two-level cycle on this coarse space for a field’s matrix (SciPy CSR), built: a TwoLevelCycle.

Parameters:

field – The field.
matrix – Its (scaled) matrix.
context (SetUpContext) – Its context: the node coordinates, the midside corners, the scaling.
cycle (dict) – The cycle’s parameters, see cycleParameters().

check(field, context, cycle)[source]

Rigid-body modes are those of a 2D or 3D displacement field with node coordinates (carriesRigidBodyModes()), on a block backend (the scalar backend takes them as the near null-space of its hierarchy itself, AMG()).

Parameters:

context (SetUpContext)
cycle (dict)

Return type:

None

describe()[source]

The settings of the coarse space, for the views.

Return type:: dict

sizesOf(cycle)[source]

What a built cycle reports about its coarse level (for the run statistics).

Return type:: dict

name = 'rigid-body modes per aggregate': The name of the coarse space in the views.

class edelweissfe.linsolve.composition.twolevel.TwoLevel(coarseSpace, amgclTree=None, backendBlockSize=None, backendPrecision=None, dropTolerance=0.0, lumpDroppedEntries=False)[source]

A two-level cycle on a field’s matrix with the coarse space coarseSpace; see the module documentation.

Reused while only the values of the matrix change, like AMG.

Parameters:

coarseSpace (CoarseSpace) – A CoarseSpace: RigidBodyModesPerAggregate or CornerNodeInterpolation.
amgclTree (dict) – The field’s AMGCL tree: its relax (a Chebyshev polynomial), npre and npost are the fine smoother; with RigidBodyModesPerAggregate its coarsening also builds the coarse level, and the whole tree the AMG hierarchy below. None: the coarse space’s default (CoarseSpace.defaultTree()), the default tree of the field’s dimension (defaultFieldTree()) for RigidBodyModesPerAggregate, and for CornerNodeInterpolation on a 3D displacement field the smoother of the corner-node two-level cycle (defaultCornerNodeTwoLevelTree()).
backendBlockSize (int) – The storage of the fine operator and of \(P\), \(P^T\), as for AMG; None: the field’s default (backendOf()).
backendPrecision (str) – The storage of the fine operator and of \(P\), \(P^T\), as for AMG; None: the field’s default (backendOf()).
dropTolerance (float) – Build the cycle (its smoother, coarse level and hierarchy) on the filtered matrix, see droppedSmallEntries(); 0 (the default) builds it on the matrix itself.
lumpDroppedEntries (bool) – Build the cycle (its smoother, coarse level and hierarchy) on the filtered matrix, see droppedSmallEntries(); 0 (the default) builds it on the matrix itself.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

children()[source]

The (role, child) pairs of this component, in the order the children are used (empty for a leaf). A component with children overrides this: the tree is walked through it (setModel(), finishSolve(), isFixedLinearOperator, …, and the views), and the role is the placeholder of the child in the component’s formula (Description.formula).

Return type:: tuple

describe()[source]

What this component is, for the views.

Return type:: Description

Condensing hard dofs exactly: SchurComplement, a block LU factorization around a child, and the condensation it is built from, SchurComplementCondensation.

Some dofs make a system hard for AMG: the dofs of embedded bars tied to their host elements by bond springs (\(u_{bar} \approx \sum_i N_i u_i\), which an aggregation coarse space cannot represent), or the dofs of a softening zone. We call these dofs \(b\) (the condensed dofs: “bars”) and all others \(k\) (the kept, or remaining, dofs). Ordered this way, the system reads

\[\begin{split}\begin{bmatrix} A_{bb} & A_{bk} \\ A_{kb} & A_{kk} \end{bmatrix} \begin{bmatrix} x_b \\ x_k \end{bmatrix} = \begin{bmatrix} r_b \\ r_k \end{bmatrix}.\end{split}\]

The block \(A_{bb}\) is factorized exactly (a direct solver). Eliminating \(x_b\) leaves the remaining dofs with the Schur complement \(S = A_{kk} - A_{kb} A_{bb}^{-1} A_{bk}\) – for bars, the concrete stiffened by the bars as seen through their bond – on which the child of SchurComplement is built (e.g. a field-split AMG). Gaussian elimination by blocks is the block LU factorization of \(A\), and its inverse is a preconditioner of the outer Krylov solve:

\[y_b = A_{bb}^{-1} r_b, \quad z_k = M_{\tilde S}^{-1} (r_k - A_{kb} y_b), \quad z_b = A_{bb}^{-1} (r_b - A_{bk} z_k),\]

with \(M_{\tilde S}^{-1}\) the child (schurComplementSolve) on an approximation \(\tilde S\) of \(S\). With \(\tilde S = S\) and an exact sweep it is \(A^{-1}\).

It solves with \(A_{bb}\) twice. The default (blockFactors = "upper") inverts 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\), with one solve:

\[z_k = M_{\tilde S}^{-1} r_k, \quad z_b = A_{bb}^{-1} (r_b - A_{bk} z_k):\]

the remaining dofs first, then the condensed ones follow them (for bars: the bars follow the concrete). Its last step solves the rows of \(b\) exactly, so after the first outer iteration the residual of these rows is zero, and the outer Krylov solve iterates on the rows of \(k\) alone, with the operator \(S \tilde S^{-1}\): the same operator as with the full block LU factorization. With exact \(S\), GMRES converges in two iterations instead of one; with \(\tilde S\) and AMG, the outer iterations hardly grow, and every one saves a solve with \(A_{bb}\).

\(S\) is dense wherever \(A_{bb}^{-1}\) is, e.g. along a whole bar. Its entries decay quickly away from the diagonal, though (a bar on its bond is a beam on an elastic foundation: a load is felt only a few host elements away), so \(\tilde S\) keeps only the significant ones:

a banded component of \(A_{bb}\) (a bar, ordered along the bar) is cut into overlapping windows of nodes along it; for the dofs of a window’s core, \(A_{bb}^{-1}\) is replaced by the inverse of the window’s block (restricted additive Schwarz); a short component is a single, exact window,
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, in batches of right-hand sides,
entries \(|\tilde S_{ij}| < \tau \sqrt{|\tilde S_{ii} \tilde S_{jj}|}\) are dropped (dropTolerance).

Which dofs are condensed is decided by markers, the same ones that mark elements for adaptive refinement (edelweissfe.adaptivity.marking): \(b\) holds all dofs of the interior nodes of the marked elements, i.e., of the nodes that have at least one element, all of them marked:

bar node:   bar, bond                -> all marked -> condensed
host node:  bond, concrete, concrete -> not all    -> kept

For embedded bars, mark the bar and the bond elements. For a softening zone, mark the elements of a perElement fieldOutput above a threshold, with a halo so that the zone’s own nodes become interior. Only the dofs of nodes can be condensed, and always all fields of a node together.

class edelweissfe.linsolve.composition.schurcomplement.SchurComplement(condensedDofMarkers=(), schurComplementSolve=None, condensedBlockSolve=None, approximation='local', dropTolerance=0.0001, windowNodes=50, haloNodes=None, blockFactors='upper')[source]

The block LU preconditioner around the Schur complement of the condensed dofs; see the module documentation.

Every solve selects the condensed dofs again and sets up condensedBlockSolve on \(A_{bb}\) (factorizes it); \(\tilde S\) and the child are rebuilt when a rebuild is asked for, and when the condensed dofs change. If \(A_{bb}\) cannot be factorized (a bond that lost all stiffness), that solve does without the condensation: the child is built on the whole \(A\).

Parameters:

condensedDofMarkers (list) – The markers whose marked elements’ interior nodes are condensed (e.g. an element set marker); the dofs they select are condensedDofs().
schurComplementSolve (Preconditioner) – \(M_{\tilde S}^{-1}\), built on \(\tilde S\).
condensedBlockSolve (DirectSolver) – The exact solve with \(A_{bb}\), a DirectSolver (default: SuperLU). The zones corrected exactly in \(\tilde S\) are factorized by a second one of the same package and options.
approximation (str) – How \(\tilde S\) is formed, see SchurComplementCondensation.
dropTolerance (float) – How \(\tilde S\) is formed, see SchurComplementCondensation.
windowNodes (int) – How \(\tilde S\) is formed, see SchurComplementCondensation.
haloNodes (int) – How \(\tilde S\) is formed, see SchurComplementCondensation.
blockFactors (str) – "upper" (the default): only the block upper triangular factor \(U\), one solve with \(A_{bb}\); "both": both factors of the block LU factorization, two solves with \(A_{bb}\).

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

builtMatrices()[source]

(name, scipy matrix) pairs of the matrices this component formed and keeps (for the matrix view); empty for most.

Return type:: tuple

children()[source]

The (role, child) pairs of this component, in the order the children are used (empty for a leaf). A component with children overrides this: the tree is walked through it (setModel(), finishSolve(), isFixedLinearOperator, …, and the views), and the role is the placeholder of the child in the component’s formula (Description.formula).

Return type:: tuple

condensedDofs()[source]

The dofs a Schur complement in this component condensed in the last solve, or None (for a degradation dump, which records them for an offline replay without the model).

Return type:: ndarray | None

describe()[source]

What this component is, for the views.

Return type:: Description

finishSolve()[source]

The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children. What the next solve may reuse (a preconditioner) is kept. Without this, a solver would hold a copy of the last matrix between solves – while the nonlinear solver assembles the next one.

Return type:: None

setModel(model, dofManager)[source]

Take a new model (the first, or one after a change of the mesh), and hand it to the children. A component that derives something from the model (e.g. which dofs a marker selects) overrides this.

Return type:: None

setUpNotes()[source]

Short facts about what this component and its children built in the last set-up or update, one string each, for the “linear solve” column of the nonlinear solver’s table at verbosity "debug": at most 24 characters, the width of that column. Most components have none.

Return type:: tuple

class edelweissfe.linsolve.composition.schurcomplement.SchurComplementCondensation(condensedDofMarkers, condensedBlockSolve=None, approximation='local', dropTolerance=0.0001, windowNodes=50, haloNodes=None)[source]

The condensed dofs \(b\), the factorization of \(A_{bb}\), the approximate Schur complement \(\tilde S\) of all other dofs \(k\), and the preconditioner; see the module documentation.

SchurComplement uses it in three steps: update() on every solve (select \(b\), split \(A\); the Schur complement then factorizes \(A_{bb}\) with its condensedBlockSolve), approximateSchurComplement() whenever it rebuilds its child (built on \(\tilde S\) instead of \(A\)), and blockLU() in every outer iteration.

Parameters:

condensedDofMarkers (list) – The markers whose marked elements’ interior nodes are condensed: they select condensedDofs.
condensedBlockSolve (DirectSolver) – The exact solve with (a part of) \(A_{bb}\), a DirectSolver: here the one of the zones corrected exactly (see _exactCorrection()), set up on their block of \(A_{bb}\) and applied to batches of right-hand sides. Default: SuperLU.
approximation (str) – "local": \(\tilde S\) as described in the module documentation; "none": \(\tilde S = A_{kk}\) (the coupling through \(b\) is ignored).
dropTolerance (float) – The relative size \(\tau\) below which entries of \(\tilde S\) are dropped.
windowNodes (int) – The nodes of a window’s core along a chain, and of its halo on either side (None: as many as the core).
haloNodes (int | None) – The nodes of a window’s core along a chain, and of its halo on either side (None: as many as the core).

Attributes:

approximation: How \(\tilde S\) is formed: "local" or "none" (\(\tilde S = A_{kk}\)).

Methods

`approximateSchurComplement`(A)	\(\tilde S\), see the module documentation.
`forget`()	Drop what was formed from the last matrix: its blocks \(A_{bb}\), \(A_{bk}\), \(A_{kb}\), and what `approximateSchurComplement()` formed from them (the factorization of the zones).
`hasCondensedDofs`()	Whether anything is condensed (after the last `update()`).
`remainingFieldBlocks`(blocks)	The field blocks of the remaining dofs: every field without its condensed dofs, in the same order (and contiguous, as the remaining dofs keep their order).
`remainingNodes`(block)	Which nodes of the field block `block` (of the whole system) remain, are not condensed: a boolean mask over the field's nodes, in its dof order.
`setModel`(model, dofManager)	Take the model and its dofs: check that what the markers read exists, and find the dofs of each node.
`update`(A)	Select the condensed dofs and cut the blocks \(A_{bb}\), \(A_{bk}\), \(A_{kb}\) out of `A`.

approximateSchurComplement(A)[source]

\(\tilde S\), see the module documentation. Also counts its diagonal entries that are not positive (nonPositiveDiagonalEntries): for a symmetric positive definite \(A\) there are none, so any means that the condensed zone did not take all the softening.

Parameters:: A (csr_matrix) – The system matrix of the last update().
Return type:: csr_matrix

forget()[source]

Drop what was formed from the last matrix: its blocks \(A_{bb}\), \(A_{bk}\), \(A_{kb}\), and what approximateSchurComplement() formed from them (the factorization of the zones). What is reused for the next matrix stays: the selection of the condensed dofs, and where the entries of the blocks are in a matrix of the same pattern.

Return type:: None

hasCondensedDofs()[source]

Whether anything is condensed (after the last update()).

Return type:: bool

remainingFieldBlocks(blocks)[source]

The field blocks of the remaining dofs: every field without its condensed dofs, in the same order (and contiguous, as the remaining dofs keep their order). A field whose dofs are all condensed (e.g. the rotations of embedded beams) has no block.

Parameters:: blocks (list)
Return type:: list

remainingNodes(block)[source]

Which nodes of the field block block (of the whole system) remain, are not condensed: a boolean mask over the field’s nodes, in its dof order. Condensation takes whole nodes, and a field’s dofs are node-major.

Parameters:: block (FieldBlock)
Return type:: ndarray

setModel(model, dofManager)[source]

Take the model and its dofs: check that what the markers read exists, and find the dofs of each node.

Return type:: None

update(A)[source]

Select the condensed dofs and cut the blocks \(A_{bb}\), \(A_{bk}\), \(A_{kb}\) out of A.

Parameters:: A (csr_matrix) – The system matrix, CSR in canonical format.
Return type:: None

property approximation: str

"local" or "none" (\(\tilde S = A_{kk}\)).

Type:: How \(\tilde S\) is formed

condensedDofs: The condensed dofs b and all other (remaining) dofs k, sorted.

nonPositiveDiagonalEntries: The number of diagonal entries of the last \(\\tilde S\) that are not positive (see approximateSchurComplement()).

selectedAllDofs: Whether the markers selected every dof at the last update(), which leaves nothing for the field split: then nothing is condensed (a direct solver is the better choice for such a system).

class edelweissfe.linsolve.composition.schurcomplement.Window(dofs, core)[source]

Create new instance of Window(dofs, core)

Methods

`count`(value, /)	Return number of occurrences of value.
`index`(value[, start, stop])	Return first index of value.

Parameters:

dofs (ndarray)
core (ndarray)

core: ndarray: The positions of its core within dofs.

dofs: ndarray

its core and the halos on either side, in order along the chain.

Type:: The dofs of the window, positions in b

edelweissfe.linsolve.composition.schurcomplement.blockLU(r, condensation, schurComplementSolve, condensedBlockSolve, blockFactors)[source]

One application of the block LU preconditioner, see the module documentation:

\[y_b = A_{bb}^{-1} r_b, \qquad z_k = M_{\tilde S}^{-1} (r_k - A_{kb} y_b), \qquad z_b = A_{bb}^{-1} (r_b - A_{bk} z_k)\]

or, with blockFactors="upper", \(z_k = M_{\tilde S}^{-1} r_k\) (the remaining dofs first, then the condensed ones follow them: the bars follow the concrete). The application of SchurComplement.

Parameters:

r (ndarray) – The vector to precondition (all dofs).
condensation – The SchurComplementCondensation: the dofs \(b, k\) and the blocks \(A_{bk}, A_{kb}\).
schurComplementSolve – \(M_{\tilde S}^{-1}\), applied to a vector of the remaining dofs \(k\).
condensedBlockSolve – \(A_{bb}^{-1}\), applied to a vector of the condensed dofs \(b\).
blockFactors (str) – "upper" or "both", see the module documentation.

Return type:

ndarray

edelweissfe.linsolve.composition.schurcomplement.BAND_STORAGE_MAX_FRACTION = 0.25: A window’s block is solved in LAPACK band storage while its band (lower + upper + 1) holds at most this fraction of its columns (_solveTransposed()): a window along a bar has a band of a few dofs, so the banded solve costs a small fraction of the dense one, while for a short window, whose band reaches across it, the band storage no longer pays.

edelweissfe.linsolve.composition.schurcomplement.CHAIN_BANDWIDTH_DIVISOR = 4: A component of \(A_{bb}\) is a chain (a bar) if its bandwidth, in the reverse Cuthill-McKee ordering, is at most a window’s core divided by this: a chain couples each node only to a few neighbours along it, while a zone of a 2D or 3D mesh couples it across the zone’s width.

edelweissfe.linsolve.composition.schurcomplement.EXACT_CORRECTION_BATCH = 256: The number of right-hand sides per batch of the exact correction of a zone (_exactCorrection()) …

edelweissfe.linsolve.composition.schurcomplement.EXACT_CORRECTION_MAX_ENTRIES = 33554432: … at most, and fewer when the zones are so large that a batch would exceed this many entries.

edelweissfe.linsolve.composition.schurcomplement.PARALLEL_BLOCK_ENTRIES = 1048576: The split of \(A\) into its blocks and the drop of the small entries of \(\tilde S\) work on blocks of about this many entries in parallel threads, so that the temporaries of each block stay small. The threads write into arrays of the calling thread: in the free-threaded interpreter, an array made by a worker thread and dropped by another is freed only when that worker runs again.

edelweissfe.linsolve.composition.schurcomplement.THINNING_FRACTION = 0.1: The correction of every window, and of every batch of a zone, is thinned against the diagonal of \(A_{kk}\) at this fraction of the drop tolerance before the parts are summed, which bounds their memory; the final drop is against the diagonal of \(\tilde S\) itself. A fraction well below 1 keeps what the final drop keeps, also an entry to which several parts contribute.

An exact solve as a component: DirectSolver, \(z = A^{-1} r\) by a sparse LU factorization.

A direct solver alone is the composed solver of type=pardiso (superlu, klu, …): one leaf, set up (i.e. factorized) on every solve and applied once. Inside a larger tree it solves a block exactly, e.g. the condensed dofs of a Schur complement, factorized once per set-up and applied in every outer iteration: every backend factorizes in factorize() and only back-substitutes afterwards.

class edelweissfe.linsolve.composition.direct.DirectSolver(package='superlu', **options)[source]

\(z = A^{-1} r\), exactly: one of the registered direct solver packages (superlu, umfpack, pardiso, klu, …; see DirectSolverBackend), factorized once per set-up – a singular matrix raises numpy.linalg.LinAlgError there – and applied by back substitution.

The backend is created here, at construction, and lives as long as this component: a backend that keeps its symbolic factorization while only the values of the matrix change (PARDISO’s reuseSymbolicFactorization) keeps it from one set-up to the next, while the numeric factorization of the last matrix is released before the next one is computed (_forget()). A backend whose optional extension is not installed raises ImportError here.

Parameters:

package (str) – The registered name of the direct solver package (registry category directsolver).
**options – The options of the package, e.g. reuseSymbolicFactorization=True for pardiso.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

describe()[source]

What this component is, for the views.

Return type:: Description

property containsDirectSolver: bool: Whether this component is, or has among its children, a direct solver (a sparse LU factorization). Such a factorization wants the explicitly stored zeros pruned out of the system matrix (needsPrunedSystemMatrix).

solver: The factorizing backend of this component.

Combining preconditioners: Multiplicative, Additive, OnField, CoarseCorrection, and the simplest smoother, Jacobi.

They are the vocabulary of Schwarz methods, and methods which look unrelated share their structure:

a correction on one field is \(R_i^T M_i^{-1} R_i\), with \(R_i\) the restriction to the field’s dofs and \(M_i^{-1}\) a preconditioner of its diagonal block \(A_{ii} = R_i A R_i^T\): OnField(name, M_i);
block Gauss–Seidel over the fields is their multiplicative combination, Multiplicative([OnField("displacement", M_u), OnField("nonlocal damage", M_d)]), and block Jacobi their additive one, Additive([...]); a symmetric sweep visits the fields forward and back, see gaussSeidelOrder();
a two-level cycle is a multiplicative combination of a smoother, a coarse correction and the smoother again, Multiplicative([Jacobi(), CoarseCorrection(P, coarseSolve), Jacobi()]).

This is also the fast way. A combination asks each correction which field it acts on (fieldIndex), and computes the residual of a correction on field \(i\) only on that field’s rows, from the fields corrected so far:

\[r_i - \sum_{j:\, z_j \neq 0} A_{ij} z_j\]

with the blocks \(A_{ij}\) formed once per matrix in one threaded pass (fieldSplit()) and read only as products. A correction on all dofs gets the residual of all fields.

class edelweissfe.linsolve.composition.combinators.Additive(corrections)[source]

Corrections of the same residual, summed: \(z = \sum_k M_k^{-1} r\). Over corrections on single fields, this is block Jacobi.

Parameters:

corrections (list | Callable) – The preconditioners \(M_k^{-1}\), all set up for the same matrix; or a function context -> list that makes them at the first set-up.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	`("M1", correction)`, `("M2", ...)`: each distinct correction once, in the order of first use (none while the corrections are still to be made per field).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

describe()[source]

What this component is, for the views.

Return type:: Description

class edelweissfe.linsolve.composition.combinators.CoarseCorrection(prolongation, coarseSolve)[source]

The coarse-grid correction \(z = P A_c^{-1} P^T r\) with the Galerkin coarse operator \(A_c = P^T A P\). With an exact coarse solve, \(P A_c^{-1} P^T A\) is the \(A\)-orthogonal projection onto the range of \(P\).

\(P\) is built at a set-up and kept until the next rebuild (see requiresRebuild). On a change of the values alone it recomputes \(A_c = P^T A P\) and hands it to the coarse solve as an update, which keeps what it can (an AMG its hierarchy). TwoLevel keeps its whole cycle instead, \(A_c\) included, until a rebuild.

The coarse solve is a component like any other, set up on \(A_c\) with a context that says what the coarse dofs are. A callable that builds \(P\) knows that, and returns it with \(P\): for the P1 interpolant, context.ofCoarseNodes(field, isCorner) – the field on its corner nodes, their coordinates – on which DiagonalScaling(TwoLevel(coarseSpace=RigidBodyModesPerAggregate())) works as on any displacement field; for coarse dofs that are no nodes, context.ofCoarseDofs(nCoarse). A fixed \(P\) says nothing about its columns: its coarse dofs are one field of scalar dofs, context.ofCoarseDofs.

\(P\) acts on the matrix the correction gets. Under a DiagonalScaling that is \(\hat A = D A D\), whose variables are \(D^{-1} u\): an interpolation of displacements is then \(D^{-1} P_{geo} D_c\) (the scaling is context.scaling), as CornerNodeInterpolation builds it. A \(P\) written out keeps the Dirichlet rows, which the native cycle decouples instead: a \(P\) with zero columns for the Dirichlet dofs would leave empty rows in \(A_c\), which a DiagonalScaling of the coarse solve cannot scale.

\(P^T r\), \(P z_c\) and \(P^T A P\) are AMGCL’s threaded products (PyAMGCLProlongation), the same as SciPy’s up to the order of the sums.

With smoothers around it, Multiplicative([smoother, CoarseCorrection(P, coarseSolve), smoother]) is a two-level cycle spelled out: any \(P\), any coarse solve, but slower than the C++ cycle TwoLevel, whose smoother reads a compact single-precision copy of the matrix; both are kept on purpose (see the “two-level cycle: two ways” and “which form to use” sections of doc/source/documentation/composedlinsolvers.rst). For production runs, TwoLevel.

Parameters:

prolongation – \(P\), a sparse matrix; or a callable (A, context) -> (P, coarseContext) that builds it for the matrix and says what its coarse dofs are.
coarseSolve (Preconditioner) – \(A_c^{-1}\), e.g. a DirectSolver, an AMG, or a two-level cycle of its own.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

children()[source]

The (role, child) pairs of this component, in the order the children are used (empty for a leaf). A component with children overrides this: the tree is walked through it (setModel(), finishSolve(), isFixedLinearOperator, …, and the views), and the role is the placeholder of the child in the component’s formula (Description.formula).

Return type:: tuple

describe()[source]

What this component is, for the views.

Return type:: Description

class edelweissfe.linsolve.composition.combinators.Jacobi(omega=0.6666666666666666)[source]

Damped Jacobi, the simplest smoother: \(z = \omega \operatorname{diag}(A)^{-1} r\).

A smoother for learning and for small problems; for production runs the Chebyshev smoother of TwoLevel and AMG (the “which form to use” section of doc/source/documentation/composedlinsolvers.rst).

Parameters:

omega (float) – The damping. Damped Jacobi damps the fast error modes only for \(0 < \omega < 2 / \lambda_{max}\), \(\lambda_{max}\) the largest eigenvalue of \(\operatorname{diag}(A)^{-1} A\). The default 2/3 is the textbook value for the Laplacian; quadratic elements have a larger \(\lambda_{max}\) and need a smaller damping (testfile TwoLevelSpelledOut).

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

describe()[source]

What this component is, for the views.

Return type:: Description

class edelweissfe.linsolve.composition.combinators.Multiplicative(corrections)[source]

Successive corrections, each one of the residual the previous ones left:

\[z \leftarrow M_1^{-1} r, \qquad z \leftarrow z + M_k^{-1} (r - A z), \quad k = 2, \ldots\]

A correction on one field (OnField) gets the residual of that field alone, see the module documentation.

Parameters:

corrections (list | Callable) – The preconditioners \(M_k^{-1}\), all set up for the same matrix; one may appear several times (it is set up once). Or a function context -> list that makes them at the first set-up, when the fields are known.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	`("M1", correction)`, `("M2", ...)`: each distinct correction once, in the order of first use (none while the corrections are still to be made per field).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

describe()[source]

What this component is, for the views.

Return type:: Description

finishSolve()[source]

The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children. What the next solve may reuse (a preconditioner) is kept. Without this, a solver would hold a copy of the last matrix between solves – while the nonlinear solver assembles the next one.

Return type:: None

class edelweissfe.linsolve.composition.combinators.OnField(fieldName, inner)[source]

A correction on the dofs of one field: \(z = R^T M^{-1} R r\), with \(R\) the restriction to the field and \(M\) built on its diagonal block \(A_{ff} = R A R^T\).

The block comes from the field split of the matrix (fieldSplit(), one threaded pass for all fields), and is cut out only if inner builds on it: a preconditioner reused for the new values never needs it.

Parameters:

fieldName (str) – The field.
inner (Preconditioner) – \(M^{-1}\), the preconditioner of the field’s block.

Attributes:

containsDirectSolver: Whether this component is, or has among its children, a direct solver (a sparse LU factorization).
containsKrylovSolver: Whether this component is, or has among its children, a Krylov solver.
fieldIndex: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs.
isFixedLinearOperator: Whether \(z = M^{-1} r\) is one fixed linear operator between two set-ups: true for a component whose children all are; a Krylov solve is not (its result depends on \(r\) nonlinearly, through the tolerance and the subspace it builds).
reusesPreconditioner: Whether this component is, or has among its children, a Krylov solver that keeps its preconditioner across solves (RebuildPolicy with reuse).

Methods

`__call__`(r)	Apply this component, \(z = M^{-1} r\): `_apply()`, recorded.
`applyRestricted`(r_f)	\(M^{-1} r_f\) on the dofs of `fieldIndex` alone: `_applyRestricted()`, recorded as an application.
`builtMatrices`()	`(name, scipy matrix)` pairs of the matrices this component formed and keeps (for the matrix view); empty for most.
`children`()	The `(role, child)` pairs of this component, in the order the children are used (empty for a leaf).
`condensedDofs`()	The dofs a Schur complement in this component condensed in the last solve, or `None` (for a degradation dump, which records them for an offline replay without the model).
`describe`()	What this component is, for the views.
`finishSolve`()	The solve is done: release what only this solve needed (e.g. the field blocks of its matrix, which the next solve forms anew from its own matrix), and tell the children.
`forget`()	Drop what this component and its children built (`_forget()`), so that nothing is held twice while it is built again.
`setModel`(model, dofManager)	Take a new model (the first, or one after a change of the mesh), and hand it to the children.
`setUp`(A, context)	Build this component for the matrix `A` (a `SystemMatrix`), and its children for their sub-problems: `forget()` what was built before, then `_setUp()`, recorded.
`setUpNotes`()	Short facts about what this component and its children built in the last set-up or update, one string each, for the "linear solve" column of the nonlinear solver's table at verbosity `"debug"`: at most 24 characters, the width of that column.
`solveRecord`()	What the last Krylov solve in this component did, a `SolveRecord`: a Krylov solver's own, else that of the first child that has one (e.g. the Krylov solve under a diagonal scaling), else `None`.
`update`(A, change)	Take the changed matrix `A`: `_update()`, recorded.
`withRebuildReasons`(change)	`change`, with the reasons for a rebuild that this component and its children add on their own (a Krylov solver's `RebuildPolicy`: a new increment, a stale preconditioner).

children()[source]

The (role, child) pairs of this component, in the order the children are used (empty for a leaf). A component with children overrides this: the tree is walked through it (setModel(), finishSolve(), isFixedLinearOperator, …, and the views), and the role is the placeholder of the child in the component’s formula (Description.formula).

Return type:: tuple

describe()[source]

What this component is, for the views.

Return type:: Description

property fieldIndex: int | None: The field this component acts on, \(R\) in \(z = R^T M^{-1} R r\), as the index of the field in its context (known after setUp()), or None for all dofs. A combination (combinators) then forms the residual of that field alone and calls applyRestricted().

edelweissfe.linsolve.composition.combinators.gaussSeidelOrder(fieldCorrections, sweeps=1, symmetric=False)[source]

The corrections of block Gauss–Seidel sweeps over fields, in the order a Multiplicative applies them: sweeps forward sweeps, each followed by a backward one if symmetric – which starts one field before the last, since the last field was just corrected with everything else current. For two fields u, d: [u, d], symmetric [u, d, u], two symmetric sweeps [u, d, u, u, d, u] (a field corrected twice in a row is corrected twice: the second correction is one of the residual the first left).

Parameters:

fieldCorrections (list)
sweeps (int)
symmetric (bool)

Return type:

list

edelweissfe.linsolve.composition.combinators.multiplicativeCorrections(blocks, corrections, r)[source]

Successive corrections, each of the residual the previous ones left: \(z \leftarrow z + R_k^T M_k^{-1} R_k (r - A z)\), with the residual of correction \(k\) formed only on its field’s rows (see the module documentation). The application of Multiplicative.

Parameters:

blocks (FieldSplitBlocks) – The field blocks \(A_{ij}\) (their products with vectors).
corrections (list) – (i, apply) per correction, in order: the index i of its field (None: all dofs) and its \(M_k^{-1}\), a callable on the residual of its dofs.
r (ndarray) – The vector to precondition.

Returns:

\(z\).

Return type:

np.ndarray

Presets: plain functions that build a composed solver.

A preset is nothing but a function returning a tree, so there is no preset machinery: read it to see what it builds, copy it into a solver script to change it.

edelweissfe.linsolve.composition.presets.blockAMG(rigidBodyModesPerAggregate=True, cornerNodeInterpolation=True, tolerance=None, maxIterations=100, sweeps=1, symmetric=False, rebuild=None, krylovOnScaledSystem=False, trueResidualRestarts=2, residualGapCompensation=None, dropTolerance=0.0, lumpDroppedEntries=False)[source]

The blockamg linear solver (formerly the monolithic BlockAMGSolver, whose numbers it reproduces): block Gauss–Seidel over the fields of the diagonally scaled matrix, inside LGMRES:

LGMRES(
    DiagonalScaling(
        Multiplicative([OnField(f, M_f) for every field f])   # M_f: blockAMGFieldPreconditioner(f)
    )
)

with each field’s preconditioner chosen for it by measurement (blockAMGFieldPreconditioner()): the two-level cycle on the corner nodes for a quadratic 3D displacement field, on the rigid-body modes for a linear one, one AMG hierarchy for any other field.

The corrections are made at the first set-up, for the fields of the model (Multiplicative takes a function of the context). A tree, not a solver: a solver script binds it, linearSolver = blockAMG().

Parameters:

rigidBodyModesPerAggregate (bool) – Whether a linear 3D displacement field, and the coarse field of the corner-node interpolation, get the two-level cycle on their rigid-body modes (else one AMG hierarchy, on the translations).
cornerNodeInterpolation (bool) – Whether a quadratic 3D displacement field gets the two-level cycle on its corner nodes (else the case of a linear one).
tolerance – The tolerance per solve (default: the Eisenstat–Walker forcing tolerance).
maxIterations (int) – The outer iteration cap.
sweeps (int) – The block Gauss–Seidel sweeps, see gaussSeidelOrder().
symmetric (bool) – The block Gauss–Seidel sweeps, see gaussSeidelOrder().
rebuild (RebuildPolicy) – When the preconditioner is rebuilt.
krylovOnScaledSystem (bool) – Run LGMRES on the scaled matrix, DiagonalScaling(LGMRES(Multiplicative(...))), with its true-residual trueResidualRestarts and residualGapCompensation (see KrylovSolver).
trueResidualRestarts (int) – See KrylovSolver; only with krylovOnScaledSystem.
residualGapCompensation (ResidualGapCompensation) – See KrylovSolver; only with krylovOnScaledSystem.
dropTolerance (float) – Build every field’s preconditioner on its filtered block, see droppedSmallEntries().
lumpDroppedEntries (bool) – Build every field’s preconditioner on its filtered block, see droppedSmallEntries().

Return type:

Preconditioner

edelweissfe.linsolve.composition.presets.blockAMGFieldPreconditioner(field, context, rigidBodyModesPerAggregate=True, cornerNodeInterpolation=True, dropTolerance=0.0, lumpDroppedEntries=False, amgclTree=None)[source]

The preconditioner of a field in the blockamg preset:

a quadratic 3D displacement field (with node coordinates): the two-level cycle on its corner nodes, TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve)), whose coarse field (linear) gets the next case, under a diagonal scaling of its own and with the smoother of the corner-node two-level cycle (defaultCornerNodeTwoLevelTree());
a linear 3D displacement field (with node coordinates): the two-level cycle on its rigid-body modes, TwoLevel(coarseSpace=RigidBodyModesPerAggregate()) (linearDisplacementFieldPreconditioner());
any other field: one AMG hierarchy, AMG().

amgclTree is the field’s AMGCL tree (None: the component’s default); dropTolerance and lumpDroppedEntries reach every level. A field is quadratic if it has midside nodes (isQuadratic(), asked only with cornerNodeInterpolation); an element whose topology the corner-node interpolation does not know (e.g. a Tet10) is an error then, cornerNodeInterpolation=False gives its field the rigid-body cycle.

Parameters:

rigidBodyModesPerAggregate (bool)
cornerNodeInterpolation (bool)
dropTolerance (float)
lumpDroppedEntries (bool)
amgclTree (dict)

edelweissfe.linsolve.composition.presets.isQuadratic(field, context)[source]

Whether a field has midside nodes (midsideCorners()); an element whose topology the corner-node interpolation does not know raises.

Return type:: bool

edelweissfe.linsolve.composition.presets.linearDisplacementFieldPreconditioner(rigidBodyModesPerAggregate, amgclTreeOfTwoLevel=None, amgclTreeOfAMG=None, dropTolerance=0.0, lumpDroppedEntries=False)[source]

The preconditioner of a linear 3D displacement field in the blockamg preset, a field of its own or the coarse field of the corner-node interpolation: the two-level cycle on its rigid-body modes, TwoLevel(coarseSpace=RigidBodyModesPerAggregate(), amgclTree=amgclTreeOfTwoLevel), or with rigidBodyModesPerAggregate false one AMG hierarchy on the translations, AMG(amgclTree=amgclTreeOfAMG). Each tree None: the component’s default.

Parameters:

rigidBodyModesPerAggregate (bool)
amgclTreeOfTwoLevel (dict)
amgclTreeOfAMG (dict)
dropTolerance (float)
lumpDroppedEntries (bool)

Views of a composed solver, all drawn from the components’ describe() and the tree of their children():

structure(): the tree as text, for the terminal and the journal (also print(solver));
formula(): the composed approximate inverse, in unicode or LaTeX;
dataFlow(): one application as a graphviz graph (DOT text, or SVG if the graphviz package is installed);
runStatistics(): the tree annotated with what each component recorded while it ran;
spy(): the sparsity patterns of the matrices the components were built on (matplotlib);
html(): all of the above in one self-contained page (also shown inline in Jupyter).

edelweissfe.linsolve.composition.views.dataFlow(component, svg=False)[source]

One application as a graphviz graph: each component a box around its steps, a step that delegates leading into its child’s box and back out. DOT text; with svg, rendered by the graphviz package (raises ImportError without it).

Parameters:: svg (bool)
Return type:: str

edelweissfe.linsolve.composition.views.formula(component, latex=False)[source]

The composed approximate inverse, each child’s formula put in the place of its role, e.g. LGMRES(A, D · (Rᵀ AMG(nonlocal damage)⁻¹ R ∘ Rᵀ TwoLevel(displacement; rigid-body modes)⁻¹ R) · D)⁻¹.

Parameters:: latex (bool)
Return type:: str

edelweissfe.linsolve.composition.views.html(component, standalone=True)[source]

The structure, formula, run statistics and data flow in one self-contained HTML page (the data flow as SVG if the graphviz package is installed, else as DOT text).

Parameters:: standalone (bool)
Return type:: str

edelweissfe.linsolve.composition.views.runStatistics(component)[source]

The tree annotated with what each component recorded (inclusive times: a component’s contain its children’s).

Return type:: str

edelweissfe.linsolve.composition.views.spy(component, A=None, markersize=0.05)[source]

The sparsity patterns of A (if given) and of the matrices the components were built on, side by side (matplotlib). Returns the figure.

Parameters:: markersize (float)

edelweissfe.linsolve.composition.views.structure(component)[source]

The tree as text: one line per component, its parameters, and the role of each child.

Example:

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)

Return type:: str

edelweissfe.linsolve.composition.views.writeHtml(component, path)[source]

Write html() to path.

Parameters:: path (str)
Return type:: None

The machinery

The adapter, the matrix the tree is set up on, the dumps, and the native cycles behind TwoLevel.

The matrix handed down a composed solver: SystemMatrix.

A component gets its matrix as a SystemMatrix, not as a bare SciPy matrix, for one reason: speed. A diagonal scaling \(\hat A = D A D\) followed by a split into field blocks \(\hat A_{ij}\) would cost a chain of serial SciPy passes over \(A\) – a scaled copy, then a row and column slice per block. Instead, SystemMatrix.scaled() only records \(D\), and SystemMatrix.fieldSplit() forms \(\hat A\) and all its blocks in one threaded native pass (PyAMGCLEquilibratedSystem). The same pass also keeps \(A\) itself as the operator of an outer Krylov solve (SystemMatrix.krylovOperator()), so a field split under a diagonal scaling under a Krylov solver reads \(A\) exactly once per solve.

A leaf that needs the matrix itself (a direct solver, the Schur complement) asks for SystemMatrix.matrix().

class edelweissfe.linsolve.composition.systemmatrix.FieldBlockMatrix(blocksOf, i)[source]

The diagonal block \(\hat A_{ii}\) of a field split as a SystemMatrix, cut out only when it is first used: a child that reuses its preconditioner never needs it, and the split itself is formed only when the block or a product with it is first asked for.

Parameters:

blocksOf – () -> FieldSplitBlocks, the field split of the whole matrix.
i (int) – The index of the field.

Attributes:

A
nnz: The stored entries of \(A\) (a proxy for its sparsity pattern).
shape

Methods

`absoluteDiagonal`()	\(\lvert \hat a_{ii} \rvert\), threaded.
`diagonal`()	\(\hat a_{ii}\): the diagonal of \(A\), scaled, without forming \(\hat A = D A D\).
`fieldSplit`(bounds)	The block as a split of one field: its products come from the split of the whole matrix, not from a copy.
`krylovOperator`()	\(\hat A\) as the operator of AMGCL's Krylov solvers: \(A\) itself (applied with `onUnscaledOperator`), or \(D A D\) for a Krylov solve on a scaled matrix, `DiagonalScaling(LGMRES(...))`.
`matrix`()	\(\hat A\) itself, for a component built on the whole matrix (AMG, a direct solver): \(A\), or \(D A D\) formed in the threaded pass of a split into one block (`fieldSplit()`), which also keeps \(A\) as the Krylov operator.
`releaseFormedMatrices`()	Drop what was formed from this matrix in a solve, the native field split and the operator of a Krylov solve (formed again when next asked for), and keep only \(A\): for a matrix a component keeps between solves (the \(\tilde S\) of a Schur complement).
`scaled`(D)	\(D \hat A D\), without forming it (see `fieldSplit()`).
`useAsKrylovOperator`()	Say that a Krylov solve runs on this matrix: a field split formed afterwards keeps \(\hat A\) itself as the operator of that solve, in its one threaded pass (see `krylovOperator()`).

fieldSplit(bounds)[source]

The block as a split of one field: its products come from the split of the whole matrix, not from a copy.

Parameters:: bounds (list)
Return type:: FieldSplitBlocks

class edelweissfe.linsolve.composition.systemmatrix.FieldSplitBlocks(native, bounds)[source]

The field blocks \(\hat A_{ij}\) of a (scaled) matrix, formed in one threaded pass: the diagonal blocks, and products with the coupling blocks, as a block Gauss–Seidel sweep needs them.

Parameters:

native – The PyAMGCLEquilibratedSystem that holds them.
bounds (list) – The (start, stop) dof range of each field.

Methods

`couplingTimes`(i, j, x)	\(\hat A_{ij} x\), for \(i \neq j\).
`diagonalBlock`(i)	\(\hat A_{ii}\) as a new SciPy CSR matrix.
`diagonalBlockTimes`(i, x)	\(\hat A_{ii} x\).

couplingTimes(i, j, x)[source]

\(\hat A_{ij} x\), for \(i \neq j\).

Parameters:

i (int)
j (int)
x (ndarray)

Return type:

ndarray

diagonalBlock(i)[source]

\(\hat A_{ii}\) as a new SciPy CSR matrix.

Parameters:: i (int)
Return type:: csr_matrix

diagonalBlockTimes(i, x)[source]

\(\hat A_{ii} x\).

Parameters:

i (int)
x (ndarray)

Return type:

ndarray

class edelweissfe.linsolve.composition.systemmatrix.SystemMatrix(A, scaling=None, _unscaled=None)[source]

A matrix \(A\), as handed to a component: optionally with a diagonal scaling, \(\hat A = D A D\).

See the module documentation: \(\hat A\) is formed on demand, together with its field blocks, in one threaded pass.

Parameters:

A – The matrix (SciPy sparse, square). Made canonical (CSR, sorted columns, no duplicates) if needed.
scaling (np.ndarray | None) – \(D\) as a vector, or None for \(D = I\).
_unscaled (SystemMatrix | None)

Attributes:

nnz: The stored entries of \(A\) (a proxy for its sparsity pattern).
shape

Methods

`absoluteDiagonal`()	\(\lvert \hat a_{ii} \rvert\), threaded.
`diagonal`()	\(\hat a_{ii}\): the diagonal of \(A\), scaled, without forming \(\hat A = D A D\).
`fieldSplit`(bounds)	The field blocks \(\hat A_{ij}\) for the field ranges `bounds`, formed with \(\hat A\) in one threaded pass at the first request and kept with this matrix: every component handed this matrix (the corrections of a combination, each on its field) gets the same blocks.
`krylovOperator`()	\(\hat A\) as the operator of AMGCL's Krylov solvers: \(A\) itself (applied with `onUnscaledOperator`), or \(D A D\) for a Krylov solve on a scaled matrix, `DiagonalScaling(LGMRES(...))`.
`matrix`()	\(\hat A\) itself, for a component built on the whole matrix (AMG, a direct solver): \(A\), or \(D A D\) formed in the threaded pass of a split into one block (`fieldSplit()`), which also keeps \(A\) as the Krylov operator.
`releaseFormedMatrices`()	Drop what was formed from this matrix in a solve, the native field split and the operator of a Krylov solve (formed again when next asked for), and keep only \(A\): for a matrix a component keeps between solves (the \(\tilde S\) of a Schur complement).
`scaled`(D)	\(D \hat A D\), without forming it (see `fieldSplit()`).
`useAsKrylovOperator`()	Say that a Krylov solve runs on this matrix: a field split formed afterwards keeps \(\hat A\) itself as the operator of that solve, in its one threaded pass (see `krylovOperator()`).

absoluteDiagonal()[source]

\(\lvert \hat a_{ii} \rvert\), threaded.

Return type:: ndarray

diagonal()[source]

\(\hat a_{ii}\): the diagonal of \(A\), scaled, without forming \(\hat A = D A D\). Of a field block (FieldBlockMatrix) it is the diagonal of the block, which is cut out of the field split once and kept.

Return type:: ndarray

fieldSplit(bounds)[source]

The field blocks \(\hat A_{ij}\) for the field ranges bounds, formed with \(\hat A\) in one threaded pass at the first request and kept with this matrix: every component handed this matrix (the corrections of a combination, each on its field) gets the same blocks. That pass also keeps the operator of the Krylov solve (krylovOperator()): \(D A D\) if the solve runs on this scaled matrix (useAsKrylovOperator()), else \(A\), for a solve on the unscaled matrix this one is a view of.

Parameters:: bounds (list)
Return type:: FieldSplitBlocks

krylovOperator()[source]: \(\hat A\) as the operator of AMGCL’s Krylov solvers: \(A\) itself (applied with onUnscaledOperator), or \(D A D\) for a Krylov solve on a scaled matrix, DiagonalScaling(LGMRES(...)). Kept from the field split of the matrix if there was one (see useAsKrylovOperator()), else built in a threaded pass.

matrix()[source]

\(\hat A\) itself, for a component built on the whole matrix (AMG, a direct solver): \(A\), or \(D A D\) formed in the threaded pass of a split into one block (fieldSplit()), which also keeps \(A\) as the Krylov operator.

Return type:: csr_matrix

releaseFormedMatrices()[source]

Drop what was formed from this matrix in a solve, the native field split and the operator of a Krylov solve (formed again when next asked for), and keep only \(A\): for a matrix a component keeps between solves (the \(\tilde S\) of a Schur complement).

Return type:: None

scaled(D)[source]

\(D \hat A D\), without forming it (see fieldSplit()).

Parameters:: D (ndarray)
Return type:: SystemMatrix

useAsKrylovOperator()[source]

Say that a Krylov solve runs on this matrix: a field split formed afterwards keeps \(\hat A\) itself as the operator of that solve, in its one threaded pass (see krylovOperator()).

Return type:: None

property nnz: int: The stored entries of \(A\) (a proxy for its sparsity pattern).

edelweissfe.linsolve.composition.systemmatrix.canonicalMatrix(A)[source]

A as CSR with sorted columns and no duplicate entries: AMGCL’s block backends read the columns as sorted, and no backend sums duplicate entries. A canonical copy if needed, leaving the caller’s matrix as it is.

Return type:: csr_matrix

The one adapter between a nonlinear solver and a composed solver: ComposedLinearSolver, and the solver script of *linearSolver, type=custom (loadSolverScript()).

The nonlinear solver calls its linear solver as x = solver(A, b) once per Newton iteration. ComposedLinearSolver turns that call into the life cycle of a composed solver: it finds what changed since the last call (a new pattern, a resized system, a jump of the right-hand side that marks a new increment), hands the changed matrix down the tree (update()), and applies the root to \(b\). It is also what reports on the solves: the journal, the warnings, the dumps (solvedumps), and the hot reload of a solver script. Every type of linear solver is a tree behind this adapter (edelweissfe.linsolve.types); a direct solver is a tree of one leaf.

In an input file, a solver script is the custom type:

*linearSolver, name=mySolver, type=custom
script=solver.py

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

solver.py builds the tree with plain constructors and binds it to the global linearSolver; it holds only the math, the adapter’s settings (verbosity, dumps, hot reload) are data lines of the block:

linearSolver = LGMRES(DiagonalScaling(...))

class edelweissfe.linsolve.composition.solver.ComposedLinearSolver(root, residualGrowthFactor=4.0, verbosity='warning', warningOuterIterations=100, dumps=None, nodeCoordinates=None, midsideCorners=None)[source]

A composed solver (the root of a tree of preconditioners, usually a Krylov solver) as a LinearSolver, callable as x = solver(A, b): the driver of the tree.

Per call it finds what changed and hands it down the tree as a MatrixChange: "first" (nothing built yet), "model" (after setModel()), "reload" (a new tree from a hot reload), "amr" (another field layout or size), "pattern" (another number of stored entries, an O(1) proxy for the sparsity pattern), or only new values. Whether a right-hand side jumped by more than residualGrowthFactor (a new increment) is passed along too; the Krylov solver decides what it means.

Parameters:

root (Preconditioner) – The composed solver.
residualGrowthFactor (float) – A right-hand side this many times larger than the previous one marks a new increment (a cutback after a diverged attempt starts below its last residual, and is not one).
verbosity (str) – "silent", "warning" (default: nothing on a normal solve; a warning for a slow or unconverged one), "info" (the structure of the solver at the first solve, in a frame, and the “linear solve” column of the nonlinear solver’s table) or "debug" (also each true-residual restart of a Krylov solve, and why the preconditioner was rebuilt or reused, in that column).
warningOuterIterations (int) – A solve with more outer iterations than this gets a warning (a sign that the preconditioner no longer fits). A solve that misses its tolerance gets a warning of its own, whatever this says.
dumps (SolveDumps) – A SolveDumps: write the systems of chosen solves to disk, for an offline replay. None (default): none.
nodeCoordinates (dict) – {fieldName: (nNodes, dimension) array}, coordinates to use instead of the model’s (an offline replay).
midsideCorners (dict) – {fieldName: (isCorner, edgeEndpoints)}, midside corners to use instead of those derived from the model.

Attributes:

lastSolveSummary
needsPrunedSystemMatrix: Whether the tree factorizes (a block of) the system matrix directly (containsDirectSolver), whose reordering is sensitive to explicitly stored zeros, and no Krylov solver of the tree keeps its preconditioner across solves (reusesPreconditioner), which a pattern pruned anew in every iteration would defeat.
reportsSolveSummary: The “linear solve” column of the nonlinear solver: the outer iterations of a Krylov solve, at verbosity "info" or "debug".

Methods

`__call__`(A, b)	Solve \(A x = b\), for one right-hand side `b` (a vector) or several (the columns of a matrix, e.g. the two of an arc-length method), each with the same set-up of the tree.
`coordinatesOf`(fieldName)	The node coordinates of a field in its dof order: given ones first, else the model's, else `None`.
`fromScript`(script[, hotReload])	The composed solver of a solver script (see `loadSolverScript()`) behind the adapter.
`midsideCornersOf`(fieldName)	The midside corners of a field (which nodes are corners, and the corners of each midside node): a given one first, else derived from the model, else `None`.
`setFieldStructure`(fields)	Receive the ordered field blocks of the DOF vector directly (in DOF order), bypassing `setModel()`.
`setJournal`(journal)	Receive the shared `Journal` instance.
`setModel`(model, dofManager)	Take the model: the field layout (base class), and for the tree what is derived from it.

classmethod fromScript(script, hotReload=False, **settings)[source]

The composed solver of a solver script (see loadSolverScript()) behind the adapter.

Parameters:

script (str) – The solver script.
hotReload (bool) – Build the tree again whenever the script changed on disk, before a solve, and continue with it (see _hotReload()).
**settings – The settings of the adapter, see ComposedLinearSolver.

Return type:

ComposedLinearSolver

coordinatesOf(fieldName)[source]

The node coordinates of a field in its dof order: given ones first, else the model’s, else None.

Parameters:: fieldName (str)
Return type:: ndarray | None

midsideCornersOf(fieldName)[source]

The midside corners of a field (which nodes are corners, and the corners of each midside node): a given one first, else derived from the model, else None. An element whose topology the corner-node interpolation does not know is an error (UnsupportedElementTopology).

Parameters:: fieldName (str)

setModel(model, dofManager)[source]

Take the model: the field layout (base class), and for the tree what is derived from it. What the tree built for the previous model is forgotten right away: the tree is built anew at the next solve, and the nonlinear solver builds the new equation system before that.

Return type:: None

hotReload: Whether the tree is built again whenever the solver script changes on disk (see _hotReload()).

identification = 'ComposedSolver': Display name for this solver’s own log lines; overridden by solvers with their own name.

lastNnz

a new count is a new sparsity pattern.

Type:: The stored entries of the last matrix

lastResidualNorm: The right-hand side norm of the last solve, None before the first one and after a new model.

property lastSolveSummary: LinearSolveSummary | None: The type of the None singleton.

property needsPrunedSystemMatrix: bool: Whether the tree factorizes (a block of) the system matrix directly (containsDirectSolver), whose reordering is sensitive to explicitly stored zeros, and no Krylov solver of the tree keeps its preconditioner across solves (reusesPreconditioner), which a pattern pruned anew in every iteration would defeat. A tree of Krylov solvers and AMG, the default, does not need it.

property reportsSolveSummary: bool

the outer iterations of a Krylov solve, at verbosity "info" or "debug".

Type:: The “linear solve” column of the nonlinear solver

script: The solver script the tree was built by (see fromScript()), None for a tree built in Python.

edelweissfe.linsolve.composition.solver.loadSolverScript(path)[source]

The composed solver a solver script defines: the tree of components it binds to its global linearSolver.

The script holds only the math, e.g. linearSolver = LGMRES(DiagonalScaling(...)); the settings of the adapter (verbosity, dumps, hot reload) are data lines of its *linearSolver, type=custom block.

Parameters:: path (str) – The script.
Return type:: Preconditioner

edelweissfe.linsolve.composition.solver.tileFieldBlocks(fieldStructure, n, owner)[source]

The field blocks tiling [0, n), in dof order, with any trailing dofs not covered by a node field (e.g. scalar variables) folded into a final scalar block.

Parameters:

fieldStructure (list | None) – The field blocks of the model, as setModel() derived them.
n (int) – The number of dofs.
owner (str) – The name error messages start with.

Return type:

list

Writing the systems of chosen solves to disk, for an offline replay: SolveDumps.

Comparing solver variants by rerunning the simulation is slow and confounded: every run repeats the assembly, and a variant that changes the Newton iterates changes the sequence of matrices being compared. Dumping authentic systems \((A, b)\) once and replaying them offline – through a ComposedLinearSolver with setFieldStructure and the dumped nodeCoordinates and midsideCorners – makes the comparison fast and controlled.

Which solves are dumped is the trigger:

"every": every solve (a sequence: whether a symbolic factorization or a preconditioner can be reused across Newton iterations is a question about a sequence, not about one matrix);
solve numbers, e.g. {12, 13, 14}: these solves of each solver (counted from 1);
"slow": a solve that needs more outer iterations than a threshold, or stops at the iteration cap short of its tolerance. A preconditioner that degrades – a reused hierarchy that no longer fits the Jacobian, a damaged operator AMG cannot handle – typically shows only tens of hours into an analysis, at a solve nobody can name in advance; this trigger decides by the solve’s own outcome.

One snapshot cannot tell an intrinsically hard operator from a solve made slow by its own state (a stale preconditioner): a fresh solver replaying the matrix alone may converge in a fraction of the live iterations. So a dump can also hold the solves that preceded a dumped one (contextSolves), for a replay through one persistent solver in the same order; and every dumped solve records why it was (not) rebuilt.

maxDumps caps the systems written to one directory in a process: a coupled 3D system of a million dofs takes gigabytes on disk. The solvers of all analysis steps that dump to the same directory share that cap.

On disk, in directory: A_<instance>_<solve>.npz (SciPy, uncompressed), b_<instance>_<solve>.npy, schurDofs_<instance>_<solve>.npy (the condensed dofs of a Schur complement, if any), geometry_<instance>_<solve>.npz (the node coordinates of the vector fields, <field>_coordinates, and the midside corners the tree used, <field>_isCorner and <field>_edgeEndpoints, if any) and one line per dumped solve in manifest.jsonl: the field blocks and file names, and the bookkeeping of a Krylov solve in the names of its SolveRecord (iterations, residual, eta, capHit, trueResidualRestarts, rebuildReason, scaledSystem), with patternChanged, newIncrement and previousIterations. <instance> is the ordinal of the solver in the process (one per analysis step).

class edelweissfe.linsolve.composition.solvedumps.SolveDumps(directory, trigger='slow', threshold=None, maxDumps=10, contextSolves=0)[source]

The dumps of a composed solver, see the module documentation; the data lines dumps, dumpDirectory, maxDumps and dumpContextSolves of every type of linear solver.

Parameters:

directory (str) – Where the dumps go; created if missing.
trigger – Which solves are dumped: "slow" (default), "every", or a collection of solve numbers.
threshold (int) – For "slow": a solve with more outer iterations than this is dumped; None: the warning threshold of the solver (warningOuterIterations), so that a dump and its warning fire on the same solves. A solve that stops at the iteration cap short of its tolerance is dumped whatever the threshold.
maxDumps (int) – The most systems written to directory in this process, across all solvers that dump there (e.g. one per analysis step), context solves included: a disk-space guard for a run that degrades on every remaining solve.
contextSolves (int) – How many of the solves preceding a dumped one are dumped with it (0: only the slow solve itself). They are kept in memory, as copies, only while this is positive.

Attributes:

contextSolves
dumpsWritten: The systems written to directory in this process, by every solver that dumps there.

Methods

`afterSolve`(instance, solveCount, A, b, ...)	Take a finished solve: dump it, with its context window, if the trigger selects it and the budget allows, then keep it (copied) as context for the next ones.
`keepContextSolves`(contextSolves)	Keep this many preceding solves from now on (the newest of those kept so far stay).
`moveTo`(directory)	Write the dumps to `directory` from now on (created if missing; `None`: write none).
`selects`(solveCount, slow)	Whether the trigger selects this solve, given whether it was slow.
`thresholdOr`(warningThreshold)	The outer iterations above which a solve is slow: `threshold`, or `warningThreshold`.

afterSolve(instance, solveCount, A, b, blocks, state, condensedDofs, geometry, slow, log)[source]

Take a finished solve: dump it, with its context window, if the trigger selects it and the budget allows, then keep it (copied) as context for the next ones.

Parameters:

instance (int) – The ordinal of the solver in this process and of the solve: the file names.
solveCount (int) – The ordinal of the solver in this process and of the solve: the file names.
A (csr_matrix) – The raw system as handed to the solver.
b (ndarray) – The raw system as handed to the solver.
blocks (list) – Its FieldBlock layout.
state (dict) – The solve’s bookkeeping, recorded in the manifest (see the module documentation).
condensedDofs (ndarray | None) – The dofs a Schur complement condensed in this solve, or None.
geometry (dict) – {name: array}: the node coordinates and midside corners of the fields, written to geometry_*.npz (nothing if empty).
slow (bool) – Whether the solve was slow (see thresholdOr()), for the trigger "slow".
log – (level, message) -> None.

Return type:

None

keepContextSolves(contextSolves)[source]

Keep this many preceding solves from now on (the newest of those kept so far stay).

Parameters:: contextSolves (int)
Return type:: None

moveTo(directory)[source]

Write the dumps to directory from now on (created if missing; None: write none).

Parameters:: directory (str | None)
Return type:: None

selects(solveCount, slow)[source]

Whether the trigger selects this solve, given whether it was slow.

Parameters:

solveCount (int)
slow (bool)

Return type:

bool

thresholdOr(warningThreshold)[source]

The outer iterations above which a solve is slow: threshold, or warningThreshold.

Parameters:: warningThreshold (int)
Return type:: int

dumpedSolves: The solves already written, so that overlapping context windows never write (or count) one twice.

property dumpsWritten: int: The systems written to directory in this process, by every solver that dumps there.

dumpsWrittenTo: dict = {}

what maxDumps bounds.

Type:: The systems written to each directory (by absolute path) in this process, by every solver

history

{"solveCount", "A", "b", "blocks", "state", "condensedDofs", "geometry"}.

Type:: The last contextSolves solves, oldest first

What the two-level cycles of a vector field (TwoLevel) share: the corner-node two-level cycle (cornernodetwolevelcycle) and the rigid-body coarse level (rigidbodymodestwolevelcycle).

Both are one two-level cycle – fine Chebyshev sweeps on the field’s scaled block, and the coarse correction on the Galerkin coarse operator – that differ in their coarse level and coarse solve: the rigid-body coarse level solves with an AMGCL hierarchy of its own, the corner-node two-level cycle with a component of the composed solver. Both are built and applied in C++: the shared cycle is TwoLevelCycle in amgcl-wrapper.hpp, wrapped as _PyAMGCLTwoLevelCycle in amgcl.pyx. TwoLevelCycle mirrors those bases here: it creates the C++ object, which validates every parameter, and passes the cycle, the fine spectral radius a build estimates, the timings and the report through. A derived class chooses the parameters and builds its coarse level, in its build, from what that coarse level needs.

class edelweissfe.linsolve.composition.twolevelcycle.TwoLevelCycle(cppClass, parameters)[source]

A two-level preconditioner of one vector field’s equilibrated diagonal block, built once per hierarchy refresh and applied many times within the outer Krylov solve, with the build / applyPreconditioner life cycle of PyAMGCLSolver.

Parameters:

cppClass – The Cython class of the C++ preconditioner, created here from parameters.
parameters (dict) – Its parameters, which its constructor checks – the one place they are checked.

Raises:

ValueError – For parameters the C++ class cannot use, before anything is built.

Attributes:

fineSpectralRadius: The estimate of the fine spectral radius the smoother was built with: relax.spectralRadiusSteps Arnoldi steps from a start vector seeded per row, zero on the Dirichlet rows (amgcl-spectral-radius.hpp).

Methods

`applyPreconditioner`(r)	One two-level cycle: fine Chebyshev sweeps and the coarse correction on the free dofs; the Dirichlet rows are solved exactly, `x[i] = r[i] / A[i, i]`.
`cycleParameters`(amgclTree, backendBlockSize, ...)	The parameters of the cycle, in the one spelling both C++ classes read (`TwoLevelCycleParameters` in `amgcl-wrapper.hpp`).
`report`()	The coarse level in figures, and the report of its coarse solve (diagnostics, not a convergence gate).
`timings`()	Wall time [s] of the last build's stages and the cumulative apply (fine vs. coarse) -- kept measured, not assumed, since the split decides where the next optimization goes.

static cycleParameters(amgclTree, backendBlockSize, backendPrecision)[source]

The parameters of the cycle, in the one spelling both C++ classes read (TwoLevelCycleParameters in amgcl-wrapper.hpp).

Parameters:

amgclTree (dict) – The field’s AMGCL tree. Its relax (a Chebyshev polynomial, without Jacobi scaling), npre and npost are the fine smoother, with AMGCL’s defaults for what it leaves out; the default trees, the one place of EdelweissFE’s defaults, are defaultFieldTree() and, for the corner-node two-level cycle, defaultCornerNodeTwoLevelTree().
backendBlockSize (int) – The storage of the fine operator: 1 (scalar entries), or 2/3 (node-wise blocks of the field’s components: less memory traffic per product, the same cycle up to the order of summation).
backendPrecision (str) – The value type of the stored fine operator and of \(P\), \(P^T\): "double", or "float" (not for 2x2 blocks) for half the bytes per entry; vectors and sums stay double, and the spectral radius is estimated on the double operator.

applyPreconditioner(r)[source]

One two-level cycle: fine Chebyshev sweeps and the coarse correction on the free dofs; the Dirichlet rows are solved exactly, x[i] = r[i] / A[i, i].

Parameters:: r (ndarray)
Return type:: ndarray

report()[source]

The coarse level in figures, and the report of its coarse solve (diagnostics, not a convergence gate).

Return type:: str

timings()[source]

Wall time [s] of the last build’s stages and the cumulative apply (fine vs. coarse) – kept measured, not assumed, since the split decides where the next optimization goes.

Return type:: dict

property fineSpectralRadius: float

relax.spectralRadiusSteps Arnoldi steps from a start vector seeded per row, zero on the Dirichlet rows (amgcl-spectral-radius.hpp).

Type:: The estimate of the fine spectral radius the smoother was built with

The rigid-body coarse level of a vector field: the cycle of TwoLevel(coarseSpace=RigidBodyModesPerAggregate()) (TwoLevel).

Why

The near null-space of linear elasticity – the displacements that cost (almost) no strain energy – are the rigid-body modes: three translations and three rotations in 3D, two translations and one rotation in 2D. Smoothed aggregation needs them in its coarse space. An error the Chebyshev smoother leaves alone – its eigenvalues lie below the smoother’s window, which starts at a fraction of the spectral radius by construction – is removed by the coarse-grid correction only as far as the coarse space can represent it on every aggregate.

The block backends of AMG, the default for 2D and 3D vector fields, take the translations alone (AMGCL’s hierarchy on node-wise blocks cannot take more): their coarse space holds displacements that are constant on every aggregate. In a bending-dominated structure – a slender beam, a tunnel lining, a plate – the smooth error is bending, which rotates every aggregate almost as a rigid body. The translations represent it poorly, the smoother does not touch it, and the outer Krylov solve pays with iterations.

What

The rigid-body coarse level keeps the fine level of the block backend – the field’s equilibrated block \(\hat A = D A D\) in node-wise blocks, smoothed by the Chebyshev polynomial of the field’s tree – and builds the first coarse level itself, from all rigid-body modes:

aggregates: AMGCL’s plain aggregation of the nodes, as AMGCL’s own first level on the block backend;
tentative prolongation: the field’s rigid-body modes about its centroid, divided by the equilibration (the near null-space of \(\hat A\) is \(D^{-1}\) times that of \(A\)) – built in Python by rigidBodyNullspace(), the one definition of the modes, and handed to C++ – restricted to every aggregate, zero on the Dirichlet rows (which the cycle solves exactly), and made orthonormal by a QR factorization, as AMGCL’s tentative prolongation does with a near null-space. A mode that depends linearly on the ones before it – a rotation of a single-node aggregate, the rotation about the line through a collinear aggregate, a mode whose components are all constrained – is dropped;
smoothed prolongation: \(P = (I - \omega\, \mathrm{diag}_B(\hat A_F)^{-1} \hat A_F) P_{tent}\), AMGCL’s, with \(\mathrm{diag}_B\) the node-wise \(B \times B\) diagonal blocks of the filtered operator;
Galerkin coarse operator: \(A_c = P^T \hat A P\), 6 coarse dofs per aggregate in 3D (3 in 2D), on which an AMGCL hierarchy continues the field’s tree, with the R factors of the QR factorizations – the coefficients of the rigid-body modes in the coarse basis – as its near null-space.

One application is the two-level cycle the corner-node interpolation (cornernodetwolevelcycle) uses too: restriction of the residual, one V-cycle of the AMGCL hierarchy, prolongation, and the fine Chebyshev sweeps. Everything runs in C++: the steps in amgcl-rigid-body-coarse-level.hpp, the cycle in amgcl-wrapper.hpp (RigidBodyModesTwoLevelCycleT).

When

It is the default of the blockAMG() preset for the 3D displacement field (of any element type), and pays most where the solves would need many outer iterations although the operator is (nearly) elastic: typically a bending-dominated structure. The rigid-body modes are the near null-space of elasticity on the displacement field; a rotation field, or any other vector field, has another, so a solver script gives it the coarse level only by name (OnField("displacement", TwoLevel(coarseSpace=RigidBodyModesPerAggregate())), e.g. for a 2D displacement field). It needs the field’s node coordinates (read from the model), and the field’s tree configures its smoother, its aggregation and the AMGCL hierarchy below it. Its setup costs more than AMGCL’s hierarchy on the translations (the Galerkin product with 6 coarse dofs per aggregate, the larger hierarchy below), so on a model where that converges in a few iterations it gains little. The measurements are in The rigid-body coarse level.

class edelweissfe.linsolve.composition.rigidbodymodestwolevelcycle.RigidBodyModesTwoLevelCycle(amgclTree, backendBlockSize, backendPrecision='double', coarseHierarchyPrecision='double', prolongation='smoothed')[source]

The preconditioner of one vector field with a rigid-body coarse level (see the module docstring): a two-level preconditioner (TwoLevelCycle) whose coarse level is spanned by the rigid-body modes of the aggregates. The C++ class takes any near null-space; this class hands it the rigid-body modes (nearNullspace()).

Parameters:

amgclTree (dict) – The parameters of the cycle, see cycleParameters(); backendBlockSize is the field’s components per node, 3 (6 rigid-body modes per aggregate) or 2 (3). The AMGCL tree means what it means to AMGCL’s own hierarchy, with AMGCL’s defaults for what it leaves out: relax, npre and npost the fine smoother; coarsening (smoothed aggregation) the aggregation and the prolongation smoothing of the rigid-body coarse level; and the whole tree the AMGCL hierarchy below it, with the strong-connection threshold halved, as AMGCL halves it from one level to the next.
backendBlockSize (int) – The parameters of the cycle, see cycleParameters(); backendBlockSize is the field’s components per node, 3 (6 rigid-body modes per aggregate) or 2 (3). The AMGCL tree means what it means to AMGCL’s own hierarchy, with AMGCL’s defaults for what it leaves out: relax, npre and npost the fine smoother; coarsening (smoothed aggregation) the aggregation and the prolongation smoothing of the rigid-body coarse level; and the whole tree the AMGCL hierarchy below it, with the strong-connection threshold halved, as AMGCL halves it from one level to the next.
backendPrecision (str) – The parameters of the cycle, see cycleParameters(); backendBlockSize is the field’s components per node, 3 (6 rigid-body modes per aggregate) or 2 (3). The AMGCL tree means what it means to AMGCL’s own hierarchy, with AMGCL’s defaults for what it leaves out: relax, npre and npost the fine smoother; coarsening (smoothed aggregation) the aggregation and the prolongation smoothing of the rigid-body coarse level; and the whole tree the AMGCL hierarchy below it, with the strong-connection threshold halved, as AMGCL halves it from one level to the next.
coarseHierarchyPrecision (str) – The value type of the AMGCL hierarchy below the rigid-body coarse level (its whole setup): "double" or "float".
prolongation (str) – "smoothed" (default) or "tentative": the tentative prolongation alone, a cheaper setup and a sparser coarse operator at the price of more outer iterations.

Raises:

ValueError – For settings that cannot be used, before anything is built.

Attributes:

aggregates: The aggregate of every node, < 0 for a node in none (e.g.
droppedModes: The number of modes dropped as linearly dependent.
fineSpectralRadius: The estimate of the fine spectral radius the smoother was built with: relax.spectralRadiusSteps Arnoldi steps from a start vector seeded per row, zero on the Dirichlet rows (amgcl-spectral-radius.hpp).
nAggregates: The number of aggregates.
nCoarse: The coarse dofs: one per mode of the near null-space on every aggregate.
prolongationSmoothingWeight: The damping weight omega of the prolongation smoothing (0 for the tentative prolongation).

Methods

`applyPreconditioner`(r)	One two-level cycle: fine Chebyshev sweeps and the coarse correction on the free dofs; the Dirichlet rows are solved exactly, `x[i] = r[i] / A[i, i]`.
`build`(A, scaling, coordinates)	Build the rigid-body coarse level, the AMGCL hierarchy below it and the fine smoother, all in C++, on the near null-space of `nearNullspace()`.
`coarseNullspace`()	The near null-space of the hierarchy below: the coefficients of the near null-space -- the field's rigid-body modes, about its centroid -- in the coarse basis, per aggregate the R factor of its QR factorization (for inspection).
`coarseOperator`()	The Galerkin coarse operator `P^T \hat A P` as the hierarchy below holds it, a scipy CSR matrix (for inspection); a dropped mode has a zero row and column but for its diagonal entry.
`cycleParameters`(amgclTree, backendBlockSize, ...)	The parameters of the cycle, in the one spelling both C++ classes read (`TwoLevelCycleParameters` in `amgcl-wrapper.hpp`).
`nearNullspace`(block, coordinates, scaling)	The near null-space whose vectors span the coarse space of every aggregate: the field's rigid-body modes, divided by the equilibration -- `rigidBodyNullspace()`, the one definition of the modes (6 in 3D, 3 in 2D).
`prolongationMatrix`()	The prolongation `P` as the cycle uses it, a scipy CSR matrix (for inspection).
`report`()	The coarse level in figures, and the report of its coarse solve (diagnostics, not a convergence gate).
`timings`()	Wall time [s] of the last build's stages and the cumulative apply (fine vs. coarse) -- kept measured, not assumed, since the split decides where the next optimization goes.

static nearNullspace(block, coordinates, scaling)[source]

The near null-space whose vectors span the coarse space of every aggregate: the field’s rigid-body modes, divided by the equilibration – rigidBodyNullspace(), the one definition of the modes (6 in 3D, 3 in 2D). The C++ class takes any near null-space of B or B (B + 1) / 2 columns; a subclass may hand it another (the tests’ control takes the translations, AMGCL’s own coarse space).

Parameters:

block (FieldBlock)
coordinates (ndarray)
scaling (ndarray)

Return type:

ndarray

build(A, scaling, coordinates)[source]

Build the rigid-body coarse level, the AMGCL hierarchy below it and the fine smoother, all in C++, on the near null-space of nearNullspace().

Parameters:

A – The field’s equilibrated diagonal block \(\hat A = D A D\) (scipy CSR), node-major, with its Dirichlet identity rows.
scaling – The field’s slice of the equilibration \(D\): the modes are divided by it.
coordinates – The field’s node coordinates in its dof order, (nNodes, >= blockSize) (columns beyond the field’s dimension, e.g. z of a 2D field, are ignored): the modes are built from them.

Return type:

None

coarseNullspace()[source]

The near null-space of the hierarchy below: the coefficients of the near null-space – the field’s rigid-body modes, about its centroid – in the coarse basis, per aggregate the R factor of its QR factorization (for inspection).

Return type:: ndarray

coarseOperator()[source]: The Galerkin coarse operator P^T \hat A P as the hierarchy below holds it, a scipy CSR matrix (for inspection); a dropped mode has a zero row and column but for its diagonal entry.

prolongationMatrix()[source]: The prolongation P as the cycle uses it, a scipy CSR matrix (for inspection).

property aggregates: ndarray: The aggregate of every node, < 0 for a node in none (e.g. a fully constrained one).

property droppedModes: int: The number of modes dropped as linearly dependent.

property nAggregates: int: The number of aggregates.

property nCoarse: int

one per mode of the near null-space on every aggregate.

Type:: The coarse dofs

property prolongationSmoothingWeight: float: The damping weight omega of the prolongation smoothing (0 for the tentative prolongation).

p-multigrid (Galerkin P1 corner-node) preconditioner for one field’s diagonal block: the cycle of TwoLevel(coarseSpace=CornerNodeInterpolation(coarseSolve)) (TwoLevel).

Precondition the quadratic operator through its linear (P1) part: \(\nu\) Chebyshev sweeps on the field’s own (scaled) block, restrict the residual through \(P^T\), solve on the Galerkin-projected \(A_c = P^T A P\) with the coarse component, prolong through \(P\), \(\nu\) more Chebyshev sweeps.

The prolongation interpolates the scaled variables

The block is \(\hat A = D A D\) (a DiagonalScaling above), and a displacement \(u\) of the physical problem is \(\hat u = D^{-1} u\) in its variables. The P1 interpolant acts on \(u\): \(P_{geo}\) keeps a corner node’s displacement and gives a midside node the mean of its two edge-endpoint corners. In the scaled variables the same interpolation is

\[P = D^{-1} P_{geo} D_c: \quad P_{cc} = 1, \qquad P_{ma} = \tfrac{1}{2} \frac{s_a}{s_m}, \quad P_{mb} = \tfrac{1}{2} \frac{s_b}{s_m}\]

for a midside dof \(m\) on the edge \((a, b)\), with \(s\) the diagonal of \(D\) and \(D_c\) its restriction to the corners. The plain weights \(\tfrac{1}{2}, \tfrac{1}{2}\) would interpolate \(\hat u\) as if it were \(u\): a scaled translation \(D^{-1} t\), the smoothest error there is, would be off by the ratios \(s_a / s_m\) at every midside node, and the coarse level could not remove it. With \(P\) the translations are reproduced exactly, and

\[A_c = P^T \hat A P = D_c \, P_{geo}^T A P_{geo} \, D_c\]

is the P1 operator of the physical problem in the corners’ scaled variables: the operator of the corner-node field, scaled by \(D_c\). Without a scaling above (\(D = I\)), \(P = P_{geo}\).

The coarse problem is a field of its own

The coarse solve is a component of the composed solver, set up on \(A_c\) with the context of the corner-node field (ofCoarseNodes(): the same field on the coarse nodes, their coordinates, the scaling \(D_c\)), so that every component works on it unchanged. Its default, DiagonalScaling(AMG()), works on any field; for a 3D displacement field DiagonalScaling(TwoLevel(coarseSpace=RigidBodyModesPerAggregate())) is faster (The corner-node two-level cycle with scaled weights and a composed coarse solve), and the blockamg preset passes it, with the smoother of the corner-node two-level cycle (defaultCornerNodeTwoLevelTree()). Both scale \(A_c\) to a unit diagonal once more: where a midside diagonal is tiny against its corners’, the weights \(s_a / s_m\) are large, \(A_c\) inherits entries of their square, and a Chebyshev smoother on the unscaled \(A_c\) would work with a spectral radius far above that of a unit-diagonal operator.

The cycle runs in C++ (CornerNodeTwoLevelCycleT in edelweissfe/linsolve/amgcl/amgcl-wrapper.hpp, wrapped as PyAMGCLCornerNodeTwoLevelCycle), OpenMP-threaded throughout, and calls the coarse component back once per application.

Dirichlet rows. A field’s diagonal block still carries its Dirichlet identity rows (the Dirichlet elimination runs before the linear solver), and the identity rows of multi-point-constraint slaves, which are treated alike. The fine smoothing and the Galerkin projection act on the free DOFs only, with the Dirichlet DOFs decoupled entirely: the cycle is that of the free submatrix, and a Dirichlet DOF that is only masked after a step still couples into the free DOFs within it. A free midside node whose edge-endpoint corner is Dirichlet-constrained keeps only the surviving weight of its free endpoint – no renormalization (the constrained corner contributes exactly zero to a homogeneous Newton correction, which the dropped weight already encodes). A free midside dof whose both endpoint dofs are constrained (a symmetry plane and a tie on one edge, issue #223) would get an empty row of \(P\) – no coarse correction at all – so its node becomes a coarse node of its own (CornerNodeTwoLevelCycle.coarseNodes), as buildMidsideCorners() makes a node a corner on an AMR edge conflict; it is decided at the build, where the Dirichlet rows are known, not in the topology map. A Dirichlet coarse dof is a zero column of \(P\) and an identity row of \(A_c\): the coarse field keeps its Dirichlet rows, as the fine one does.

The C++ implementation does not extract the free submatrix, it decouples the Dirichlet DOFs exactly inside the full block layout instead: the fine right-hand side and iterate are kept zero on them (a Dirichlet row holds only its diagonal, so every residual and every Chebyshev update stays zero there, whatever the columns hold), and P has no rows for them. That is the free-submatrix cycle up to floating-point summation order, and it is what allows a node-wise 3x3 block storage of the fine operator (fineBlockSize) on partially constrained nodes.

class edelweissfe.linsolve.composition.cornernodetwolevelcycle.CornerNodeTwoLevelCycle(isCorner, edgeEndpoints, amgclTree, backendBlockSize=1, backendPrecision='double')[source]

A corner-node two-level cycle for one vector field’s diagonal block, a two-level preconditioner (TwoLevelCycle) whose coarse level is the P1 corner nodes, and whose coarse solve is a callable, set with setCoarseSolve().

Parameters:

isCorner (ndarray) – This field’s midside corners (edelweissfe.numerics.p1topology.buildMidsideCorners()), in the field’s own node order.
edgeEndpoints (ndarray) – This field’s midside corners (edelweissfe.numerics.p1topology.buildMidsideCorners()), in the field’s own node order.
amgclTree (dict) – The parameters of the cycle, see cycleParameters(): the fine Chebyshev smoother is the relax, npre and npost of the field’s AMGCL tree. With npre 0 the cycle starts with the coarse-grid correction of the right-hand side itself and ends with npost sweeps: one smoother application fewer per cycle than a symmetric one, at the price of a nonsymmetric preconditioner – of no concern to the right-preconditioned GMRES around it, which needs a fixed linear operator, not a symmetric one. Ending on the smoother is the right way round for that one sweep: it damps the high-frequency error the prolongated correction leaves on the fine level.
backendBlockSize (int) – The parameters of the cycle, see cycleParameters(): the fine Chebyshev smoother is the relax, npre and npost of the field’s AMGCL tree. With npre 0 the cycle starts with the coarse-grid correction of the right-hand side itself and ends with npost sweeps: one smoother application fewer per cycle than a symmetric one, at the price of a nonsymmetric preconditioner – of no concern to the right-preconditioned GMRES around it, which needs a fixed linear operator, not a symmetric one. Ending on the smoother is the right way round for that one sweep: it damps the high-frequency error the prolongated correction leaves on the fine level.
backendPrecision (str) – The parameters of the cycle, see cycleParameters(): the fine Chebyshev smoother is the relax, npre and npost of the field’s AMGCL tree. With npre 0 the cycle starts with the coarse-grid correction of the right-hand side itself and ends with npost sweeps: one smoother application fewer per cycle than a symmetric one, at the price of a nonsymmetric preconditioner – of no concern to the right-preconditioned GMRES around it, which needs a fixed linear operator, not a symmetric one. Ending on the smoother is the right way round for that one sweep: it damps the high-frequency error the prolongated correction leaves on the fine level.

Attributes:

coarseNodes: Per node, whether the last build made it a coarse node: the corner nodes, and every midside node with a free dof whose two edge-endpoint dofs are both constrained (see the module documentation).
fineSpectralRadius: The estimate of the fine spectral radius the smoother was built with: relax.spectralRadiusSteps Arnoldi steps from a start vector seeded per row, zero on the Dirichlet rows (amgcl-spectral-radius.hpp).
nCoarse: The coarse dofs: all dofs of the coarseNodes.

Methods

`applyPreconditioner`(r)	One two-level cycle: fine Chebyshev sweeps and the coarse correction on the free dofs; the Dirichlet rows are solved exactly, `x[i] = r[i] / A[i, i]`.
`build`(A, scaling)	Build the prolongation, the Galerkin coarse operator and the fine Chebyshev smoother, all in C++, and return the coarse operator, on which the caller sets up the coarse solve before the first application.
`cycleParameters`(amgclTree, backendBlockSize, ...)	The parameters of the cycle, in the one spelling both C++ classes read (`TwoLevelCycleParameters` in `amgcl-wrapper.hpp`).
`prolongationMatrix`()	The prolongation `P` as the cycle uses it, a scipy CSR matrix (for inspection).
`report`()	The coarse level in figures, and the report of its coarse solve (diagnostics, not a convergence gate).
`setCoarseSolve`(coarseSolve)	The coarse solve, a callable `r_c -> z_c` with \(z_c \approx A_c^{-1} r_c\), e.g. a component of the composed solver once it is set up on the coarse operator `build()` returns; called once per application.
`timings`()	Wall time [s] of the last build's stages and the cumulative apply (fine vs. coarse) -- kept measured, not assumed, since the split decides where the next optimization goes.

build(A, scaling)[source]

Build the prolongation, the Galerkin coarse operator and the fine Chebyshev smoother, all in C++, and return the coarse operator, on which the caller sets up the coarse solve before the first application.

Parameters:

A – This field’s scaled diagonal block \(\hat A = D A D\) (scipy CSR), still carrying its Dirichlet identity rows.
scaling (ndarray) – This field’s slice of the scaling \(D\): the weights of the prolongation (see the module documentation).

Returns:

\(A_c = P^T \hat A P\) on all dofs of the coarseNodes, in node order and then by component, with an identity row for each Dirichlet coarse dof; 0 x 0 for a field without a free dof (the cycle is then the exact diagonal inverse and needs no coarse solve).

Return type:

sp.csr_matrix

prolongationMatrix()[source]

The prolongation P as the cycle uses it, a scipy CSR matrix (for inspection).

Return type:: csr_matrix

setCoarseSolve(coarseSolve)[source]

The coarse solve, a callable r_c -> z_c with \(z_c \approx A_c^{-1} r_c\), e.g. a component of the composed solver once it is set up on the coarse operator build() returns; called once per application. Set before or after the build, before the first application; the cycle calls the one set last.

Return type:: None

property coarseNodes: ndarray

the corner nodes, and every midside node with a free dof whose two edge-endpoint dofs are both constrained (see the module documentation).

Type:: Per node, whether the last build made it a coarse node

property nCoarse: int

all dofs of the coarseNodes.

Type:: The coarse dofs