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

2.1.2 Building blocks for programming preconditioners

In \(\S\)2.1.1, we saw the Jacobi method given as a local preconditioner. We now delve deeper into such local preconditioners, which are often useful ingredients/smoothers in more complicated preconditioning strategies. In this tutorial we will see

• Jacobi and Gauss-Seidel smoothers/preconditioners,

• their block versions,

• how to select your own smoothing blocks, and

• how to additively combine preconditioners.

from ngsolve import *
from ngsolve.webgui import Draw
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(grad(u)*grad(v)*dx + u*v*dx)
    f = LinearForm(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.0144918
 0.0565868
 0.0913563
 0.112048
 0.141222
 0.195761
 0.247646
 0.351665
 0.446235
 0.524795
 0.632006
 0.756244
 0.835742
 0.927184
 1.03588
 1.16687
 1.25179
 1.36664
 1.47109
 1.57336
 1.68741
 1.81492
 1.93739
 2.01627
 2.13865
 2.22417
 2.31655
 2.37388
  2.4388
 2.49019
 2.50489
 2.53176
 2.80787

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

max(lams)/min(lams)

193.75502520680772

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)

1561.8273803217185

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)

CG iteration 1, residual = 0.05514173105997197
CG iteration 2, residual = 0.09362431014556011
CG iteration 3, residual = 0.10183798655456601
CG iteration 4, residual = 0.08946401515399124
CG iteration 5, residual = 0.09864154662595712
CG iteration 6, residual = 0.08752084132820948
CG iteration 7, residual = 0.09196141307218668
CG iteration 8, residual = 0.07931974777210128
CG iteration 9, residual = 0.05262078540101382
CG iteration 10, residual = 0.04323706125676102
CG iteration 11, residual = 0.035206451492069774
CG iteration 12, residual = 0.027795575409704183
CG iteration 13, residual = 0.022311587802186686
CG iteration 14, residual = 0.019516985933755487
CG iteration 15, residual = 0.013880031216222442
CG iteration 16, residual = 0.0108194142234518
CG iteration 17, residual = 0.009432534548762434
CG iteration 18, residual = 0.00701882909110787
CG iteration 19, residual = 0.004666218187002316
CG iteration 20, residual = 0.003164824338234743
CG iteration 21, residual = 0.0017730431581364215
CG iteration 22, residual = 0.0011842359412064284
CG iteration 23, residual = 0.0010076069006142196
CG iteration 24, residual = 0.0007466659458352507
CG iteration 25, residual = 0.0004697350490174966
CG iteration 26, residual = 0.0003348113460042972
CG iteration 27, residual = 0.0002829748057927384
CG iteration 28, residual = 0.00021828093905188454
CG iteration 29, residual = 0.00015449917737738025
CG iteration 30, residual = 0.00011300774305439106
CG iteration 31, residual = 8.206625685090374e-05
CG iteration 32, residual = 6.899528154493053e-05
CG iteration 33, residual = 5.795236525379512e-05
CG iteration 34, residual = 4.3057180145713543e-05
CG iteration 35, residual = 3.120366863521651e-05
CG iteration 36, residual = 2.345226320419527e-05
CG iteration 37, residual = 1.7861388530143805e-05
CG iteration 38, residual = 1.416656738033425e-05
CG iteration 39, residual = 1.0900081367054178e-05
CG iteration 40, residual = 7.9139094892054e-06
CG iteration 41, residual = 5.132078419950257e-06
CG iteration 42, residual = 3.953008836659825e-06
CG iteration 43, residual = 2.781054144862857e-06
CG iteration 44, residual = 1.977974718431345e-06
CG iteration 45, residual = 1.3918628904040142e-06
CG iteration 46, residual = 9.50396087518867e-07
CG iteration 47, residual = 6.381955827846508e-07
CG iteration 48, residual = 4.5071157676475914e-07
CG iteration 49, residual = 3.2397903934175013e-07
CG iteration 50, residual = 1.9978680980298134e-07
CG iteration 51, residual = 1.378203022117175e-07
CG iteration 52, residual = 1.0164595645020918e-07
CG iteration 53, residual = 7.211289007062626e-08
CG iteration 54, residual = 5.185357563205716e-08
CG iteration 55, residual = 3.731338299178322e-08
CG iteration 56, residual = 2.7804750047049427e-08
CG iteration 57, residual = 1.8525785318317213e-08
CG iteration 58, residual = 1.1844642126790734e-08
CG iteration 59, residual = 7.155287249925312e-09
CG iteration 60, residual = 4.760016905917839e-09
CG iteration 61, residual = 3.630196148823671e-09
CG iteration 62, residual = 2.627479145645254e-09
CG iteration 63, residual = 1.7306647971566386e-09
CG iteration 64, residual = 1.0486752192572573e-09
CG iteration 65, residual = 7.032650198799314e-10
CG iteration 66, residual = 5.364281702399671e-10
CG iteration 67, residual = 4.198031896535958e-10
CG iteration 68, residual = 2.910771579991341e-10
CG iteration 69, residual = 1.7349690893068683e-10
CG iteration 70, residual = 1.0534699620096689e-10
CG iteration 71, residual = 7.302316532601926e-11
CG iteration 72, residual = 5.603304209058976e-11
CG iteration 73, residual = 4.089963029578208e-11
CG iteration 74, residual = 2.4897862393951943e-11
CG iteration 75, residual = 1.4508864644591884e-11
CG iteration 76, residual = 1.0177902501246067e-11
CG iteration 77, residual = 7.705336097670046e-12
CG iteration 78, residual = 5.342321883071467e-12
CG iteration 79, residual = 3.49502049072479e-12
CG iteration 80, residual = 2.0168020920534605e-12
CG iteration 81, residual = 1.2027998593776186e-12
CG iteration 82, residual = 8.409546049149749e-13
CG iteration 83, residual = 6.23777637864478e-13
CG iteration 84, residual = 4.2057902491579505e-13
CG iteration 85, residual = 2.698437985703656e-13
CG iteration 86, residual = 1.7604807657639508e-13
CG iteration 87, residual = 1.1316954419940909e-13
CG iteration 88, residual = 7.890045691006462e-14
CG iteration 89, residual = 5.31393441329676e-14

WebGLScene

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(). It is 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.08154902419641553
it# 1 , res = 0.07698831978567876
it# 2 , res = 0.07343571088587378
it# 3 , res = 0.07022987452159993
it# 4 , res = 0.06734784593987088
it# 5 , res = 0.06472557072763736
it# 6 , res = 0.06231335823672829
it# 7 , res = 0.060074879731911224
it# 8 , res = 0.05798285966901513
it# 9 , res = 0.056016316674589804
it# 10 , res = 0.05415877995370276
it# 11 , res = 0.05239709685483901
it# 12 , res = 0.05072061170350224
it# 13 , res = 0.04912058411072688
it# 14 , res = 0.047589766385055604
it# 15 , res = 0.04612208963244292
it# 16 , res = 0.044712426078080604
it# 17 , res = 0.04335640617752859
it# 18 , res = 0.04205027600466386
it# 19 , res = 0.040790784835190934
it# 20 , res = 0.03957509574733402
it# 21 , res = 0.03840071401068155
it# 22 , res = 0.037265429377108035
it# 23 , res = 0.03616726933615246
it# 24 , res = 0.035104461082522555
it# 25 , res = 0.03407540044859191
it# 26 , res = 0.033078626433638486
it# 27 , res = 0.03211280024986055
it# 28 , res = 0.031176688027209725
it# 29 , res = 0.030269146491763342
it# 30 , res = 0.029389111067791993
it# 31 , res = 0.028535585960640954
it# 32 , res = 0.027707635862510547
it# 33 , res = 0.026904378991044666
it# 34 , res = 0.026124981225003154
it# 35 , res = 0.025368651145027645
it# 36 , res = 0.024634635822794893
it# 37 , res = 0.023922217230401895
it# 38 , res = 0.023230709164974103
it# 39 , res = 0.022559454602299422
it# 40 , res = 0.021907823408608803
it# 41 , res = 0.02127521035211056
it# 42 , res = 0.020661033366095424
it# 43 , res = 0.020064732023773787
it# 44 , res = 0.01948576619185389
it# 45 , res = 0.018923614835481054
it# 46 , res = 0.018377774951776855
it# 47 , res = 0.01784776061301663
it# 48 , res = 0.017333102103614936
it# 49 , res = 0.01683334513767716
it# 50 , res = 0.016348050146014176
it# 51 , res = 0.01587679162328914
it# 52 , res = 0.015419157527437848
it# 53 , res = 0.01497474872472229
it# 54 , res = 0.014543178474801948
it# 55 , res = 0.014124071951048874
it# 56 , res = 0.013717065792048452
it# 57 , res = 0.013321807680817812
it# 58 , res = 0.012937955948776165
it# 59 , res = 0.012565179201917739
it# 60 , res = 0.012203155966991142
it# 61 , res = 0.01185157435579114
it# 62 , res = 0.011510131745914623
it# 63 , res = 0.011178534476553453
it# 64 , res = 0.010856497558072513
it# 65 , res = 0.010543744394280822
it# 66 , res = 0.010240006516432833
it# 67 , res = 0.009945023328111762
it# 68 , res = 0.009658541860242528
it# 69 , res = 0.009380316535567288
it# 70 , res = 0.009110108941986281
it# 71 , res = 0.008847687614229343
it# 72 , res = 0.008592827823380405
it# 73 , res = 0.008345311373817485
it# 74 , res = 0.00810492640717913
it# 75 , res = 0.00787146721299802
it# 76 , res = 0.007644734045675342
it# 77 , res = 0.007424532947499533
it# 78 , res = 0.007210675577431656
it# 79 , res = 0.007002979045406306
it# 80 , res = 0.006801265751912552
it# 81 , res = 0.006605363232636963
it# 82 , res = 0.006415104007966018
it# 83 , res = 0.006230325437158684
it# 84 , res = 0.006050869577009658
it# 85 , res = 0.005876583044839283
it# 86 , res = 0.005707316885649518
it# 87 , res = 0.005542926443300168
it# 88 , res = 0.005383271235562786
it# 89 , res = 0.005228214832920553
it# 90 , res = 0.005077624740985864
it# 91 , res = 0.004931372286416426
it# 92 , res = 0.004789332506213437
it# 93 , res = 0.004651384040293219
it# 94 , res = 0.004517409027226086
it# 95 , res = 0.0043872930030428095
it# 96 , res = 0.004260924803010143
it# 97 , res = 0.004138196466285642
it# 98 , res = 0.004019003143359983
it# 99 , res = 0.003903243006201826
it# 100 , res = 0.0037908171610226615
it# 101 , res = 0.00368162956358284
it# 102 , res = 0.0035755869369603617
it# 103 , res = 0.0034725986917093217
it# 104 , res = 0.0033725768483364436
it# 105 , res = 0.0032754359620263823
it# 106 , res = 0.003181093049548214
it# 107 , res = 0.0030894675182804514
it# 108 , res = 0.003000481097289107
it# 109 , res = 0.0029140577704012166
it# 110 , res = 0.002830123711212951
it# 111 , res = 0.002748607219976816
it# 112 , res = 0.002669438662312448
it# 113 , res = 0.0025925504096869245
it# 114 , res = 0.0025178767816158233
it# 115 , res = 0.0024453539895301123
it# 116 , res = 0.0023749200822653473
it# 117 , res = 0.0023065148931222033
it# 118 , res = 0.002240079988455254
it# 119 , res = 0.0021755586177435978
it# 120 , res = 0.0021128956651026596
it# 121 , res = 0.002052037602192795
it# 122 , res = 0.0019929324424863516
it# 123 , res = 0.0019355296968528689
it# 124 , res = 0.0018797803304244933
it# 125 , res = 0.001825636720705095
it# 126 , res = 0.0017730526168851677
it# 127 , res = 0.0017219831003313932
it# 128 , res = 0.0016723845462129087
it# 129 , res = 0.0016242145862339235
it# 130 , res = 0.0015774320724402159
it# 131 , res = 0.0015319970420685778
it# 132 , res = 0.0014878706834079106
it# 133 , res = 0.0014450153026454877
it# 134 , res = 0.0014033942916666334
it# 135 , res = 0.0013629720967835567
it# 136 , res = 0.0013237141883635087
it# 137 , res = 0.0012855870313333638
it# 138 , res = 0.0012485580565323208
it# 139 , res = 0.0012125956328912278
it# 140 , res = 0.001177669040412115
it# 141 , res = 0.0011437484439267946
it# 142 , res = 0.0011108048676108143
it# 143 , res = 0.0010788101702321788
it# 144 , res = 0.0010477370211121422
it# 145 , res = 0.0010175588767789867
it# 146 , res = 0.000988249958294322
it# 147 , res = 0.0009597852292320095
it# 148 , res = 0.0009321403742918709
it# 149 , res = 0.000905291778528643
it# 150 , res = 0.0008792165071802817
it# 151 , res = 0.0008538922860765539
it# 152 , res = 0.0008292974826116835
it# 153 , res = 0.0008054110872656734
it# 154 , res = 0.0007822126956577888
it# 155 , res = 0.0007596824911162986
it# 156 , res = 0.0007378012277509575
it# 157 , res = 0.0007165502140129796
it# 158 , res = 0.000695911296727994
it# 159 , res = 0.0006758668455896285
it# 160 , res = 0.000656399738098956
it# 161 , res = 0.0006374933449383716
it# 162 , res = 0.0006191315157664529
it# 163 , res = 0.0006012985654220138
it# 164 , res = 0.0005839792605254928
it# 165 , res = 0.0005671588064660554
it# 166 , res = 0.0005508228347643104
it# 167 , res = 0.0005349573907977138
it# 168 , res = 0.0005195489218806117
it# 169 , res = 0.000504584265687435
it# 170 , res = 0.0004900506390086511
it# 171 , res = 0.000475935626831143
it# 172 , res = 0.00046222717173357994
it# 173 , res = 0.0004489135635859625
it# 174 , res = 0.0004359834295468412
it# 175 , res = 0.0004234257243487381
it# 176 , res = 0.00041122972086220873
it# 177 , res = 0.00039938500093305185
it# 178 , res = 0.000387881446482983
it# 179 , res = 0.0003767092308660087
it# 180 , res = 0.00036585881047451155
it# 181 , res = 0.00035532091658712733
it# 182 , res = 0.0003450865474506564
it# 183 , res = 0.0003351469605908471
it# 184 , res = 0.00032549366534458633
it# 185 , res = 0.00031611841560624923
it# 186 , res = 0.0003070132027845393
it# 187 , res = 0.00029817024896050824
it# 188 , res = 0.00028958200024407465
it# 189 , res = 0.000281241120321089
it# 190 , res = 0.00027314048418644967
it# 191 , res = 0.00026527317205760363
it# 192 , res = 0.0002576324634637743
it# 193 , res = 0.0002502118315049194
it# 194 , res = 0.00024300493727648365
it# 195 , res = 0.00023600562445420966
it# 196 , res = 0.00022920791403576346
it# 197 , res = 0.00022260599923271026
it# 198 , res = 0.0002161942405107193
it# 199 , res = 0.0002099671607718633
it# 200 , res = 0.0002039194406758857
it# 201 , res = 0.0001980459140964335
it# 202 , res = 0.00019234156370812223
it# 203 , res = 0.00018680151670017608
it# 204 , res = 0.00018142104061430646
it# 205 , res = 0.00017619553930205915
it# 206 , res = 0.00017112054899858179
it# 207 , res = 0.0001661917345095037
it# 208 , res = 0.00016140488550806106
it# 209 , res = 0.0001567559129379277
it# 210 , res = 0.0001522408455207379
it# 211 , res = 0.0001478558263635837
it# 212 , res = 0.00014359710966402791
it# 213 , res = 0.0001394610575109252
it# 214 , res = 0.0001354441367763822
it# 215 , res = 0.00013154291609766453
it# 216 , res = 0.00012775406294641866
it# 217 , res = 0.0001240743407814057
it# 218 , res = 0.00012050060628416776
it# 219 , res = 0.00011702980667384177
it# 220 , res = 0.00011365897709931606
it# 221 , res = 0.00011038523810645488
it# 222 , res = 0.00010720579317858816
it# 223 , res = 0.00010411792634770729
it# 224 , res = 0.00010111899987406696
it# 225 , res = 9.820645199350301e-05
it# 226 , res = 9.537779472854354e-05
it# 227 , res = 9.263061176338161e-05
it# 228 , res = 8.996255637981363e-05
it# 229 , res = 8.737134945277001e-05
it# 230 , res = 8.485477750259047e-05
it# 231 , res = 8.241069080561367e-05
it# 232 , res = 8.003700155650622e-05
it# 233 , res = 7.773168208549986e-05
it# 234 , res = 7.549276312629825e-05
it# 235 , res = 7.331833213349592e-05
it# 236 , res = 7.120653164911058e-05
it# 237 , res = 6.915555771595481e-05
it# 238 , res = 6.716365833647299e-05
it# 239 , res = 6.522913197610796e-05
it# 240 , res = 6.335032610998995e-05
it# 241 , res = 6.152563581113964e-05
it# 242 , res = 5.975350237961398e-05
it# 243 , res = 5.803241201112682e-05
it# 244 , res = 5.636089450338161e-05
it# 245 , res = 5.473752200068699e-05
it# 246 , res = 5.316090777452601e-05
it# 247 , res = 5.162970503797986e-05
it# 248 , res = 5.014260579640749e-05
it# 249 , res = 4.869833972882769e-05
it# 250 , res = 4.7295673104551094e-05
it# 251 , res = 4.5933407727390914e-05
it# 252 , res = 4.4610379913235245e-05
it# 253 , res = 4.3325459495958504e-05
it# 254 , res = 4.2077548861625546e-05
it# 255 , res = 4.0865582011135517e-05
it# 256 , res = 3.968852364959525e-05
it# 257 , res = 3.854536830170849e-05
it# 258 , res = 3.743513945320208e-05
it# 259 , res = 3.6356888716436335e-05
it# 260 , res = 3.530969502041789e-05
it# 261 , res = 3.429266382389116e-05
it# 262 , res = 3.330492635116165e-05
it# 263 , res = 3.23456388502091e-05
it# 264 , res = 3.141398187156026e-05
it# 265 , res = 3.050915956835785e-05
it# 266 , res = 2.963039901714786e-05
it# 267 , res = 2.8776949556656693e-05
it# 268 , res = 2.794808214719798e-05
it# 269 , res = 2.7143088747695537e-05
it# 270 , res = 2.63612817112007e-05
it# 271 , res = 2.5601993196551614e-05
it# 272 , res = 2.486457459929837e-05
it# 273 , res = 2.414839599626255e-05
it# 274 , res = 2.3452845608110377e-05
it# 275 , res = 2.277732927693732e-05
it# 276 , res = 2.2121269958588006e-05
it# 277 , res = 2.1484107229194958e-05
it# 278 , res = 2.086529680704279e-05
it# 279 , res = 2.0264310087514962e-05
it# 280 , res = 1.968063369136689e-05
it# 281 , res = 1.9113769026571297e-05
it# 282 , res = 1.8563231861903217e-05
it# 283 , res = 1.8028551913460823e-05
it# 284 , res = 1.7509272443282183e-05
it# 285 , res = 1.7004949868614867e-05
it# 286 , res = 1.6515153383426352e-05
it# 287 , res = 1.6039464590477228e-05
it# 288 , res = 1.5577477143153493e-05
it# 289 , res = 1.5128796399608339e-05
it# 290 , res = 1.4693039084363772e-05
it# 291 , res = 1.4269832961791608e-05
it# 292 , res = 1.3858816517546784e-05
it# 293 , res = 1.3459638650571483e-05
it# 294 , res = 1.307195837190308e-05
it# 295 , res = 1.2695444514873804e-05
it# 296 , res = 1.2329775450876588e-05
it# 297 , res = 1.1974638815730643e-05
it# 298 , res = 1.162973124187441e-05
it# 299 , res = 1.12947580999763e-05
it# 300 , res = 1.0969433246774319e-05
it# 301 , res = 1.0653478781253917e-05
it# 302 , res = 1.0346624806273206e-05
it# 303 , res = 1.0048609198932356e-05
it# 304 , res = 9.759177386124823e-06
it# 305 , res = 9.478082127537557e-06
it# 306 , res = 9.205083303797985e-06
it# 307 , res = 8.939947711952543e-06
it# 308 , res = 8.682448866040939e-06
it# 309 , res = 8.432366803476527e-06
it# 310 , res = 8.189487897365318e-06
it# 311 , res = 7.953604674129256e-06
it# 312 , res = 7.724515635724234e-06
it# 313 , res = 7.502025088163965e-06
it# 314 , res = 7.285942973898772e-06
it# 315 , res = 7.07608470993066e-06
it# 316 , res = 6.872271029413329e-06
it# 317 , res = 6.674327829294091e-06
it# 318 , res = 6.482086020993533e-06
it# 319 , res = 6.295381386556745e-06
it# 320 , res = 6.1140544376999e-06
it# 321 , res = 5.937950279905525e-06
it# 322 , res = 5.766918480499939e-06
it# 323 , res = 5.60081293924513e-06
it# 324 , res = 5.439491764420611e-06
it# 325 , res = 5.282817151021825e-06
it# 326 , res = 5.130655263257956e-06
it# 327 , res = 4.982876120420412e-06
it# 328 , res = 4.839353485659886e-06
it# 329 , res = 4.6999647579132005e-06
it# 330 , res = 4.564590867607231e-06
it# 331 , res = 4.433116174644921e-06
it# 332 , res = 4.305428369699458e-06
it# 333 , res = 4.181418378623103e-06
it# 334 , res = 4.0609802685233475e-06
it# 335 , res = 3.944011158031917e-06
it# 336 , res = 3.830411128673361e-06
it# 337 , res = 3.72008314056042e-06
it# 338 , res = 3.612932948353923e-06
it# 339 , res = 3.508869021553855e-06
it# 340 , res = 3.4078024657587257e-06
it# 341 , res = 3.3096469474056794e-06
it# 342 , res = 3.2143186192523254e-06
it# 343 , res = 3.121736049134269e-06
it# 344 , res = 3.0318201507800805e-06
it# 345 , res = 2.9444941154522837e-06
it# 346 , res = 2.8596833468532936e-06
it# 347 , res = 2.7773153973529115e-06
it# 348 , res = 2.6973199061001595e-06
it# 349 , res = 2.6196285385272846e-06
it# 350 , res = 2.544174928878084e-06
it# 351 , res = 2.470894622396847e-06
it# 352 , res = 2.3997250211568448e-06
it# 353 , res = 2.3306053301730986e-06
it# 354 , res = 2.263476505417967e-06
it# 355 , res = 2.198281203733331e-06
it# 356 , res = 2.1349637335319523e-06
it# 357 , res = 2.073470007259759e-06
it# 358 , res = 2.0137474953240006e-06
it# 359 , res = 1.9557451811584393e-06
it# 360 , res = 1.8994135175146881e-06
it# 361 , res = 1.8447043844505433e-06
it# 362 , res = 1.7915710479349495e-06
it# 363 , res = 1.73996812000261e-06
it# 364 , res = 1.6898515200674736e-06
it# 365 , res = 1.6411784372189632e-06
it# 366 , res = 1.5939072933705362e-06
it# 367 , res = 1.5479977085433242e-06
it# 368 , res = 1.5034104653202554e-06
it# 369 , res = 1.4601074759521327e-06
it# 370 , res = 1.418051750179263e-06
it# 371 , res = 1.3772073625921952e-06
it# 372 , res = 1.3375394229042957e-06
it# 373 , res = 1.2990140456352713e-06
it# 374 , res = 1.2615983215689688e-06
it# 375 , res = 1.2252602888234533e-06
it# 376 , res = 1.189968906786175e-06
it# 377 , res = 1.1556940284814694e-06
it# 378 , res = 1.1224063753201257e-06
it# 379 , res = 1.0900775122469149e-06
it# 380 , res = 1.0586798227662105e-06
it# 381 , res = 1.0281864864032394e-06
it# 382 , res = 9.985714548756043e-07
it# 383 , res = 9.698094301513545e-07
it# 384 , res = 9.418758428638932e-07
it# 385 , res = 9.147468315834136e-07
it# 386 , res = 8.883992218533449e-07
it# 387 , res = 8.628105068173562e-07
it# 388 , res = 8.379588278829127e-07
it# 389 , res = 8.138229561761493e-07
it# 390 , res = 7.903822741126206e-07
it# 391 , res = 7.67616758021396e-07
it# 392 , res = 7.455069610624926e-07
it# 393 , res = 7.240339963435305e-07
it# 394 , res = 7.031795210377193e-07
it# 395 , res = 6.829257208162232e-07
it# 396 , res = 6.632552942086831e-07
it# 397 , res = 6.441514383504427e-07
it# 398 , res = 6.255978340289856e-07
it# 399 , res = 6.075786324828107e-07
it# 400 , res = 5.900784408825796e-07
it# 401 , res = 5.730823104701048e-07
it# 402 , res = 5.565757225576183e-07
it# 403 , res = 5.405445766353661e-07
it# 404 , res = 5.249751785330479e-07
it# 405 , res = 5.098542285134347e-07
it# 406 , res = 4.951688097972318e-07
it# 407 , res = 4.809063776786613e-07
it# 408 , res = 4.670547489215618e-07
it# 409 , res = 4.5360209090883067e-07
it# 410 , res = 4.405369121622175e-07
it# 411 , res = 4.2784805197373333e-07
it# 412 , res = 4.1552467128732793e-07
it# 413 , res = 4.0355624291314913e-07
it# 414 , res = 3.9193254333167524e-07
it# 415 , res = 3.806436431084319e-07
it# 416 , res = 3.6967989903746164e-07
it# 417 , res = 3.5903194551811313e-07
it# 418 , res = 3.4869068688601037e-07
it# 419 , res = 3.3864728924563425e-07
it# 420 , res = 3.28893173308512e-07
it# 421 , res = 3.1942000696038433e-07
it# 422 , res = 3.1021969770142184e-07
it# 423 , res = 3.012843865375272e-07
it# 424 , res = 2.926064407798554e-07
it# 425 , res = 2.84178447397882e-07
it# 426 , res = 2.759932068587635e-07
it# 427 , res = 2.6804372724914127e-07
it# 428 , res = 2.6032321776647044e-07
it# 429 , res = 2.5282508354480654e-07
it# 430 , res = 2.4554291932050166e-07
it# 431 , res = 2.384705045718945e-07
it# 432 , res = 2.3160179769254755e-07
it# 433 , res = 2.249309313956211e-07
it# 434 , res = 2.1845220726953672e-07
it# 435 , res = 2.1216009094746265e-07
it# 436 , res = 2.0604920751978884e-07
it# 437 , res = 2.0011433686775713e-07
it# 438 , res = 1.9435040934648096e-07
it# 439 , res = 1.8875250117010408e-07
it# 440 , res = 1.8331583049674858e-07
it# 441 , res = 1.7803575328498423e-07
it# 442 , res = 1.7290775892230725e-07
it# 443 , res = 1.679274671100223e-07
it# 444 , res = 1.63090623468206e-07
it# 445 , res = 1.5839309635601212e-07
it# 446 , res = 1.538308727869663e-07
it# 447 , res = 1.4940005581631076e-07
it# 448 , res = 1.4509686047357325e-07
it# 449 , res = 1.4091761067543774e-07
it# 450 , res = 1.3685873666760947e-07
it# 451 , res = 1.32916770936485e-07
it# 452 , res = 1.290883464343434e-07
it# 453 , res = 1.2537019275174507e-07
it# 454 , res = 1.2175913346189742e-07
it# 455 , res = 1.1825208438579894e-07
it# 456 , res = 1.1484604930311197e-07
it# 457 , res = 1.1153811889608629e-07
it# 458 , res = 1.0832546730243801e-07
it# 459 , res = 1.0520535029686134e-07
it# 460 , res = 1.0217510256359273e-07
it# 461 , res = 9.923213567106059e-08
it# 462 , res = 9.637393543659313e-08
it# 463 , res = 9.359806039407237e-08
it# 464 , res = 9.090213940937919e-08
it# 465 , res = 8.828386956686662e-08
it# 466 , res = 8.574101414547337e-08
it# 467 , res = 8.327140103116006e-08
it# 468 , res = 8.087292047905739e-08
it# 469 , res = 7.85435240012563e-08
it# 470 , res = 7.62812213554447e-08
it# 471 , res = 7.408408030466951e-08
it# 472 , res = 7.195022384530282e-08
it# 473 , res = 6.987782915735972e-08
it# 474 , res = 6.786512611295113e-08
it# 475 , res = 6.591039524280722e-08
it# 476 , res = 6.401196689041484e-08
it# 477 , res = 6.216821915649564e-08
it# 478 , res = 6.037757731894831e-08
it# 479 , res = 5.863851167680967e-08
it# 480 , res = 5.694953649515649e-08
it# 481 , res = 5.5309209178042494e-08
it# 482 , res = 5.3716128635881e-08
it# 483 , res = 5.216893383954753e-08
it# 484 , res = 5.066630304735927e-08
it# 485 , res = 4.9206952876958935e-08
it# 486 , res = 4.7789636478784324e-08
it# 487 , res = 4.641314344576823e-08
it# 488 , res = 4.507629772423715e-08
it# 489 , res = 4.377795745146812e-08
it# 490 , res = 4.2517013318117004e-08
it# 491 , res = 4.1292388484611016e-08
it# 492 , res = 4.0103036864817407e-08
it# 493 , res = 3.894794214905974e-08
it# 494 , res = 3.782611791405422e-08
it# 495 , res = 3.673660576630237e-08
it# 496 , res = 3.567847494563905e-08
it# 497 , res = 3.465082183630406e-08
it# 498 , res = 3.365276823877756e-08
it# 499 , res = 3.2683461797417596e-08

WebGLScene

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)

CG iteration 1, residual = 0.09414079549704508
CG iteration 2, residual = 0.12203689059859807
CG iteration 3, residual = 0.08630349626448042
CG iteration 4, residual = 0.05476506137338164
CG iteration 5, residual = 0.027039871611641356
CG iteration 6, residual = 0.014941028668260052
CG iteration 7, residual = 0.0066946309143105535
CG iteration 8, residual = 0.00284839806274339
CG iteration 9, residual = 0.0009877015081514713
CG iteration 10, residual = 0.0003968124927947862
CG iteration 11, residual = 0.0001240754457456034
CG iteration 12, residual = 6.330209216165876e-05
CG iteration 13, residual = 2.9785257662869855e-05
CG iteration 14, residual = 1.3390306292216037e-05
CG iteration 15, residual = 5.363365741544913e-06
CG iteration 16, residual = 2.0462661924717113e-06
CG iteration 17, residual = 8.231874195477868e-07
CG iteration 18, residual = 3.4544567760902384e-07
CG iteration 19, residual = 1.295123016019344e-07
CG iteration 20, residual = 5.3553595524320744e-08
CG iteration 21, residual = 1.7611013168353518e-08
CG iteration 22, residual = 5.970141704993225e-09
CG iteration 23, residual = 2.4674252558032007e-09
CG iteration 24, residual = 7.865934453395216e-10
CG iteration 25, residual = 3.1722300643176425e-10
CG iteration 26, residual = 1.1006534206433601e-10
CG iteration 27, residual = 4.095450344654651e-11
CG iteration 28, residual = 1.5996345010343813e-11
CG iteration 29, residual = 5.326653784450914e-12
CG iteration 30, residual = 1.9887046841316273e-12
CG iteration 31, residual = 5.843515407392285e-13
CG iteration 32, residual = 2.267856969958932e-13
CG iteration 33, residual = 6.065405903933637e-14

WebGLScene

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

20.311617625163535

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 \(\S\)1.8.

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

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

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

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

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

35.05232462837895

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

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.993158
 0.997949
 0.999962
 1.35426
 1.82282
 1.90685
 1.96916
 1.98517
  1.9997

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)
Draw(gfu)
mg = Preconditioner(a, 'multigrid')  # register mg to a
a.Assemble()                         # assemble on coarsest mesh

<ngsolve.comp.BilinearForm at 0x7fa59ee8fab0>

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.826237
 0.849126
 0.897069
 0.933508
 1.00501
 1.11965
 1.36946
 1.43186
 1.54355
 1.65858
 1.76409
 1.85156
 1.90736
 1.97186
 1.99987

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