This page was generated from unit-2.11-matrixfree/matrixfree.ipynb.
2.11 Matrix free operator application¶
Usually, we assemble matrices in sparse matrix format. Certain methods allow to improve performance and reduce memory requirements considerably by avoiding building and storing the system matrix. The counterpart of this approach is that it is now difficult to build preconditioners
Hybrid mixed method¶
consider the hybrid mixed method
\[\begin{split}\begin{array}{ccccl}
A \sigma & + & B^T u & = & f \\
B u & & & = & g
\end{array}\end{split}\]
where the matrix \(A\) comes from
\[\int \sigma \tau dx\]
and the matrix \(B\) from
\[\sum_T \int_T \operatorname{div} \, \sigma \, v_T - \int_{\partial T} \sigma_n v_F\]
For \(\sigma\) we use a VectorL2 space, with a Piola-mapping to the physical elements
[1]:
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import unit_square
from ngsolve.fem import MixedFE
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
order=1
Sigma = VectorL2(mesh, order=order, piola=True)
Vt = L2(mesh, order=order-1)
Vf = FacetFESpace(mesh, order=order)
V = FESpace([Vt,Vf])
sigma,tau = Sigma.TnT()
(ut,uf), (vt,vf) = V.TnT()
Generate Mesh from spline geometry
Boundary mesh done, np = 12
CalcLocalH: 12 Points 0 Elements 0 Surface Elements
Meshing domain 1 / 1
load internal triangle rules
Surface meshing done
Edgeswapping, topological
Smoothing
Split improve
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Split improve
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Split improve
Combine improve
Smoothing
Update mesh topology
Update clusters
[2]:
a = BilinearForm(Sigma)
a += sigma*tau * dx
b = BilinearForm(trialspace=Sigma, testspace=V)
b += div(sigma) * vt * dx
n = specialcf.normal(mesh.dim)
dS = dx(element_boundary=True)
b += -sigma*n*vf * dS
b.Assemble()
assemble mixed bilinearform
[2]:
<ngsolve.comp.BilinearForm at 0x7f70cd7cff30>
compute element matrices …
we should observe that the matrices from the \(A\) operator are block-diagonal with \(2\times 2\) blocks, and each element matrix from \(B\) belongs to one of two equivalence classes.
[3]:
for el in mesh.Elements():
felS = Sigma.GetFE(el)
trafo = mesh.GetTrafo(el)
elmat = a.integrators[0].CalcElementMatrix(felS, trafo)
print (elmat)
0.5 0 0 0.5 0 0
0 0.25 0 0 0.25 0
0 0 0.0833333 0 0 0.0833333
0.5 0 0 1 0 0
0 0.25 0 0 0.5 0
0 0 0.0833333 0 0 0.166667
0.417528 0 0 0.286279 0 0
0 0.208764 0 0 0.143139 0
0 0 0.0695879 0 0 0.0477131
0.286279 0 0 0.795051 0 0
0 0.143139 0 0 0.397525 0
0 0 0.0477131 0 0 0.132508
0.57872 0 0 0.0785934 0 0
0 0.28936 0 0 0.0392967 0
0 0 0.0964533 0 0 0.0130989
0.0785934 0 0 0.442661 0 0
0 0.0392967 0 0 0.221331 0
0 0 0.0130989 0 0 0.0737769
0.783348 0 0 0.425796 0 0
0 0.391674 0 0 0.212898 0
0 0 0.130558 0 0 0.070966
0.425796 0 0 0.550588 0 0
0 0.212898 0 0 0.275294 0
0 0 0.070966 0 0 0.0917646
0.863976 0 0 0.364102 0 0
0 0.431988 0 0 0.182051 0
0 0 0.143996 0 0 0.0606837
0.364102 0 0 0.442802 0 0
0 0.182051 0 0 0.221401 0
0 0 0.0606837 0 0 0.0738004
0.794554 0 0 0.286399 0 0
0 0.397277 0 0 0.1432 0
0 0 0.132426 0 0 0.0477332
0.286399 0 0 0.417875 0 0
0 0.1432 0 0 0.208938 0
0 0 0.0477332 0 0 0.0696458
0.482624 0 0 0.12477 0 0
0 0.241312 0 0 0.062385 0
0 0 0.0804373 0 0 0.020795
0.12477 0 0 0.550258 0 0
0 0.062385 0 0 0.275129 0
0 0 0.020795 0 0 0.0917097
0.5 0 0 0.5 0 0
0 0.25 0 0 0.25 0
0 0 0.0833333 0 0 0.0833333
0.5 0 0 1 0 0
0 0.25 0 0 0.5 0
0 0 0.0833333 0 0 0.166667
0.513945 0 0 0.13755 0 0
0 0.256973 0 0 0.0687751 0
0 0 0.0856575 0 0 0.022925
0.13755 0 0 0.523247 0 0
0 0.0687751 0 0 0.261623 0
0 0 0.022925 0 0 0.0872078
0.664746 0 0 0.164849 0 0
0 0.332373 0 0 0.0824246 0
0 0 0.110791 0 0 0.0274749
0.164849 0 0 0.416964 0 0
0 0.0824246 0 0 0.208482 0
0 0 0.0274749 0 0 0.069494
0.631785 0 0 0.293178 0 0
0 0.315892 0 0 0.146589 0
0 0 0.105297 0 0 0.0488631
0.293178 0 0 0.531753 0 0
0 0.146589 0 0 0.265877 0
0 0 0.0488631 0 0 0.0886255
0.752167 0 0 0.252064 0 0
0 0.376083 0 0 0.126032 0
0 0 0.125361 0 0 0.0420106
0.252064 0 0 0.416844 0 0
0 0.126032 0 0 0.208422 0
0 0 0.0420106 0 0 0.069474
1.0139 0 0 0.721815 0 0
0 0.506948 0 0 0.360908 0
0 0 0.168983 0 0 0.120303
0.721815 0 0 0.76045 0 0
0 0.360908 0 0 0.380225 0
0 0 0.120303 0 0 0.126742
0.523407 0 0 0.137311 0 0
0 0.261703 0 0 0.0686554 0
0 0 0.0872345 0 0 0.0228851
0.137311 0 0 0.513662 0 0
0 0.0686554 0 0 0.256831 0
0 0 0.0228851 0 0 0.0856104
0.576938 0 0 0.238606 0 0
0 0.288469 0 0 0.119303 0
0 0 0.0961564 0 0 0.0397677
0.238606 0 0 0.532003 0 0
0 0.119303 0 0 0.266002 0
0 0 0.0397677 0 0 0.0886672
0.330872 0 0 0.0390758 0 0
0 0.165436 0 0 0.0195379 0
0 0 0.0551453 0 0 0.00651264
0.0390758 0 0 0.760195 0 0
0 0.0195379 0 0 0.380097 0
0 0 0.00651264 0 0 0.126699
1.01606 0 0 0.553634 0 0
0 0.508029 0 0 0.276817 0
0 0 0.169343 0 0 0.0922723
0.553634 0 0 0.547715 0 0
0 0.276817 0 0 0.273858 0
0 0 0.0922723 0 0 0.0912859
0.570517 0 0 0.197191 0 0
0 0.285258 0 0 0.0985953 0
0 0 0.0950861 0 0 0.0328651
0.197191 0 0 0.506355 0 0
0 0.0985953 0 0 0.253178 0
0 0 0.0328651 0 0 0.0843925
0.700217 0 0 0.281867 0 0
0 0.350109 0 0 0.140933 0
0 0 0.116703 0 0 0.0469778
0.281867 0 0 0.470495 0 0
0 0.140933 0 0 0.235247 0
0 0 0.0469778 0 0 0.0784158
0.5234 0 0 0.406298 0 0
0 0.2617 0 0 0.203149 0
0 0 0.0872334 0 0 0.0677164
0.406298 0 0 0.793042 0 0
0 0.203149 0 0 0.396521 0
0 0 0.0677164 0 0 0.132174
0.470317 0 0 0.281589 0 0
0 0.235158 0 0 0.140794 0
0 0 0.0783861 0 0 0.0469315
0.281589 0 0 0.70015 0 0
0 0.140794 0 0 0.350075 0
0 0 0.0469315 0 0 0.116692
0.45625 0 0 -0.00580578 0 0
0 0.228125 0 0 -0.00290289 0
0 0 0.0760417 0 0 -0.00096763
-0.00580578 0 0 0.548019 0 0
0 -0.00290289 0 0 0.274009 0
0 0 -0.00096763 0 0 0.0913365
0.793155 0 0 0.406059 0 0
0 0.396578 0 0 0.203029 0
0 0 0.132193 0 0 0.0676764
0.406059 0 0 0.52308 0 0
0 0.203029 0 0 0.26154 0
0 0 0.0676764 0 0 0.08718
0.682802 0 0 0.309198 0 0
0 0.341401 0 0 0.154599 0
0 0 0.1138 0 0 0.051533
0.309198 0 0 0.506155 0 0
0 0.154599 0 0 0.253077 0
0 0 0.051533 0 0 0.0843591
[4]:
for el in mesh.Elements():
felS = Sigma.GetFE(el)
felV = V.GetFE(el)
fel = MixedFE(felS, felV)
trafo = mesh.GetTrafo(el)
# try integratos 0 and 1 ...
elmat = b.integrators[0].CalcElementMatrix(fel, trafo)
print (elmat)
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 0 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 0 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 4.26422e-17 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 3.85716e-17 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 -1.95619e-17 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 0 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 -7.42287e-17 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 0 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 -7.42134e-17 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 0 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 2.681e-17 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 0 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 7.09354e-17 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 -1.41884e-16 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 0.5 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 1.5 -0.5 0 2.92391e-17 -1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
Geometry-free matrix multiplication¶
NGSolve provides now the possibility to benefit from the geometry independent element matrix. Just one matrix per equivalence class is stored, and the matrix vector product is performed for all elements simultaneously using dense matrix algebra.
[5]:
from netgen.csg import *
geom = CSGeometry()
geom.Add (Sphere(Pnt(50,50,50),80) \
-Cylinder(Pnt(-100,0,0),Pnt(200,0,0), 40) \
-Cylinder(Pnt(100,-100,100),Pnt(100,200,100),40)
-Cylinder(Pnt(0,100,-100), Pnt(0,100,200),40)
-Sphere(Pnt(50,50,50),50))
# geom.Draw()
mesh = Mesh(geom.GenerateMesh(maxh=25))
mesh.Curve(5)
# Draw (mesh)
Start Findpoints
main-solids: 1
Found points 36
Analyze spec points
Find edges
6 edges found
Check intersecting edges
CalcLocalH: 96 Points 0 Elements 0 Surface Elements
Start Findpoints
Analyze spec points
Find edges
6 edges found
Check intersecting edges
Start Findpoints
Analyze spec points
Find edges
6 edges found
Check intersecting edges
Surface 1 / 5
Surface meshing done
0 illegal triangles
Optimize Surface
Edgeswapping, topological
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
199 elements, 163 points
Surface 2 / 5
Surface meshing done
0 illegal triangles
Optimize Surface
Edgeswapping, topological
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
46 elements, 170 points
Surface 3 / 5
Surface meshing done
0 illegal triangles
Optimize Surface
Edgeswapping, topological
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
44 elements, 176 points
Surface 4 / 5
Surface meshing done
0 illegal triangles
Optimize Surface
Edgeswapping, topological
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
46 elements, 183 points
Surface 5 / 5
Surface meshing done
0 illegal triangles
Optimize Surface
Edgeswapping, topological
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
Edgeswapping, metric
Smoothing
Combine improve
Smoothing
93 elements, 212 points
SplitSeparateFaces
Check subdomain 1 / 1
Meshing subdomain 1 of 1
Use internal rules
Delaunay meshing
number of points: 206
blockfill local h
number of points: 206
Points: 206
Elements: 1142
Tree data entries per element: 1.57618
Tree nodes per element: 0.030648
SwapImprove
SwapImprove
SwapImprove
SwapImprove
0 degenerated elements removed
Remove intersecting
Remove outer
206 points, 498 elements
start tetmeshing
Use internal rules
Success !
209 points, 514 elements
Remove Illegal Elements
Volume Optimization
CombineImprove
ImproveMesh
SplitImprove
ImproveMesh
SwapImprove
SwapImprove2
ImproveMesh
CombineImprove
ImproveMesh
SplitImprove
ImproveMesh
SwapImprove
SwapImprove2
ImproveMesh
CombineImprove
ImproveMesh
SplitImprove
ImproveMesh
SwapImprove
SwapImprove2
ImproveMesh
Update mesh topology
Update clusters
Curve elements, order = 5
Update clusters
Curving edges
Curving faces
Same as above, with the geom_free flag for the BilinearForm:
[6]:
order=5
SetHeapSize(100*1000*1000)
Sigma = VectorL2(mesh, order=order, piola=True)
Vt = L2(mesh, order=order-1)
Vf = FacetFESpace(mesh, order=order, dirichlet=[2])
V = FESpace([Vt,Vf])
print ("Sigma.ndof =", Sigma.ndof, ", V.ndof =", V.ndof)
sigma,tau = Sigma.TnT()
(ut,uf), (vt,vf) = V.TnT()
b = BilinearForm(trialspace=Sigma, testspace=V, geom_free=True)
b += div(sigma) * vt * dx
n = specialcf.normal(mesh.dim)
dS = dx(element_boundary=True)
b += -sigma*n*vf * dS
b.Assemble()
Sigma.ndof = 82320 , V.ndof = 42098
[6]:
<ngsolve.comp.BilinearForm at 0x7f70cd7edaf0>
Check timing with geom_free=True/False:
[7]:
x = b.mat.CreateRowVector()
y = b.mat.CreateColVector()
x[:] = 1
%timeit y.data = b.mat * x
1.97 ms ± 153 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
L2 finite element spaces have a Mass method generating a linear operator for multiplication with the diagonal mass matrix. In case of variable coefficients or curved elements it uses numerical integration. It also provide an inverse operator:
[8]:
ainv = Sigma.Mass(rho=1).Inverse()
Combine the operators to form the Schur complement, which is actually a discretization of the Laplace operator:
[9]:
Laplace = b.mat @ ainv @ b.mat.T
[10]:
f = LinearForm (V)
f += 1*vt * dx
f.Assemble()
proj = Projector(V.FreeDofs(), True)
gfu = GridFunction (V)
from time import time
start = time()
with TaskManager():
solvers.CG(mat=Laplace, pre=proj, rhs=f.vec, sol=gfu.vec, maxsteps=5000, printrates=False)
print ("time =", time()-start)
assemble VOL element 490/490
time = 42.366188526153564
The iteration count is very high since we used only the most trivial preconditioner. We can add a low order coarse-grid preconditioner to get rid of the mesh dependency, and also local inverses to improve the \(p\)-dependency.
[11]:
Draw (gfu.components[1], mesh, "gfu1", draw_vol=False)
[11]:
BaseWebGuiScene
[12]:
help (solvers.CG)
Help on function CG in module ngsolve.krylovspace:
CG(mat, rhs, pre=None, sol=None, tol=1e-12, maxsteps=100, printrates=True, initialize=True, conjugate=False, callback=None)
preconditioned conjugate gradient method
Parameters
----------
mat : Matrix
The left hand side of the equation to solve. The matrix has to be spd o hermitsch.
rhs : Vector
The right hand side of the equation.
pre : Preconditioner
If provided the preconditioner is used.
sol : Vector
Start vector for CG method, if initialize is set False. Gets overwritten by the solution vector. If sol = None then a new vector is created.
tol : double
Tolerance of the residuum. CG stops if tolerance is reached.
maxsteps : int
Number of maximal steps for CG. If the maximal number is reached before the tolerance is reached CG stops.
printrates : bool
If set to True then the error of the iterations is displayed.
initialize : bool
If set to True then the initial guess for the CG method is set to zero. Otherwise the values of the vector sol, if provided, is used.
conjugate : bool
If set to True, then the complex inner product is used.
Returns
-------
(vector)
Solution vector of the CG method.