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.0146422
 0.0569646
 0.0940129
 0.111788
 0.150881
 0.204399
 0.255688
 0.355463
 0.444993
 0.530835
 0.648619
 0.755733
 0.84438
 0.938056
 1.04363
 1.17467
 1.26055
 1.36943
 1.49386
 1.59296
 1.73416
 1.84185
 1.93522
 2.04299
 2.14565
 2.21453
 2.31143
 2.40492
 2.45389
 2.51245
 2.53508
 2.59211
 2.80962

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

max(lams)/min(lams)

191.88512696418846

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)

1642.1489948259557

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.055268476235376174
iteration 1 error = 0.09382120513440023
iteration 2 error = 0.10877588296622372
iteration 3 error = 0.08914338928898918
iteration 4 error = 0.09978536222727061
iteration 5 error = 0.08655140715101732
iteration 6 error = 0.08819881897359114
iteration 7 error = 0.0780914158709465
iteration 8 error = 0.058372392825907414
iteration 9 error = 0.045405152293075235
iteration 10 error = 0.035977010634211166
iteration 11 error = 0.028485753643071918
iteration 12 error = 0.02272922001659292
iteration 13 error = 0.02074672348238179
iteration 14 error = 0.015395157671022725
iteration 15 error = 0.011470599703926708
iteration 16 error = 0.009683614939889418
iteration 17 error = 0.007588767248025651
iteration 18 error = 0.005154735556970687
iteration 19 error = 0.00346973056209065
iteration 20 error = 0.0021606808972590963
iteration 21 error = 0.0014099027081058996
iteration 22 error = 0.001072719751387697
iteration 23 error = 0.0007843516628696807
iteration 24 error = 0.0004914082413368085
iteration 25 error = 0.0003508933574579758
iteration 26 error = 0.00026125701601065687
iteration 27 error = 0.00019982976817496307
iteration 28 error = 0.00014533440122529127
iteration 29 error = 0.00010363931587641238
iteration 30 error = 6.750008193327995e-05
iteration 31 error = 4.8156564723316196e-05
iteration 32 error = 3.6071745924109976e-05
iteration 33 error = 2.7393049371856237e-05
iteration 34 error = 2.0481166276782385e-05
iteration 35 error = 1.39798719208204e-05
iteration 36 error = 9.19803775832077e-06
iteration 37 error = 6.373309101789777e-06
iteration 38 error = 5.250342597533671e-06
iteration 39 error = 4.593834900387931e-06
iteration 40 error = 3.3196875424214727e-06
iteration 41 error = 2.3625425056115605e-06
iteration 42 error = 1.9109484031156676e-06
iteration 43 error = 1.657668134139862e-06
iteration 44 error = 1.3644814408819299e-06
iteration 45 error = 9.795170530012462e-07
iteration 46 error = 6.560574209501156e-07
iteration 47 error = 4.542666511267736e-07
iteration 48 error = 3.754159339010875e-07
iteration 49 error = 2.626593991994966e-07
iteration 50 error = 1.6770849068299545e-07
iteration 51 error = 1.0732847020427294e-07
iteration 52 error = 8.216337916002824e-08
iteration 53 error = 6.543169481035481e-08
iteration 54 error = 4.779373611111277e-08
iteration 55 error = 3.065008044965062e-08
iteration 56 error = 1.8714955266966125e-08
iteration 57 error = 1.4082717183446533e-08
iteration 58 error = 1.0189035368512845e-08
iteration 59 error = 6.576458149482347e-09
iteration 60 error = 4.144888619072247e-09
iteration 61 error = 2.750764930997574e-09
iteration 62 error = 2.043121950091988e-09
iteration 63 error = 1.392974103940359e-09
iteration 64 error = 9.308009514199542e-10
iteration 65 error = 6.897768149208763e-10
iteration 66 error = 4.657739641876209e-10
iteration 67 error = 3.3160617077672524e-10
iteration 68 error = 2.3811697810351095e-10
iteration 69 error = 1.5537382691243538e-10
iteration 70 error = 1.0126239363372569e-10
iteration 71 error = 7.094849986187965e-11
iteration 72 error = 5.206447757962135e-11
iteration 73 error = 3.498319984639181e-11
iteration 74 error = 2.0349906467466357e-11
iteration 75 error = 1.3508851408212583e-11
iteration 76 error = 9.218490988072588e-12
iteration 77 error = 6.780321458922074e-12
iteration 78 error = 4.700831776438064e-12
iteration 79 error = 2.7424147474774167e-12
iteration 80 error = 1.563232225319326e-12
iteration 81 error = 1.0379361277767889e-12
iteration 82 error = 7.726462874829239e-13
iteration 83 error = 5.462441363428167e-13
iteration 84 error = 3.170907566949878e-13
iteration 85 error = 1.8626204487820541e-13
iteration 86 error = 1.2565688713653832e-13
iteration 87 error = 1.0082288120266588e-13
iteration 88 error = 7.093469005091741e-14
iteration 89 error = 4.2433220344958173e-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.08199396652265978
it# 1 , res = 0.07731725740499301
it# 2 , res = 0.07375945762656286
it# 3 , res = 0.07058028411183849
it# 4 , res = 0.06771422592318167
it# 5 , res = 0.06509485073752364
it# 6 , res = 0.06267525154173857
it# 7 , res = 0.060422254822258284
it# 8 , res = 0.05831113697252155
it# 9 , res = 0.05632277687663107
it# 10 , res = 0.054442004881023895
it# 11 , res = 0.052656559707675334
it# 12 , res = 0.050956386378755825
it# 13 , res = 0.04933314170210069
it# 14 , res = 0.047779833628449
it# 15 , res = 0.046290550838356405
it# 16 , res = 0.044860255462346225
it# 17 , res = 0.04348462142718085
it# 18 , res = 0.042159906689899974
it# 19 , res = 0.04088285121042283
it# 20 , res = 0.03965059482369328
it# 21 , res = 0.03846061071164193
it# 22 , res = 0.03731065123608623
it# 23 , res = 0.03619870364788998
it# 24 , res = 0.035122953738665544
it# 25 , res = 0.03408175591305893
it# 26 , res = 0.03307360847313405
it# 27 , res = 0.03209713314859922
it# 28 , res = 0.031151058096000883
it# 29 , res = 0.030234203739481557
it# 30 , res = 0.029345470944546637
it# 31 , res = 0.028483831111360963
it# 32 , res = 0.02764831785052608
it# 33 , res = 0.02683801996597428
it# 34 , res = 0.02605207551956936
it# 35 , res = 0.025289666792567382
it# 36 , res = 0.02455001599210774
it# 37 , res = 0.023832381577833097
it# 38 , res = 0.02313605510573882
it# 39 , res = 0.022460358504356317
it# 40 , res = 0.021804641713130575
it# 41 , res = 0.021168280624956216
it# 42 , res = 0.020550675284785894
it# 43 , res = 0.01995124830440415
it# 44 , res = 0.019369443460199115
it# 45 , res = 0.01880472444631869
it# 46 , res = 0.01825657376018289
it# 47 , res = 0.01772449170111331
it# 48 , res = 0.017207995465977878
it# 49 , res = 0.016706618328344167
it# 50 , res = 0.01621990888979102
it# 51 , res = 0.015747430393820774
it# 52 , res = 0.015288760094297965
it# 53 , res = 0.014843488671588953
it# 54 , res = 0.014411219690607099
it# 55 , res = 0.013991569095836244
it# 56 , res = 0.013584164739130503
it# 57 , res = 0.013188645936694793
it# 58 , res = 0.012804663052164822
it# 59 , res = 0.012431877103132444
it# 60 , res = 0.012069959388827186
it# 61 , res = 0.011718591136972262
it# 62 , res = 0.011377463168093113
it# 63 , res = 0.011046275575776041
it# 64 , res = 0.010724737421564946
it# 65 , res = 0.010412566443345319
it# 66 , res = 0.010109488776196367
it# 67 , res = 0.009815238684817445
it# 68 , res = 0.009529558306730095
it# 69 , res = 0.00925219740554639
it# 70 , res = 0.008982913133670868
it# 71 , res = 0.008721469803864243
it# 72 , res = 0.008467638669158364
it# 73 , res = 0.008221197710658483
it# 74 , res = 0.0079819314328115
it# 75 , res = 0.0077496306657581235
it# 76 , res = 0.007524092374418903
it# 77 , res = 0.007305119473992156
it# 78 , res = 0.007092520651570713
it# 79 , res = 0.006886110193602168
it# 80 , res = 0.006685707818942339
it# 81 , res = 0.0064911385172659155
it# 82 , res = 0.006302232392616393
it# 83 , res = 0.00611882451189141
it# 84 , res = 0.0059407547580709475
it# 85 , res = 0.005767867688009802
it# 86 , res = 0.005600012394624387
it# 87 , res = 0.005437042373314604
it# 88 , res = 0.005278815392468837
it# 89 , res = 0.005125193367909578
it# 90 , res = 0.004976042241143591
it# 91 , res = 0.004831231861286871
it# 92 , res = 0.0046906358705410146
it# 93 , res = 0.004554131593103938
it# 94 , res = 0.004421599927401679
it# 95 , res = 0.004292925241534117
it# 96 , res = 0.004167995271830852
it# 97 , res = 0.004046701024419974
it# 98 , res = 0.0039289366797117175
it# 99 , res = 0.0038145994997084327
it# 100 , res = 0.0037035897380515766
it# 101 , res = 0.0035958105527214034
it# 102 , res = 0.003491167921307608
it# 103 , res = 0.003389570558772414
it# 104 , res = 0.0032909298376309144
it# 105 , res = 0.0031951597104734344
it# 106 , res = 0.0031021766347616234
it# 107 , res = 0.003011899499828273
it# 108 , res = 0.002924249556014332
it# 109 , res = 0.002839150345882201
it# 110 , res = 0.002756527637439114
it# 111 , res = 0.002676309359313116
it# 112 , res = 0.0025984255378251517
it# 113 , res = 0.0025228082358972116
it# 114 , res = 0.0024493914937454703
it# 115 , res = 0.0023781112713047374
it# 116 , res = 0.002308905392332217
it# 117 , res = 0.0022417134901425453
it# 118 , res = 0.0021764769549254557
it# 119 , res = 0.002113138882599064
it# 120 , res = 0.002051644025154375
it# 121 , res = 0.001991938742448066
it# 122 , res = 0.001933970955398641
it# 123 , res = 0.001877690100547108
it# 124 , res = 0.0018230470859415598
it# 125 , res = 0.0017699942483066927
it# 126 , res = 0.0017184853114610517
it# 127 , res = 0.001668475345944911
it# 128 , res = 0.0016199207298241971
it# 129 , res = 0.001572779110635775
it# 130 , res = 0.0015270093684409869
it# 131 , res = 0.0014825715799541468
it# 132 , res = 0.0014394269837163288
it# 133 , res = 0.0013975379462827344
it# 134 , res = 0.0013568679293938085
it# 135 , res = 0.0013173814581028712
it# 136 , res = 0.0012790440898312546
it# 137 , res = 0.001241822384323691
it# 138 , res = 0.001205683874478823
it# 139 , res = 0.0011705970380286395
it# 140 , res = 0.0011365312700416317
it# 141 , res = 0.0011034568562274599
it# 142 , res = 0.0010713449470170623
it# 143 , res = 0.0010401675323984077
it# 144 , res = 0.0010098974174837362
it# 145 , res = 0.0009805081987888814
it# 146 , res = 0.0009519742412014318
it# 147 , res = 0.0009242706556211303
it# 148 , res = 0.0008973732772491384
it# 149 , res = 0.0008712586445101768
it# 150 , res = 0.0008459039785879808
it# 151 , res = 0.0008212871635560823
it# 152 , res = 0.00079738672708707
it# 153 , res = 0.0007741818217226824
it# 154 , res = 0.000751652206690006
it# 155 , res = 0.0007297782302458191
it# 156 , res = 0.0007085408125351023
it# 157 , res = 0.000687921428948629
it# 158 , res = 0.0006679020939649223
it# 159 , res = 0.0006484653454615427
it# 160 , res = 0.0006295942294843894
it# 161 , res = 0.0006112722854585709
it# 162 , res = 0.0005934835318313558
it# 163 , res = 0.0005762124521314797
it# 164 , res = 0.0005594439814351318
it# 165 , res = 0.0005431634932252416
it# 166 , res = 0.0005273567866333878
it# 167 , res = 0.0005120100740527534
it# 168 , res = 0.0004971099691124781
it# 169 , res = 0.000482643475000023
it# 170 , res = 0.0004685979731259661
it# 171 , res = 0.0004549612121161387
it# 172 , res = 0.00044172129712579555
it# 173 , res = 0.00042886667946417245
it# 174 , res = 0.0004163861465210031
it# 175 , res = 0.0004042688119862393
it# 176 , res = 0.0003925041063538846
it# 177 , res = 0.0003810817677033673
it# 178 , res = 0.0003699918327479813
it# 179 , res = 0.000359224628144358
it# 180 , res = 0.0003487707620550286
it# 181 , res = 0.0003386211159558031
it# 182 , res = 0.0003287668366824013
it# 183 , res = 0.00031919932870817195
it# 184 , res = 0.0003099102466466672
it# 185 , res = 0.00030089148797168386
it# 186 , res = 0.0002921351859502598
it# 187 , res = 0.00028363370278081167
it# 188 , res = 0.00027537962293064784
it# 189 , res = 0.0002673657466680166
it# 190 , res = 0.0002595850837818965
it# 191 , res = 0.0002520308474847589
it# 192 , res = 0.00024469644849252516
it# 193 , res = 0.00023757548927710257
it# 194 , res = 0.0002306617584863183
it# 195 , res = 0.00022394922552545428
it# 196 , res = 0.0002174320352974221
it# 197 , res = 0.0002111045030953945
it# 198 , res = 0.00020496110964420515
it# 199 , res = 0.00019899649628651468
it# 200 , res = 0.00019320546030786358
it# 201 , res = 0.00018758295039912723
it# 202 , res = 0.0001821240622502779
it# 203 , res = 0.00017682403427280333
it# 204 , res = 0.0001716782434457796
it# 205 , res = 0.00016668220128384465
it# 206 , res = 0.00016183154992193183
it# 207 , res = 0.0001571220583138212
it# 208 , res = 0.0001525496185421695
it# 209 , res = 0.00014811024223477119
it# 210 , res = 0.00014380005708565806
it# 211 , res = 0.00013961530347783046
it# 212 , res = 0.00013555233120333225
it# 213 , res = 0.0001316075962799637
it# 214 , res = 0.00012777765785943902
it# 215 , res = 0.0001240591752262794
it# 216 , res = 0.00012044890488370944
it# 217 , res = 0.00011694369772498202
it# 218 , res = 0.00011354049628559912
it# 219 , res = 0.0001102363320773642
it# 220 , res = 0.00010702832299837968
it# 221 , res = 0.00010391367081932578
it# 222 , res = 0.0001008896587429546
it# 223 , res = 9.795364903389919e-05
it# 224 , res = 9.510308071809582e-05
it# 225 , res = 9.233546734860298e-05
it# 226 , res = 8.964839483763948e-05
it# 227 , res = 8.703951934966807e-05
it# 228 , res = 8.450656525792764e-05
it# 229 , res = 8.204732315890951e-05
it# 230 , res = 7.965964794571668e-05
it# 231 , res = 7.734145693626838e-05
it# 232 , res = 7.509072805711113e-05
it# 233 , res = 7.290549807962978e-05
it# 234 , res = 7.078386090751606e-05
it# 235 , res = 6.872396591403433e-05
it# 236 , res = 6.672401632798244e-05
it# 237 , res = 6.478226766636268e-05
it# 238 , res = 6.289702621268093e-05
it# 239 , res = 6.10666475393615e-05
it# 240 , res = 5.928953507443867e-05
it# 241 , res = 5.7564138707116403e-05
it# 242 , res = 5.5888953437252723e-05
it# 243 , res = 5.426251806178346e-05
it# 244 , res = 5.268341390068297e-05
it# 245 , res = 5.1150263558853625e-05
it# 246 , res = 4.96617297252513e-05
it# 247 , res = 4.821651400611386e-05
it# 248 , res = 4.681335579261834e-05
it# 249 , res = 4.545103116096362e-05
it# 250 , res = 4.412835180505032e-05
it# 251 , res = 4.284416400008339e-05
it# 252 , res = 4.159734759579847e-05
it# 253 , res = 4.0386815039926855e-05
it# 254 , res = 3.9211510428683486e-05
it# 255 , res = 3.807040858722379e-05
it# 256 , res = 3.696251417375631e-05
it# 257 , res = 3.5886860812497565e-05
it# 258 , res = 3.48425102502679e-05
it# 259 , res = 3.382855153824706e-05
it# 260 , res = 3.284410023727345e-05
it# 261 , res = 3.188829764626996e-05
it# 262 , res = 3.096031005360359e-05
it# 263 , res = 3.005932800943092e-05
it# 264 , res = 2.918456562032126e-05
it# 265 , res = 2.8335259862581566e-05
it# 266 , res = 2.751066991796155e-05
it# 267 , res = 2.671007652694327e-05
it# 268 , res = 2.5932781360723084e-05
it# 269 , res = 2.5178106413420863e-05
it# 270 , res = 2.4445393409430762e-05
it# 271 , res = 2.3734003230188554e-05
it# 272 , res = 2.3043315355594082e-05
it# 273 , res = 2.237272732419575e-05
it# 274 , res = 2.1721654206306987e-05
it# 275 , res = 2.1089528094718563e-05
it# 276 , res = 2.047579760885736e-05
it# 277 , res = 1.9879927414041956e-05
it# 278 , res = 1.930139775438914e-05
it# 279 , res = 1.873970399968147e-05
it# 280 , res = 1.8194356204732325e-05
it# 281 , res = 1.7664878682859914e-05
it# 282 , res = 1.7150809589990343e-05
it# 283 , res = 1.665170052242772e-05
it# 284 , res = 1.616711612566449e-05
it# 285 , res = 1.569663371409862e-05
it# 286 , res = 1.5239842903356072e-05
it# 287 , res = 1.4796345251561002e-05
it# 288 , res = 1.4365753911625136e-05
it# 289 , res = 1.3947693294535714e-05
it# 290 , res = 1.3541798741404767e-05
it# 291 , res = 1.3147716205063178e-05
it# 292 , res = 1.2765101941873399e-05
it# 293 , res = 1.2393622211144193e-05
it# 294 , res = 1.2032952984882649e-05
it# 295 , res = 1.1682779664291678e-05
it# 296 , res = 1.1342796806137763e-05
it# 297 , res = 1.101270785562064e-05
it# 298 , res = 1.069222488833597e-05
it# 299 , res = 1.0381068358594983e-05
it# 300 , res = 1.0078966855942043e-05
it# 301 , res = 9.785656868256144e-06
it# 302 , res = 9.500882551887051e-06
it# 303 , res = 9.22439550876982e-06
it# 304 , res = 8.955954569213672e-06
it# 305 , res = 8.695325581869288e-06
it# 306 , res = 8.442281209817224e-06
it# 307 , res = 8.196600731595376e-06
it# 308 , res = 7.958069848991073e-06
it# 309 , res = 7.726480500359148e-06
it# 310 , res = 7.5016306786672095e-06
it# 311 , res = 7.283324255644664e-06
it# 312 , res = 7.071370810608158e-06
it# 313 , res = 6.865585464213489e-06
it# 314 , res = 6.665788717627288e-06
it# 315 , res = 6.4718062950551606e-06
it# 316 , res = 6.283468993068674e-06
it# 317 , res = 6.100612531888193e-06
it# 318 , res = 5.923077412204406e-06
it# 319 , res = 5.750708776853639e-06
it# 320 , res = 5.583356274861587e-06
it# 321 , res = 5.420873930863177e-06
it# 322 , res = 5.263120017354191e-06
it# 323 , res = 5.109956931301146e-06
it# 324 , res = 4.961251074243699e-06
it# 325 , res = 4.81687273545772e-06
it# 326 , res = 4.6766959789478486e-06
it# 327 , res = 4.540598533663811e-06
it# 328 , res = 4.408461686839927e-06
it# 329 , res = 4.280170180225982e-06
it# 330 , res = 4.15561210983039e-06
it# 331 , res = 4.0346788283394116e-06
it# 332 , res = 3.917264849826195e-06
it# 333 , res = 3.8032677585993065e-06
it# 334 , res = 3.6925881189076508e-06
it# 335 , res = 3.5851293891525383e-06
it# 336 , res = 3.480797836891751e-06
it# 337 , res = 3.379502457538117e-06
it# 338 , res = 3.2811548950345725e-06
it# 339 , res = 3.1856693639969216e-06
it# 340 , res = 3.0929625761941997e-06
it# 341 , res = 3.002953666685411e-06
it# 342 , res = 2.9155641242091964e-06
it# 343 , res = 2.8307177217457707e-06
it# 344 , res = 2.7483404512101958e-06
it# 345 , res = 2.668360457703739e-06
it# 346 , res = 2.590707977610022e-06
it# 347 , res = 2.515315277535422e-06
it# 348 , res = 2.442116595264478e-06
it# 349 , res = 2.371048082090965e-06
it# 350 , res = 2.302047747675439e-06
it# 351 , res = 2.2350554055017598e-06
it# 352 , res = 2.1700126205236743e-06
it# 353 , res = 2.106862658411901e-06
it# 354 , res = 2.0455504357193963e-06
it# 355 , res = 1.986022472071452e-06
it# 356 , res = 1.9282268432623813e-06
it# 357 , res = 1.8721131364555543e-06
it# 358 , res = 1.81763240557013e-06
it# 359 , res = 1.7647371291272442e-06
it# 360 , res = 1.7133811684660818e-06
it# 361 , res = 1.6635197276141179e-06
it# 362 , res = 1.61510931424864e-06
it# 363 , res = 1.5681077019330307e-06
it# 364 , res = 1.5224738927142603e-06
it# 365 , res = 1.478168081896708e-06
it# 366 , res = 1.4351516231327402e-06
it# 367 , res = 1.3933869947390823e-06
it# 368 , res = 1.3528377671761933e-06
it# 369 , res = 1.3134685704243746e-06
it# 370 , res = 1.2752450642842654e-06
it# 371 , res = 1.2381339081393516e-06
it# 372 , res = 1.2021027308623645e-06
it# 373 , res = 1.167120103979903e-06
it# 374 , res = 1.133155513229285e-06
it# 375 , res = 1.100179332879383e-06
it# 376 , res = 1.0681627987026839e-06
it# 377 , res = 1.0370779840025082e-06
it# 378 , res = 1.0068977745774457e-06
it# 379 , res = 9.77595845225149e-07
it# 380 , res = 9.49146637190966e-07
it# 381 , res = 9.215253350463405e-07
it# 382 , res = 8.94707845872964e-07
it# 383 , res = 8.686707776926894e-07
it# 384 , res = 8.433914194122403e-07
it# 385 , res = 8.188477204912424e-07
it# 386 , res = 7.9501827271491e-07
it# 387 , res = 7.718822903450391e-07
it# 388 , res = 7.494195926596353e-07
it# 389 , res = 7.276105864511644e-07
it# 390 , res = 7.064362482381901e-07
it# 391 , res = 6.858781087378007e-07
it# 392 , res = 6.659182356139318e-07
it# 393 , res = 6.465392186637866e-07
it# 394 , res = 6.277241542488565e-07
it# 395 , res = 6.094566307476893e-07
it# 396 , res = 5.917207140002562e-07
it# 397 , res = 5.74500933723008e-07
it# 398 , res = 5.577822695529143e-07
it# 399 , res = 5.415501385570616e-07
it# 400 , res = 5.257903819689547e-07
it# 401 , res = 5.104892531140999e-07
it# 402 , res = 4.95633405376482e-07
it# 403 , res = 4.812098805625705e-07
it# 404 , res = 4.6720609759666646e-07
it# 405 , res = 4.536098415148426e-07
it# 406 , res = 4.404092528184786e-07
it# 407 , res = 4.275928168548029e-07
it# 408 , res = 4.151493545975895e-07
it# 409 , res = 4.0306801194597477e-07
it# 410 , res = 3.91338250820969e-07
it# 411 , res = 3.7994983971236213e-07
it# 412 , res = 3.6889284498964117e-07
it# 413 , res = 3.5815762227901083e-07
it# 414 , res = 3.477348070997161e-07
it# 415 , res = 3.3761530839279655e-07
it# 416 , res = 3.277902992728228e-07
it# 417 , res = 3.1825120954890035e-07
it# 418 , res = 3.08989718792838e-07
it# 419 , res = 2.9999774852556904e-07
it# 420 , res = 2.9126745527559135e-07
it# 421 , res = 2.8279122408491574e-07
it# 422 , res = 2.745616613206929e-07
it# 423 , res = 2.66571588613773e-07
it# 424 , res = 2.5881403661896573e-07
it# 425 , res = 2.5128223861513653e-07
it# 426 , res = 2.439696250761983e-07
it# 427 , res = 2.3686981720791803e-07
it# 428 , res = 2.2997662230087905e-07
it# 429 , res = 2.232840276673475e-07
it# 430 , res = 2.1678619540883574e-07
it# 431 , res = 2.1047745789177792e-07
it# 432 , res = 2.0435231215561485e-07
it# 433 , res = 1.98405415549691e-07
it# 434 , res = 1.9263158065462254e-07
it# 435 , res = 1.8702577138558587e-07
it# 436 , res = 1.8158309772489977e-07
it# 437 , res = 1.7629881250585378e-07
it# 438 , res = 1.711683061945983e-07
it# 439 , res = 1.661871038396135e-07
it# 440 , res = 1.6135086042137514e-07
it# 441 , res = 1.5665535747291945e-07
it# 442 , res = 1.5209649926273006e-07
it# 443 , res = 1.476703091932269e-07
it# 444 , res = 1.4337292662509268e-07
it# 445 , res = 1.3920060295234216e-07
it# 446 , res = 1.3514969899125614e-07
it# 447 , res = 1.312166811365594e-07
it# 448 , res = 1.2739811887139178e-07
it# 449 , res = 1.23690681219578e-07
it# 450 , res = 1.2009113452116902e-07
it# 451 , res = 1.1659633890291611e-07
it# 452 , res = 1.1320324604749718e-07
it# 453 , res = 1.099088961863434e-07
it# 454 , res = 1.0671041592109074e-07
it# 455 , res = 1.036050151795172e-07
it# 456 , res = 1.005899853350563e-07
it# 457 , res = 9.766269645445286e-08
it# 458 , res = 9.482059523508558e-08
it# 459 , res = 9.206120258082928e-08
it# 460 , res = 8.938211147267951e-08
it# 461 , res = 8.678098511026308e-08
it# 462 , res = 8.42555546547635e-08
it# 463 , res = 8.180361731745376e-08
it# 464 , res = 7.942303425487311e-08
it# 465 , res = 7.711172895061227e-08
it# 466 , res = 7.48676854651636e-08
it# 467 , res = 7.268894626927355e-08
it# 468 , res = 7.057361097309903e-08
it# 469 , res = 6.851983464747854e-08
it# 470 , res = 6.652582545455414e-08
it# 471 , res = 6.458984428669503e-08
it# 472 , res = 6.271020260347482e-08
it# 473 , res = 6.088526076940412e-08
it# 474 , res = 5.911342685509827e-08
it# 475 , res = 5.739315546375917e-08
it# 476 , res = 5.57229459733297e-08
it# 477 , res = 5.410134166149575e-08
it# 478 , res = 5.252692792369822e-08
it# 479 , res = 5.0998331498715014e-08
it# 480 , res = 4.9514218993658135e-08
it# 481 , res = 4.8073296108415345e-08
it# 482 , res = 4.6674305667397673e-08
it# 483 , res = 4.531602763329963e-08
it# 484 , res = 4.399727699118459e-08
it# 485 , res = 4.271690364457906e-08
it# 486 , res = 4.147379062463229e-08
it# 487 , res = 4.026685373335395e-08
it# 488 , res = 3.9095040158823025e-08
it# 489 , res = 3.795732770252726e-08
it# 490 , res = 3.68527240729449e-08
it# 491 , res = 3.5780265877967195e-08
it# 492 , res = 3.47390172709674e-08
it# 493 , res = 3.372807030178621e-08
it# 494 , res = 3.274654317672153e-08
it# 495 , res = 3.179357958055258e-08
it# 496 , res = 3.0868348366140785e-08
it# 497 , res = 2.9970042509991136e-08
it# 498 , res = 2.9097878488331975e-08
it# 499 , res = 2.8251095420155788e-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.09420141514100706
iteration 1 error = 0.12232793581312436
iteration 2 error = 0.08623239093524722
iteration 3 error = 0.05426269195969581
iteration 4 error = 0.027115353602136988
iteration 5 error = 0.014554766540388403
iteration 6 error = 0.006741759003015277
iteration 7 error = 0.0029652505195927746
iteration 8 error = 0.0010202948609567514
iteration 9 error = 0.0003933967952848954
iteration 10 error = 0.0001259416341172326
iteration 11 error = 5.0505266724961675e-05
iteration 12 error = 2.051021951612557e-05
iteration 13 error = 9.927367343961224e-06
iteration 14 error = 4.737718575764489e-06
iteration 15 error = 2.098844389517805e-06
iteration 16 error = 8.058310454478196e-07
iteration 17 error = 3.6891509958089896e-07
iteration 18 error = 1.2466359079563973e-07
iteration 19 error = 5.8855326711849246e-08
iteration 20 error = 1.9412108061419848e-08
iteration 21 error = 7.810537561963098e-09
iteration 22 error = 2.664515971735314e-09
iteration 23 error = 8.723735164966984e-10
iteration 24 error = 3.5611763235218816e-10
iteration 25 error = 9.29598896871641e-11
iteration 26 error = 3.653928697965554e-11
iteration 27 error = 1.4521570443047326e-11
iteration 28 error = 5.115685035026149e-12
iteration 29 error = 2.025751023436859e-12
iteration 30 error = 7.234574847827328e-13
iteration 31 error = 2.361855198521576e-13
iteration 32 error = 8.159369624610737e-14

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

20.22190236497056

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, 358, 359, 40, 935, 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, 356, 357, 932, 933}, {352, 353, 356, 357, 932, 934, 72, 361, 360}, {72, 158, 159, 160, 161, 866, 356, 357, 934, 358, 360, 361, 359, 40, 935, 364, 365, 938}, {547, 72, 74, 158, 159, 160, 161, 546, 867, 164, 165, 166, 167, 40, 41, 362, 363, 869, 358, 359, 935, 361, 360, 364, 365, 938, 940, 366, 367, 372, 373, 999}, {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, 586, 587, 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}, {540, 541, 544, 545, 546, 547, 548, 549, 932, 933, 934, 40, 938, 558, 559, 71, 72, 74, 350, 351, 992, 352, 353, 99, 356, 357, 998, 358, 360, 361, 359, 999, 364, 365, 366, 367}, {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, 584, 585, 76, 77, 78, 101, 103, 1000, 1018}, {576, 577, 580, 581, 582, 583, 584, 585, 588, 77, 78, 398, 399, 400, 402, 403, 401, 79, 406, 407, 408, 409, 589, 103, 1003, 46, 47, 947, 948, 950, 1018}, {404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 416, 417, 418, 419, 47, 48, 49, 949, 950, 951, 954, 582, 583, 584, 585, 586, 587, 588, 589, 78, 79, 81, 598, 599, 103, 1003, 1004}, {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}, {586, 587, 588, 589, 590, 79, 591, 81, 80, 594, 595, 1038, 598, 599, 414, 415, 416, 417, 418, 419, 422, 423, 424, 425, 103, 1004, 49, 50, 953, 954, 955}, {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, 772, 524, 525, 526, 527, 1044, 1045, 664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 678, 679, 773, 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}, {714, 1037, 1038, 715, 1040, 1046, 700, 701, 702, 703, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 78, 79, 77, 81, 590, 591, 80, 594, 598, 599, 595, 596, 597, 101, 103, 106, 1003, 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, 770, 771, 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, 128, 1023, 1055, 770, 771, 1067, 1069, 824, 825, 90, 91, 738, 739, 742, 743, 108, 113, 114, 124, 766, 767}, {768, 1025, 1026, 769, 772, 773, 774, 775, 789, 650, 651, 652, 653, 654, 655, 658, 659, 660, 661, 1042, 1044, 664, 665, 1045, 788, 668, 669, 670, 671, 672, 673, 1055, 1056, 91, 92, 94, 95, 113, 114, 117, 124, 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}, {772, 773, 774, 1031, 1032, 775, 778, 779, 780, 781, 788, 1045, 789, 666, 667, 668, 669, 94, 670, 96, 671, 1056, 676, 677, 678, 679, 1066, 114, 115, 117, 124}, {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, 128, 770, 772, 773, 774, 775, 771, 131, 778, 779, 780, 781, 782, 783, 767, 788, 789, 1055, 1056, 1066, 1069, 1072, 1075, 822, 823, 824, 825, 826, 827, 1086, 836, 837, 840, 841, 113, 114, 115, 117, 124, 766, 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, 768, 770, 771, 769, 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.949091191376844

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

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.992901
 0.997685
 0.999956
 1.34076
 1.80275
 1.84065
 1.96021
 1.98344
 1.99982

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.