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()

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.0147174
 0.0567491
 0.0945183
 0.112112
 0.151298
 0.198636
 0.24877
 0.348224
 0.444189
 0.527129
 0.637083
 0.750586
 0.833432
 0.924852
  1.0317
 1.16536
 1.24096
 1.35922
 1.47459
 1.58078
  1.7088
 1.81197
 1.90832
 2.02129
 2.13615
 2.19683
 2.29152
 2.37404
 2.41815
 2.46677
 2.50584
 2.53477
 2.80933

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

max(lams)/min(lams)

190.88563050697397

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)

1579.1517245003301

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)

it =  0  err =  0.0552429538794627
it =  1  err =  0.09376127597496152
it =  2  err =  0.10563655334471139
it =  3  err =  0.09047524204572518
it =  4  err =  0.10156639303380757
it =  5  err =  0.08511608120027855
it =  6  err =  0.0893829909061473
it =  7  err =  0.07614379382605639
it =  8  err =  0.056030859412492806
it =  9  err =  0.044600040694044295
it =  10  err =  0.03429954565483888
it =  11  err =  0.026636010440679703
it =  12  err =  0.02202117072115868
it =  13  err =  0.019822026120637968
it =  14  err =  0.01416646426422339
it =  15  err =  0.01070060863461819
it =  16  err =  0.009571375613913156
it =  17  err =  0.007208975913139408
it =  18  err =  0.0047717359322376205
it =  19  err =  0.003242954481368875
it =  20  err =  0.0018977551986906563
it =  21  err =  0.001296818593136872
it =  22  err =  0.0010093628325403131
it =  23  err =  0.0006896194692057572
it =  24  err =  0.00043705829093208933
it =  25  err =  0.0003116769183444989
it =  26  err =  0.00024542906858424215
it =  27  err =  0.00017976018096536563
it =  28  err =  0.00012434796924257836
it =  29  err =  8.95498455335112e-05
it =  30  err =  5.8870435359964907e-05
it =  31  err =  3.999679326968226e-05
it =  32  err =  3.104991407486517e-05
it =  33  err =  2.2218022282812195e-05
it =  34  err =  1.626576046078482e-05
it =  35  err =  1.1311320610596612e-05
it =  36  err =  6.793738789486321e-06
it =  37  err =  4.689149557252623e-06
it =  38  err =  3.6035551986113253e-06
it =  39  err =  2.8242577045720175e-06
it =  40  err =  1.8262611771623447e-06
it =  41  err =  1.266147960682338e-06
it =  42  err =  9.337678125450523e-07
it =  43  err =  7.738628687396238e-07
it =  44  err =  6.20763692868977e-07
it =  45  err =  4.769547705750374e-07
it =  46  err =  3.553501060242699e-07
it =  47  err =  2.721308595988448e-07
it =  48  err =  2.3142701012470213e-07
it =  49  err =  1.724857792460835e-07
it =  50  err =  1.235367428705939e-07
it =  51  err =  8.71051252292543e-08
it =  52  err =  7.106716139603171e-08
it =  53  err =  5.350104755530211e-08
it =  54  err =  3.716949778368705e-08
it =  55  err =  2.538326342638807e-08
it =  56  err =  1.775005919187837e-08
it =  57  err =  1.2518978809226337e-08
it =  58  err =  8.631363162732257e-09
it =  59  err =  5.782085263588922e-09
it =  60  err =  4.152921569492078e-09
it =  61  err =  2.871994988919029e-09
it =  62  err =  1.9516058000901353e-09
it =  63  err =  1.293673967020714e-09
it =  64  err =  8.568986693801478e-10
it =  65  err =  5.990020557525543e-10
it =  66  err =  4.5339036780984956e-10
it =  67  err =  3.2162605713968247e-10
it =  68  err =  2.07886516133772e-10
it =  69  err =  1.2894495975029143e-10
it =  70  err =  9.095291166339456e-11
it =  71  err =  6.790800669954515e-11
it =  72  err =  4.762123539324343e-11
it =  73  err =  2.744739646934626e-11
it =  74  err =  1.6851259114909124e-11
it =  75  err =  1.1659311697226216e-11
it =  76  err =  8.515477989532025e-12
it =  77  err =  6.137230241301649e-12
it =  78  err =  3.750437390366042e-12
it =  79  err =  2.1322299725577025e-12
it =  80  err =  1.3650751853220057e-12
it =  81  err =  9.827922929329267e-13
it =  82  err =  7.157964303441344e-13
it =  83  err =  4.4742629744929425e-13
it =  84  err =  2.781675338996337e-13
it =  85  err =  1.8128595210009506e-13
it =  86  err =  1.3511328451708182e-13
it =  87  err =  9.400002686212036e-14
it =  88  err =  5.971248533660745e-14
it =  89  err =  3.800343865160828e-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.08188746699328223
it# 1 , res = 0.07722434690357707
it# 2 , res = 0.07372671533944081
it# 3 , res = 0.070583627991871
it# 4 , res = 0.06774390313157987
it# 5 , res = 0.06514233211838658
it# 6 , res = 0.0627335665084685
it# 7 , res = 0.06048606960312826
it# 8 , res = 0.0583765116098597
it# 9 , res = 0.05638686344051164
it# 10 , res = 0.05450278376824646
it# 11 , res = 0.052712631494810536
it# 12 , res = 0.05100681340471882
it# 13 , res = 0.04937732916040655
it# 14 , res = 0.04781744070093716
it# 15 , res = 0.04632142425784601
it# 16 , res = 0.04488437971976207
it# 17 , res = 0.04350208133966725
it# 18 , res = 0.04217085918722511
it# 19 , res = 0.040887504033854834
it# 20 , res = 0.03964919043713668
it# 21 , res = 0.03845341416237436
it# 22 , res = 0.037297941020162484
it# 23 , res = 0.03618076486792295
it# 24 , res = 0.035100073013665106
it# 25 , res = 0.034054217628338476
it# 26 , res = 0.03304169205488234
it# 27 , res = 0.03206111112095301
it# 28 , res = 0.03111119473432636
it# 29 , res = 0.030190754176397806
it# 30 , res = 0.02929868061817242
it# 31 , res = 0.028433935470679173
it# 32 , res = 0.027595542252392156
it# 33 , res = 0.026782579713482975
it# 34 , res = 0.025994176003240298
it# 35 , res = 0.025229503704902026
it# 36 , res = 0.024487775593094237
it# 37 , res = 0.023768240994406047
it# 38 , res = 0.023070182652386414
it# 39 , res = 0.022392914015290238
it# 40 , res = 0.021735776878911962
it# 41 , res = 0.02109813932836695
it# 42 , res = 0.020479393932185508
it# 43 , res = 0.0198789561499152
it# 44 , res = 0.01929626292089699
it# 45 , res = 0.018730771407234385
it# 46 , res = 0.018181957868394546
it# 47 , res = 0.017649316648553456
it# 48 , res = 0.017132359260836703
it# 49 , res = 0.016630613555136786
it# 50 , res = 0.016143622958285746
it# 51 , res = 0.015670945777114878
it# 52 , res = 0.015212154556391027
it# 53 , res = 0.014766835484837103
it# 54 , res = 0.014334587843463373
it# 55 , res = 0.013915023491290732
it# 56 , res = 0.013507766384262652
it# 57 , res = 0.013112452123743283
it# 58 , res = 0.01272872753150956
it# 59 , res = 0.012356250248569333
it# 60 , res = 0.011994688355500042
it# 61 , res = 0.011643720012311
it# 62 , res = 0.011303033116088448
it# 63 , res = 0.010972324974906866
it# 64 , res = 0.010651301996676518
it# 65 , res = 0.01033967939175966
it# 66 , res = 0.010037180888323654
it# 67 , res = 0.009743538459521264
it# 68 , res = 0.009458492061685116
it# 69 , res = 0.009181789382816496
it# 70 , res = 0.008913185600720451
it# 71 , res = 0.00865244315020789
it# 72 , res = 0.008399331498840557
it# 73 , res = 0.008153626930745528
it# 74 , res = 0.007915112338070542
it# 75 , res = 0.007683577019686464
it# 76 , res = 0.007458816486780842
it# 77 , res = 0.0072406322750133
it# 78 , res = 0.007028831762928834
it# 79 , res = 0.00682322799635261
it# 80 , res = 0.006623639518504259
it# 81 , res = 0.0064298902055936625
it# 82 , res = 0.00624180910767266
it# 83 , res = 0.006059230294533262
it# 84 , res = 0.005881992706456748
it# 85 , res = 0.00570994000962835
it# 86 , res = 0.00554292045604443
it# 87 , res = 0.005380786747747555
it# 88 , res = 0.005223395905235011
it# 89 , res = 0.005070609139892996
it# 90 , res = 0.00492229173031736
it# 91 , res = 0.004778312902388466
it# 92 , res = 0.004638545712971615
it# 93 , res = 0.004502866937125571
it# 94 , res = 0.004371156958700316
it# 95 , res = 0.004243299664215574
it# 96 , res = 0.004119182339914103
it# 97 , res = 0.0039986955718870655
it# 98 , res = 0.0038817331491755903
it# 99 , res = 0.0037681919697531863
it# 100 , res = 0.0036579719493001154
it# 101 , res = 0.0035509759326828212
it# 102 , res = 0.003447109608053369
it# 103 , res = 0.003346281423490175
it# 104 , res = 0.0032484025061012586
it# 105 , res = 0.0031533865835154558
it# 106 , res = 0.0030611499076874105
it# 107 , res = 0.0029716111809495395
it# 108 , res = 0.0028846914842398648
it# 109 , res = 0.0028003142074421346
it# 110 , res = 0.0027184049817751458
it# 111 , res = 0.00263889161416779
it# 112 , res = 0.0025617040235636898
it# 113 , res = 0.0024867741790949483
it# 114 , res = 0.002414036040071473
it# 115 , res = 0.002343425497730811
it# 116 , res = 0.0022748803186968852
it# 117 , res = 0.002208340090096674
it# 118 , res = 0.0021437461662861066
it# 119 , res = 0.002081041617137274
it# 120 , res = 0.0020201711778412837
it# 121 , res = 0.0019610812001815857
it# 122 , res = 0.0019037196052349337
it# 123 , res = 0.0018480358374574006
it# 124 , res = 0.0017939808201157185
it# 125 , res = 0.0017415069120220806
it# 126 , res = 0.0016905678655379076
it# 127 , res = 0.0016411187858054966
it# 128 , res = 0.0015931160911728955
it# 129 , res = 0.0015465174747786331
it# 130 , res = 0.0015012818672589137
it# 131 , res = 0.0014573694005482858
it# 132 , res = 0.0014147413727373802
it# 133 , res = 0.0013733602139612208
it# 134 , res = 0.0013331894532842084
it# 135 , res = 0.0012941936865542164
it# 136 , res = 0.0012563385451980532
it# 137 , res = 0.0012195906659291861
it# 138 , res = 0.0011839176613418485
it# 139 , res = 0.0011492880913655961
it# 140 , res = 0.0011156714355551331
it# 141 , res = 0.0010830380661901314
it# 142 , res = 0.001051359222162593
it# 143 , res = 0.0010206069836284299
it# 144 , res = 0.000990754247399467
it# 145 , res = 0.0009617747030562004
it# 146 , res = 0.0009336428097588429
it# 147 , res = 0.0009063337737374186
it# 148 , res = 0.0008798235264392089
it# 149 , res = 0.0008540887033164393
it# 150 , res = 0.0008291066232334186
it# 151 , res = 0.0008048552684776956
it# 152 , res = 0.0007813132653537901
it# 153 , res = 0.0007584598653464741
it# 154 , res = 0.0007362749268335144
it# 155 , res = 0.0007147388973346553
it# 156 , res = 0.0006938327962787361
it# 157 , res = 0.0006735381982759462
it# 158 , res = 0.0006538372168779915
it# 159 , res = 0.0006347124888146239
it# 160 , res = 0.0006161471586908728
it# 161 , res = 0.0005981248641309025
it# 162 , res = 0.0005806297213580052
it# 163 , res = 0.0005636463111954288
it# 164 , res = 0.0005471596654762811
it# 165 , res = 0.0005311552538525713
it# 166 , res = 0.0005156189709881925
it# 167 , res = 0.0005005371241274187
it# 168 , res = 0.0004858964210273502
it# 169 , res = 0.0004716839582427327
it# 170 , res = 0.0004578872097535058
it# 171 , res = 0.00044449401592534777
it# 172 , res = 0.0004314925727930603
it# 173 , res = 0.00041887142165663877
it# 174 , res = 0.0004066194389830285
it# 175 , res = 0.00039472582660195183
it# 176 , res = 0.0003831801021893964
it# 177 , res = 0.0003719720900289774
it# 178 , res = 0.0003610919120438523
it# 179 , res = 0.0003505299790906294
it# 180 , res = 0.00034027698250789864
it# 181 , res = 0.0003303238859125612
it# 182 , res = 0.0003206619172353159
it# 183 , res = 0.0003112825609899892
it# 184 , res = 0.0003021775507679126
it# 185 , res = 0.00029333886195255515
it# 186 , res = 0.0002847587046474964
it# 187 , res = 0.0002764295168106401
it# 188 , res = 0.0002683439575891355
it# 189 , res = 0.00026049490085015283
it# 190 , res = 0.0002528754289001816
it# 191 , res = 0.00024547882638775807
it# 192 , res = 0.00023829857438565385
it# 193 , res = 0.0002313283446449708
it# 194 , res = 0.0002245619940179679
it# 195 , res = 0.00021799355904395507
it# 196 , res = 0.00021161725069323184
it# 197 , res = 0.0002054274492649516
it# 198 , res = 0.00019941869943498848
it# 199 , res = 0.00019358570544661507
it# 200 , res = 0.00018792332644448632
it# 201 , res = 0.00018242657194280574
it# 202 , res = 0.00017709059742761224
it# 203 , res = 0.00017191070008671232
it# 204 , res = 0.00016688231466526791
it# 205 , res = 0.0001620010094418838
it# 206 , res = 0.0001572624823231584
it# 207 , res = 0.00015266255705173652
it# 208 , res = 0.0001481971795256785
it# 209 , res = 0.00014386241422535776
it# 210 , res = 0.000139654440745373
it# 211 , res = 0.00013556955042671692
it# 212 , res = 0.00013160414308897544
it# 213 , res = 0.00012775472385679446
it# 214 , res = 0.00012401790007968077
it# 215 , res = 0.00012039037834241112
it# 216 , res = 0.00011686896156184269
it# 217 , res = 0.00011345054616969878
it# 218 , res = 0.00011013211937687863
it# 219 , res = 0.00010691075651853734
it# 220 , res = 0.00010378361847564718
it# 221 , res = 0.00010074794917417807
it# 222 , res = 9.780107315451159e-05
it# 223 , res = 9.494039321463996e-05
it# 224 , res = 9.216338812063321e-05
it# 225 , res = 8.946761038422856e-05
it# 226 , res = 8.685068410658038e-05
it# 227 , res = 8.43103028838818e-05
it# 228 , res = 8.18442277741219e-05
it# 229 , res = 7.945028532470486e-05
it# 230 , res = 7.712636565632926e-05
it# 231 , res = 7.48704206038908e-05
it# 232 , res = 7.268046191072803e-05
it# 233 , res = 7.055455947688054e-05
it# 234 , res = 6.84908396576898e-05
it# 235 , res = 6.648748361246788e-05
it# 236 , res = 6.454272570200841e-05
it# 237 , res = 6.265485193148067e-05
it# 238 , res = 6.0822198440911274e-05
it# 239 , res = 5.904315003758845e-05
it# 240 , res = 5.731613877374046e-05
it# 241 , res = 5.563964256388931e-05
it# 242 , res = 5.401218384315243e-05
it# 243 , res = 5.243232826592873e-05
it# 244 , res = 5.089868344091972e-05
it# 245 , res = 4.940989770425869e-05
it# 246 , res = 4.796465892848132e-05
it# 247 , res = 4.656169336529855e-05
it# 248 , res = 4.5199764524172606e-05
it# 249 , res = 4.3877672081498534e-05
it# 250 , res = 4.259425082312461e-05
it# 251 , res = 4.134836961760845e-05
it# 252 , res = 4.013893041870703e-05
it# 253 , res = 3.896486729856328e-05
it# 254 , res = 3.78251455072447e-05
it# 255 , res = 3.671876056153479e-05
it# 256 , res = 3.564473735900218e-05
it# 257 , res = 3.460212931928312e-05
it# 258 , res = 3.359001754922468e-05
it# 259 , res = 3.260751003347081e-05
it# 260 , res = 3.165374084805117e-05
it# 261 , res = 3.0727869397128364e-05
it# 262 , res = 2.9829079672429483e-05
it# 263 , res = 2.8956579533915393e-05
it# 264 , res = 2.8109600011486396e-05
it# 265 , res = 2.7287394627475064e-05
it# 266 , res = 2.6489238738755307e-05
it# 267 , res = 2.5714428898179984e-05
it# 268 , res = 2.496228223383795e-05
it# 269 , res = 2.4232135848204618e-05
it# 270 , res = 2.3523346233538318e-05
it# 271 , res = 2.2835288704665166e-05
it# 272 , res = 2.21673568483086e-05
it# 273 , res = 2.1518961988821297e-05
it# 274 , res = 2.0889532669423523e-05
it# 275 , res = 2.02785141482473e-05
it# 276 , res = 1.9685367909793856e-05
it# 277 , res = 1.9109571190099023e-05
it# 278 , res = 1.8550616515816507e-05
it# 279 , res = 1.8008011257515832e-05
it# 280 , res = 1.7481277194949456e-05
it# 281 , res = 1.69699500958973e-05
it# 282 , res = 1.6473579307005876e-05
it# 283 , res = 1.599172735632133e-05
it# 284 , res = 1.5523969567979177e-05
it# 285 , res = 1.5069893688067786e-05
it# 286 , res = 1.4629099520931172e-05
it# 287 , res = 1.4201198576666792e-05
it# 288 , res = 1.3785813728605795e-05
it# 289 , res = 1.3382578881169277e-05
it# 290 , res = 1.2991138647189104e-05
it# 291 , res = 1.2611148034202663e-05
it# 292 , res = 1.2242272140982288e-05
it# 293 , res = 1.1884185862230419e-05
it# 294 , res = 1.1536573601775814e-05
it# 295 , res = 1.1199128994795986e-05
it# 296 , res = 1.0871554637607727e-05
it# 297 , res = 1.0553561825465693e-05
it# 298 , res = 1.0244870298398076e-05
it# 299 , res = 9.94520799393071e-06
it# 300 , res = 9.654310807421864e-06
it# 301 , res = 9.371922359310339e-06
it# 302 , res = 9.09779376915432e-06
it# 303 , res = 8.831683436157511e-06
it# 304 , res = 8.573356826435306e-06
it# 305 , res = 8.322586266268553e-06
it# 306 , res = 8.079150741194668e-06
it# 307 , res = 7.842835701749244e-06
it# 308 , res = 7.61343287365298e-06
it# 309 , res = 7.3907400748551066e-06
it# 310 , res = 7.1745610371265864e-06
it# 311 , res = 6.964705233078211e-06
it# 312 , res = 6.7609877081550325e-06
it# 313 , res = 6.563228917892473e-06
it# 314 , res = 6.37125456916473e-06
it# 315 , res = 6.1848954674263486e-06
it# 316 , res = 6.0039873666731936e-06
it# 317 , res = 5.828370825104019e-06
it# 318 , res = 5.6578910647768515e-06
it# 319 , res = 5.492397834966206e-06
it# 320 , res = 5.331745279660969e-06
it# 321 , res = 5.175791809251664e-06
it# 322 , res = 5.024399975416623e-06
it# 323 , res = 4.877436350607186e-06
it# 324 , res = 4.734771409597086e-06
it# 325 , res = 4.596279415924091e-06
it# 326 , res = 4.461838311122528e-06
it# 327 , res = 4.331329606680353e-06
it# 328 , res = 4.20463827982624e-06
it# 329 , res = 4.081652672491875e-06
it# 330 , res = 3.962264392153359e-06
it# 331 , res = 3.846368217376085e-06
it# 332 , res = 3.7338620041057196e-06
it# 333 , res = 3.624646595723967e-06
it# 334 , res = 3.5186257364782257e-06
it# 335 , res = 3.41570598576586e-06
it# 336 , res = 3.3157966363133553e-06
it# 337 , res = 3.2188096337976836e-06
it# 338 , res = 3.124659499753784e-06
it# 339 , res = 3.0332632556795047e-06
it# 340 , res = 2.9445403502950694e-06
it# 341 , res = 2.8584125886796177e-06
it# 342 , res = 2.7748040627766155e-06
it# 343 , res = 2.6936410850708316e-06
it# 344 , res = 2.614852123292599e-06
it# 345 , res = 2.5383677373785364e-06
it# 346 , res = 2.4641205186285242e-06
it# 347 , res = 2.3920450299783966e-06
it# 348 , res = 2.322077748249715e-06
it# 349 , res = 2.254157008560629e-06
it# 350 , res = 2.188222949436737e-06
it# 351 , res = 2.1242174605447647e-06
it# 352 , res = 2.062084131273931e-06
it# 353 , res = 2.0017682010236502e-06
it# 354 , res = 1.943216510802902e-06
it# 355 , res = 1.8863774566698614e-06
it# 356 , res = 1.8312009440389662e-06
it# 357 , res = 1.7776383434551769e-06
it# 358 , res = 1.725642448325568e-06
it# 359 , res = 1.6751674322659445e-06
it# 360 , res = 1.626168809640083e-06
it# 361 , res = 1.5786033960095647e-06
it# 362 , res = 1.5324292700730822e-06
it# 363 , res = 1.4876057366979776e-06
it# 364 , res = 1.4440932908895086e-06
it# 365 , res = 1.4018535835147626e-06
it# 366 , res = 1.3608493869294268e-06
it# 367 , res = 1.3210445624918917e-06
it# 368 , res = 1.2824040286847239e-06
it# 369 , res = 1.2448937297178072e-06
it# 370 , res = 1.208480606665911e-06
it# 371 , res = 1.1731325667981551e-06
it# 372 , res = 1.1388184566293222e-06
it# 373 , res = 1.1055080336287691e-06
it# 374 , res = 1.0731719400144773e-06
it# 375 , res = 1.0417816766584586e-06
it# 376 , res = 1.0113095780255569e-06
it# 377 , res = 9.81728787903465e-07
it# 378 , res = 9.530132353376798e-07
it# 379 , res = 9.251376123838071e-07
it# 380 , res = 8.98077350882184e-07
it# 381 , res = 8.718086016851057e-07
it# 382 , res = 8.463082129412395e-07
it# 383 , res = 8.215537103492866e-07
it# 384 , res = 7.975232764237535e-07
it# 385 , res = 7.741957324250015e-07
it# 386 , res = 7.51550518735848e-07
it# 387 , res = 7.295676771787723e-07
it# 388 , res = 7.082278334392152e-07
it# 389 , res = 6.87512179801473e-07
it# 390 , res = 6.674024587516438e-07
it# 391 , res = 6.478809466088582e-07
it# 392 , res = 6.289304385830849e-07
it# 393 , res = 6.10534232581628e-07
it# 394 , res = 5.92676115306649e-07
it# 395 , res = 5.753403477421899e-07
it# 396 , res = 5.585116510374611e-07
it# 397 , res = 5.421751935806867e-07
it# 398 , res = 5.26316577229265e-07
it# 399 , res = 5.109218250891858e-07
it# 400 , res = 4.959773692980557e-07
it# 401 , res = 4.814700384602469e-07
it# 402 , res = 4.673870470608611e-07
it# 403 , res = 4.537159828261768e-07
it# 404 , res = 4.4044479712211523e-07
it# 405 , res = 4.2756179328143293e-07
it# 406 , res = 4.150556171649392e-07
it# 407 , res = 4.029152466233702e-07
it# 408 , res = 3.911299815418348e-07
it# 409 , res = 3.796894353718273e-07
it# 410 , res = 3.685835249918889e-07
it# 411 , res = 3.5780246232981713e-07
it# 412 , res = 3.473367456095796e-07
it# 413 , res = 3.371771509527446e-07
it# 414 , res = 3.2731472433563826e-07
it# 415 , res = 3.177407734559517e-07
it# 416 , res = 3.0844686049922066e-07
it# 417 , res = 2.9942479442570043e-07
it# 418 , res = 2.90666623653193e-07
it# 419 , res = 2.821646291456673e-07
it# 420 , res = 2.739113179710938e-07
it# 421 , res = 2.6589941592524137e-07
it# 422 , res = 2.5812186197899556e-07
it# 423 , res = 2.50571801363738e-07
it# 424 , res = 2.432425798993089e-07
it# 425 , res = 2.3612773801313813e-07
it# 426 , res = 2.2922100511835983e-07
it# 427 , res = 2.2251629417768254e-07
it# 428 , res = 2.1600769579936372e-07
it# 429 , res = 2.096894739650195e-07
it# 430 , res = 2.0355606000079677e-07
it# 431 , res = 1.9760204830017322e-07
it# 432 , res = 1.918221911749318e-07
it# 433 , res = 1.8621139499558315e-07
it# 434 , res = 1.807647145316534e-07
it# 435 , res = 1.7547734924351117e-07
it# 436 , res = 1.7034463935493467e-07
it# 437 , res = 1.6536206126113185e-07
it# 438 , res = 1.605252233444167e-07
it# 439 , res = 1.558298629240251e-07
it# 440 , res = 1.5127184174281003e-07
it# 441 , res = 1.4684714260665445e-07
it# 442 , res = 1.4255186586665661e-07
it# 443 , res = 1.3838222597296236e-07
it# 444 , res = 1.34334547811458e-07
it# 445 , res = 1.3040526427248405e-07
it# 446 , res = 1.265909122401279e-07
it# 447 , res = 1.2288812992216935e-07
it# 448 , res = 1.1929365395623852e-07
it# 449 , res = 1.1580431637459208e-07
it# 450 , res = 1.124170419153674e-07
it# 451 , res = 1.0912884493480585e-07
it# 452 , res = 1.0593682795225098e-07
it# 453 , res = 1.0283817727863086e-07
it# 454 , res = 9.983016207390087e-08
it# 455 , res = 9.691013127841093e-08
it# 456 , res = 9.407551142898596e-08
it# 457 , res = 9.132380413686277e-08
it# 458 , res = 8.86525842307432e-08
it# 459 , res = 8.60594974365478e-08
it# 460 , res = 8.354225833623986e-08
it# 461 , res = 8.1098648503562e-08
it# 462 , res = 7.872651426105001e-08
it# 463 , res = 7.6423764931776e-08
it# 464 , res = 7.418837096155935e-08
it# 465 , res = 7.201836225429941e-08
it# 466 , res = 6.991182619042975e-08
it# 467 , res = 6.78669063101455e-08
it# 468 , res = 6.588180021798733e-08
it# 469 , res = 6.395475856220575e-08
it# 470 , res = 6.208408288947246e-08
it# 471 , res = 6.026812434121364e-08
it# 472 , res = 5.850528261021156e-08
it# 473 , res = 5.6794003973223614e-08
it# 474 , res = 5.513278016553754e-08
it# 475 , res = 5.352014715025261e-08
it# 476 , res = 5.195468374719275e-08
it# 477 , res = 5.043500994915756e-08
it# 478 , res = 4.8959786614117284e-08
it# 479 , res = 4.7527713590018824e-08
it# 480 , res = 4.613752869026745e-08
it# 481 , res = 4.4788006600563246e-08
it# 482 , res = 4.347795817643553e-08
it# 483 , res = 4.220622855955817e-08
it# 484 , res = 4.097169695151601e-08
it# 485 , res = 3.9773275393897505e-08
it# 486 , res = 3.8609907693730896e-08
it# 487 , res = 3.7480568529599996e-08
it# 488 , res = 3.6384262375878236e-08
it# 489 , res = 3.532002323881674e-08
it# 490 , res = 3.428691304295612e-08
it# 491 , res = 3.3284021481142804e-08
it# 492 , res = 3.2310464258721004e-08
it# 493 , res = 3.136538369415709e-08
it# 494 , res = 3.044794659852228e-08
it# 495 , res = 2.955734464106953e-08
it# 496 , res = 2.8692792816992623e-08
it# 497 , res = 2.785352899073144e-08
it# 498 , res = 2.7038813752531944e-08
it# 499 , res = 2.6247928821865992e-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)

it =  0  err =  0.09424486770773638
it =  1  err =  0.12236774743788827
it =  2  err =  0.08652906347550857
it =  3  err =  0.054137526987692905
it =  4  err =  0.026793982074225228
it =  5  err =  0.014311745868200618
it =  6  err =  0.006816019517118323
it =  7  err =  0.0029406308760702143
it =  8  err =  0.0009933184655932467
it =  9  err =  0.0003708994087191265
it =  10  err =  0.00011720010012206816
it =  11  err =  4.70067609563941e-05
it =  12  err =  1.779930861045365e-05
it =  13  err =  7.686455404841756e-06
it =  14  err =  3.6669252677909825e-06
it =  15  err =  1.7946977068048559e-06
it =  16  err =  6.735810433048691e-07
it =  17  err =  3.2415094439455677e-07
it =  18  err =  1.0339551054266279e-07
it =  19  err =  5.166621063401919e-08
it =  20  err =  1.7279913947029172e-08
it =  21  err =  7.0893161459153975e-09
it =  22  err =  2.554186586737014e-09
it =  23  err =  7.586153011096213e-10
it =  24  err =  3.2077035238463533e-10
it =  25  err =  7.796554929308715e-11
it =  26  err =  2.8362247195196405e-11
it =  27  err =  1.0250393448983244e-11
it =  28  err =  2.925060917273775e-12
it =  29  err =  1.1336539520765186e-12
it =  30  err =  3.99000387538362e-13
it =  31  err =  1.5422278015511608e-13
it =  32  err =  5.990054914571175e-14

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

20.139558452131674

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)

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

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)

34.629152916562035

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.976713457276087

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: 11111111111111111111111111111111111100000000000000
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: 0000000000000000000000000000000000000000000000

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.992983
 0.997563
 0.99996
 1.34273
 1.80296
 1.83822
 1.95805
 1.98195
 1.99978

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

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.811825
 0.845488
 0.895515
 0.932095
 0.997157
 1.11324
  1.3694
 1.42823
 1.54131
 1.65742
 1.76352
 1.84966
  1.9073
 1.97199
 1.99987

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