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

[1]:

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

[2]:

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)
f = LinearForm(fes)
f += v*dx
gfu = GridFunction(fes)
return mesh, fes, a, f, gfu

[3]:

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

[3]:

<ngsolve.comp.LinearForm at 0x7f345a3cbdf0>


## 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:

[4]:

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.

[5]:

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

[5]:

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


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

[6]:

max(lams)/min(lams)

[6]:

194.19561167462462


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

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

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


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

[7]:

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


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

IdentityMatrix(fes.ndof, complex=False).

[8]:

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

[8]:

1678.5881226918377


Clearly the point Jacobi preconditioner has reduced the condition number.

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

[9]:

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

iteration 0 error = 0.05505568711828939
iteration 1 error = 0.09357794070949911
iteration 2 error = 0.11255359311920529
iteration 3 error = 0.0883578866157642
iteration 4 error = 0.09465752852240354
iteration 5 error = 0.08675791614352733
iteration 6 error = 0.09174564400414802
iteration 7 error = 0.07950174974020313
iteration 8 error = 0.05558110481441834
iteration 9 error = 0.044799001641920515
iteration 10 error = 0.03885061649304077
iteration 11 error = 0.030408675499253898
iteration 12 error = 0.02288741507787512
iteration 13 error = 0.02051092486797348
iteration 14 error = 0.015405453021112554
iteration 15 error = 0.011893431227990001
iteration 16 error = 0.010292148243983716
iteration 17 error = 0.007888479867611086
iteration 18 error = 0.005481703824765564
iteration 19 error = 0.003612015168020746
iteration 20 error = 0.0023840872403187567
iteration 21 error = 0.0015057839927216509
iteration 22 error = 0.0011613533060643182
iteration 23 error = 0.000869420245521821
iteration 24 error = 0.0005329860738693037
iteration 25 error = 0.00037108669875231834
iteration 26 error = 0.0002859501917016282
iteration 27 error = 0.00022783224021485763
iteration 28 error = 0.000166594323733906
iteration 29 error = 0.00011944094829848384
iteration 30 error = 8.147974841738306e-05
iteration 31 error = 5.5020653070099096e-05
iteration 32 error = 4.190089073735408e-05
iteration 33 error = 3.369706727003576e-05
iteration 34 error = 2.7766665129423223e-05
iteration 35 error = 1.9736038657175716e-05
iteration 36 error = 1.363994533705487e-05
iteration 37 error = 1.0506189808377937e-05
iteration 38 error = 8.983446326952995e-06
iteration 39 error = 7.660521208331822e-06
iteration 40 error = 5.801953313066017e-06
iteration 41 error = 4.051197313069691e-06
iteration 42 error = 2.9961283895352596e-06
iteration 43 error = 2.3762446752543294e-06
iteration 44 error = 1.8104198799132735e-06
iteration 45 error = 1.2473066060491126e-06
iteration 46 error = 8.123842440232735e-07
iteration 47 error = 5.602213594389701e-07
iteration 48 error = 4.297300852190805e-07
iteration 49 error = 3.046923636128493e-07
iteration 50 error = 1.9532592484907017e-07
iteration 51 error = 1.3346842701507282e-07
iteration 52 error = 9.387214810153153e-08
iteration 53 error = 7.042267578742155e-08
iteration 54 error = 5.0876599295554645e-08
iteration 55 error = 3.6287644199056415e-08
iteration 56 error = 2.4168652132655575e-08
iteration 57 error = 1.6475086046502934e-08
iteration 58 error = 1.0716935189887449e-08
iteration 59 error = 6.890521019254731e-09
iteration 60 error = 4.550452935048494e-09
iteration 61 error = 3.135967176609285e-09
iteration 62 error = 2.033380757491356e-09
iteration 63 error = 1.331589593370284e-09
iteration 64 error = 9.48052142950352e-10
iteration 65 error = 7.191889610955045e-10
iteration 66 error = 5.584212662146425e-10
iteration 67 error = 3.963364255887182e-10
iteration 68 error = 2.642865235172716e-10
iteration 69 error = 1.7946264132246952e-10
iteration 70 error = 1.1724755382242814e-10
iteration 71 error = 7.742596210186727e-11
iteration 72 error = 5.2131525221998164e-11
iteration 73 error = 3.62410036123939e-11
iteration 74 error = 2.3341859703488752e-11
iteration 75 error = 1.5421560226282028e-11
iteration 76 error = 1.0252664873785537e-11
iteration 77 error = 6.927119148775129e-12
iteration 78 error = 4.713860806170039e-12
iteration 79 error = 3.1545465764504343e-12
iteration 80 error = 1.9657493943078776e-12
iteration 81 error = 1.2237497611515653e-12
iteration 82 error = 8.299364158574324e-13
iteration 83 error = 5.756065313208486e-13
iteration 84 error = 4.107851399532092e-13
iteration 85 error = 2.6545116923689193e-13
iteration 86 error = 1.6948927421017946e-13
iteration 87 error = 1.1040106281723313e-13
iteration 88 error = 7.528580384049475e-14
iteration 89 error = 5.3907824314899504e-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.

[10]:

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

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

it# 0 , res = 0.08119734356602216
it# 1 , res = 0.07696557260572953
it# 2 , res = 0.07344471412858082
it# 3 , res = 0.07025007673116736
it# 4 , res = 0.06735911057304736
it# 5 , res = 0.06472122579225087
it# 6 , res = 0.062294223622381836
it# 7 , res = 0.06004424891557193
it# 8 , res = 0.05794435368653123
it# 9 , res = 0.05597310993397294
it# 10 , res = 0.05411339643696332
it# 11 , res = 0.05235143129498116
it# 12 , res = 0.05067602643031652
it# 13 , res = 0.049078019166521124
it# 14 , res = 0.047549839165887446
it# 15 , res = 0.04608517709824829
it# 16 , res = 0.04467872944058139
it# 17 , res = 0.043326000417520126
it# 18 , res = 0.04202314712044809
it# 19 , res = 0.04076685752716331
it# 20 , res = 0.03955425380652286
it# 21 , res = 0.038382815214659936
it# 22 , res = 0.037250316285503696
it# 23 , res = 0.03615477704179933
it# 24 , res = 0.035094422710862015
it# 25 , res = 0.03406765099651578
it# 26 , res = 0.03307300538722954
it# 27 , res = 0.03210915330716986
it# 28 , res = 0.03117486816799722
it# 29 , res = 0.030269014573674553
it# 30 , res = 0.029390536082149768
it# 31 , res = 0.02853844504667718
it# 32 , res = 0.02771181415333076
it# 33 , res = 0.026909769345597963
it# 34 , res = 0.02613148388613044
it# 35 , res = 0.025376173353034902
it# 36 , res = 0.02464309140603591
it# 37 , res = 0.02393152618837694
it# 38 , res = 0.023240797254965688
it# 39 , res = 0.02257025293720382
it# 40 , res = 0.02191926807110181
it# 41 , res = 0.021287242028414446
it# 42 , res = 0.02067359700122097
it# 43 , res = 0.02007777649909595
it# 44 , res = 0.019499244025139206
it# 45 , res = 0.018937481902962742
it# 46 , res = 0.01839199023151119
it# 47 , res = 0.01786228594851358
it# 48 , res = 0.01734790198658913
it# 49 , res = 0.016848386508685472
it# 50 , res = 0.016363302211716505
it# 51 , res = 0.015892225689077555
it# 52 , res = 0.015434746844209748
it# 53 , res = 0.01499046834862759
it# 54 , res = 0.014559005138852558
it# 55 , res = 0.014139983947550387
it# 56 , res = 0.01373304286488472
it# 57 , res = 0.013337830926694156
it# 58 , res = 0.012954007726595163
it# 59 , res = 0.012581243049534415
it# 60 , res = 0.012219216524660239
it# 61 , res = 0.011867617295679266
it# 62 , res = 0.01152614370711084
it# 63 , res = 0.011194503005064697
it# 64 , res = 0.010872411051340005
it# 65 , res = 0.010559592049799267
it# 66 , res = 0.010255778284096096
it# 67 , res = 0.009960709865945666
it# 68 , res = 0.009674134493220828
it# 69 , res = 0.009395807217235921
it# 70 , res = 0.009125490218650528
it# 71 , res = 0.00886295259148325
it# 72 , res = 0.008607970134778353
it# 73 , res = 0.008360325151510997
it# 74 , res = 0.008119806254357436
it# 75 , res = 0.007886208177987818
it# 76 , res = 0.007659331597570535
it# 77 , res = 0.007438982953202166
it# 78 , res = 0.007224974279999439
it# 79 , res = 0.007017123043610654
it# 80 , res = 0.006815251980920437
it# 81 , res = 0.006619188945737578
it# 82 , res = 0.006428766759272515
it# 83 , res = 0.006243823065218345
it# 84 , res = 0.006064200189265575
it# 85 , res = 0.005889745002886657
it# 86 , res = 0.00572030879124145
it# 87 , res = 0.005555747125053994
it# 88 , res = 0.0053959197363275336
it# 89 , res = 0.005240690397767275
it# 90 , res = 0.005089926805785061
it# 91 , res = 0.004943500466970432
it# 92 , res = 0.004801286587914051
it# 93 , res = 0.0046631639682756165
it# 94 , res = 0.004529014896994936
it# 95 , res = 0.004398725051543881
it# 96 , res = 0.004272183400128201
it# 97 , res = 0.004149282106743879
it# 98 , res = 0.004029916439004409
it# 99 , res = 0.003913984678649816
it# 100 , res = 0.0038013880346601174
it# 101 , res = 0.003692030558892565
it# 102 , res = 0.003585819064165283
it# 103 , res = 0.0034826630447189954
it# 104 , res = 0.0033824745989797455
it# 105 , res = 0.00328516835455907
it# 106 , res = 0.0031906613954232324
it# 107 , res = 0.003098873191168193
it# 108 , res = 0.003009725528337407
it# 109 , res = 0.002923142443723854
it# 110 , res = 0.0028390501595974916
it# 111 , res = 0.002757377020801111
it# 112 , res = 0.0026780534336611917
it# 113 , res = 0.002601011806660501
it# 114 , res = 0.0025261864928203598
it# 115 , res = 0.0024535137337449242
it# 116 , res = 0.0023829316052758204
it# 117 , res = 0.0023143799647147513
it# 118 , res = 0.0022478003995643862
it# 119 , res = 0.0021831361777465536
it# 120 , res = 0.002120332199253914
it# 121 , res = 0.002059334949193202
it# 122 , res = 0.0020000924521809734
it# 123 , res = 0.0019425542280515152
it# 124 , res = 0.00188667124884025
it# 125 , res = 0.0018323958970042887
it# 126 , res = 0.0017796819248471253
it# 127 , res = 0.001728484415108924
it# 128 , res = 0.0016787597426922676
it# 129 , res = 0.001630465537489084
it# 130 , res = 0.0015835606482767284
it# 131 , res = 0.0015380051076530215
it# 132 , res = 0.0014937600979809229
it# 133 , res = 0.0014507879183114146
it# 134 , res = 0.001409051952260031
it# 135 , res = 0.0013685166368059222
it# 136 , res = 0.0013291474319896562
it# 137 , res = 0.001290910791481944
it# 138 , res = 0.0012537741339998386
it# 139 , res = 0.0012177058155444591
it# 140 , res = 0.0011826751024382805
it# 141 , res = 0.0011486521451370662
it# 142 , res = 0.0011156079527963382
it# 143 , res = 0.0010835143685682465
it# 144 , res = 0.0010523440456102749
it# 145 , res = 0.001022070423782998
it# 146 , res = 0.0009926677070188294
it# 147 , res = 0.0009641108413425826
it# 148 , res = 0.0009363754935217705
it# 149 , res = 0.0009094380303344205
it# 150 , res = 0.0008832754984305099
it# 151 , res = 0.0008578656047740584
it# 152 , res = 0.0008331866976479217
it# 153 , res = 0.0008092177482045541
it# 154 , res = 0.000785938332547231
it# 155 , res = 0.0007633286143278866
it# 156 , res = 0.0007413693278441148
it# 157 , res = 0.0007200417616235552
it# 158 , res = 0.0006993277424800745
it# 159 , res = 0.0006792096200284055
it# 160 , res = 0.0006596702516446203
it# 161 , res = 0.0006406929878589201
it# 162 , res = 0.0006222616581690505
it# 163 , res = 0.00060436055726136
it# 164 , res = 0.0005869744316283719
it# 165 , res = 0.0005700884665716965
it# 166 , res = 0.0005536882735783817
it# 167 , res = 0.000537759878060613
it# 168 , res = 0.0005222897074479053
it# 169 , res = 0.0005072645796222377
it# 170 , res = 0.0004926716916858104
it# 171 , res = 0.0004784986090509398
it# 172 , res = 0.00046473325484572476
it# 173 , res = 0.00045136389962282847
it# 174 , res = 0.0004383791513645787
it# 175 , res = 0.00042576794577677067
it# 176 , res = 0.0004135195368599968
it# 177 , res = 0.00040162348775344846
it# 178 , res = 0.00039006966184131563
it# 179 , res = 0.0003788482141164032
it# 180 , res = 0.0003679495827897905
it# 181 , res = 0.0003573644811442805
it# 182 , res = 0.000347083889621145
it# 183 , res = 0.00033709904813420744
it# 184 , res = 0.0003274014486056376
it# 185 , res = 0.0003179828277161132
it# 186 , res = 0.0003088351598637657
it# 187 , res = 0.00029995065032554954
it# 188 , res = 0.0002913217286154457
it# 189 , res = 0.00028294104203361865
it# 190 , res = 0.00027480144940109265
it# 191 , res = 0.0002668960149748972
it# 192 , res = 0.00025921800253807676
it# 193 , res = 0.000251760869659655
it# 194 , res = 0.0002445182621203887
it# 195 , res = 0.00023748400849752467
it# 196 , res = 0.00023065211490659756
it# 197 , res = 0.0002240167598942669
it# 198 , res = 0.00021757228947727304
it# 199 , res = 0.00021131321232509877
it# 200 , res = 0.00020523419508077932
it# 201 , res = 0.000199330057816347
it# 202 , res = 0.00019359576961870548
it# 203 , res = 0.00018802644430348806
it# 204 , res = 0.00018261733625069224
it# 205 , res = 0.00017736383636165128
it# 206 , res = 0.00017226146813114382
it# 207 , res = 0.00016730588383317808
it# 208 , res = 0.00016249286081648773
it# 209 , res = 0.00015781829790597314
it# 210 , res = 0.0001532782119084258
it# 211 , res = 0.00014886873421832
it# 212 , res = 0.0001445861075214198
it# 213 , res = 0.00014042668259322085
it# 214 , res = 0.00013638691518935573
it# 215 , res = 0.00013246336302542817
it# 216 , res = 0.00012865268284443935
it# 217 , res = 0.00012495162756690073
it# 218 , res = 0.00012135704352520126
it# 219 , res = 0.0001178658677755353
it# 220 , res = 0.00011447512548853386
it# 221 , res = 0.00011118192741384818
it# 222 , res = 0.00010798346741852885
it# 223 , res = 0.00010487702009640282
it# 224 , res = 0.0001018599384444072
it# 225 , res = 9.892965160873351e-05
it# 226 , res = 9.608366269287981e-05
it# 227 , res = 9.331954663063205e-05
it# 228 , res = 9.063494812001635e-05
it# 229 , res = 8.802757961551883e-05
it# 230 , res = 8.5495219379298e-05
it# 231 , res = 8.303570958838247e-05
it# 232 , res = 8.064695449507518e-05
it# 233 , res = 7.832691864221051e-05
it# 234 , res = 7.607362512815404e-05
it# 235 , res = 7.38851539222216e-05
it# 236 , res = 7.175964022889903e-05
it# 237 , res = 6.96952728988565e-05
it# 238 , res = 6.769029288535838e-05
it# 239 , res = 6.574299174563177e-05
it# 240 , res = 6.385171018519757e-05
it# 241 , res = 6.201483664378394e-05
it# 242 , res = 6.0230805921899757e-05
it# 243 , res = 5.849809784785605e-05
it# 244 , res = 5.681523598151747e-05
it# 245 , res = 5.5180786357110256e-05
it# 246 , res = 5.3593356260585004e-05
it# 247 , res = 5.205159304330611e-05
it# 248 , res = 5.0554182969542835e-05
it# 249 , res = 4.909985009658866e-05
it# 250 , res = 4.768735518807469e-05
it# 251 , res = 4.631549465744327e-05
it# 252 , res = 4.498309954286539e-05
it# 253 , res = 4.36890345111116e-05
it# 254 , res = 4.243219688964488e-05
it# 255 , res = 4.121151572765292e-05
it# 256 , res = 4.002595088329978e-05
it# 257 , res = 3.8874492136984886e-05
it# 258 , res = 3.775615833118976e-05
it# 259 , res = 3.6669996534226535e-05
it# 260 , res = 3.561508122802995e-05
it# 261 , res = 3.45905135193489e-05
it# 262 , res = 3.359542037477192e-05
it# 263 , res = 3.262895387554767e-05
it# 264 , res = 3.1690290496133045e-05
it# 265 , res = 3.077863040160895e-05
it# 266 , res = 2.9893196766707928e-05
it# 267 , res = 2.9033235114061888e-05
it# 268 , res = 2.8198012670473635e-05
it# 269 , res = 2.7386817743212115e-05
it# 270 , res = 2.6598959113254317e-05
it# 271 , res = 2.5833765446684974e-05
it# 272 , res = 2.5090584722258722e-05
it# 273 , res = 2.436878367582749e-05
it# 274 , res = 2.3667747261157187e-05
it# 275 , res = 2.2986878125337466e-05
it# 276 , res = 2.2325596099946585e-05
it# 277 , res = 2.1683337706852214e-05
it# 278 , res = 2.1059555677893493e-05
it# 279 , res = 2.0453718488728933e-05
it# 280 , res = 1.9865309905590275e-05
it# 281 , res = 1.9293828545663746e-05
it# 282 , res = 1.873878745013851e-05
it# 283 , res = 1.8199713668585545e-05
it# 284 , res = 1.7676147856358916e-05
it# 285 , res = 1.716764388322107e-05
it# 286 , res = 1.6673768453234718e-05
it# 287 , res = 1.6194100735217993e-05
it# 288 , res = 1.572823200465937e-05
it# 289 , res = 1.5275765294734605e-05
it# 290 , res = 1.483631505894191e-05
it# 291 , res = 1.440950684171057e-05
it# 292 , res = 1.399497696002339e-05
it# 293 , res = 1.3592372193153448e-05
it# 294 , res = 1.3201349481631834e-05
it# 295 , res = 1.2821575635203542e-05
it# 296 , res = 1.2452727048636231e-05
it# 297 , res = 1.2094489426362285e-05
it# 298 , res = 1.1746557514090899e-05
it# 299 , res = 1.1408634839154064e-05
it# 300 , res = 1.1080433457791323e-05
it# 301 , res = 1.0761673709773286e-05
it# 302 , res = 1.0452083980019454e-05
it# 303 , res = 1.0151400467229204e-05
it# 304 , res = 9.859366958950267e-06
it# 305 , res = 9.575734613805014e-06
it# 306 , res = 9.300261748557444e-06
it# 307 , res = 9.032713632887942e-06
it# 308 , res = 8.772862289351362e-06
it# 309 , res = 8.52048629847047e-06
it# 310 , res = 8.275370611002951e-06
it# 311 , res = 8.037306363685937e-06
it# 312 , res = 7.806090702165665e-06
it# 313 , res = 7.581526607646251e-06
it# 314 , res = 7.3634227293118e-06
it# 315 , res = 7.1515932208259704e-06
it# 316 , res = 6.945857582274746e-06
it# 317 , res = 6.74604050629952e-06
it# 318 , res = 6.551971729182755e-06
it# 319 , res = 6.363485884694585e-06
it# 320 , res = 6.180422364198965e-06
it# 321 , res = 6.002625179413754e-06
it# 322 , res = 5.8299428292598825e-06
it# 323 , res = 5.66222817129389e-06
it# 324 , res = 5.499338296012727e-06
it# 325 , res = 5.341134404672039e-06
it# 326 , res = 5.18748169216689e-06
it# 327 , res = 5.038249230868083e-06
it# 328 , res = 4.893309859858835e-06
it# 329 , res = 4.752540076507952e-06
it# 330 , res = 4.615819930754405e-06
it# 331 , res = 4.483032923469462e-06
it# 332 , res = 4.354065906944336e-06
it# 333 , res = 4.228808988352778e-06
it# 334 , res = 4.1071554363812995e-06
it# 335 , res = 3.989001590112266e-06
it# 336 , res = 3.8742467705605776e-06
it# 337 , res = 3.762793195206251e-06
it# 338 , res = 3.6545458947066984e-06
it# 339 , res = 3.549412631384941e-06
it# 340 , res = 3.44730382126366e-06
it# 341 , res = 3.348132457516503e-06
it# 342 , res = 3.2518140361437767e-06
it# 343 , res = 3.1582664843194334e-06
it# 344 , res = 3.067410090307078e-06
it# 345 , res = 2.9791674351533595e-06
it# 346 , res = 2.893463327586205e-06
it# 347 , res = 2.8102247390729675e-06
it# 348 , res = 2.7293807420582952e-06
it# 349 , res = 2.6508624493284514e-06
it# 350 , res = 2.5746029555447406e-06
it# 351 , res = 2.5005372800047275e-06
it# 352 , res = 2.428602311397095e-06
it# 353 , res = 2.3587367537913716e-06
it# 354 , res = 2.2908810751191085e-06
it# 355 , res = 2.224977455242853e-06
it# 356 , res = 2.1609697379225345e-06
it# 357 , res = 2.0988033820302652e-06
it# 358 , res = 2.038425415820471e-06
it# 359 , res = 1.979784391192782e-06
it# 360 , res = 1.9228303400422064e-06
it# 361 , res = 1.86751473182168e-06
it# 362 , res = 1.8137904323959279e-06
it# 363 , res = 1.7616116631852195e-06
it# 364 , res = 1.710933962500022e-06
it# 365 , res = 1.6617141480405906e-06
it# 366 , res = 1.613910279576796e-06
it# 367 , res = 1.567481623597086e-06
it# 368 , res = 1.522388617989059e-06
it# 369 , res = 1.4785928391996155e-06
it# 370 , res = 1.4360569688617811e-06
it# 371 , res = 1.3947447621082154e-06
it# 372 , res = 1.3546210169659847e-06
it# 373 , res = 1.3156515438445063e-06
it# 374 , res = 1.27780313688417e-06
it# 375 , res = 1.2410435455193557e-06
it# 376 , res = 1.20534144690881e-06
it# 377 , res = 1.1706664192937224e-06
it# 378 , res = 1.1369889161115048e-06
it# 379 , res = 1.1042802408675386e-06
it# 380 , res = 1.072512522340481e-06
it# 381 , res = 1.0416586912721046e-06
it# 382 , res = 1.0116924572609038e-06
it# 383 , res = 9.825882860519578e-07
it# 384 , res = 9.543213779212043e-07
it# 385 , res = 9.268676467172414e-07
it# 386 , res = 9.002036990374693e-07
it# 387 , res = 8.743068146792112e-07
it# 388 , res = 8.49154926748328e-07
it# 389 , res = 8.247266035600695e-07
it# 390 , res = 8.010010294992781e-07
it# 391 , res = 7.779579880918734e-07
it# 392 , res = 7.555778445297398e-07
it# 393 , res = 7.338415284689586e-07
it# 394 , res = 7.127305185586999e-07
it# 395 , res = 6.922268259952688e-07
it# 396 , res = 6.723129797348643e-07
it# 397 , res = 6.52972011145823e-07
it# 398 , res = 6.341874397022707e-07
it# 399 , res = 6.159432591872167e-07
it# 400 , res = 5.982239236226157e-07
it# 401 , res = 5.810143344804237e-07
it# 402 , res = 5.642998274814528e-07
it# 403 , res = 5.48066160098755e-07
it# 404 , res = 5.322994998415674e-07
it# 405 , res = 5.169864116029918e-07
it# 406 , res = 5.021138473580336e-07
it# 407 , res = 4.876691342245933e-07
it# 408 , res = 4.736399636878286e-07
it# 409 , res = 4.6001438155557564e-07
it# 410 , res = 4.467807776072234e-07
it# 411 , res = 4.3392787531824345e-07
it# 412 , res = 4.214447228242701e-07
it# 413 , res = 4.093206832812781e-07
it# 414 , res = 3.975454257726523e-07
it# 415 , res = 3.8610891653779513e-07
it# 416 , res = 3.7500141055927573e-07
it# 417 , res = 3.642134431479294e-07
it# 418 , res = 3.5373582178470617e-07
it# 419 , res = 3.4355961871592467e-07
it# 420 , res = 3.336761625887141e-07
it# 421 , res = 3.240770319307513e-07
it# 422 , res = 3.1475404704815144e-07
it# 423 , res = 3.056992640780088e-07
it# 424 , res = 2.9690496729704544e-07
it# 425 , res = 2.8836366314048594e-07
it# 426 , res = 2.8006807353371656e-07
it# 427 , res = 2.720111297267254e-07
it# 428 , res = 2.641859666903443e-07
it# 429 , res = 2.565859162847454e-07
it# 430 , res = 2.4920450273121613e-07
it# 431 , res = 2.4203543627645496e-07
it# 432 , res = 2.3507260804371125e-07
it# 433 , res = 2.283100851126359e-07
it# 434 , res = 2.2174210510631445e-07
it# 435 , res = 2.1536307145382084e-07
it# 436 , res = 2.0916754861423567e-07
it# 437 , res = 2.031502573761954e-07
it# 438 , res = 1.9730607042775195e-07
it# 439 , res = 1.9163000789331896e-07
it# 440 , res = 1.861172332360944e-07
it# 441 , res = 1.807630489240887e-07
it# 442 , res = 1.7556289285271858e-07
it# 443 , res = 1.7051233376704244e-07
it# 444 , res = 1.6560706824401367e-07
it# 445 , res = 1.608429163935476e-07
it# 446 , res = 1.5621581875407193e-07
it# 447 , res = 1.5172183250744885e-07
it# 448 , res = 1.4735712843300102e-07
it# 449 , res = 1.4311798736849158e-07
it# 450 , res = 1.3900079692954408e-07
it# 451 , res = 1.3500204922810624e-07
it# 452 , res = 1.3111833662057514e-07
it# 453 , res = 1.2734634987053666e-07
it# 454 , res = 1.2368287493932198e-07
it# 455 , res = 1.20124790127106e-07
it# 456 , res = 1.1666906358766115e-07
it# 457 , res = 1.1331275074173177e-07
it# 458 , res = 1.100529916738983e-07
it# 459 , res = 1.0688700855493636e-07
it# 460 , res = 1.0381210399306976e-07
it# 461 , res = 1.0082565762606738e-07
it# 462 , res = 9.79251248214122e-08
it# 463 , res = 9.510803392718951e-08
it# 464 , res = 9.23719845903885e-08
it# 465 , res = 8.971464531215289e-08
it# 466 , res = 8.71337519116926e-08
it# 467 , res = 8.462710516719772e-08
it# 468 , res = 8.219256903327963e-08
it# 469 , res = 7.982806920203638e-08
it# 470 , res = 7.753159096448283e-08
it# 471 , res = 7.530117723681155e-08
it# 472 , res = 7.313492769685997e-08
it# 473 , res = 7.103099633872838e-08
it# 474 , res = 6.898759056141798e-08
it# 475 , res = 6.700296896013302e-08
it# 476 , res = 6.507544066154603e-08
it# 477 , res = 6.3203363048093e-08
it# 478 , res = 6.138514099161166e-08
it# 479 , res = 5.961922529860007e-08
it# 480 , res = 5.790411099885679e-08
it# 481 , res = 5.623833684638744e-08
it# 482 , res = 5.462048343618774e-08
it# 483 , res = 5.304917195076123e-08
it# 484 , res = 5.1523063722207655e-08
it# 485 , res = 5.0040858239050835e-08
it# 486 , res = 4.86012926604021e-08
it# 487 , res = 4.7203140117398695e-08
it# 488 , res = 4.584520943575234e-08
it# 489 , res = 4.4526343338363224e-08
it# 490 , res = 4.3245418146443594e-08
it# 491 , res = 4.2001342484177105e-08
it# 492 , res = 4.079305600114823e-08
it# 493 , res = 3.961952933918804e-08
it# 494 , res = 3.84797624760545e-08
it# 495 , res = 3.737278413472251e-08
it# 496 , res = 3.629765116928812e-08
it# 497 , res = 3.5253447403659214e-08
it# 498 , res = 3.4239283263301304e-08
it# 499 , res = 3.32542941094594e-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.

[11]:

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

[12]:

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

iteration 0 error = 0.09386193808775962
iteration 1 error = 0.12289904446040574
iteration 2 error = 0.08558846385984858
iteration 3 error = 0.053786375091255334
iteration 4 error = 0.02799513368494736
iteration 5 error = 0.014969155114208049
iteration 6 error = 0.006689317102568868
iteration 7 error = 0.0028248604314505604
iteration 8 error = 0.0010332508846612996
iteration 9 error = 0.00040318420091610234
iteration 10 error = 0.000126153484125304
iteration 11 error = 5.134068138544936e-05
iteration 12 error = 1.86267529855379e-05
iteration 13 error = 7.86571512990075e-06
iteration 14 error = 3.310916766942548e-06
iteration 15 error = 1.5461433004860725e-06
iteration 16 error = 6.174156961107636e-07
iteration 17 error = 3.2008170496573533e-07
iteration 18 error = 1.243696230710687e-07
iteration 19 error = 5.906688350434197e-08
iteration 20 error = 2.207041609152557e-08
iteration 21 error = 7.324210399036425e-09
iteration 22 error = 3.0122986220535976e-09
iteration 23 error = 8.584333680651104e-10
iteration 24 error = 3.6570806010889224e-10
iteration 25 error = 1.2692446936158052e-10
iteration 26 error = 4.461371781815046e-11
iteration 27 error = 1.975092094382798e-11
iteration 28 error = 5.5343243734924135e-12
iteration 29 error = 2.307643414710973e-12
iteration 30 error = 6.724334725524573e-13
iteration 31 error = 1.9582752398605173e-13
iteration 32 error = 7.628337311332365e-14

[13]:

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

[13]:

20.294597459170067


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

## A Block Jacobi preconditioner¶

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

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

[14]:

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

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

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


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

[15]:

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

[15]:

35.59893732672581


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

[16]:

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

[16]:

2.9828997439629763


## Add a coarse grid correction¶

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

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

[17]:

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

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

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


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

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

The matrix coarsepre below represents

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

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:

[19]:

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

[19]:

       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.

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.

[20]:

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:

[21]:

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

[21]:

 0.993322
0.998163
0.999953
1.34169
1.80171
1.82729
1.95319
1.97402
1.99976


## Combining multigrid with block smoothers¶

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

[22]:

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

[22]:

<ngsolve.comp.BilinearForm at 0x7f3420a07eb0>


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

[23]:

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.

[24]:

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.

[25]:

mgblock = mg.mat + blockgs

[26]:

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

[26]:

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


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