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)

iteration 0 error = 0.05524295387946269
iteration 1 error = 0.09376127597496148
iteration 2 error = 0.10563655334471124
iteration 3 error = 0.09047524204572506
iteration 4 error = 0.10156639303380743
iteration 5 error = 0.08511608120027846
iteration 6 error = 0.0893829909061472
iteration 7 error = 0.07614379382605627
iteration 8 error = 0.05603085941249269
iteration 9 error = 0.04460004069404419
iteration 10 error = 0.034299545654838816
iteration 11 error = 0.02663601044067967
iteration 12 error = 0.022021170721158653
iteration 13 error = 0.019822026120637937
iteration 14 error = 0.014166464264223367
iteration 15 error = 0.010700608634618176
iteration 16 error = 0.00957137561391314
iteration 17 error = 0.007208975913139397
iteration 18 error = 0.004771735932237613
iteration 19 error = 0.0032429544813688685
iteration 20 error = 0.0018977551986906515
iteration 21 error = 0.0012968185931368685
iteration 22 error = 0.0010093628325403105
iteration 23 error = 0.0006896194692057556
iteration 24 error = 0.00043705829093208836
iteration 25 error = 0.0003116769183444981
iteration 26 error = 0.00024542906858424156
iteration 27 error = 0.00017976018096536533
iteration 28 error = 0.00012434796924257817
iteration 29 error = 8.954984553351105e-05
iteration 30 error = 5.887043535996482e-05
iteration 31 error = 3.99967932696822e-05
iteration 32 error = 3.1049914074865136e-05
iteration 33 error = 2.221802228281217e-05
iteration 34 error = 1.626576046078481e-05
iteration 35 error = 1.1311320610596612e-05
iteration 36 error = 6.793738789486323e-06
iteration 37 error = 4.68914955725263e-06
iteration 38 error = 3.6035551986113257e-06
iteration 39 error = 2.824257704571959e-06
iteration 40 error = 1.8262611771621658e-06
iteration 41 error = 1.2661479606819396e-06
iteration 42 error = 9.337678125441254e-07
iteration 43 error = 7.738628687367862e-07
iteration 44 error = 6.207636928561359e-07
iteration 45 error = 4.769547705304751e-07
iteration 46 error = 3.5535010590239097e-07
iteration 47 error = 2.721308593219393e-07
iteration 48 error = 2.314270091202767e-07
iteration 49 error = 1.7248577406309877e-07
iteration 50 error = 1.235367271772596e-07
iteration 51 error = 8.710509299897504e-08
iteration 52 error = 7.10670658540855e-08
iteration 53 error = 5.350063578422047e-08
iteration 54 error = 3.716799344015083e-08
iteration 55 error = 2.5379369054693467e-08
iteration 56 error = 1.7741493006066332e-08
iteration 57 error = 1.2494645047723942e-08
iteration 58 error = 8.576379183017597e-09
iteration 59 error = 5.662817970438263e-09
iteration 60 error = 4.028594813783876e-09
iteration 61 error = 2.8404614922835592e-09
iteration 62 error = 1.950999638380831e-09
iteration 63 error = 1.2973379061827987e-09
iteration 64 error = 8.575870820395727e-10
iteration 65 error = 5.991281854653724e-10
iteration 66 error = 4.5341628216596803e-10
iteration 67 error = 3.2162951752282944e-10
iteration 68 error = 2.0788772995548513e-10
iteration 69 error = 1.2894507620446703e-10
iteration 70 error = 9.095293935957835e-11
iteration 71 error = 6.790801003583008e-11
iteration 72 error = 4.7621235730775816e-11
iteration 73 error = 2.7447396486358056e-11
iteration 74 error = 1.685125911846479e-11
iteration 75 error = 1.1659311698121712e-11
iteration 76 error = 8.515477989678494e-12
iteration 77 error = 6.13723024131398e-12
iteration 78 error = 3.750437390363605e-12
iteration 79 error = 2.132229972553307e-12
iteration 80 error = 1.3650751853148497e-12
iteration 81 error = 9.827922929217147e-13
iteration 82 error = 7.157964303282535e-13
iteration 83 error = 4.4742629742889706e-13
iteration 84 error = 2.781675338709927e-13
iteration 85 error = 1.8128595205771743e-13
iteration 86 error = 1.3511328445321168e-13
iteration 87 error = 9.400002677651537e-14
iteration 88 error = 5.971248522469686e-14
iteration 89 error = 3.800343849495923e-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.07722434690357706
it# 2 , res = 0.07372671533944081
it# 3 , res = 0.07058362799187098
it# 4 , res = 0.06774390313157985
it# 5 , res = 0.06514233211838658
it# 6 , res = 0.06273356650846847
it# 7 , res = 0.060486069603128245
it# 8 , res = 0.05837651160985971
it# 9 , res = 0.05638686344051165
it# 10 , res = 0.054502783768246424
it# 11 , res = 0.05271263149481052
it# 12 , res = 0.05100681340471881
it# 13 , res = 0.04937732916040654
it# 14 , res = 0.047817440700937144
it# 15 , res = 0.04632142425784601
it# 16 , res = 0.044884379719762044
it# 17 , res = 0.04350208133966723
it# 18 , res = 0.04217085918722511
it# 19 , res = 0.04088750403385481
it# 20 , res = 0.03964919043713667
it# 21 , res = 0.03845341416237432
it# 22 , res = 0.037297941020162456
it# 23 , res = 0.036180764867922925
it# 24 , res = 0.03510007301366511
it# 25 , res = 0.034054217628338455
it# 26 , res = 0.033041692054882325
it# 27 , res = 0.03206111112095299
it# 28 , res = 0.031111194734326386
it# 29 , res = 0.03019075417639773
it# 30 , res = 0.0292986806181724
it# 31 , res = 0.028433935470679107
it# 32 , res = 0.027595542252392167
it# 33 , res = 0.026782579713482888
it# 34 , res = 0.025994176003240305
it# 35 , res = 0.025229503704902036
it# 36 , res = 0.02448777559309414
it# 37 , res = 0.023768240994406116
it# 38 , res = 0.02307018265238637
it# 39 , res = 0.02239291401529028
it# 40 , res = 0.02173577687891186
it# 41 , res = 0.021098139328366832
it# 42 , res = 0.02047939393218556
it# 43 , res = 0.01987895614991515
it# 44 , res = 0.019296262920896937
it# 45 , res = 0.0187307714072344
it# 46 , res = 0.018181957868394615
it# 47 , res = 0.01764931664855348
it# 48 , res = 0.017132359260836703
it# 49 , res = 0.016630613555136724
it# 50 , res = 0.016143622958285718
it# 51 , res = 0.015670945777114795
it# 52 , res = 0.015212154556390963
it# 53 , res = 0.014766835484837025
it# 54 , res = 0.014334587843463264
it# 55 , res = 0.013915023491290745
it# 56 , res = 0.013507766384262583
it# 57 , res = 0.013112452123743193
it# 58 , res = 0.012728727531509582
it# 59 , res = 0.012356250248569327
it# 60 , res = 0.011994688355500013
it# 61 , res = 0.011643720012311022
it# 62 , res = 0.011303033116088452
it# 63 , res = 0.010972324974906848
it# 64 , res = 0.010651301996676554
it# 65 , res = 0.010339679391759622
it# 66 , res = 0.010037180888323687
it# 67 , res = 0.009743538459521257
it# 68 , res = 0.00945849206168522
it# 69 , res = 0.009181789382816524
it# 70 , res = 0.008913185600720436
it# 71 , res = 0.008652443150207873
it# 72 , res = 0.008399331498840572
it# 73 , res = 0.008153626930745587
it# 74 , res = 0.007915112338070544
it# 75 , res = 0.007683577019686401
it# 76 , res = 0.007458816486780851
it# 77 , res = 0.0072406322750132035
it# 78 , res = 0.0070288317629288565
it# 79 , res = 0.006823227996352576
it# 80 , res = 0.006623639518504249
it# 81 , res = 0.006429890205593719
it# 82 , res = 0.006241809107672673
it# 83 , res = 0.00605923029453324
it# 84 , res = 0.005881992706456684
it# 85 , res = 0.005709940009628347
it# 86 , res = 0.005542920456044424
it# 87 , res = 0.005380786747747469
it# 88 , res = 0.00522339590523501
it# 89 , res = 0.005070609139893011
it# 90 , res = 0.004922291730317393
it# 91 , res = 0.004778312902388384
it# 92 , res = 0.0046385457129716106
it# 93 , res = 0.0045028669371256266
it# 94 , res = 0.004371156958700312
it# 95 , res = 0.0042432996642156155
it# 96 , res = 0.004119182339914166
it# 97 , res = 0.003998695571887038
it# 98 , res = 0.0038817331491754993
it# 99 , res = 0.0037681919697531664
it# 100 , res = 0.0036579719493002078
it# 101 , res = 0.003550975932682843
it# 102 , res = 0.003447109608053307
it# 103 , res = 0.003346281423490148
it# 104 , res = 0.003248402506101294
it# 105 , res = 0.003153386583515393
it# 106 , res = 0.003061149907687406
it# 107 , res = 0.002971611180949587
it# 108 , res = 0.002884691484239744
it# 109 , res = 0.0028003142074421853
it# 110 , res = 0.002718404981775048
it# 111 , res = 0.0026388916141677848
it# 112 , res = 0.002561704023563668
it# 113 , res = 0.0024867741790948654
it# 114 , res = 0.002414036040071463
it# 115 , res = 0.002343425497730812
it# 116 , res = 0.0022748803186968237
it# 117 , res = 0.0022083400900966537
it# 118 , res = 0.002143746166286021
it# 119 , res = 0.0020810416171371804
it# 120 , res = 0.0020201711778412416
it# 121 , res = 0.0019610812001815917
it# 122 , res = 0.00190371960523491
it# 123 , res = 0.00184803583745751
it# 124 , res = 0.0017939808201156134
it# 125 , res = 0.0017415069120220714
it# 126 , res = 0.0016905678655379874
it# 127 , res = 0.0016411187858055665
it# 128 , res = 0.0015931160911729673
it# 129 , res = 0.001546517474778554
it# 130 , res = 0.0015012818672589911
it# 131 , res = 0.0014573694005482359
it# 132 , res = 0.0014147413727373832
it# 133 , res = 0.0013733602139613498
it# 134 , res = 0.0013331894532842574
it# 135 , res = 0.0012941936865542177
it# 136 , res = 0.001256338545198091
it# 137 , res = 0.0012195906659291577
it# 138 , res = 0.0011839176613417895
it# 139 , res = 0.001149288091365617
it# 140 , res = 0.0011156714355551071
it# 141 , res = 0.0010830380661901024
it# 142 , res = 0.001051359222162636
it# 143 , res = 0.0010206069836285014
it# 144 , res = 0.0009907542473994843
it# 145 , res = 0.0009617747030561617
it# 146 , res = 0.0009336428097589232
it# 147 , res = 0.0009063337737374909
it# 148 , res = 0.000879823526439194
it# 149 , res = 0.0008540887033163882
it# 150 , res = 0.0008291066232334746
it# 151 , res = 0.0008048552684777187
it# 152 , res = 0.0007813132653538532
it# 153 , res = 0.0007584598653465038
it# 154 , res = 0.0007362749268333651
it# 155 , res = 0.0007147388973345909
it# 156 , res = 0.0006938327962787634
it# 157 , res = 0.0006735381982758542
it# 158 , res = 0.0006538372168778939
it# 159 , res = 0.0006347124888146653
it# 160 , res = 0.0006161471586909099
it# 161 , res = 0.0005981248641309302
it# 162 , res = 0.0005806297213580232
it# 163 , res = 0.0005636463111953613
it# 164 , res = 0.0005471596654762852
it# 165 , res = 0.0005311552538525114
it# 166 , res = 0.0005156189709881307
it# 167 , res = 0.0005005371241273079
it# 168 , res = 0.00048589642102729347
it# 169 , res = 0.0004716839582426986
it# 170 , res = 0.0004578872097534966
it# 171 , res = 0.0004444940159253884
it# 172 , res = 0.0004314925727930988
it# 173 , res = 0.0004188714216566528
it# 174 , res = 0.00040661943898294786
it# 175 , res = 0.00039472582660186477
it# 176 , res = 0.00038318010218932916
it# 177 , res = 0.0003719720900290091
it# 178 , res = 0.00036109191204390705
it# 179 , res = 0.00035052997909066003
it# 180 , res = 0.00034027698250788563
it# 181 , res = 0.0003303238859125127
it# 182 , res = 0.00032066191723534635
it# 183 , res = 0.0003112825609899833
it# 184 , res = 0.00030217755076781747
it# 185 , res = 0.0002933388619525272
it# 186 , res = 0.0002847587046475142
it# 187 , res = 0.0002764295168105995
it# 188 , res = 0.00026834395758916106
it# 189 , res = 0.00026049490085014064
it# 190 , res = 0.0002528754289000828
it# 191 , res = 0.00024547882638772034
it# 192 , res = 0.00023829857438562005
it# 193 , res = 0.00023132834464501784
it# 194 , res = 0.00022456199401807423
it# 195 , res = 0.00021799355904394656
it# 196 , res = 0.00021161725069317853
it# 197 , res = 0.00020542744926498743
it# 198 , res = 0.00019941869943495994
it# 199 , res = 0.00019358570544658065
it# 200 , res = 0.00018792332644447182
it# 201 , res = 0.00018242657194282564
it# 202 , res = 0.00017709059742762072
it# 203 , res = 0.0001719107000867175
it# 204 , res = 0.0001668823146651631
it# 205 , res = 0.00016200100944190082
it# 206 , res = 0.0001572624823231725
it# 207 , res = 0.0001526625570516961
it# 208 , res = 0.00014819717952564376
it# 209 , res = 0.0001438624142253989
it# 210 , res = 0.00013965444074527463
it# 211 , res = 0.00013556955042670317
it# 212 , res = 0.00013160414308902428
it# 213 , res = 0.00012775472385676722
it# 214 , res = 0.0001240179000796989
it# 215 , res = 0.00012039037834232827
it# 216 , res = 0.00011686896156168353
it# 217 , res = 0.00011345054616972604
it# 218 , res = 0.00011013211937688923
it# 219 , res = 0.00010691075651840176
it# 220 , res = 0.0001037836184755832
it# 221 , res = 0.0001007479491741605
it# 222 , res = 9.780107315455316e-05
it# 223 , res = 9.494039321474211e-05
it# 224 , res = 9.216338812055096e-05
it# 225 , res = 8.946761038429759e-05
it# 226 , res = 8.685068410665026e-05
it# 227 , res = 8.431030288387458e-05
it# 228 , res = 8.184422777409835e-05
it# 229 , res = 7.945028532476943e-05
it# 230 , res = 7.712636565635817e-05
it# 231 , res = 7.487042060389033e-05
it# 232 , res = 7.268046191071427e-05
it# 233 , res = 7.055455947688313e-05
it# 234 , res = 6.849083965761614e-05
it# 235 , res = 6.648748361250343e-05
it# 236 , res = 6.454272570197884e-05
it# 237 , res = 6.265485193152187e-05
it# 238 , res = 6.082219844087528e-05
it# 239 , res = 5.904315003757068e-05
it# 240 , res = 5.7316138773749814e-05
it# 241 , res = 5.5639642563891405e-05
it# 242 , res = 5.4012183843123874e-05
it# 243 , res = 5.243232826599922e-05
it# 244 , res = 5.089868344099463e-05
it# 245 , res = 4.940989770431604e-05
it# 246 , res = 4.7964658928451316e-05
it# 247 , res = 4.6561693365391046e-05
it# 248 , res = 4.519976452420281e-05
it# 249 , res = 4.387767208155239e-05
it# 250 , res = 4.259425082311131e-05
it# 251 , res = 4.134836961752538e-05
it# 252 , res = 4.013893041875812e-05
it# 253 , res = 3.8964867298544575e-05
it# 254 , res = 3.7825145507320105e-05
it# 255 , res = 3.67187605614981e-05
it# 256 , res = 3.564473735898124e-05
it# 257 , res = 3.46021293192268e-05
it# 258 , res = 3.3590017549306215e-05
it# 259 , res = 3.260751003352665e-05
it# 260 , res = 3.16537408480087e-05
it# 261 , res = 3.072786939707079e-05
it# 262 , res = 2.9829079672480437e-05
it# 263 , res = 2.895657953388457e-05
it# 264 , res = 2.8109600011462015e-05
it# 265 , res = 2.728739462742902e-05
it# 266 , res = 2.6489238738840895e-05
it# 267 , res = 2.5714428898127996e-05
it# 268 , res = 2.496228223376348e-05
it# 269 , res = 2.4232135848119416e-05
it# 270 , res = 2.352334623356665e-05
it# 271 , res = 2.2835288704664153e-05
it# 272 , res = 2.2167356848262494e-05
it# 273 , res = 2.1518961988854427e-05
it# 274 , res = 2.0889532669413016e-05
it# 275 , res = 2.0278514148241554e-05
it# 276 , res = 1.96853679098048e-05
it# 277 , res = 1.9109571190044037e-05
it# 278 , res = 1.8550616515844327e-05
it# 279 , res = 1.8008011257503154e-05
it# 280 , res = 1.7481277195029734e-05
it# 281 , res = 1.6969950095892863e-05
it# 282 , res = 1.6473579306977232e-05
it# 283 , res = 1.599172735630305e-05
it# 284 , res = 1.5523969568028234e-05
it# 285 , res = 1.5069893688064872e-05
it# 286 , res = 1.462909952101963e-05
it# 287 , res = 1.4201198576668426e-05
it# 288 , res = 1.3785813728638868e-05
it# 289 , res = 1.3382578881192376e-05
it# 290 , res = 1.2991138647245321e-05
it# 291 , res = 1.26111480341776e-05
it# 292 , res = 1.2242272141015785e-05
it# 293 , res = 1.1884185862203153e-05
it# 294 , res = 1.1536573601788474e-05
it# 295 , res = 1.119912899482601e-05
it# 296 , res = 1.0871554637587846e-05
it# 297 , res = 1.055356182550256e-05
it# 298 , res = 1.0244870298456642e-05
it# 299 , res = 9.945207993949436e-06
it# 300 , res = 9.654310807439873e-06
it# 301 , res = 9.371922359284201e-06
it# 302 , res = 9.09779376917742e-06
it# 303 , res = 8.831683436149856e-06
it# 304 , res = 8.573356826419919e-06
it# 305 , res = 8.322586266259702e-06
it# 306 , res = 8.079150741256952e-06
it# 307 , res = 7.842835701744413e-06
it# 308 , res = 7.613432873644087e-06
it# 309 , res = 7.3907400749504876e-06
it# 310 , res = 7.174561037159982e-06
it# 311 , res = 6.964705233123946e-06
it# 312 , res = 6.760987708153711e-06
it# 313 , res = 6.563228917804502e-06
it# 314 , res = 6.3712545691564565e-06
it# 315 , res = 6.18489546744773e-06
it# 316 , res = 6.003987366643113e-06
it# 317 , res = 5.828370825122277e-06
it# 318 , res = 5.657891064782736e-06
it# 319 , res = 5.4923978348955354e-06
it# 320 , res = 5.331745279592933e-06
it# 321 , res = 5.175791809168539e-06
it# 322 , res = 5.024399975407873e-06
it# 323 , res = 4.877436350567942e-06
it# 324 , res = 4.734771409516403e-06
it# 325 , res = 4.596279416013542e-06
it# 326 , res = 4.461838311104216e-06
it# 327 , res = 4.331329606666986e-06
it# 328 , res = 4.2046382797720485e-06
it# 329 , res = 4.081652672397343e-06
it# 330 , res = 3.962264392220424e-06
it# 331 , res = 3.846368217368049e-06
it# 332 , res = 3.733862004007347e-06
it# 333 , res = 3.6246465956546562e-06
it# 334 , res = 3.518625736392811e-06
it# 335 , res = 3.4157059857089384e-06
it# 336 , res = 3.315796636327627e-06
it# 337 , res = 3.2188096338976797e-06
it# 338 , res = 3.1246594997349553e-06
it# 339 , res = 3.0332632555529122e-06
it# 340 , res = 2.944540350398397e-06
it# 341 , res = 2.8584125887162624e-06
it# 342 , res = 2.7748040628272905e-06
it# 343 , res = 2.693641085107519e-06
it# 344 , res = 2.6148521232083246e-06
it# 345 , res = 2.538367737382733e-06
it# 346 , res = 2.46412051870885e-06
it# 347 , res = 2.3920450300653624e-06
it# 348 , res = 2.322077748310852e-06
it# 349 , res = 2.254157008477709e-06
it# 350 , res = 2.188222949428747e-06
it# 351 , res = 2.1242174606068192e-06
it# 352 , res = 2.0620841313610375e-06
it# 353 , res = 2.0017682010020356e-06
it# 354 , res = 1.9432165107675366e-06
it# 355 , res = 1.8863774566260217e-06
it# 356 , res = 1.831200944002719e-06
it# 357 , res = 1.7776383434778075e-06
it# 358 , res = 1.725642448316669e-06
it# 359 , res = 1.6751674321763824e-06
it# 360 , res = 1.6261688096594287e-06
it# 361 , res = 1.5786033960364125e-06
it# 362 , res = 1.5324292701882895e-06
it# 363 , res = 1.4876057366225872e-06
it# 364 , res = 1.4440932908837875e-06
it# 365 , res = 1.401853583466103e-06
it# 366 , res = 1.3608493869929485e-06
it# 367 , res = 1.321044562593136e-06
it# 368 , res = 1.2824040286224487e-06
it# 369 , res = 1.244893729800185e-06
it# 370 , res = 1.208480606725003e-06
it# 371 , res = 1.1731325667051715e-06
it# 372 , res = 1.1388184565807493e-06
it# 373 , res = 1.1055080336778896e-06
it# 374 , res = 1.0731719400031575e-06
it# 375 , res = 1.041781676703286e-06
it# 376 , res = 1.0113095780830772e-06
it# 377 , res = 9.817287879581581e-07
it# 378 , res = 9.530132353467065e-07
it# 379 , res = 9.251376123702588e-07
it# 380 , res = 8.980773508734169e-07
it# 381 , res = 8.7180860167856e-07
it# 382 , res = 8.463082129426579e-07
it# 383 , res = 8.215537102718795e-07
it# 384 , res = 7.975232764890159e-07
it# 385 , res = 7.741957323952881e-07
it# 386 , res = 7.515505187143169e-07
it# 387 , res = 7.295676772087463e-07
it# 388 , res = 7.082278334032112e-07
it# 389 , res = 6.875121797847203e-07
it# 390 , res = 6.67402458689065e-07
it# 391 , res = 6.478809467274684e-07
it# 392 , res = 6.289304386078325e-07
it# 393 , res = 6.105342325085935e-07
it# 394 , res = 5.926761153083626e-07
it# 395 , res = 5.753403477719973e-07
it# 396 , res = 5.585116511278068e-07
it# 397 , res = 5.421751935506215e-07
it# 398 , res = 5.263165771788591e-07
it# 399 , res = 5.10921825083133e-07
it# 400 , res = 4.959773693263004e-07
it# 401 , res = 4.814700385854812e-07
it# 402 , res = 4.673870470438019e-07
it# 403 , res = 4.5371598290605653e-07
it# 404 , res = 4.4044479707762335e-07
it# 405 , res = 4.275617933734257e-07
it# 406 , res = 4.1505561709995274e-07
it# 407 , res = 4.0291524644007435e-07
it# 408 , res = 3.9112998155937857e-07
it# 409 , res = 3.796894353898043e-07
it# 410 , res = 3.685835249698386e-07
it# 411 , res = 3.578024623058505e-07
it# 412 , res = 3.473367456116609e-07
it# 413 , res = 3.371771509567545e-07
it# 414 , res = 3.273147243026346e-07
it# 415 , res = 3.177407734453116e-07
it# 416 , res = 3.08446860518488e-07
it# 417 , res = 2.9942479448627774e-07
it# 418 , res = 2.9066662361765253e-07
it# 419 , res = 2.821646291264795e-07
it# 420 , res = 2.7391131799993214e-07
it# 421 , res = 2.6589941599859617e-07
it# 422 , res = 2.581218620334006e-07
it# 423 , res = 2.5057180142395547e-07
it# 424 , res = 2.432425798182805e-07
it# 425 , res = 2.3612773801521413e-07
it# 426 , res = 2.2922100510857406e-07
it# 427 , res = 2.2251629413835374e-07
it# 428 , res = 2.1600769583557618e-07
it# 429 , res = 2.0968947396800128e-07
it# 430 , res = 2.03556060042931e-07
it# 431 , res = 1.9760204826318257e-07
it# 432 , res = 1.9182219126901564e-07
it# 433 , res = 1.862113949939069e-07
it# 434 , res = 1.8076471445718104e-07
it# 435 , res = 1.7547734925373195e-07
it# 436 , res = 1.7034463935197135e-07
it# 437 , res = 1.6536206122581804e-07
it# 438 , res = 1.6052522328802906e-07
it# 439 , res = 1.5582986291599775e-07
it# 440 , res = 1.5127184181828352e-07
it# 441 , res = 1.4684714265885206e-07
it# 442 , res = 1.4255186579601687e-07
it# 443 , res = 1.3838222589884002e-07
it# 444 , res = 1.3433454783942237e-07
it# 445 , res = 1.3040526427341992e-07
it# 446 , res = 1.2659091230179675e-07
it# 447 , res = 1.2288812995313153e-07
it# 448 , res = 1.1929365394516867e-07
it# 449 , res = 1.1580431635992562e-07
it# 450 , res = 1.1241704188375741e-07
it# 451 , res = 1.0912884500864415e-07
it# 452 , res = 1.0593682792250536e-07
it# 453 , res = 1.0283817718510378e-07
it# 454 , res = 9.983016211458954e-08
it# 455 , res = 9.691013137813251e-08
it# 456 , res = 9.407551148732647e-08
it# 457 , res = 9.132380410060052e-08
it# 458 , res = 8.865258421219522e-08
it# 459 , res = 8.60594974435908e-08
it# 460 , res = 8.354225837398092e-08
it# 461 , res = 8.109864859009565e-08
it# 462 , res = 7.872651425066996e-08
it# 463 , res = 7.642376496338523e-08
it# 464 , res = 7.418837101295464e-08
it# 465 , res = 7.201836212610933e-08
it# 466 , res = 6.99118261961644e-08
it# 467 , res = 6.786690627352963e-08
it# 468 , res = 6.588180022898236e-08
it# 469 , res = 6.395475861177617e-08
it# 470 , res = 6.208408285901248e-08
it# 471 , res = 6.02681243001086e-08
it# 472 , res = 5.850528257718286e-08
it# 473 , res = 5.6794003944493336e-08
it# 474 , res = 5.5132780176713715e-08
it# 475 , res = 5.3520147212873365e-08
it# 476 , res = 5.195468361744761e-08
it# 477 , res = 5.043500987967098e-08
it# 478 , res = 4.89597866053002e-08
it# 479 , res = 4.752771363749689e-08
it# 480 , res = 4.613752872055795e-08
it# 481 , res = 4.4788006660357823e-08
it# 482 , res = 4.347795818260758e-08
it# 483 , res = 4.2206228505797764e-08
it# 484 , res = 4.097169690848634e-08
it# 485 , res = 3.977327540255158e-08
it# 486 , res = 3.860990765130689e-08
it# 487 , res = 3.74805684752609e-08
it# 488 , res = 3.6384262446282424e-08
it# 489 , res = 3.5320023274594164e-08
it# 490 , res = 3.4286913113784354e-08
it# 491 , res = 3.3284021466337866e-08
it# 492 , res = 3.2310464276205306e-08
it# 493 , res = 3.136538369494838e-08
it# 494 , res = 3.0447946681557026e-08
it# 495 , res = 2.955734466989283e-08
it# 496 , res = 2.8692792771418512e-08
it# 497 , res = 2.7853529053565008e-08
it# 498 , res = 2.7038813661485597e-08
it# 499 , res = 2.6247928762558323e-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.09424486770773638
iteration 1 error = 0.12236774743788821
iteration 2 error = 0.08652906347550848
iteration 3 error = 0.054137526987692794
iteration 4 error = 0.02679398207422516
iteration 5 error = 0.01431174586820057
iteration 6 error = 0.006816019517118301
iteration 7 error = 0.0029406308760702034
iteration 8 error = 0.0009933184655932428
iteration 9 error = 0.00037089940871912517
iteration 10 error = 0.00011720010012206778
iteration 11 error = 4.7006760956393956e-05
iteration 12 error = 1.7799308610453605e-05
iteration 13 error = 7.68645540484174e-06
iteration 14 error = 3.6669252677909787e-06
iteration 15 error = 1.7946977068048542e-06
iteration 16 error = 6.735810433048677e-07
iteration 17 error = 3.241509443945552e-07
iteration 18 error = 1.0339551054266204e-07
iteration 19 error = 5.166621063401873e-08
iteration 20 error = 1.7279913947029037e-08
iteration 21 error = 7.089316145915357e-09
iteration 22 error = 2.5541865867370014e-09
iteration 23 error = 7.586153011096176e-10
iteration 24 error = 3.207703523846336e-10
iteration 25 error = 7.796554929308681e-11
iteration 26 error = 2.836224719519651e-11
iteration 27 error = 1.0250393448983554e-11
iteration 28 error = 2.9250609172746305e-12
iteration 29 error = 1.1336539520784918e-12
iteration 30 error = 3.990003875426433e-13
iteration 31 error = 1.5422278016429751e-13
iteration 32 error = 5.990054916435788e-14

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

20.13955845213168

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, 217, 216, 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, 248, 244, 245, 55, 56, 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, 280, 24, 281, 25, 465, 58, 59}, {471, 905, 907, 908, 272, 273, 276, 277, 278, 279, 280, 281, 26, 283, 24, 25, 282, 286, 287, 59, 60, 470}, {476, 907, 909, 910, 278, 279, 25, 282, 27, 284, 26, 283, 285, 288, 289, 286, 287, 292, 293, 477, 60, 61}, {909, 911, 912, 26, 27, 28, 285, 284, 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}, {73, 74, 550, 367, 869, 551, 868, 164, 166, 167, 165, 41, 170, 171, 40, 362, 363, 871, 42, 172, 173, 368, 369, 939, 370, 371, 376, 377, 940, 366, 372, 373, 941}, {384, 385, 73, 377, 997, 552, 375, 553, 994, 376, 100, 870, 871, 378, 41, 42, 170, 172, 173, 171, 873, 176, 177, 368, 369, 43, 178, 179, 374, 939, 370, 371, 379}, {384, 385, 388, 389, 75, 184, 375, 994, 995, 100, 872, 873, 42, 43, 44, 875, 942, 176, 177, 178, 179, 564, 565, 182, 183, 374, 185, 378, 379, 380, 381, 382, 383}, {386, 387, 388, 389, 382, 383, 390, 391, 394, 75, 395, 76, 189, 380, 874, 43, 44, 875, 45, 877, 942, 943, 944, 562, 563, 182, 183, 184, 185, 188, 381, 190, 191}, {194, 195, 386, 387, 196, 197, 392, 393, 390, 391, 394, 395, 76, 396, 77, 397, 400, 401, 876, 45, 44, 877, 46, 879, 943, 945, 946, 570, 571, 188, 189, 190, 191}, {576, 577, 194, 195, 196, 197, 200, 201, 392, 393, 202, 203, 398, 399, 396, 77, 397, 400, 401, 78, 402, 403, 406, 407, 45, 878, 46, 879, 47, 881, 945, 947, 948}, {582, 79, 200, 201, 202, 203, 204, 205, 398, 399, 206, 207, 78, 402, 404, 405, 403, 406, 407, 408, 409, 412, 413, 46, 47, 880, 881, 48, 883, 947, 949, 583, 950}, {79, 204, 205, 206, 207, 144, 145, 13, 142, 404, 405, 143, 210, 211, 212, 213, 410, 411, 408, 409, 412, 413, 416, 417, 47, 48, 49, 882, 883, 884, 949, 886, 951}, {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, 440, 441, 437, 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, 443, 436, 437, 446, 447, 960, 958, 442, 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, 460, 461, 270, 271, 457, 458, 459, 462, 22, 23, 463, 86, 970, 466, 467, 937, 456, 56, 57, 58, 966}, {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, 88, 474, 475, 478, 479, 626, 627}, {907, 908, 910, 278, 279, 280, 281, 282, 25, 26, 283, 286, 287, 288, 289, 59, 60, 61, 967, 971, 972, 470, 471, 472, 89, 473, 87, 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}, {64, 65, 66, 509, 512, 513, 981, 979, 507, 917, 150, 151, 918, 919, 154, 155, 921, 311, 30, 93, 504, 508, 506, 308, 501, 310, 309, 500, 505, 314, 315, 316, 317}, {320, 321, 66, 65, 67, 322, 323, 512, 509, 513, 514, 515, 520, 521, 981, 662, 983, 920, 921, 663, 923, 984, 93, 95, 314, 508, 507, 506, 315, 316, 317, 510, 511}, {320, 321, 322, 323, 67, 66, 326, 327, 68, 516, 517, 328, 329, 518, 519, 524, 525, 520, 521, 982, 983, 664, 665, 922, 923, 514, 925, 94, 95, 515, 985, 510, 511}, {67, 68, 516, 326, 327, 328, 329, 517, 69, 332, 333, 522, 523, 334, 335, 524, 525, 526, 527, 982, 530, 531, 986, 987, 924, 925, 94, 927, 96, 666, 667, 518, 519}, {68, 69, 70, 537, 522, 523, 332, 333, 334, 335, 528, 529, 338, 339, 340, 341, 526, 527, 530, 531, 986, 532, 533, 536, 926, 927, 96, 929, 97, 988, 674, 675, 989}, {69, 70, 71, 528, 529, 338, 339, 340, 341, 534, 535, 344, 345, 346, 347, 532, 533, 536, 537, 928, 929, 97, 931, 988, 98, 990, 538, 542, 680, 681, 543, 539, 991}, {992, 993, 70, 71, 72, 538, 534, 535, 344, 345, 346, 347, 540, 541, 350, 351, 352, 353, 930, 931, 98, 933, 990, 542, 543, 544, 545, 99, 548, 549, 686, 687, 539}, {992, 547, 548, 71, 72, 74, 540, 541, 350, 351, 352, 353, 544, 545, 932, 933, 358, 359, 356, 357, 40, 938, 364, 365, 549, 998, 558, 559, 999, 99, 366, 367, 546}, {550, 551, 552, 41, 42, 939, 556, 553, 941, 554, 560, 561, 555, 557, 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, 560, 556, 557, 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, 696, 568, 569, 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}, {582, 583, 584, 585, 588, 589, 78, 79, 81, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 416, 417, 418, 419, 1003, 47, 48, 49, 949, 950, 951, 954}, {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, 88, 87, 474, 472, 473, 475, 478, 479, 480, 481, 89, 744, 745, 109, 110, 633, 626, 627, 628, 629, 630, 631, 632, 1014, 1017, 636, 637}, {1034, 1041, 107, 58, 59, 968, 969, 464, 465, 466, 467, 468, 469, 86, 87, 472, 473, 474, 475, 88, 85, 732, 733, 996, 618, 619, 620, 621, 110, 632, 624, 625, 626, 627, 1012, 622, 623, 1016, 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, 108, 109, 1015, 628, 629, 1013, 631, 630, 1014, 638, 639}, {640, 641, 1024, 644, 645, 646, 647, 648, 649, 652, 653, 639, 61, 62, 973, 974, 976, 89, 90, 91, 482, 483, 484, 485, 486, 487, 738, 739, 490, 491, 492, 493, 108, 113, 1013, 638, 1023}, {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}, {64, 65, 512, 66, 513, 514, 515, 1027, 656, 657, 658, 979, 980, 981, 662, 663, 984, 659, 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}, {71, 72, 1033, 74, 540, 541, 542, 543, 544, 545, 992, 99, 548, 549, 98, 993, 546, 547, 998, 686, 687, 558, 559, 112, 560, 561, 692, 693, 689, 688, 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, 570, 571, 700, 701, 702, 703, 704, 705, 578, 579, 580, 581, 699, 572, 573, 574, 714, 715, 76, 77, 575, 712, 713, 1084, 730, 731, 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, 1040, 1046, 700, 701, 702, 703, 576, 577, 578, 579, 580, 581, 584, 585, 586, 715, 587, 714, 77, 78, 590, 81, 591, 80, 594, 595, 598, 599, 596, 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, 1017, 87, 89, 736, 737, 740, 741, 744, 745, 746, 747, 108, 109, 110, 748, 749, 752, 753, 1015, 628, 629, 630, 631, 632, 1014, 633, 123, 125}, {637, 1041, 1048, 632, 1057, 802, 803, 1058, 87, 88, 732, 733, 120, 734, 735, 744, 745, 107, 748, 109, 110, 750, 751, 753, 626, 627, 749, 752, 630, 631, 1016, 633, 634, 635, 636, 1017, 125}, {1039, 790, 791, 1047, 1022, 794, 795, 1060, 1062, 552, 553, 554, 555, 556, 557, 694, 695, 696, 697, 73, 100, 759, 111, 112, 754, 755, 756, 757, 118, 119, 760, 758, 761, 1019, 764, 765, 126}, {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, 834, 835, 134, 848, 849, 790, 791, 1047, 792, 794, 795, 793, 796, 797, 798, 799, 1080, 1062, 759, 111, 1071, 756, 757, 758, 119, 760, 118, 761, 126, 1081}, {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, 1089, 1091, 135, 842, 843, 1054, 736, 737, 1059, 740, 741, 742, 743, 746, 747, 108, 109, 748, 749, 816, 817, 1068, 818, 819, 820, 821, 123, 125, 830, 831}, {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}, {704, 129, 130, 706, 707, 133, 710, 711, 712, 713, 1082, 844, 845, 846, 847, 850, 851, 705, 101, 102, 122, 812, 813, 814, 815, 1073, 1074, 698, 699, 1084, 1085}, {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, 1089, 1090, 1091, 132, 131, 1094, 135, 840, 841, 842, 843, 1095, 852, 853, 858, 859, 862, 863, 816, 817, 818, 819, 820, 821, 125, 123, 828, 829, 830, 831}]

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

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

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.