This page was generated from unit-2.1.2-blockjacobi/blockjacobi.ipynb.

2.1.2 Building blocks for programming preconditioners

In this tutorial we will discuss

• Jacobi and Gauss-Seidel smoothers/preconditioners,

• their block versions,

• how to select your own smoothing blocks, and

• how to additively combine preconditioners.

import netgen.gui
from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

def Setup(h=0.1, p=3):
    mesh = Mesh(unit_square.GenerateMesh(maxh=h))
    fes = H1(mesh, order=p, dirichlet="left|bottom")
    u, v = fes.TnT()
    a = BilinearForm(fes)
    a += grad(u)*grad(v)*dx + u*v*dx
    f = LinearForm(fes)
    f += v*dx
    gfu = GridFunction(fes)
    return mesh, fes, a, f, gfu

mesh, fes, a, f, gfu = Setup(h=0.1, p=3)
a.Assemble()
f.Assemble()

<ngsolve.comp.LinearForm at 0x7f659826b0b0>

Create a Jacobi-preconditioner

Let \(A=\) a.mat be the assembled matrix, which can be decomposed based on FreeDofs (\(F\)) the remainder (\(D\)), as in \(\S\)1.3,

\[\begin{split}A = \left( \begin{array}{cc} A_{FF} & A_{FD} \\ A_{DF} & A_{DD} \end{array} \right).\end{split}\]

Then the matrix form of the point Jacobi preconditioner is

\[\begin{split}J = \left( \begin{array}{cc} \text{diag}(A_{FF})^{-1} & 0 \\ 0 & 0 \end{array} \right),\end{split}\]

which can be obtained in NGSolve using CreateSmoother:

preJpoint = a.mat.CreateSmoother(fes.FreeDofs())

NGSolve also gives us a facility to quickly check an estimate of the condition number of the preconditioned matrix by applying the Lanczos algorithm on the preconditioned system.

lams = EigenValues_Preconditioner(mat=a.mat, pre=preJpoint)
lams

 0.0144743
 0.0565951
 0.0940742
 0.111584
 0.148079
 0.204136
 0.279083
 0.379019
 0.45319
 0.53846
 0.662266
 0.764709
 0.854406
 0.966403
 1.07186
 1.18706
 1.28597
 1.40124
 1.50818
 1.62481
 1.74747
 1.88402
  1.9919
 2.07749
 2.18043
 2.29565
 2.38514
 2.45388
 2.50384
 2.56324
 2.58875
  2.6348
 2.81085

An estimate of the condition number \(\kappa = \lambda_{\text{max}} / \lambda_{\text{min}}\) is therefore given as follows:

max(lams)/min(lams)

194.19561167462462

One might wonder if we have gained anything by this point Jacobi preconditioning. What if we did not precondition at all?

Not preconditioning is the same as preconditioning by the identity operator on \(F\)-dofs. One way to realize this identity operator in NGSolve is through the projection into the space of free dofs (i.e., the space spanned by the shape functions corresponding to free dofs). NGSolve provides

Projector(mask, range)   # mask: bit array; range: bool

which projects into the space spanned by the shape functions of the degrees of freedom marked as range in mask.

preI = Projector(mask=fes.FreeDofs(), range=True)

Note that another way to obtain the identity matrix in NGSolve is

IdentityMatrix(fes.ndof, complex=False).

lams = EigenValues_Preconditioner(mat=a.mat, pre=preI)
max(lams)/min(lams)

1678.5881226918377

Clearly the point Jacobi preconditioner has reduced the condition number.

We can use preconditioners within iterative solvers provided by NGSolve’s solvers module (which has MinRes, QMR etc.) Here is an illustration of its use within CG, or the conjugate gradient solver:

solvers.CG(mat=a.mat, pre=preJpoint, rhs=f.vec, sol=gfu.vec)
Draw(gfu)

iteration 0 error = 0.055055687118289384
iteration 1 error = 0.09357794070949911
iteration 2 error = 0.1125535931192053
iteration 3 error = 0.08835788661576421
iteration 4 error = 0.09465752852240356
iteration 5 error = 0.08675791614352736
iteration 6 error = 0.09174564400414806
iteration 7 error = 0.07950174974020316
iteration 8 error = 0.055581104814418365
iteration 9 error = 0.044799001641920536
iteration 10 error = 0.038850616493040765
iteration 11 error = 0.03040867549925388
iteration 12 error = 0.022887415077875117
iteration 13 error = 0.020510924867973467
iteration 14 error = 0.015405453021112544
iteration 15 error = 0.011893431227989994
iteration 16 error = 0.010292148243983711
iteration 17 error = 0.00788847986761108
iteration 18 error = 0.005481703824765562
iteration 19 error = 0.0036120151680207454
iteration 20 error = 0.0023840872403187567
iteration 21 error = 0.0015057839927216498
iteration 22 error = 0.001161353306064317
iteration 23 error = 0.0008694202455218208
iteration 24 error = 0.000532986073869304
iteration 25 error = 0.00037108669875231867
iteration 26 error = 0.0002859501917016284
iteration 27 error = 0.00022783224021485782
iteration 28 error = 0.00016659432373390605
iteration 29 error = 0.00011944094829848387
iteration 30 error = 8.147974841738305e-05
iteration 31 error = 5.502065307009907e-05
iteration 32 error = 4.190089073735409e-05
iteration 33 error = 3.369706727003579e-05
iteration 34 error = 2.7766665129423247e-05
iteration 35 error = 1.973603865717574e-05
iteration 36 error = 1.363994533705489e-05
iteration 37 error = 1.0506189808377953e-05
iteration 38 error = 8.983446326953014e-06
iteration 39 error = 7.66052120833184e-06
iteration 40 error = 5.801953313066026e-06
iteration 41 error = 4.0511973130697e-06
iteration 42 error = 2.9961283895352735e-06
iteration 43 error = 2.376244675254369e-06
iteration 44 error = 1.8104198799133953e-06
iteration 45 error = 1.247306606049404e-06
iteration 46 error = 8.123842440240046e-07
iteration 47 error = 5.602213594412261e-07
iteration 48 error = 4.2973008522608137e-07
iteration 49 error = 3.046923636321491e-07
iteration 50 error = 1.9532592489801808e-07
iteration 51 error = 1.3346842716412627e-07
iteration 52 error = 9.387214842741348e-08
iteration 53 error = 7.042267663035101e-08
iteration 54 error = 5.087660221401061e-08
iteration 55 error = 3.628765551820581e-08
iteration 56 error = 2.4168689891085906e-08
iteration 57 error = 1.6475168255150912e-08
iteration 58 error = 1.0717196210734519e-08
iteration 59 error = 6.891402341379671e-09
iteration 60 error = 4.5529888574774655e-09
iteration 61 error = 3.143607507984258e-09
iteration 62 error = 2.0550328408099652e-09
iteration 63 error = 1.3865747340215791e-09
iteration 64 error = 1.021543719719666e-09
iteration 65 error = 7.443399528840255e-10
iteration 66 error = 5.504334989079559e-10
iteration 67 error = 3.929608302906005e-10
iteration 68 error = 2.63606489195036e-10
iteration 69 error = 1.7938452265190134e-10
iteration 70 error = 1.1723394107561655e-10
iteration 71 error = 7.74247604431213e-11
iteration 72 error = 5.2131329238562905e-11
iteration 73 error = 3.624096808459988e-11
iteration 74 error = 2.334185557902436e-11
iteration 75 error = 1.5421559600265853e-11
iteration 76 error = 1.0252664746513643e-11
iteration 77 error = 6.927119090603895e-12
iteration 78 error = 4.713860694246031e-12
iteration 79 error = 3.154546343263564e-12
iteration 80 error = 1.965748862576144e-12
iteration 81 error = 1.2237482757836687e-12
iteration 82 error = 8.299332439534429e-13
iteration 83 error = 5.756006946531209e-13
iteration 84 error = 4.1076273233521287e-13
iteration 85 error = 2.6538607399531165e-13
iteration 86 error = 1.6936824425901149e-13
iteration 87 error = 1.1027310317902923e-13
iteration 88 error = 7.510563570013987e-14
iteration 89 error = 5.341700604287314e-14

Gauss-Seidel smoothing

The same point Jacobi smoother object can also used to perform point Gauss-Seidel smoothing. One step of the classical Gauss-Seidel iteration is realized by the method preJpoint.Smooth(). Its well known that this iteration converges for matrices like \(A\). Below we show how to use it as a linear iterative solver.

gfu.vec[:] = 0
res = f.vec.CreateVector()              # residual
projres = f.vec.CreateVector()          # residual projected to freedofs
proj = Projector(fes.FreeDofs(), True)

for i in range(500):
    preJpoint.Smooth(gfu.vec, f.vec)    # one step of point Gauss-Seidel
    res.data = f.vec - a.mat*gfu.vec
    projres.data = proj * res
    print ("it#", i, ", res =", Norm(projres))
Draw (gfu)

it# 0 , res = 0.08119734356602216
it# 1 , res = 0.07696557260572953
it# 2 , res = 0.07344471412858082
it# 3 , res = 0.07025007673116736
it# 4 , res = 0.06735911057304737
it# 5 , res = 0.06472122579225086
it# 6 , res = 0.06229422362238182
it# 7 , res = 0.06004424891557192
it# 8 , res = 0.057944353686531225
it# 9 , res = 0.055973109933972955
it# 10 , res = 0.054113396436963314
it# 11 , res = 0.052351431294981175
it# 12 , res = 0.05067602643031652
it# 13 , res = 0.04907801916652114
it# 14 , res = 0.047549839165887446
it# 15 , res = 0.0460851770982483
it# 16 , res = 0.04467872944058138
it# 17 , res = 0.043326000417520126
it# 18 , res = 0.04202314712044809
it# 19 , res = 0.04076685752716332
it# 20 , res = 0.039554253806522846
it# 21 , res = 0.03838281521465994
it# 22 , res = 0.037250316285503716
it# 23 , res = 0.036154777041799285
it# 24 , res = 0.03509442271086202
it# 25 , res = 0.034067650996515805
it# 26 , res = 0.03307300538722952
it# 27 , res = 0.03210915330716984
it# 28 , res = 0.031174868167997208
it# 29 , res = 0.030269014573674532
it# 30 , res = 0.029390536082149733
it# 31 , res = 0.028538445046677165
it# 32 , res = 0.027711814153330697
it# 33 , res = 0.026909769345598
it# 34 , res = 0.02613148388613044
it# 35 , res = 0.02537617335303486
it# 36 , res = 0.024643091406035895
it# 37 , res = 0.023931526188376996
it# 38 , res = 0.023240797254965737
it# 39 , res = 0.02257025293720385
it# 40 , res = 0.021919268071101814
it# 41 , res = 0.021287242028414432
it# 42 , res = 0.020673597001220963
it# 43 , res = 0.020077776499096003
it# 44 , res = 0.019499244025139226
it# 45 , res = 0.018937481902962683
it# 46 , res = 0.018391990231511154
it# 47 , res = 0.017862285948513607
it# 48 , res = 0.017347901986589184
it# 49 , res = 0.016848386508685476
it# 50 , res = 0.016363302211716526
it# 51 , res = 0.01589222568907755
it# 52 , res = 0.015434746844209805
it# 53 , res = 0.014990468348627651
it# 54 , res = 0.014559005138852527
it# 55 , res = 0.014139983947550387
it# 56 , res = 0.013733042864884788
it# 57 , res = 0.013337830926694137
it# 58 , res = 0.012954007726595104
it# 59 , res = 0.012581243049534374
it# 60 , res = 0.012219216524660272
it# 61 , res = 0.01186761729567923
it# 62 , res = 0.011526143707110838
it# 63 , res = 0.011194503005064704
it# 64 , res = 0.010872411051339972
it# 65 , res = 0.010559592049799227
it# 66 , res = 0.010255778284096104
it# 67 , res = 0.00996070986594569
it# 68 , res = 0.00967413449322085
it# 69 , res = 0.009395807217235968
it# 70 , res = 0.009125490218650605
it# 71 , res = 0.00886295259148328
it# 72 , res = 0.008607970134778417
it# 73 , res = 0.008360325151510988
it# 74 , res = 0.008119806254357473
it# 75 , res = 0.007886208177987825
it# 76 , res = 0.007659331597570521
it# 77 , res = 0.007438982953202125
it# 78 , res = 0.007224974279999399
it# 79 , res = 0.007017123043610688
it# 80 , res = 0.00681525198092036
it# 81 , res = 0.006619188945737575
it# 82 , res = 0.006428766759272549
it# 83 , res = 0.006243823065218441
it# 84 , res = 0.006064200189265577
it# 85 , res = 0.005889745002886688
it# 86 , res = 0.005720308791241462
it# 87 , res = 0.005555747125054004
it# 88 , res = 0.005395919736327513
it# 89 , res = 0.005240690397767337
it# 90 , res = 0.0050899268057850505
it# 91 , res = 0.004943500466970448
it# 92 , res = 0.0048012865879140205
it# 93 , res = 0.004663163968275585
it# 94 , res = 0.004529014896994875
it# 95 , res = 0.00439872505154392
it# 96 , res = 0.004272183400128148
it# 97 , res = 0.004149282106743925
it# 98 , res = 0.004029916439004412
it# 99 , res = 0.0039139846786498305
it# 100 , res = 0.003801388034660207
it# 101 , res = 0.0036920305588925405
it# 102 , res = 0.0035858190641653146
it# 103 , res = 0.003482663044718992
it# 104 , res = 0.003382474598979751
it# 105 , res = 0.0032851683545590846
it# 106 , res = 0.003190661395423251
it# 107 , res = 0.0030988731911681896
it# 108 , res = 0.003009725528337382
it# 109 , res = 0.0029231424437238173
it# 110 , res = 0.002839050159597507
it# 111 , res = 0.002757377020801113
it# 112 , res = 0.0026780534336612307
it# 113 , res = 0.002601011806660502
it# 114 , res = 0.0025261864928203637
it# 115 , res = 0.002453513733744944
it# 116 , res = 0.0023829316052758494
it# 117 , res = 0.0023143799647148042
it# 118 , res = 0.002247800399564422
it# 119 , res = 0.0021831361777466013
it# 120 , res = 0.0021203321992539622
it# 121 , res = 0.002059334949193195
it# 122 , res = 0.0020000924521809313
it# 123 , res = 0.0019425542280515418
it# 124 , res = 0.001886671248840185
it# 125 , res = 0.0018323958970043153
it# 126 , res = 0.0017796819248470654
it# 127 , res = 0.0017284844151089314
it# 128 , res = 0.0016787597426922202
it# 129 , res = 0.0016304655374890543
it# 130 , res = 0.0015835606482767016
it# 131 , res = 0.0015380051076530712
it# 132 , res = 0.0014937600979808916
it# 133 , res = 0.001450787918311334
it# 134 , res = 0.0014090519522599629
it# 135 , res = 0.0013685166368059551
it# 136 , res = 0.0013291474319896584
it# 137 , res = 0.0012909107914819628
it# 138 , res = 0.0012537741339998221
it# 139 , res = 0.0012177058155444208
it# 140 , res = 0.0011826751024382766
it# 141 , res = 0.0011486521451371055
it# 142 , res = 0.001115607952796332
it# 143 , res = 0.0010835143685682792
it# 144 , res = 0.001052344045610326
it# 145 , res = 0.0010220704237829882
it# 146 , res = 0.0009926677070189497
it# 147 , res = 0.0009641108413425693
it# 148 , res = 0.000936375493521789
it# 149 , res = 0.0009094380303344124
it# 150 , res = 0.00088327549843047
it# 151 , res = 0.0008578656047741365
it# 152 , res = 0.0008331866976479747
it# 153 , res = 0.000809217748204604
it# 154 , res = 0.0007859383325472232
it# 155 , res = 0.0007633286143278829
it# 156 , res = 0.0007413693278440734
it# 157 , res = 0.00072004176162352
it# 158 , res = 0.0006993277424801499
it# 159 , res = 0.0006792096200283831
it# 160 , res = 0.0006596702516445981
it# 161 , res = 0.000640692987858939
it# 162 , res = 0.0006222616581690231
it# 163 , res = 0.0006043605572613286
it# 164 , res = 0.0005869744316283162
it# 165 , res = 0.0005700884665716905
it# 166 , res = 0.0005536882735783805
it# 167 , res = 0.0005377598780606356
it# 168 , res = 0.0005222897074479154
it# 169 , res = 0.0005072645796222704
it# 170 , res = 0.0004926716916857934
it# 171 , res = 0.00047849860905098944
it# 172 , res = 0.0004647332548458016
it# 173 , res = 0.0004513638996228422
it# 174 , res = 0.00043837915136459997
it# 175 , res = 0.00042576794577676053
it# 176 , res = 0.00041351953686000957
it# 177 , res = 0.0004016234877534652
it# 178 , res = 0.00039006966184135564
it# 179 , res = 0.00037884821411634677
it# 180 , res = 0.0003679495827898038
it# 181 , res = 0.00035736448114428727
it# 182 , res = 0.00034708388962110544
it# 183 , res = 0.0003370990481342521
it# 184 , res = 0.00032740144860569144
it# 185 , res = 0.0003179828277160959
it# 186 , res = 0.00030883515986375145
it# 187 , res = 0.00029995065032549723
it# 188 , res = 0.0002913217286154386
it# 189 , res = 0.0002829410420336342
it# 190 , res = 0.0002748014494010721
it# 191 , res = 0.0002668960149749203
it# 192 , res = 0.0002592180025380226
it# 193 , res = 0.000251760869659717
it# 194 , res = 0.0002445182621203749
it# 195 , res = 0.00023748400849751703
it# 196 , res = 0.0002306521149066934
it# 197 , res = 0.0002240167598942401
it# 198 , res = 0.0002175722894772837
it# 199 , res = 0.00021131321232501228
it# 200 , res = 0.00020523419508074997
it# 201 , res = 0.00019933005781636346
it# 202 , res = 0.00019359576961872258
it# 203 , res = 0.0001880264443035157
it# 204 , res = 0.00018261733625069224
it# 205 , res = 0.00017736383636158545
it# 206 , res = 0.0001722614681311297
it# 207 , res = 0.00016730588383320972
it# 208 , res = 0.00016249286081649724
it# 209 , res = 0.00015781829790598563
it# 210 , res = 0.00015327821190840607
it# 211 , res = 0.0001488687342183249
it# 212 , res = 0.0001445861075214779
it# 213 , res = 0.0001404266825932195
it# 214 , res = 0.00013638691518935554
it# 215 , res = 0.00013246336302550932
it# 216 , res = 0.00012865268284443593
it# 217 , res = 0.00012495162756689338
it# 218 , res = 0.00012135704352518684
it# 219 , res = 0.00011786586777547874
it# 220 , res = 0.00011447512548852358
it# 221 , res = 0.0001111819274138385
it# 222 , res = 0.0001079834674186035
it# 223 , res = 0.00010487702009636327
it# 224 , res = 0.00010185993844440758
it# 225 , res = 9.892965160877539e-05
it# 226 , res = 9.608366269284962e-05
it# 227 , res = 9.331954663062235e-05
it# 228 , res = 9.063494811998718e-05
it# 229 , res = 8.802757961552227e-05
it# 230 , res = 8.549521937938734e-05
it# 231 , res = 8.303570958838159e-05
it# 232 , res = 8.064695449506837e-05
it# 233 , res = 7.83269186422367e-05
it# 234 , res = 7.607362512815241e-05
it# 235 , res = 7.388515392216386e-05
it# 236 , res = 7.17596402289278e-05
it# 237 , res = 6.969527289892685e-05
it# 238 , res = 6.769029288538315e-05
it# 239 , res = 6.574299174562836e-05
it# 240 , res = 6.385171018522758e-05
it# 241 , res = 6.201483664379827e-05
it# 242 , res = 6.023080592190492e-05
it# 243 , res = 5.849809784782394e-05
it# 244 , res = 5.681523598150952e-05
it# 245 , res = 5.518078635709628e-05
it# 246 , res = 5.35933562605996e-05
it# 247 , res = 5.2051593043333585e-05
it# 248 , res = 5.055418296951271e-05
it# 249 , res = 4.909985009665584e-05
it# 250 , res = 4.768735518808334e-05
it# 251 , res = 4.631549465742178e-05
it# 252 , res = 4.498309954279398e-05
it# 253 , res = 4.368903451108626e-05
it# 254 , res = 4.243219688961322e-05
it# 255 , res = 4.1211515727676516e-05
it# 256 , res = 4.002595088327588e-05
it# 257 , res = 3.887449213692731e-05
it# 258 , res = 3.775615833135097e-05
it# 259 , res = 3.666999653430887e-05
it# 260 , res = 3.561508122796114e-05
it# 261 , res = 3.459051351934268e-05
it# 262 , res = 3.359542037478209e-05
it# 263 , res = 3.262895387559145e-05
it# 264 , res = 3.1690290496062226e-05
it# 265 , res = 3.077863040161718e-05
it# 266 , res = 2.9893196766781647e-05
it# 267 , res = 2.903323511412632e-05
it# 268 , res = 2.8198012670515743e-05
it# 269 , res = 2.738681774318998e-05
it# 270 , res = 2.6598959113182042e-05
it# 271 , res = 2.5833765446644516e-05
it# 272 , res = 2.5090584722274026e-05
it# 273 , res = 2.4368783675830648e-05
it# 274 , res = 2.3667747261198902e-05
it# 275 , res = 2.2986878125309097e-05
it# 276 , res = 2.232559609994721e-05
it# 277 , res = 2.1683337706855653e-05
it# 278 , res = 2.1059555677894848e-05
it# 279 , res = 2.045371848871139e-05
it# 280 , res = 1.9865309905505643e-05
it# 281 , res = 1.929382854567591e-05
it# 282 , res = 1.873878745014132e-05
it# 283 , res = 1.819971366853517e-05
it# 284 , res = 1.7676147856397842e-05
it# 285 , res = 1.716764388320948e-05
it# 286 , res = 1.6673768453249704e-05
it# 287 , res = 1.619410073523324e-05
it# 288 , res = 1.5728232004674416e-05
it# 289 , res = 1.5275765294731234e-05
it# 290 , res = 1.483631505894444e-05
it# 291 , res = 1.4409506841717038e-05
it# 292 , res = 1.3994976960021093e-05
it# 293 , res = 1.3592372193119796e-05
it# 294 , res = 1.320134948164082e-05
it# 295 , res = 1.2821575635219436e-05
it# 296 , res = 1.2452727048639514e-05
it# 297 , res = 1.2094489426316882e-05
it# 298 , res = 1.1746557514059799e-05
it# 299 , res = 1.1408634839074057e-05
it# 300 , res = 1.1080433457799627e-05
it# 301 , res = 1.076167370972917e-05
it# 302 , res = 1.0452083980066552e-05
it# 303 , res = 1.0151400467235138e-05
it# 304 , res = 9.859366958945784e-06
it# 305 , res = 9.57573461378288e-06
it# 306 , res = 9.300261748523745e-06
it# 307 , res = 9.03271363290797e-06
it# 308 , res = 8.772862289409516e-06
it# 309 , res = 8.520486298492814e-06
it# 310 , res = 8.275370610927065e-06
it# 311 , res = 8.037306363628625e-06
it# 312 , res = 7.80609070210676e-06
it# 313 , res = 7.581526607684569e-06
it# 314 , res = 7.3634227294293486e-06
it# 315 , res = 7.151593220831914e-06
it# 316 , res = 6.945857582272965e-06
it# 317 , res = 6.746040506335811e-06
it# 318 , res = 6.551971729211386e-06
it# 319 , res = 6.3634858846987185e-06
it# 320 , res = 6.180422364230604e-06
it# 321 , res = 6.002625179414401e-06
it# 322 , res = 5.829942829251844e-06
it# 323 , res = 5.662228171337829e-06
it# 324 , res = 5.499338296010946e-06
it# 325 , res = 5.341134404691619e-06
it# 326 , res = 5.18748169213308e-06
it# 327 , res = 5.0382492308510415e-06
it# 328 , res = 4.893309859950636e-06
it# 329 , res = 4.752540076463343e-06
it# 330 , res = 4.615819930722382e-06
it# 331 , res = 4.483032923527409e-06
it# 332 , res = 4.3540659069707406e-06
it# 333 , res = 4.2288089883893184e-06
it# 334 , res = 4.107155436379231e-06
it# 335 , res = 3.989001590061696e-06
it# 336 , res = 3.874246770541752e-06
it# 337 , res = 3.762793195240183e-06
it# 338 , res = 3.654545894717841e-06
it# 339 , res = 3.5494126313937357e-06
it# 340 , res = 3.4473038212552832e-06
it# 341 , res = 3.3481324574600897e-06
it# 342 , res = 3.2518140361648322e-06
it# 343 , res = 3.158266484317209e-06
it# 344 , res = 3.067410090281997e-06
it# 345 , res = 2.9791674351496838e-06
it# 346 , res = 2.89346332757956e-06
it# 347 , res = 2.810224739077168e-06
it# 348 , res = 2.7293807420502645e-06
it# 349 , res = 2.6508624493427535e-06
it# 350 , res = 2.5746029555152176e-06
it# 351 , res = 2.5005372800512605e-06
it# 352 , res = 2.428602311338707e-06
it# 353 , res = 2.3587367538395157e-06
it# 354 , res = 2.290881075135971e-06
it# 355 , res = 2.224977455256337e-06
it# 356 , res = 2.1609697378851587e-06
it# 357 , res = 2.09880338203221e-06
it# 358 , res = 2.0384254158352543e-06
it# 359 , res = 1.9797843912149633e-06
it# 360 , res = 1.922830340027735e-06
it# 361 , res = 1.8675147318644292e-06
it# 362 , res = 1.8137904324192994e-06
it# 363 , res = 1.7616116631103953e-06
it# 364 , res = 1.710933962531527e-06
it# 365 , res = 1.6617141480427942e-06
it# 366 , res = 1.613910279542295e-06
it# 367 , res = 1.5674816235004455e-06
it# 368 , res = 1.5223886180371287e-06
it# 369 , res = 1.4785928392486356e-06
it# 370 , res = 1.4360569688455717e-06
it# 371 , res = 1.3947447621351398e-06
it# 372 , res = 1.354621016951224e-06
it# 373 , res = 1.3156515438365656e-06
it# 374 , res = 1.27780313688044e-06
it# 375 , res = 1.2410435455207698e-06
it# 376 , res = 1.205341446928076e-06
it# 377 , res = 1.1706664192604899e-06
it# 378 , res = 1.136988916100811e-06
it# 379 , res = 1.1042802409227474e-06
it# 380 , res = 1.0725125222997412e-06
it# 381 , res = 1.0416586913264574e-06
it# 382 , res = 1.0116924572853684e-06
it# 383 , res = 9.8258828602897e-07
it# 384 , res = 9.543213779265529e-07
it# 385 , res = 9.268676466974556e-07
it# 386 , res = 9.002036991007446e-07
it# 387 , res = 8.743068147159219e-07
it# 388 , res = 8.491549267811265e-07
it# 389 , res = 8.247266034991804e-07
it# 390 , res = 8.010010294761169e-07
it# 391 , res = 7.779579881130243e-07
it# 392 , res = 7.555778444798169e-07
it# 393 , res = 7.3384152850465e-07
it# 394 , res = 7.127305185337442e-07
it# 395 , res = 6.922268260496095e-07
it# 396 , res = 6.723129797802287e-07
it# 397 , res = 6.529720111683081e-07
it# 398 , res = 6.341874397303743e-07
it# 399 , res = 6.159432591347765e-07
it# 400 , res = 5.982239236260788e-07
it# 401 , res = 5.810143344228642e-07
it# 402 , res = 5.642998274385035e-07
it# 403 , res = 5.480661601137316e-07
it# 404 , res = 5.322994997562004e-07
it# 405 , res = 5.169864115400374e-07
it# 406 , res = 5.021138473923472e-07
it# 407 , res = 4.876691341834604e-07
it# 408 , res = 4.736399637616994e-07
it# 409 , res = 4.6001438156872165e-07
it# 410 , res = 4.4678077759507473e-07
it# 411 , res = 4.3392787533707315e-07
it# 412 , res = 4.214447228712117e-07
it# 413 , res = 4.093206832598425e-07
it# 414 , res = 3.975454257588739e-07
it# 415 , res = 3.86108916528486e-07
it# 416 , res = 3.7500141056865286e-07
it# 417 , res = 3.64213443151104e-07
it# 418 , res = 3.5373582183549464e-07
it# 419 , res = 3.4355961870799337e-07
it# 420 , res = 3.33676162690721e-07
it# 421 , res = 3.240770319704296e-07
it# 422 , res = 3.147540470982013e-07
it# 423 , res = 3.0569926406741627e-07
it# 424 , res = 2.9690496723295485e-07
it# 425 , res = 2.883636630262676e-07
it# 426 , res = 2.800680735357045e-07
it# 427 , res = 2.7201112970909123e-07
it# 428 , res = 2.641859666775173e-07
it# 429 , res = 2.565859162548356e-07
it# 430 , res = 2.4920450271428235e-07
it# 431 , res = 2.420354362896398e-07
it# 432 , res = 2.350726080120727e-07
it# 433 , res = 2.283100850899965e-07
it# 434 , res = 2.2174210506310934e-07
it# 435 , res = 2.1536307146633887e-07
it# 436 , res = 2.0916754860267494e-07
it# 437 , res = 2.0315025732306493e-07
it# 438 , res = 1.9730607041994718e-07
it# 439 , res = 1.916300078641151e-07
it# 440 , res = 1.8611723317234976e-07
it# 441 , res = 1.807630490609844e-07
it# 442 , res = 1.7556289285130043e-07
it# 443 , res = 1.7051233378054055e-07
it# 444 , res = 1.6560706817395565e-07
it# 445 , res = 1.6084291635080218e-07
it# 446 , res = 1.5621581872113273e-07
it# 447 , res = 1.5172183253947863e-07
it# 448 , res = 1.4735712849807978e-07
it# 449 , res = 1.4311798737663376e-07
it# 450 , res = 1.3900079700571098e-07
it# 451 , res = 1.3500204919236473e-07
it# 452 , res = 1.3111833657294248e-07
it# 453 , res = 1.2734634989198147e-07
it# 454 , res = 1.236828748658876e-07
it# 455 , res = 1.2012479007271662e-07
it# 456 , res = 1.1666906353064313e-07
it# 457 , res = 1.1331275071756048e-07
it# 458 , res = 1.1005299158742253e-07
it# 459 , res = 1.068870085409945e-07
it# 460 , res = 1.0381210399766806e-07
it# 461 , res = 1.0082565763267403e-07
it# 462 , res = 9.792512481974481e-08
it# 463 , res = 9.510803400955188e-08
it# 464 , res = 9.237198461100925e-08
it# 465 , res = 8.971464533593739e-08
it# 466 , res = 8.713375192852615e-08
it# 467 , res = 8.462710516784603e-08
it# 468 , res = 8.21925689797035e-08
it# 469 , res = 7.982806924717622e-08
it# 470 , res = 7.753159093898629e-08
it# 471 , res = 7.530117731114372e-08
it# 472 , res = 7.313492768513725e-08
it# 473 , res = 7.103099635939458e-08
it# 474 , res = 6.898759056916919e-08
it# 475 , res = 6.700296899119797e-08
it# 476 , res = 6.507544065674031e-08
it# 477 , res = 6.32033630646422e-08
it# 478 , res = 6.138514104116617e-08
it# 479 , res = 5.961922525144198e-08
it# 480 , res = 5.790411099528519e-08
it# 481 , res = 5.6238336905067064e-08
it# 482 , res = 5.462048328716683e-08
it# 483 , res = 5.304917192087867e-08
it# 484 , res = 5.152306366387732e-08
it# 485 , res = 5.0040858255584435e-08
it# 486 , res = 4.8601292618442114e-08
it# 487 , res = 4.720314016707623e-08
it# 488 , res = 4.5845209401926625e-08
it# 489 , res = 4.452634332717489e-08
it# 490 , res = 4.324541814437533e-08
it# 491 , res = 4.200134248162821e-08
it# 492 , res = 4.0793056035516903e-08
it# 493 , res = 3.961952930255926e-08
it# 494 , res = 3.8479762458281354e-08
it# 495 , res = 3.737278414058256e-08
it# 496 , res = 3.6297651188265186e-08
it# 497 , res = 3.5253447406999905e-08
it# 498 , res = 3.4239283255189594e-08
it# 499 , res = 3.325429409712161e-08

Implement a forward-backward GS preconditioner

The same point Jacobi smoother object is also able to perform a Gauss-Seidel iteration after reversing the ordering of the points, i.e., a backward Gauss-Seidel sweep. One can combine the forward and backward sweeps to construct a symmetric preconditioner, often called the symmetrized Gauss-Seidel preconditioner. This offers a good simple illustration of how to construct NGSolve preconditioners from within python.

class SymmetricGS(BaseMatrix):
    def __init__ (self, smoother):
        super(SymmetricGS, self).__init__()
        self.smoother = smoother
    def Mult (self, x, y):
        y[:] = 0.0
        self.smoother.Smooth(y, x)
        self.smoother.SmoothBack(y,x)
    def Height (self):
        return self.smoother.height
    def Width (self):
        return self.smoother.height

preGS = SymmetricGS(preJpoint)
solvers.CG(mat=a.mat, pre=preGS, rhs=f.vec, sol=gfu.vec)
Draw(gfu)

iteration 0 error = 0.09386193808775962
iteration 1 error = 0.12289904446040575
iteration 2 error = 0.08558846385984852
iteration 3 error = 0.05378637509125533
iteration 4 error = 0.027995133684947375
iteration 5 error = 0.014969155114208058
iteration 6 error = 0.006689317102568868
iteration 7 error = 0.0028248604314505604
iteration 8 error = 0.0010332508846613
iteration 9 error = 0.00040318420091610223
iteration 10 error = 0.00012615348412530394
iteration 11 error = 5.134068138544936e-05
iteration 12 error = 1.8626752985537923e-05
iteration 13 error = 7.865715129900787e-06
iteration 14 error = 3.310916766942577e-06
iteration 15 error = 1.5461433004860907e-06
iteration 16 error = 6.174156961107723e-07
iteration 17 error = 3.200817049657394e-07
iteration 18 error = 1.2436962307106972e-07
iteration 19 error = 5.906688350434222e-08
iteration 20 error = 2.2070416091525614e-08
iteration 21 error = 7.3242103990364256e-09
iteration 22 error = 3.0122986220535967e-09
iteration 23 error = 8.584333680651098e-10
iteration 24 error = 3.657080601088919e-10
iteration 25 error = 1.2692446936158058e-10
iteration 26 error = 4.461371781815044e-11
iteration 27 error = 1.9750920943828053e-11
iteration 28 error = 5.534324373492611e-12
iteration 29 error = 2.307643414711425e-12
iteration 30 error = 6.724334725535219e-13
iteration 31 error = 1.9582752398903902e-13
iteration 32 error = 7.628337312031245e-14

lams = EigenValues_Preconditioner(mat=a.mat, pre=preGS)
max(lams)/min(lams)

20.294597459170067

Note that the condition number now is better than that of the system preconditioned by point Jacobi.

A Block Jacobi preconditioner

The point Jacobi preconditioner is based on inverses of 1 x 1 diagonal blocks. Condition numbers can be improved by using larger blocks. It is possible to group dofs into blocks within python and construct an NGSolve preconditioner based on the blocks.

Here is an example that constructs vertex-based blocks, using the mesh querying techniques we learnt in 1.8.

def VertexPatchBlocks(mesh, fes):
    blocks = []
    freedofs = fes.FreeDofs()
    for v in mesh.vertices:
        vdofs = set()
        for el in mesh[v].elements:
            vdofs |= set(d for d in fes.GetDofNrs(el)
                         if freedofs[d])
        blocks.append(vdofs)
    return blocks

blocks = VertexPatchBlocks(mesh, fes)
print(blocks)

[{160, 873, 159}, {13, 207, 48, 145, 146, 208, 143, 211, 212, 144, 889, 891}, {257, 2, 258, 908, 147, 148, 21, 22, 149, 150}, {65, 312, 924, 315, 926, 316, 311, 152, 151, 155, 156, 30}, {160, 161, 162, 165, 166, 359, 40, 873, 874, 360, 942, 159}, {161, 162, 165, 166, 167, 40, 168, 874, 41, 171, 172, 875, 363, 364, 876}, {167, 168, 41, 42, 875, 171, 172, 173, 174, 877, 177, 178, 878, 369, 370}, {375, 376, 42, 43, 877, 173, 174, 879, 177, 178, 179, 180, 880, 183, 184}, {382, 43, 44, 879, 881, 882, 179, 180, 381, 183, 184, 185, 186, 189, 190}, {192, 195, 196, 387, 388, 44, 45, 881, 883, 884, 185, 186, 189, 190, 191}, {192, 195, 196, 197, 198, 201, 202, 393, 394, 45, 46, 883, 885, 886, 191}, {197, 198, 201, 202, 203, 204, 205, 206, 399, 400, 46, 47, 885, 887, 888}, {203, 204, 205, 206, 207, 208, 145, 146, 405, 406, 47, 48, 887, 889, 890}, {13, 14, 143, 144, 145, 146, 211, 212, 209, 210, 213, 214, 217, 218, 411, 412, 48, 49, 891, 892, 893}, {13, 14, 15, 209, 210, 213, 214, 215, 216, 217, 218, 219, 220, 223, 224, 416, 415, 49, 50, 892, 894, 895}, {896, 897, 14, 15, 16, 215, 216, 219, 220, 221, 222, 223, 224, 225, 226, 229, 230, 421, 422, 50, 51, 894}, {896, 898, 899, 15, 16, 17, 221, 222, 225, 226, 227, 228, 229, 230, 231, 232, 235, 236, 427, 428, 51, 52}, {898, 900, 901, 16, 17, 18, 227, 228, 231, 232, 233, 234, 235, 236, 237, 238, 241, 242, 433, 52, 53, 434}, {900, 902, 903, 17, 18, 19, 439, 440, 233, 234, 237, 238, 239, 240, 241, 242, 243, 244, 53, 54, 247, 248}, {446, 902, 904, 905, 18, 19, 20, 55, 249, 239, 240, 243, 244, 245, 54, 247, 248, 246, 250, 253, 254, 445}, {256, 259, 260, 451, 452, 904, 906, 907, 19, 20, 21, 245, 246, 55, 56, 249, 250, 251, 252, 253, 254, 255}, {256, 257, 2, 259, 260, 258, 263, 264, 906, 908, 147, 20, 21, 148, 22, 149, 150, 943, 56, 251, 252, 255}, {257, 2, 258, 259, 261, 262, 260, 263, 265, 266, 264, 908, 269, 270, 909, 457, 458, 147, 148, 21, 22, 149, 150, 23, 943, 944, 56, 57}, {261, 262, 265, 266, 267, 268, 269, 270, 909, 910, 271, 272, 275, 276, 461, 22, 23, 24, 462, 911, 57, 58}, {267, 268, 910, 271, 272, 912, 273, 275, 276, 274, 277, 23, 24, 281, 278, 25, 282, 465, 466, 913, 58, 59}, {471, 472, 912, 273, 274, 914, 915, 277, 278, 279, 24, 281, 25, 282, 284, 280, 26, 283, 288, 287, 59, 60}, {914, 916, 917, 279, 280, 25, 26, 283, 284, 285, 27, 287, 288, 289, 290, 286, 477, 293, 294, 478, 60, 61}, {916, 918, 919, 26, 27, 28, 285, 286, 289, 290, 291, 292, 293, 294, 295, 296, 483, 484, 299, 300, 61, 62}, {918, 920, 921, 27, 28, 29, 291, 292, 295, 296, 297, 298, 299, 300, 301, 302, 489, 490, 305, 306, 62, 63}, {64, 920, 922, 923, 28, 29, 30, 297, 298, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 495, 496, 63}, {64, 65, 151, 152, 922, 155, 924, 29, 30, 156, 925, 303, 304, 307, 308, 309, 310, 311, 312, 501, 502}, {65, 321, 66, 322, 155, 156, 926, 927, 928, 508, 507, 315, 316, 317, 318}, {512, 321, 66, 322, 323, 67, 324, 327, 328, 927, 929, 930, 317, 318, 511}, {323, 67, 324, 68, 327, 328, 329, 330, 517, 518, 333, 334, 929, 931, 932}, {68, 69, 329, 330, 523, 524, 333, 334, 335, 336, 339, 340, 931, 933, 934}, {69, 70, 335, 336, 529, 530, 339, 340, 341, 342, 345, 346, 933, 935, 936}, {70, 71, 341, 342, 535, 536, 345, 346, 347, 348, 351, 352, 935, 937, 938}, {71, 72, 347, 348, 541, 542, 351, 352, 353, 354, 357, 358, 937, 939, 940}, {353, 354, 357, 358, 72, 361, 362, 939, 941}, {72, 159, 160, 161, 162, 357, 358, 359, 40, 873, 361, 362, 360, 941, 942, 365, 366, 945}, {159, 160, 161, 162, 547, 548, 165, 166, 167, 40, 41, 168, 942, 945, 947, 72, 74, 359, 360, 361, 874, 363, 364, 876, 362, 365, 366, 367, 368, 1006, 373, 374}, {165, 166, 167, 168, 41, 40, 171, 172, 42, 173, 174, 551, 552, 946, 947, 948, 73, 74, 875, 363, 364, 876, 367, 878, 369, 370, 371, 372, 368, 373, 374, 377, 378}, {385, 386, 41, 42, 43, 171, 172, 173, 174, 553, 177, 178, 179, 180, 946, 73, 100, 554, 1001, 1002, 877, 878, 880, 369, 370, 371, 372, 375, 376, 377, 378, 379, 380}, {384, 385, 386, 389, 390, 42, 43, 44, 177, 178, 179, 180, 949, 565, 183, 184, 185, 186, 566, 75, 100, 1001, 1003, 879, 880, 882, 375, 376, 379, 380, 381, 382, 383}, {384, 387, 388, 389, 390, 391, 392, 395, 396, 43, 44, 45, 563, 564, 949, 950, 183, 184, 185, 186, 951, 189, 190, 191, 192, 75, 76, 881, 882, 884, 381, 382, 383}, {387, 388, 391, 392, 393, 394, 395, 396, 397, 398, 401, 402, 44, 45, 46, 950, 952, 953, 571, 572, 189, 190, 191, 192, 195, 196, 197, 198, 76, 77, 883, 884, 886}, {393, 394, 397, 398, 399, 400, 401, 402, 403, 404, 407, 408, 45, 46, 47, 952, 954, 955, 577, 578, 195, 196, 197, 198, 201, 202, 203, 204, 77, 78, 885, 886, 888}, {399, 400, 403, 404, 405, 406, 407, 408, 409, 410, 413, 414, 46, 47, 48, 954, 956, 957, 583, 584, 201, 202, 203, 204, 205, 206, 79, 207, 208, 78, 887, 888, 890}, {13, 143, 144, 145, 146, 405, 406, 409, 410, 411, 412, 413, 414, 417, 418, 47, 48, 49, 956, 958, 205, 206, 79, 207, 208, 211, 212, 213, 214, 889, 890, 891, 893}, {13, 14, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 425, 426, 48, 49, 50, 958, 960, 961, 587, 588, 79, 209, 210, 211, 212, 213, 214, 81, 217, 218, 219, 220, 892, 893, 895}, {897, 14, 15, 415, 416, 419, 420, 421, 422, 423, 424, 425, 426, 429, 430, 49, 50, 51, 959, 960, 962, 591, 80, 81, 592, 215, 216, 217, 218, 219, 220, 223, 224, 225, 226, 894, 895}, {896, 897, 899, 15, 16, 421, 422, 423, 424, 427, 428, 429, 430, 431, 432, 50, 51, 52, 435, 436, 959, 963, 964, 80, 593, 82, 594, 221, 222, 223, 224, 225, 226, 229, 230, 231, 232}, {898, 899, 901, 16, 17, 427, 428, 431, 432, 433, 434, 51, 52, 53, 435, 436, 438, 437, 441, 442, 963, 965, 966, 82, 83, 601, 602, 227, 228, 229, 230, 231, 232, 235, 236, 237, 238}, {900, 901, 903, 17, 18, 433, 434, 52, 53, 54, 439, 440, 441, 442, 438, 437, 443, 444, 447, 448, 965, 967, 968, 83, 84, 607, 608, 233, 234, 235, 236, 237, 238, 241, 242, 243, 244}, {902, 903, 905, 18, 19, 53, 54, 55, 439, 440, 443, 444, 445, 446, 447, 448, 449, 450, 453, 454, 967, 969, 970, 84, 85, 613, 614, 239, 240, 241, 242, 243, 244, 247, 248, 249, 250}, {256, 904, 905, 907, 19, 20, 54, 55, 56, 445, 446, 449, 450, 451, 452, 453, 454, 455, 456, 969, 459, 460, 971, 972, 85, 86, 619, 620, 245, 246, 247, 248, 249, 250, 253, 254, 255}, {256, 257, 258, 259, 260, 263, 264, 265, 906, 907, 266, 20, 21, 22, 943, 944, 55, 56, 57, 451, 452, 455, 456, 457, 458, 459, 460, 971, 973, 463, 464, 86, 251, 252, 253, 254, 255}, {261, 262, 263, 264, 265, 266, 269, 270, 909, 271, 272, 911, 22, 23, 944, 56, 57, 58, 457, 458, 459, 460, 461, 462, 463, 464, 973, 977, 467, 468, 86}, {267, 268, 269, 910, 271, 272, 270, 911, 275, 276, 277, 278, 23, 24, 913, 57, 58, 59, 461, 462, 975, 463, 465, 466, 464, 977, 469, 470, 467, 88, 468, 86, 475, 476, 1004, 625, 626}, {912, 273, 274, 275, 276, 277, 278, 913, 24, 281, 25, 282, 283, 284, 58, 59, 60, 974, 975, 976, 465, 466, 469, 470, 471, 472, 473, 474, 87, 475, 88, 476, 479, 480, 627, 628, 915}, {914, 915, 917, 279, 280, 25, 26, 283, 284, 281, 282, 287, 288, 289, 290, 59, 60, 61, 974, 978, 979, 471, 472, 473, 474, 87, 89, 477, 478, 479, 480, 481, 482, 485, 486, 629, 630}, {640, 916, 917, 919, 26, 27, 285, 286, 287, 288, 289, 290, 293, 294, 295, 296, 60, 61, 62, 978, 980, 981, 89, 90, 477, 478, 481, 482, 483, 484, 485, 486, 487, 488, 491, 492, 639}, {645, 646, 918, 919, 921, 27, 28, 291, 292, 293, 294, 295, 296, 299, 300, 301, 302, 61, 62, 63, 980, 982, 983, 90, 91, 483, 484, 487, 488, 489, 490, 491, 492, 493, 494, 497, 498}, {651, 652, 920, 921, 923, 28, 29, 297, 298, 299, 300, 301, 302, 305, 306, 307, 308, 62, 63, 64, 982, 984, 985, 91, 92, 489, 490, 493, 494, 495, 496, 497, 498, 499, 500, 503, 504}, {657, 658, 922, 923, 29, 30, 925, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 63, 64, 65, 984, 986, 987, 92, 93, 495, 496, 499, 500, 501, 502, 503, 504, 505, 506, 509, 510}, {513, 514, 151, 152, 155, 924, 156, 30, 925, 926, 928, 309, 310, 311, 312, 315, 316, 317, 318, 64, 65, 66, 986, 988, 93, 501, 502, 505, 506, 507, 508, 509, 510}, {512, 513, 514, 515, 516, 521, 522, 663, 664, 927, 928, 930, 315, 316, 317, 318, 321, 66, 322, 65, 67, 323, 324, 988, 93, 990, 95, 991, 507, 508, 509, 510, 511}, {512, 515, 516, 517, 518, 519, 520, 521, 522, 525, 526, 665, 666, 929, 930, 932, 321, 66, 323, 67, 324, 322, 327, 328, 68, 329, 330, 989, 94, 990, 95, 992, 511}, {517, 518, 519, 520, 523, 524, 525, 526, 527, 528, 531, 532, 667, 668, 931, 932, 934, 67, 68, 69, 327, 328, 329, 330, 333, 334, 335, 336, 989, 94, 96, 993, 994}, {523, 524, 527, 528, 529, 530, 531, 532, 533, 534, 537, 538, 675, 676, 933, 934, 936, 68, 69, 70, 333, 334, 335, 336, 339, 340, 341, 342, 96, 993, 97, 995, 996}, {529, 530, 533, 534, 535, 536, 537, 538, 539, 540, 543, 544, 935, 936, 681, 938, 682, 69, 70, 71, 339, 340, 341, 342, 345, 346, 347, 348, 97, 98, 995, 997, 998}, {535, 536, 539, 540, 541, 542, 543, 544, 545, 546, 549, 550, 937, 938, 940, 687, 688, 70, 71, 72, 345, 346, 347, 348, 351, 352, 353, 354, 98, 99, 997, 999, 1000}, {541, 542, 545, 546, 547, 548, 549, 550, 40, 939, 940, 941, 559, 560, 945, 71, 72, 74, 351, 352, 353, 354, 99, 357, 358, 359, 360, 361, 362, 999, 365, 366, 1005, 367, 368, 1006}, {1027, 1046, 1048, 551, 552, 41, 42, 553, 554, 555, 556, 557, 558, 561, 562, 946, 948, 699, 700, 73, 74, 119, 100, 1002, 111, 369, 370, 371, 372, 373, 374, 759, 760, 377, 378, 379, 380}, {1026, 1027, 547, 548, 549, 550, 551, 40, 41, 552, 555, 556, 559, 560, 561, 562, 947, 948, 693, 694, 72, 73, 74, 99, 363, 364, 365, 1006, 1005, 367, 368, 366, 371, 372, 373, 374, 111}, {384, 385, 386, 389, 390, 391, 392, 1043, 1044, 43, 44, 563, 564, 949, 565, 567, 568, 951, 566, 569, 570, 697, 698, 575, 576, 709, 710, 75, 76, 100, 102, 1003, 1008, 118, 381, 382, 383}, {387, 388, 389, 390, 391, 392, 395, 396, 397, 398, 44, 45, 563, 564, 950, 951, 567, 953, 568, 571, 572, 573, 574, 575, 576, 701, 702, 579, 580, 75, 76, 77, 101, 102, 1007, 1008, 1009}, {1025, 393, 394, 395, 396, 397, 398, 401, 402, 403, 404, 1042, 45, 46, 952, 953, 955, 571, 572, 573, 574, 704, 577, 578, 579, 580, 581, 582, 703, 585, 586, 76, 77, 78, 101, 103, 1007}, {1025, 399, 400, 401, 402, 403, 404, 407, 408, 409, 410, 46, 47, 954, 955, 957, 577, 578, 581, 582, 583, 584, 585, 586, 77, 78, 79, 589, 590, 103, 1010}, {405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 417, 418, 419, 420, 47, 48, 49, 956, 957, 958, 961, 583, 584, 585, 586, 587, 588, 589, 78, 79, 590, 81, 599, 600, 103, 1010, 1011}, {1045, 1052, 421, 422, 423, 424, 425, 426, 429, 430, 431, 432, 50, 51, 959, 962, 964, 715, 716, 591, 80, 81, 592, 593, 82, 594, 597, 598, 595, 596, 599, 600, 605, 606, 103, 106, 1015}, {1045, 415, 416, 417, 418, 419, 420, 423, 424, 425, 426, 49, 50, 960, 961, 962, 587, 588, 589, 590, 79, 591, 81, 80, 592, 595, 596, 599, 600, 103, 1011}, {427, 428, 429, 430, 431, 432, 51, 52, 435, 436, 437, 438, 963, 964, 966, 719, 80, 593, 82, 594, 83, 597, 598, 720, 601, 602, 603, 604, 605, 606, 609, 610, 104, 106, 1012, 1015, 1016}, {433, 434, 435, 52, 53, 438, 437, 436, 441, 442, 443, 444, 965, 966, 968, 717, 718, 82, 83, 84, 601, 602, 603, 604, 607, 608, 609, 610, 611, 612, 615, 104, 616, 105, 1012, 1013, 1014}, {53, 54, 439, 440, 441, 442, 443, 444, 447, 448, 449, 450, 967, 968, 970, 83, 84, 85, 725, 726, 607, 608, 611, 612, 613, 614, 615, 616, 105, 617, 618, 107, 623, 624, 1013, 1017, 1018}, {1041, 54, 55, 445, 446, 447, 448, 449, 450, 453, 454, 455, 456, 969, 970, 972, 84, 85, 86, 88, 613, 614, 635, 617, 618, 619, 620, 107, 621, 623, 624, 622, 625, 626, 1017, 1019, 636}, {55, 56, 57, 58, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 971, 972, 973, 463, 464, 461, 462, 977, 85, 86, 467, 468, 469, 470, 88, 619, 620, 1004, 621, 622, 625, 626, 1019}, {1024, 643, 644, 637, 59, 60, 974, 976, 979, 471, 87, 89, 472, 473, 474, 475, 88, 476, 479, 480, 481, 482, 747, 748, 109, 110, 627, 628, 629, 630, 631, 632, 633, 634, 1021, 638, 1023}, {1041, 1049, 107, 58, 59, 975, 976, 465, 466, 467, 468, 469, 470, 87, 88, 473, 474, 475, 476, 86, 85, 735, 736, 635, 619, 1004, 620, 621, 622, 110, 625, 626, 627, 628, 623, 624, 633, 634, 1019, 636, 637, 638, 1023}, {640, 641, 642, 643, 644, 647, 648, 60, 61, 978, 979, 981, 87, 89, 90, 477, 478, 479, 480, 481, 482, 739, 740, 485, 486, 487, 488, 108, 109, 629, 630, 631, 632, 1020, 1021, 1022, 639}, {640, 641, 642, 1028, 645, 646, 647, 648, 649, 650, 1029, 653, 654, 61, 62, 980, 981, 983, 89, 90, 91, 483, 484, 485, 486, 487, 488, 741, 742, 491, 492, 493, 494, 108, 112, 1020, 639}, {1028, 645, 646, 1030, 1031, 649, 650, 651, 652, 653, 654, 655, 656, 661, 662, 62, 63, 982, 983, 985, 90, 91, 92, 489, 490, 491, 492, 493, 494, 112, 497, 498, 499, 500, 113, 761, 762}, {1030, 1032, 651, 652, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 1050, 673, 674, 63, 64, 984, 985, 91, 92, 987, 93, 95, 495, 496, 497, 498, 499, 500, 113, 503, 504, 505, 506}, {513, 514, 515, 516, 1032, 657, 658, 659, 660, 663, 664, 64, 65, 66, 986, 987, 92, 93, 988, 991, 95, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510}, {768, 517, 518, 519, 520, 521, 522, 1036, 525, 526, 527, 528, 665, 666, 667, 668, 1051, 670, 671, 672, 673, 674, 669, 1053, 679, 680, 67, 68, 989, 94, 95, 992, 96, 994, 113, 116, 767}, {512, 513, 514, 515, 516, 519, 520, 521, 522, 1032, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 1050, 1051, 669, 670, 673, 674, 66, 67, 92, 93, 990, 95, 991, 992, 94, 113, 511}, {773, 774, 1033, 523, 524, 525, 526, 527, 528, 1036, 1037, 531, 532, 533, 534, 667, 668, 671, 672, 675, 676, 677, 678, 679, 680, 683, 684, 68, 69, 94, 96, 993, 994, 97, 996, 114, 116}, {771, 772, 1033, 1034, 1035, 529, 530, 531, 532, 533, 534, 537, 538, 539, 540, 675, 676, 677, 678, 681, 682, 683, 684, 685, 686, 689, 690, 69, 70, 96, 97, 98, 995, 996, 998, 114, 115}, {1034, 779, 780, 1038, 1039, 535, 536, 537, 538, 539, 540, 543, 544, 545, 546, 681, 682, 685, 686, 687, 688, 689, 690, 691, 692, 695, 696, 70, 71, 97, 98, 99, 997, 998, 1000, 115, 117}, {1026, 1038, 1040, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 687, 688, 559, 560, 561, 562, 693, 694, 691, 692, 695, 696, 71, 72, 74, 98, 99, 999, 1000, 1005, 111, 117, 757, 758}, {384, 385, 386, 1044, 1048, 793, 794, 553, 42, 43, 554, 557, 558, 565, 566, 1079, 569, 698, 697, 570, 699, 700, 73, 75, 119, 100, 1001, 1002, 1003, 118, 375, 376, 377, 378, 379, 380, 383}, {129, 1042, 1062, 1074, 1075, 571, 572, 701, 702, 703, 576, 704, 705, 579, 580, 581, 582, 706, 573, 574, 707, 715, 76, 77, 716, 575, 711, 712, 733, 734, 101, 102, 103, 106, 708, 1007, 1009}, {129, 130, 1043, 795, 796, 563, 564, 1075, 1077, 567, 568, 569, 570, 1078, 701, 702, 573, 576, 574, 575, 707, 708, 709, 710, 711, 712, 713, 714, 75, 76, 847, 848, 101, 102, 1008, 1009, 118}, {1025, 715, 1042, 716, 1045, 1052, 1062, 703, 704, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 78, 79, 590, 81, 77, 591, 80, 592, 595, 599, 600, 596, 597, 598, 705, 706, 101, 103, 106, 1010, 1011}, {1056, 1057, 1058, 807, 808, 717, 718, 719, 720, 721, 82, 83, 722, 723, 724, 727, 728, 601, 602, 603, 604, 605, 606, 731, 732, 609, 610, 611, 612, 104, 105, 106, 1012, 1014, 1016, 121, 122}, {1056, 1059, 1060, 809, 810, 717, 718, 721, 722, 83, 84, 725, 726, 727, 728, 729, 730, 607, 608, 609, 610, 611, 612, 737, 738, 615, 616, 105, 104, 617, 107, 618, 1013, 1014, 121, 1018, 123}, {129, 1052, 1057, 1062, 815, 816, 1074, 1076, 703, 704, 705, 706, 707, 708, 715, 716, 719, 80, 593, 82, 594, 720, 597, 598, 595, 596, 723, 724, 603, 604, 605, 606, 731, 732, 733, 734, 101, 103, 104, 106, 1015, 1016, 122}, {1041, 1049, 1059, 1061, 84, 85, 725, 726, 88, 729, 730, 735, 736, 737, 738, 613, 614, 615, 616, 617, 618, 107, 105, 621, 622, 623, 624, 110, 123, 753, 754, 1017, 1018, 635, 636, 637, 638}, {640, 641, 642, 643, 644, 1029, 647, 648, 649, 650, 1020, 766, 1063, 1067, 1068, 825, 826, 89, 90, 739, 740, 741, 742, 743, 744, 745, 746, 108, 109, 749, 750, 112, 126, 124, 765, 1022, 639}, {1024, 641, 642, 643, 644, 1063, 1064, 1066, 823, 824, 125, 87, 89, 739, 740, 743, 744, 747, 108, 109, 748, 110, 749, 750, 751, 752, 755, 629, 630, 631, 632, 633, 634, 756, 124, 1021, 1022}, {1024, 1049, 1061, 1064, 1065, 819, 820, 125, 87, 88, 756, 735, 736, 737, 738, 635, 747, 748, 109, 110, 107, 751, 753, 754, 627, 628, 752, 755, 631, 632, 633, 634, 123, 636, 637, 638, 1023}, {1026, 1027, 1040, 789, 1046, 1047, 790, 551, 552, 555, 556, 557, 558, 559, 560, 561, 562, 693, 694, 695, 696, 73, 74, 99, 759, 111, 760, 117, 757, 119, 758}, {769, 770, 1028, 645, 646, 647, 648, 649, 650, 1029, 1031, 653, 654, 655, 656, 1055, 801, 802, 1067, 1069, 90, 91, 126, 741, 742, 745, 746, 108, 112, 113, 120, 761, 762, 763, 764, 765, 766}, {768, 769, 770, 1030, 1031, 651, 652, 653, 654, 655, 656, 785, 786, 659, 660, 661, 662, 665, 1050, 666, 1051, 669, 670, 671, 672, 673, 674, 1053, 1054, 1055, 91, 92, 94, 95, 112, 113, 116, 120, 761, 762, 763, 764, 767}, {128, 771, 772, 773, 774, 775, 776, 1033, 777, 1035, 778, 1037, 781, 782, 787, 788, 675, 676, 677, 678, 679, 680, 683, 684, 685, 686, 1070, 1071, 1073, 837, 838, 96, 97, 114, 115, 116, 127}, {771, 772, 131, 775, 776, 1034, 1035, 779, 780, 781, 1039, 782, 783, 784, 791, 792, 681, 682, 683, 684, 685, 686, 1070, 689, 690, 691, 692, 1080, 1081, 839, 840, 97, 98, 114, 115, 117, 127}, {768, 769, 770, 128, 773, 774, 777, 778, 1036, 1037, 785, 786, 787, 788, 667, 668, 1053, 669, 671, 672, 670, 1054, 803, 804, 677, 678, 679, 680, 1071, 1072, 94, 96, 113, 114, 116, 120, 767}, {131, 779, 780, 1038, 1039, 1040, 783, 784, 789, 790, 1047, 791, 792, 799, 800, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 1081, 1082, 98, 99, 759, 111, 760, 115, 117, 757, 119, 758}, {130, 131, 1043, 1044, 793, 794, 795, 796, 797, 798, 799, 800, 1077, 565, 566, 567, 568, 697, 698, 699, 700, 569, 570, 1079, 1087, 1091, 709, 710, 713, 714, 75, 851, 852, 100, 102, 118, 119}, {131, 789, 1046, 1047, 790, 1048, 793, 794, 791, 792, 797, 799, 800, 798, 553, 554, 555, 556, 557, 558, 1079, 697, 698, 699, 700, 1082, 1087, 73, 100, 759, 111, 757, 117, 119, 760, 758, 118}, {768, 769, 770, 128, 135, 766, 785, 786, 787, 788, 1054, 1055, 801, 802, 803, 804, 805, 806, 1069, 1072, 835, 836, 1092, 1093, 843, 844, 112, 113, 116, 120, 761, 762, 763, 764, 765, 126, 767}, {132, 133, 1056, 1058, 1060, 807, 808, 809, 810, 811, 812, 813, 814, 817, 818, 821, 822, 1083, 1088, 1094, 717, 718, 721, 722, 723, 724, 727, 728, 729, 730, 859, 860, 104, 105, 121, 122, 123}, {129, 133, 1057, 1058, 807, 808, 813, 814, 815, 816, 817, 818, 1076, 1085, 1088, 719, 720, 721, 722, 723, 724, 849, 850, 731, 732, 733, 734, 104, 106, 121, 122}, {132, 1059, 1060, 1061, 809, 810, 1065, 811, 812, 819, 820, 821, 822, 1083, 1084, 829, 830, 725, 726, 727, 728, 729, 730, 735, 736, 737, 738, 105, 107, 110, 753, 754, 755, 756, 121, 123, 125}, {134, 1063, 1066, 1068, 823, 824, 825, 826, 827, 828, 831, 832, 1089, 833, 834, 1098, 739, 740, 743, 744, 745, 746, 108, 109, 749, 750, 751, 752, 124, 125, 126}, {132, 134, 1064, 1065, 1066, 819, 820, 821, 822, 823, 824, 827, 1084, 829, 830, 828, 832, 1089, 831, 1090, 861, 862, 747, 748, 109, 110, 751, 752, 753, 754, 755, 756, 750, 749, 123, 124, 125}, {134, 135, 766, 801, 802, 805, 806, 1067, 1068, 1069, 825, 826, 827, 828, 833, 834, 835, 836, 1092, 1098, 1099, 867, 868, 741, 742, 743, 744, 745, 746, 108, 112, 124, 120, 763, 764, 765, 126}, {128, 771, 772, 131, 775, 776, 777, 778, 136, 781, 782, 783, 784, 1070, 1073, 1080, 837, 838, 839, 840, 841, 842, 1095, 1096, 845, 846, 857, 858, 114, 115, 127}, {128, 773, 774, 775, 776, 777, 778, 135, 136, 785, 786, 787, 788, 803, 804, 805, 806, 1071, 1072, 1073, 837, 838, 1093, 1096, 841, 842, 843, 844, 845, 846, 1103, 871, 872, 114, 116, 120, 127}, {129, 130, 133, 815, 816, 817, 1074, 1075, 1076, 818, 1078, 701, 702, 1085, 1086, 705, 706, 707, 708, 711, 712, 713, 714, 847, 848, 849, 850, 853, 854, 731, 732, 733, 734, 101, 102, 106, 122}, {129, 130, 131, 133, 136, 795, 796, 797, 798, 1077, 1078, 1086, 1091, 709, 710, 711, 712, 713, 714, 1097, 1100, 847, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 865, 866, 102, 118}, {130, 131, 136, 779, 780, 781, 782, 783, 784, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799, 800, 1080, 1081, 1082, 1087, 1091, 839, 840, 841, 842, 1095, 1097, 851, 852, 855, 856, 857, 858, 115, 117, 118, 119, 127}, {132, 133, 134, 809, 810, 811, 812, 813, 814, 819, 820, 821, 822, 1083, 1084, 829, 830, 831, 832, 1090, 1094, 1101, 859, 860, 861, 862, 863, 864, 121, 123, 125}, {129, 130, 132, 133, 134, 136, 807, 808, 811, 812, 813, 814, 815, 816, 817, 818, 1085, 1086, 1088, 1094, 1100, 1101, 1102, 847, 848, 849, 850, 853, 854, 855, 856, 859, 860, 861, 862, 863, 864, 865, 866, 869, 870, 121, 122}, {132, 133, 134, 135, 136, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 1089, 1090, 833, 834, 835, 836, 1098, 1099, 1101, 1102, 1104, 859, 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 124, 125, 126}, {128, 134, 135, 136, 801, 802, 803, 804, 805, 806, 833, 834, 835, 836, 1092, 1093, 843, 844, 1099, 845, 846, 1103, 1104, 867, 868, 869, 870, 871, 872, 120, 126}, {128, 130, 131, 133, 134, 135, 136, 844, 1103, 1104, 837, 838, 839, 840, 841, 842, 1095, 1096, 845, 846, 1097, 1100, 1102, 843, 851, 852, 853, 854, 855, 856, 857, 858, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 127}]

CreateBlockSmoother can now take these blocks and construct a block Jacobi preconditioner.

blockjac = a.mat.CreateBlockSmoother(blocks)
lams = EigenValues_Preconditioner(mat=a.mat, pre=blockjac)
max(lams)/min(lams)

35.59893732672581

Multiplicative smoothers and its symmetrized version often yield better condition numbers in practice. We can apply the same code we wrote above for symmetrization (SymmetricGS) to the block smoother:

blockgs = SymmetricGS(blockjac)
lams = EigenValues_Preconditioner(mat=a.mat, pre=blockgs)
max(lams)/min(lams)

2.9828997439629763

Add a coarse grid correction

Dependence of the condition number on degrees of freedom can often be reduced by preconditioners that appropriately use a coarse grid correction. We can experiment with coarse grid corrections using NGSolve’s python interface. Here is an example on how to precondition with a coarse grid correction made using the lowest order subspace of fes.

In the example below, note that we use fes.GetDofNrs again. Previously we used it with argument el of type ElementId, while now we use it with an argument v of type MeshNode.

def VertexDofs(mesh, fes):
    vertexdofs = BitArray(fes.ndof)
    vertexdofs[:] = False
    for v in mesh.vertices:
        for d in fes.GetDofNrs(v):
            vertexdofs[d] = True
    vertexdofs &= fes.FreeDofs()
    return vertexdofs

vertexdofs = VertexDofs(mesh, fes)
print(vertexdofs)   # bit array, printed 50 chars/line

0: 00100000000001111111111111111110000000001111111111
50: 11111111111111111111111111111111111111111111111111
100: 11111111111111111111111111111111111110000000000000
150: 00000000000000000000000000000000000000000000000000
200: 00000000000000000000000000000000000000000000000000
250: 00000000000000000000000000000000000000000000000000
300: 00000000000000000000000000000000000000000000000000
350: 00000000000000000000000000000000000000000000000000
400: 00000000000000000000000000000000000000000000000000
450: 00000000000000000000000000000000000000000000000000
500: 00000000000000000000000000000000000000000000000000
550: 00000000000000000000000000000000000000000000000000
600: 00000000000000000000000000000000000000000000000000
650: 00000000000000000000000000000000000000000000000000
700: 00000000000000000000000000000000000000000000000000
750: 00000000000000000000000000000000000000000000000000
800: 00000000000000000000000000000000000000000000000000
850: 00000000000000000000000000000000000000000000000000
900: 00000000000000000000000000000000000000000000000000
950: 00000000000000000000000000000000000000000000000000
1000: 00000000000000000000000000000000000000000000000000
1050: 00000000000000000000000000000000000000000000000000
1100: 00000

Thus we have made a mask vertexdofs which reveals all free dofs associated to vertices. If these are labeled \(c\) (and the remainder is labeled \(f\)), then the matrix \(A\) can partitioned into

\[\begin{split}A = \left( \begin{array}{cc} A_{cc} & A_{cf} \\ A_{fc} & A_{ff} \end{array} \right).\end{split}\]

The matrix coarsepre below represents

\[\begin{split}\left( \begin{array}{cc} A_{cc}^{-1} & 0 \\ 0 & 0 \end{array} \right).\end{split}\]

coarsepre = a.mat.Inverse(vertexdofs)

This matrix can be used for coarse grid correction.

Pitfall!

Note that coarsepre is not appropriate as a preconditioner by itself as it has a large null space. You might get the wrong idea from the results of a Lanczos eigenvalue estimation:

EigenValues_Preconditioner(mat=a.mat, pre=coarsepre)

       1

But this result only gives the Laczos eigenvalue estimates on the range of the preconditioner. The preconditioned operator in this case is simply

\[\begin{split}\left( \begin{array}{cc} A_{cc}^{-1} & 0 \\ 0 & 0 \end{array} \right) \left( \begin{array}{cc} A_{cc} & A_{cf} \\ A_{fc} & A_{ff} \end{array} \right) = \left( \begin{array}{cc} I & A_{cc}^{-1} A_{cf} \\ 0 & 0 \end{array} \right),\end{split}\]

which is a projection into the \(c\)-dofs. Hence its no surprise that Lanczos estimated the eigenvalues of this operator (on its range) to be just one. But this does not suggest that coarsepre has any utility as a preconditioner by itself.

Additive two-grid preconditioner

One well-known correct way to combine the coarse grid correction with one of the previous smoothers is to combine them additively, to get an additive two-grid preconditioner as follows.

twogrid = coarsepre + blockgs

This addition of two operators (of type BaseMatrix) results in another operator, which is stored as an expression, to be evaluated only when needed. The 2-grid preconditioner has a very good condition number:

lam = EigenValues_Preconditioner(mat=a.mat, pre=twogrid)
lam

 0.993322
 0.998163
 0.999953
 1.34169
 1.80171
 1.82729
 1.95319
 1.97402
 1.99976

Combining multigrid with block smoothers

The twogrid preconditioner becomes more expensive on finer meshes due to the large matrix inversion required for the computation of coarsepre. This can be avoided by replacing coarsepre by the multigrid preconditioner we saw in 2.1.1. It is easy to combine your own block smoothers on the finest grid with the built-in multigrid on coarser levels, as the next example shows.

mesh, fes, a, f, gfu = Setup(h=0.5, p=3)
mg = Preconditioner(a, 'multigrid')  # register mg to a
a.Assemble()                         # assemble on coarsest mesh

<ngsolve.comp.BilinearForm at 0x7f6594e7a370>

Let us refine a few times to make the problem size larger.

for i in range(6):   # finer mesh updates & assembly
    mesh.Refine()
    fes.Update()
    a.Assemble()
print('NDOF = ', fes.ndof)

NDOF =  111361

Next, reset the block smoother to use the blocks on the finest grid.

blocks = VertexPatchBlocks(mesh, fes)
blockjac = a.mat.CreateBlockSmoother(blocks)
blockgs = SymmetricGS(blockjac)

Finally, add the multigrid and symmetrized block Gauss-Seidel preconditioners together to form a new preconditioner.

mgblock = mg.mat + blockgs

lam = EigenValues_Preconditioner(mat=a.mat, pre=mgblock)
lam

 0.81817
 0.846572
 0.896024
 0.932504
 1.00024
  1.1144
 1.36948
 1.43102
 1.54233
 1.65822
 1.76373
  1.8514
 1.90737
 1.97204
 1.99987

Although mgblock is similarly conditioned as twogrid, it is much more efficient.