5. Iterative Solvers#
So far we have used direct solvers to solve the linear system of equations. Although a direct solver can profit from the sparse matrix, it’s arithmetic complexity is sub-optimal. For large-scale problems iterative solvers are a must.
The conjugate gradient (cg) method is the standard method for symmetric and positive definite matrices. It’s convergence rate depends on a preconditioner, what is a cheap approximative inverse to the matrix.
from ngsolve import *
from ngsolve.webgui import Draw
We generate a 3D geometry and mesh using the OCC constructive solid geometry (CSG) modeler:
from netgen.occ import *
cube = Box((0,0,0),(1,1,1))
cyl = Cylinder((0,0.5,0.5),X, r=0.2, h=1)
cube.faces.name = "outer"
cyl.faces.name = "cyl"
shape = cube-cyl
ea = { "euler_angles" : [-180, -80, 165] }
Draw(shape, **ea);
Generate a mesh. To obtain finer meshes, we perform hierarchical mesh refinement:
ngmesh = OCCGeometry(shape).GenerateMesh(maxh=0.1)
mesh = Mesh(ngmesh).Curve(3)
for l in range(0):
mesh.Refine()
Draw (mesh, **ea);
We try different preconditioners:
Localis a Jacobi or block-Jacobi preconditionerBDDCis an element-level domain decomposition preconditionerMultiGridbenefits from hierarchically refined meshes
fes = H1(mesh, order=3, dirichlet="outer", wb_withedges=False)
print ("we have", fes.ndof, "unknowns")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)
f = LinearForm(1*v*dx)
pre = preconditioners.Local(a)
# pre = preconditioners.BDDC(a)
# pre = preconditioners.MultiGrid(a)
gfu = GridFunction(fes)
we have 23783 unknowns
assemble system and setup preconditioner in parallel:
with TaskManager():
a.Assemble()
f.Assemble()
solve the system using the preconditioned conjugate gradient method:
from ngsolve.krylovspace import CGSolver
from time import time
ts = time()
with TaskManager():
# inv = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())
inv = solvers.CGSolver(mat=a.mat, pre=pre, printrates=True, plotrates=False, maxiter=1000)
gfu.vec.data = inv * f.vec
te = time()
print ("needed", te-ts, "seconds for", fes.ndof, "unknowns")
CG iteration 1, residual = 0.04767218279616599
CG iteration 2, residual = 0.05351997384858849
CG iteration 3, residual = 0.05547110488845035
CG iteration 4, residual = 0.045650611795690375
CG iteration 5, residual = 0.030119809357029216
CG iteration 6, residual = 0.02356790155193162
CG iteration 7, residual = 0.01845379391732243
CG iteration 8, residual = 0.012791616931561138
CG iteration 9, residual = 0.009140044746927672
CG iteration 10, residual = 0.007236156169219033
CG iteration 11, residual = 0.006530384543816025
CG iteration 12, residual = 0.006311000334322975
CG iteration 13, residual = 0.006533383309504386
CG iteration 14, residual = 0.0068526086855284175
CG iteration 15, residual = 0.0066677062776169865
CG iteration 16, residual = 0.005552600994308546
CG iteration 17, residual = 0.004204894733746069
CG iteration 18, residual = 0.0031156825675337954
CG iteration 19, residual = 0.0022847737919975105
CG iteration 20, residual = 0.001749365949534994
CG iteration 21, residual = 0.001414406100657776
CG iteration 22, residual = 0.0012972124397142145
CG iteration 23, residual = 0.00123831105554679
CG iteration 24, residual = 0.0011543087077594983
CG iteration 25, residual = 0.0010376512900760835
CG iteration 26, residual = 0.000872242295741223
CG iteration 27, residual = 0.0006795873617928119
CG iteration 28, residual = 0.0005162389432181474
CG iteration 29, residual = 0.00038886105714954503
CG iteration 30, residual = 0.0003104401031959345
CG iteration 31, residual = 0.00025377040218641715
CG iteration 32, residual = 0.0002092337700717176
CG iteration 33, residual = 0.00017233958284846744
CG iteration 34, residual = 0.0001417861697254287
CG iteration 35, residual = 0.00011569049448405617
CG iteration 36, residual = 9.657739592630568e-05
CG iteration 37, residual = 8.14842118899809e-05
CG iteration 38, residual = 7.104601974209938e-05
CG iteration 39, residual = 6.505764680544796e-05
CG iteration 40, residual = 5.8623468548756195e-05
CG iteration 41, residual = 5.013631378348798e-05
CG iteration 42, residual = 4.217247926828296e-05
CG iteration 43, residual = 3.2767730254554715e-05
CG iteration 44, residual = 2.454408316296603e-05
CG iteration 45, residual = 1.8726582302721968e-05
CG iteration 46, residual = 1.4504232832724984e-05
CG iteration 47, residual = 1.1545469940745142e-05
CG iteration 48, residual = 9.81609237186518e-06
CG iteration 49, residual = 8.708010640820944e-06
CG iteration 50, residual = 8.064298158390331e-06
CG iteration 51, residual = 7.570574444695999e-06
CG iteration 52, residual = 6.837364667576531e-06
CG iteration 53, residual = 5.888398167249054e-06
CG iteration 54, residual = 4.6861081701326035e-06
CG iteration 55, residual = 3.6323507742212955e-06
CG iteration 56, residual = 2.714256496535724e-06
CG iteration 57, residual = 2.0223776620583106e-06
CG iteration 58, residual = 1.5450752025279479e-06
CG iteration 59, residual = 1.2451134362748749e-06
CG iteration 60, residual = 1.0689023371297396e-06
CG iteration 61, residual = 9.710550764242155e-07
CG iteration 62, residual = 9.392457092918174e-07
CG iteration 63, residual = 8.869746309841282e-07
CG iteration 64, residual = 7.839113224051142e-07
CG iteration 65, residual = 6.45533221930883e-07
CG iteration 66, residual = 4.997255073903688e-07
CG iteration 67, residual = 3.728389565951669e-07
CG iteration 68, residual = 2.86798270080072e-07
CG iteration 69, residual = 2.2526257937416442e-07
CG iteration 70, residual = 1.8728018395169988e-07
CG iteration 71, residual = 1.6227255103386358e-07
CG iteration 72, residual = 1.457642220913908e-07
CG iteration 73, residual = 1.342630412564839e-07
CG iteration 74, residual = 1.2162816907268396e-07
CG iteration 75, residual = 1.0760756091535105e-07
CG iteration 76, residual = 9.085270910122731e-08
CG iteration 77, residual = 7.441606056269484e-08
CG iteration 78, residual = 5.9553351438221224e-08
CG iteration 79, residual = 4.686203730964407e-08
CG iteration 80, residual = 3.6118845715638345e-08
CG iteration 81, residual = 2.8572320709168047e-08
CG iteration 82, residual = 2.3045631406389168e-08
CG iteration 83, residual = 1.9036014569497077e-08
CG iteration 84, residual = 1.67390328068067e-08
CG iteration 85, residual = 1.520884233456327e-08
CG iteration 86, residual = 1.4080442409717375e-08
CG iteration 87, residual = 1.2533820428721961e-08
CG iteration 88, residual = 1.0644100907490584e-08
CG iteration 89, residual = 8.570783911421492e-09
CG iteration 90, residual = 6.639060920100317e-09
CG iteration 91, residual = 5.107952458386973e-09
CG iteration 92, residual = 4.04678114285079e-09
CG iteration 93, residual = 3.30841188456239e-09
CG iteration 94, residual = 2.8100150801170045e-09
CG iteration 95, residual = 2.418845340780697e-09
CG iteration 96, residual = 2.1020435078689607e-09
CG iteration 97, residual = 1.8444093248317163e-09
CG iteration 98, residual = 1.5844601303401638e-09
CG iteration 99, residual = 1.3779531233990642e-09
CG iteration 100, residual = 1.1772487793122262e-09
CG iteration 101, residual = 1.015302473598604e-09
CG iteration 102, residual = 8.725529128097508e-10
CG iteration 103, residual = 7.215157426398298e-10
CG iteration 104, residual = 5.814792041083531e-10
CG iteration 105, residual = 4.491435364329812e-10
CG iteration 106, residual = 3.482963960948355e-10
CG iteration 107, residual = 2.6941582204377447e-10
CG iteration 108, residual = 2.212627701529477e-10
CG iteration 109, residual = 1.9727418626931502e-10
CG iteration 110, residual = 1.8212750878300682e-10
CG iteration 111, residual = 1.7413342403516025e-10
CG iteration 112, residual = 1.6103147567798684e-10
CG iteration 113, residual = 1.3824308256987428e-10
CG iteration 114, residual = 1.125493417224692e-10
CG iteration 115, residual = 8.696314290247476e-11
CG iteration 116, residual = 6.407650473933682e-11
CG iteration 117, residual = 4.706683929469554e-11
CG iteration 118, residual = 3.6378454066742264e-11
CG iteration 119, residual = 2.988792796823192e-11
CG iteration 120, residual = 2.643518068477906e-11
CG iteration 121, residual = 2.5300436014099583e-11
CG iteration 122, residual = 2.4548671569449216e-11
CG iteration 123, residual = 2.2804125480883428e-11
CG iteration 124, residual = 1.9592728796802015e-11
CG iteration 125, residual = 1.549603108680883e-11
CG iteration 126, residual = 1.1530124956562722e-11
CG iteration 127, residual = 8.5426668391718e-12
CG iteration 128, residual = 6.330323985015678e-12
CG iteration 129, residual = 4.871813974213172e-12
CG iteration 130, residual = 4.167529826169715e-12
CG iteration 131, residual = 3.828247606062055e-12
CG iteration 132, residual = 3.683268746055522e-12
CG iteration 133, residual = 3.571692524139383e-12
CG iteration 134, residual = 3.3219380589574273e-12
CG iteration 135, residual = 2.838391782780697e-12
CG iteration 136, residual = 2.25596502262534e-12
CG iteration 137, residual = 1.6715543135350264e-12
CG iteration 138, residual = 1.1941006083397006e-12
CG iteration 139, residual = 8.822400192073492e-13
CG iteration 140, residual = 6.755777497907131e-13
CG iteration 141, residual = 5.676559860631217e-13
CG iteration 142, residual = 5.211380616415192e-13
CG iteration 143, residual = 5.050400284804928e-13
CG iteration 144, residual = 4.973407775803987e-13
CG iteration 145, residual = 4.4969325294095336e-13
CG iteration 146, residual = 3.7488546604506905e-13
CG iteration 147, residual = 2.892248893259312e-13
CG iteration 148, residual = 2.0950247755088933e-13
CG iteration 149, residual = 1.5117406969982347e-13
CG iteration 150, residual = 1.1207775072203375e-13
CG iteration 151, residual = 8.673316200942853e-14
CG iteration 152, residual = 7.529349990129952e-14
CG iteration 153, residual = 6.998795512877412e-14
CG iteration 154, residual = 6.793922880125646e-14
CG iteration 155, residual = 6.542830331502609e-14
CG iteration 156, residual = 5.933642866346914e-14
CG iteration 157, residual = 4.973172290720412e-14
CG iteration 158, residual = 3.781707630415108e-14
needed 0.020255088806152344 seconds for 23783 unknowns
Draw (gfu, **ea);
help (solvers.CGSolver)
Help on class CGSolver in module ngsolve.krylovspace:
class CGSolver(LinearSolver)
| CGSolver(
| *args,
| conjugate: bool = False,
| abstol: float = None,
| maxsteps: int = None,
| printing: bool = False,
| **kwargs
| )
|
| Method resolution order:
| CGSolver
| LinearSolver
| ngsolve.la.BaseMatrix
| pybind11_builtins.pybind11_object
| builtins.object
|
| Methods defined here:
|
| __init__(
| self,
| *args,
| conjugate: bool = False,
| abstol: float = None,
| maxsteps: int = None,
| printing: bool = False,
| **kwargs
| )
| __init__(*args, **kwargs)
| Overloaded function.
|
| 1. __init__(self: ngsolve.la.BaseMatrix) -> None
|
| 2. __init__(self: ngsolve.la.BaseMatrix, arg0: ngsolve.la.BaseVector) -> None
|
| 3. __init__(self: ngsolve.la.BaseMatrix, arg0: ngsolve.la.MultiVector) -> None
|
| 4. __init__(self: ngsolve.la.BaseMatrix, arg0: ngsolve.bla.MatrixD) -> None
|
| 5. __init__(self: ngsolve.la.BaseMatrix, arg0: object) -> None
|
| ----------------------------------------------------------------------
| Readonly properties defined here:
|
| errors
|
| ----------------------------------------------------------------------
| Data and other attributes defined here:
|
| name = 'CG'
|
| ----------------------------------------------------------------------
| Methods inherited from LinearSolver:
|
| CheckResidual(self, residual)
|
| CreateVector(self, col)
| CreateVector(self: ngsolve.la.BaseMatrix, colvector: bool = False) -> ngsolve.la.BaseVector
|
| Height(self) -> int
|
| IsComplex(self) -> bool
|
| Mult(self, x: BaseVector, y: BaseVector) -> None
| Mult(self: ngsolve.la.BaseMatrix, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| Solve = retfunc(*args, **kwargs) from netgen.TimeFunction.<locals>
|
| Update(self)
| Update(self: ngsolve.la.BaseMatrix) -> None
|
| Update matrix
|
| Width(self) -> int
|
| ----------------------------------------------------------------------
| Class methods inherited from LinearSolver:
|
| Creator(**kwargs)
|
| ----------------------------------------------------------------------
| Data descriptors inherited from LinearSolver:
|
| __dict__
| dictionary for instance variables
|
| ----------------------------------------------------------------------
| Methods inherited from ngsolve.la.BaseMatrix:
|
| AsVector(...)
| AsVector(self: ngsolve.la.BaseMatrix) -> ngsolve.la.BaseVector
|
| Interprets the matrix values as a vector
|
| CreateColVector(...)
| CreateColVector(self: ngsolve.la.BaseMatrix) -> ngsolve.la.BaseVector
|
| CreateDeviceMatrix(...)
| CreateDeviceMatrix(self: ngsolve.la.BaseMatrix) -> BaseMatrix
|
| CreateMatrix(...)
| CreateMatrix(self: ngsolve.la.BaseMatrix) -> BaseMatrix
|
| Create matrix of same dimension and same sparsestructure
|
| CreateRowVector(...)
| CreateRowVector(self: ngsolve.la.BaseMatrix) -> ngsolve.la.BaseVector
|
| CreateSparseMatrix(...)
| CreateSparseMatrix(self: ngsolve.la.BaseMatrix) -> ngla::BaseSparseMatrix
|
| GetInverseType(...)
| GetInverseType(self: ngsolve.la.BaseMatrix) -> str
|
| GetOperatorInfo(...)
| GetOperatorInfo(self: ngsolve.la.BaseMatrix) -> str
|
| Inverse(...)
| Inverse(self: ngsolve.la.BaseMatrix, freedofs: pyngcore.pyngcore.BitArray = None, inverse: None | str | object = None, flags: pyngcore.pyngcore.Flags = <pyngcore.pyngcore.Flags object at 0x1107f8cf0>) -> BaseMatrix
|
| Calculate inverse of sparse matrix
| Parameters:
|
| freedofs : BitArray
| If set, invert only the rows/columns the matrix defined by the bit array, otherwise invert the whole matrix
|
| inverse : string
| Solver to use, allowed values are:
| sparsecholesky - internal solver of NGSolve for symmetric matrices
| umfpack - solver by Suitesparse/UMFPACK (if NGSolve was configured with USE_UMFPACK=ON)
| pardiso - PARDISO, either provided by libpardiso (USE_PARDISO=ON) or Intel MKL (USE_MKL=ON).
| If neither Pardiso nor Intel MKL was linked at compile-time, NGSolve will look
| for libmkl_rt in LD_LIBRARY_PATH (Unix) or PATH (Windows) at run-time.
|
| MultAdd(...)
| MultAdd(*args, **kwargs)
| Overloaded function.
|
| 1. MultAdd(self: ngsolve.la.BaseMatrix, value: typing.SupportsFloat | typing.SupportsIndex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| 2. MultAdd(self: ngsolve.la.BaseMatrix, value: typing.SupportsComplex | typing.SupportsFloat | typing.SupportsIndex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| MultScale(...)
| MultScale(*args, **kwargs)
| Overloaded function.
|
| 1. MultScale(self: ngsolve.la.BaseMatrix, value: typing.SupportsFloat | typing.SupportsIndex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| 2. MultScale(self: ngsolve.la.BaseMatrix, value: typing.SupportsComplex | typing.SupportsFloat | typing.SupportsIndex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| MultTrans(...)
| MultTrans(*args, **kwargs)
| Overloaded function.
|
| 1. MultTrans(self: ngsolve.la.BaseMatrix, value: typing.SupportsFloat | typing.SupportsIndex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| 2. MultTrans(self: ngsolve.la.BaseMatrix, value: typing.SupportsComplex | typing.SupportsFloat | typing.SupportsIndex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| MultTransAdd(...)
| MultTransAdd(*args, **kwargs)
| Overloaded function.
|
| 1. MultTransAdd(self: ngsolve.la.BaseMatrix, value: typing.SupportsFloat | typing.SupportsIndex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| 2. MultTransAdd(self: ngsolve.la.BaseMatrix, value: typing.SupportsComplex | typing.SupportsFloat | typing.SupportsIndex, x: ngsolve.la.BaseVector, y: ngsolve.la.BaseVector) -> None
|
| SetInverseType(...)
| SetInverseType(self: ngsolve.la.BaseMatrix, arg0: str) -> str
|
| ToDense(...)
| ToDense(self: ngsolve.la.BaseMatrix) -> object
|
| __add__(...)
| __add__(self: BaseMatrix, mat: BaseMatrix) -> BaseMatrix
|
| __iadd__(...)
| __iadd__(self: ngsolve.la.BaseMatrix, mat: ngsolve.la.BaseMatrix) -> None
|
| __matmul__(...)
| __matmul__(self: BaseMatrix, mat: BaseMatrix) -> BaseMatrix
|
| __mul__(...)
| __mul__(*args, **kwargs)
| Overloaded function.
|
| 1. __mul__(self: BaseMatrix, arg0: ngsolve.la.BaseVector) -> ngla::DynamicVectorExpression
|
| 2. __mul__(self: BaseMatrix, arg0: ngsolve.la.MultiVector) -> ngsolve.la.MultiVectorExpr
|
| __neg__(...)
| __neg__(self: BaseMatrix) -> BaseMatrix
|
| __radd__(...)
| __radd__(self: BaseMatrix, arg0: typing.SupportsInt | typing.SupportsIndex) -> BaseMatrix
|
| __rmul__(...)
| __rmul__(*args, **kwargs)
| Overloaded function.
|
| 1. __rmul__(self: BaseMatrix, value: typing.SupportsFloat | typing.SupportsIndex) -> BaseMatrix
|
| 2. __rmul__(self: BaseMatrix, value: typing.SupportsComplex | typing.SupportsFloat | typing.SupportsIndex) -> BaseMatrix
|
| __str__(...)
| __str__(self: ngsolve.la.BaseMatrix) -> str
|
| __sub__(...)
| __sub__(self: BaseMatrix, mat: BaseMatrix) -> BaseMatrix
|
| __timing__(...)
| __timing__(self: ngsolve.la.BaseMatrix, runs: typing.SupportsInt | typing.SupportsIndex = 10) -> float
|
| matvec(...)
| matvec(self: BaseMatrix, arg0: ngsolve.la.BaseVector) -> ngsolve.la.BaseVector
|
| rmatvec(...)
| rmatvec(self: BaseMatrix, arg0: ngsolve.la.BaseVector) -> ngsolve.la.BaseVector
|
| ----------------------------------------------------------------------
| Static methods inherited from ngsolve.la.BaseMatrix:
|
| GetDefaultInverseType(...) method of pybind11_builtins.pybind11_detail_function_record_v1_system_libcpp_abi1 instance
| GetDefaultInverseType() -> str
|
| SetDefaultInverseType(...) method of pybind11_builtins.pybind11_detail_function_record_v1_system_libcpp_abi1 instance
| SetDefaultInverseType(type: str) -> None
|
| ----------------------------------------------------------------------
| Readonly properties inherited from ngsolve.la.BaseMatrix:
|
| H
| Return conjugate transpose of matrix (WIP, only partially supported)
|
| T
| Return transpose of matrix
|
| comm
|
| dtype
|
| height
| Height of the matrix
|
| is_complex
| is the matrix complex-valued ?
|
| is_symmetric
|
| local_mat
|
| nze
| number of non-zero elements
|
| shape
|
| width
| Width of the matrix
|
| ----------------------------------------------------------------------
| Static methods inherited from pybind11_builtins.pybind11_object:
|
| __new__(*args, **kwargs) class method of pybind11_builtins.pybind11_object
| Create and return a new object. See help(type) for accurate signature.