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.10183798655456605
CG iteration 4, residual = 0.08946401515399127
CG iteration 5, residual = 0.09864154662595712
CG iteration 6, residual = 0.0875208413282095
CG iteration 7, residual = 0.09196141307218676
CG iteration 8, residual = 0.07931974777210134
CG iteration 9, residual = 0.05262078540101384
CG iteration 10, residual = 0.04323706125676106
CG iteration 11, residual = 0.035206451492069794
CG iteration 12, residual = 0.02779557540970419
CG iteration 13, residual = 0.022311587802186693
CG iteration 14, residual = 0.01951698593375548
CG iteration 15, residual = 0.013880031216222435
CG iteration 16, residual = 0.010819414223451791
CG iteration 17, residual = 0.00943253454876243
CG iteration 18, residual = 0.007018829091107869
CG iteration 19, residual = 0.004666218187002315
CG iteration 20, residual = 0.0031648243382347433
CG iteration 21, residual = 0.001773043158136422
CG iteration 22, residual = 0.0011842359412064286
CG iteration 23, residual = 0.0010076069006142196
CG iteration 24, residual = 0.0007466659458352504
CG iteration 25, residual = 0.0004697350490174962
CG iteration 26, residual = 0.0003348113460042969
CG iteration 27, residual = 0.000282974805792738
CG iteration 28, residual = 0.0002182809390518843
CG iteration 29, residual = 0.0001544991773773801
CG iteration 30, residual = 0.00011300774305439091
CG iteration 31, residual = 8.206625685090361e-05
CG iteration 32, residual = 6.899528154493041e-05
CG iteration 33, residual = 5.7952365253795013e-05
CG iteration 34, residual = 4.305718014571346e-05
CG iteration 35, residual = 3.1203668635216434e-05
CG iteration 36, residual = 2.345226320419521e-05
CG iteration 37, residual = 1.786138853014374e-05
CG iteration 38, residual = 1.4166567380334187e-05
CG iteration 39, residual = 1.0900081367054066e-05
CG iteration 40, residual = 7.9139094892051e-06
CG iteration 41, residual = 5.132078419949403e-06
CG iteration 42, residual = 3.953008836657823e-06
CG iteration 43, residual = 2.7810541448574595e-06
CG iteration 44, residual = 1.9779747184157013e-06
CG iteration 45, residual = 1.391862890344353e-06
CG iteration 46, residual = 9.503960873233038e-07
CG iteration 47, residual = 6.381955822738259e-07
CG iteration 48, residual = 4.507115754128724e-07
CG iteration 49, residual = 3.239790356686484e-07
CG iteration 50, residual = 1.9978679812336823e-07
CG iteration 51, residual = 1.378202768364813e-07
CG iteration 52, residual = 1.0164590077440121e-07
CG iteration 53, residual = 7.211270765536797e-08
CG iteration 54, residual = 5.1852818249002047e-08
CG iteration 55, residual = 3.731126990468823e-08
CG iteration 56, residual = 2.7798775363977926e-08
CG iteration 57, residual = 1.8510038987759495e-08
CG iteration 58, residual = 1.1791679130716936e-08
CG iteration 59, residual = 7.036360860955813e-09
CG iteration 60, residual = 4.593590997803127e-09
CG iteration 61, residual = 3.502563765507461e-09
CG iteration 62, residual = 2.600049741741244e-09
CG iteration 63, residual = 1.7207356307939052e-09
CG iteration 64, residual = 1.052691683948583e-09
CG iteration 65, residual = 7.041916808499999e-10
CG iteration 66, residual = 5.366191436277449e-10
CG iteration 67, residual = 4.1984587853582006e-10
CG iteration 68, residual = 2.9108052413149205e-10
CG iteration 69, residual = 1.734972964160835e-10
CG iteration 70, residual = 1.0534704999457344e-10
CG iteration 71, residual = 7.302318446738366e-11
CG iteration 72, residual = 5.6033046151164665e-11
CG iteration 73, residual = 4.089963069254367e-11
CG iteration 74, residual = 2.4897862445365535e-11
CG iteration 75, residual = 1.4508864653548972e-11
CG iteration 76, residual = 1.0177902503131967e-11
CG iteration 77, residual = 7.705336098080402e-12
CG iteration 78, residual = 5.342321883115746e-12
CG iteration 79, residual = 3.4950204907252572e-12
CG iteration 80, residual = 2.0168020920500237e-12
CG iteration 81, residual = 1.2027998593712077e-12
CG iteration 82, residual = 8.409546049045092e-13
CG iteration 83, residual = 6.237776378490207e-13
CG iteration 84, residual = 4.205790248954814e-13
CG iteration 85, residual = 2.698437985431949e-13
CG iteration 86, residual = 1.7604807653796455e-13
CG iteration 87, residual = 1.1316954414469497e-13
CG iteration 88, residual = 7.89004568294809e-14
CG iteration 89, residual = 5.313934402099496e-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.07343571088587376
it# 3 , res = 0.07022987452159993
it# 4 , res = 0.06734784593987088
it# 5 , res = 0.06472557072763736
it# 6 , res = 0.06231335823672828
it# 7 , res = 0.06007487973191122
it# 8 , res = 0.05798285966901515
it# 9 , res = 0.0560163166745898
it# 10 , res = 0.05415877995370277
it# 11 , res = 0.052397096854838995
it# 12 , res = 0.05072061170350224
it# 13 , res = 0.04912058411072687
it# 14 , res = 0.04758976638505559
it# 15 , res = 0.046122089632442874
it# 16 , res = 0.04471242607808054
it# 17 , res = 0.04335640617752859
it# 18 , res = 0.042050276004663906
it# 19 , res = 0.04079078483519091
it# 20 , res = 0.03957509574733402
it# 21 , res = 0.038400714010681514
it# 22 , res = 0.037265429377108
it# 23 , res = 0.03616726933615244
it# 24 , res = 0.035104461082522506
it# 25 , res = 0.03407540044859192
it# 26 , res = 0.03307862643363848
it# 27 , res = 0.03211280024986049
it# 28 , res = 0.03117668802720972
it# 29 , res = 0.030269146491763346
it# 30 , res = 0.02938911106779195
it# 31 , res = 0.028535585960640916
it# 32 , res = 0.027707635862510547
it# 33 , res = 0.026904378991044624
it# 34 , res = 0.026124981225003122
it# 35 , res = 0.02536865114502756
it# 36 , res = 0.024634635822794803
it# 37 , res = 0.02392221723040188
it# 38 , res = 0.023230709164974054
it# 39 , res = 0.022559454602299456
it# 40 , res = 0.02190782340860878
it# 41 , res = 0.021275210352110597
it# 42 , res = 0.02066103336609536
it# 43 , res = 0.02006473202377378
it# 44 , res = 0.019485766191853843
it# 45 , res = 0.018923614835480995
it# 46 , res = 0.018377774951776782
it# 47 , res = 0.01784776061301662
it# 48 , res = 0.01733310210361498
it# 49 , res = 0.0168333451376772
it# 50 , res = 0.016348050146014214
it# 51 , res = 0.015876791623289133
it# 52 , res = 0.015419157527437818
it# 53 , res = 0.014974748724722192
it# 54 , res = 0.014543178474801867
it# 55 , res = 0.014124071951048981
it# 56 , res = 0.013717065792048417
it# 57 , res = 0.013321807680817796
it# 58 , res = 0.012937955948776104
it# 59 , res = 0.012565179201917688
it# 60 , res = 0.012203155966991192
it# 61 , res = 0.011851574355791103
it# 62 , res = 0.011510131745914625
it# 63 , res = 0.011178534476553445
it# 64 , res = 0.010856497558072541
it# 65 , res = 0.010543744394280774
it# 66 , res = 0.010240006516432799
it# 67 , res = 0.009945023328111709
it# 68 , res = 0.009658541860242379
it# 69 , res = 0.009380316535567357
it# 70 , res = 0.009110108941986187
it# 71 , res = 0.008847687614229383
it# 72 , res = 0.008592827823380377
it# 73 , res = 0.00834531137381751
it# 74 , res = 0.008104926407179153
it# 75 , res = 0.007871467212997921
it# 76 , res = 0.00764473404567534
it# 77 , res = 0.007424532947499521
it# 78 , res = 0.007210675577431674
it# 79 , res = 0.0070029790454063205
it# 80 , res = 0.006801265751912509
it# 81 , res = 0.006605363232636929
it# 82 , res = 0.006415104007965919
it# 83 , res = 0.006230325437158621
it# 84 , res = 0.006050869577009597
it# 85 , res = 0.005876583044839298
it# 86 , res = 0.005707316885649552
it# 87 , res = 0.005542926443300198
it# 88 , res = 0.005383271235562818
it# 89 , res = 0.005228214832920568
it# 90 , res = 0.0050776247409858545
it# 91 , res = 0.004931372286416431
it# 92 , res = 0.00478933250621343
it# 93 , res = 0.004651384040293125
it# 94 , res = 0.004517409027226128
it# 95 , res = 0.004387293003042811
it# 96 , res = 0.0042609248030101
it# 97 , res = 0.0041381964662856605
it# 98 , res = 0.004019003143360028
it# 99 , res = 0.003903243006201777
it# 100 , res = 0.0037908171610226324
it# 101 , res = 0.0036816295635828753
it# 102 , res = 0.0035755869369603816
it# 103 , res = 0.0034725986917093117
it# 104 , res = 0.0033725768483364735
it# 105 , res = 0.003275435962026333
it# 106 , res = 0.003181093049548221
it# 107 , res = 0.0030894675182804375
it# 108 , res = 0.0030004810972890815
it# 109 , res = 0.002914057770401152
it# 110 , res = 0.0028301237112129704
it# 111 , res = 0.0027486072199768546
it# 112 , res = 0.002669438662312291
it# 113 , res = 0.00259255040968693
it# 114 , res = 0.0025178767816157674
it# 115 , res = 0.0024453539895301158
it# 116 , res = 0.002374920082265268
it# 117 , res = 0.0023065148931221834
it# 118 , res = 0.002240079988455231
it# 119 , res = 0.0021755586177435657
it# 120 , res = 0.0021128956651026687
it# 121 , res = 0.002052037602192773
it# 122 , res = 0.0019929324424863694
it# 123 , res = 0.0019355296968528836
it# 124 , res = 0.0018797803304245418
it# 125 , res = 0.0018256367207049916
it# 126 , res = 0.0017730526168851044
it# 127 , res = 0.0017219831003313837
it# 128 , res = 0.0016723845462128155
it# 129 , res = 0.0016242145862339758
it# 130 , res = 0.0015774320724403052
it# 131 , res = 0.0015319970420685299
it# 132 , res = 0.0014878706834080089
it# 133 , res = 0.0014450153026454651
it# 134 , res = 0.0014033942916666217
it# 135 , res = 0.0013629720967834056
it# 136 , res = 0.001323714188363418
it# 137 , res = 0.001285587031333257
it# 138 , res = 0.0012485580565322735
it# 139 , res = 0.001212595632891133
it# 140 , res = 0.00117766904041212
it# 141 , res = 0.0011437484439267506
it# 142 , res = 0.001110804867610847
it# 143 , res = 0.0010788101702321963
it# 144 , res = 0.0010477370211122077
it# 145 , res = 0.0010175588767790374
it# 146 , res = 0.0009882499582943246
it# 147 , res = 0.0009597852292320703
it# 148 , res = 0.0009321403742917895
it# 149 , res = 0.0009052917785285848
it# 150 , res = 0.0008792165071802011
it# 151 , res = 0.0008538922860765071
it# 152 , res = 0.000829297482611704
it# 153 , res = 0.0008054110872657125
it# 154 , res = 0.0007822126956578976
it# 155 , res = 0.0007596824911163184
it# 156 , res = 0.0007378012277509669
it# 157 , res = 0.0007165502140129318
it# 158 , res = 0.0006959112967280363
it# 159 , res = 0.0006758668455897079
it# 160 , res = 0.0006563997380988923
it# 161 , res = 0.0006374933449383653
it# 162 , res = 0.0006191315157664264
it# 163 , res = 0.0006012985654220176
it# 164 , res = 0.0005839792605253681
it# 165 , res = 0.0005671588064660694
it# 166 , res = 0.0005508228347642514
it# 167 , res = 0.000534957390797723
it# 168 , res = 0.0005195489218806568
it# 169 , res = 0.0005045842656874794
it# 170 , res = 0.0004900506390086021
it# 171 , res = 0.000475935626831136
it# 172 , res = 0.00046222717173360233
it# 173 , res = 0.000448913563585921
it# 174 , res = 0.0004359834295468369
it# 175 , res = 0.00042342572434870277
it# 176 , res = 0.000411229720862241
it# 177 , res = 0.0003993850009330292
it# 178 , res = 0.0003878814464829984
it# 179 , res = 0.0003767092308660354
it# 180 , res = 0.0003658588104745886
it# 181 , res = 0.0003553209165871722
it# 182 , res = 0.00034508654745068514
it# 183 , res = 0.0003351469605908874
it# 184 , res = 0.00032549366534456015
it# 185 , res = 0.0003161184156063097
it# 186 , res = 0.0003070132027844359
it# 187 , res = 0.0002981702489604665
it# 188 , res = 0.0002895820002440241
it# 189 , res = 0.00028124112032107244
it# 190 , res = 0.00027314048418646154
it# 191 , res = 0.0002652731720575917
it# 192 , res = 0.00025763246346382053
it# 193 , res = 0.0002502118315049693
it# 194 , res = 0.0002430049372764945
it# 195 , res = 0.00023600562445425888
it# 196 , res = 0.00022920791403574262
it# 197 , res = 0.00022260599923273447
it# 198 , res = 0.00021619424051071807
it# 199 , res = 0.00020996716077184696
it# 200 , res = 0.0002039194406758542
it# 201 , res = 0.0001980459140964641
it# 202 , res = 0.00019234156370815387
it# 203 , res = 0.00018680151670015088
it# 204 , res = 0.00018142104061437403
it# 205 , res = 0.00017619553930198762
it# 206 , res = 0.0001711205489985628
it# 207 , res = 0.000166191734509535
it# 208 , res = 0.00016140488550808998
it# 209 , res = 0.0001567559129379464
it# 210 , res = 0.00015224084552079042
it# 211 , res = 0.0001478558263635987
it# 212 , res = 0.0001435971096640631
it# 213 , res = 0.00013946105751087983
it# 214 , res = 0.00013544413677633575
it# 215 , res = 0.00013154291609768467
it# 216 , res = 0.0001277540629463968
it# 217 , res = 0.00012407434078137323
it# 218 , res = 0.00012050060628417368
it# 219 , res = 0.00011702980667381509
it# 220 , res = 0.00011365897709925741
it# 221 , res = 0.00011038523810641912
it# 222 , res = 0.000107205793178504
it# 223 , res = 0.00010411792634770015
it# 224 , res = 0.00010111899987402726
it# 225 , res = 9.820645199349388e-05
it# 226 , res = 9.537779472855718e-05
it# 227 , res = 9.263061176337818e-05
it# 228 , res = 8.996255637990053e-05
it# 229 , res = 8.737134945268324e-05
it# 230 , res = 8.485477750265274e-05
it# 231 , res = 8.241069080560279e-05
it# 232 , res = 8.003700155647155e-05
it# 233 , res = 7.77316820855138e-05
it# 234 , res = 7.549276312627895e-05
it# 235 , res = 7.33183321335028e-05
it# 236 , res = 7.120653164910118e-05
it# 237 , res = 6.915555771594402e-05
it# 238 , res = 6.716365833649476e-05
it# 239 , res = 6.522913197612263e-05
it# 240 , res = 6.335032610998692e-05
it# 241 , res = 6.152563581107202e-05
it# 242 , res = 5.975350237966505e-05
it# 243 , res = 5.8032412011094095e-05
it# 244 , res = 5.636089450333515e-05
it# 245 , res = 5.473752200070755e-05
it# 246 , res = 5.316090777449845e-05
it# 247 , res = 5.162970503798831e-05
it# 248 , res = 5.014260579634797e-05
it# 249 , res = 4.8698339728960004e-05
it# 250 , res = 4.729567310464032e-05
it# 251 , res = 4.593340772732194e-05
it# 252 , res = 4.461037991320727e-05
it# 253 , res = 4.332545949600582e-05
it# 254 , res = 4.2077548861643916e-05
it# 255 , res = 4.08655820111208e-05
it# 256 , res = 3.968852364956292e-05
it# 257 , res = 3.854536830167224e-05
it# 258 , res = 3.7435139453189796e-05
it# 259 , res = 3.6356888716430046e-05
it# 260 , res = 3.530969502045383e-05
it# 261 , res = 3.429266382387839e-05
it# 262 , res = 3.330492635114455e-05
it# 263 , res = 3.2345638850267564e-05
it# 264 , res = 3.1413981871493425e-05
it# 265 , res = 3.0509159568377468e-05
it# 266 , res = 2.963039901713872e-05
it# 267 , res = 2.8776949556638008e-05
it# 268 , res = 2.7948082147119694e-05
it# 269 , res = 2.7143088747703462e-05
it# 270 , res = 2.6361281711075542e-05
it# 271 , res = 2.5601993196578055e-05
it# 272 , res = 2.4864574599332136e-05
it# 273 , res = 2.4148395996170652e-05
it# 274 , res = 2.3452845608119163e-05
it# 275 , res = 2.277732927697544e-05
it# 276 , res = 2.2121269958512952e-05
it# 277 , res = 2.148410722922485e-05
it# 278 , res = 2.0865296806944827e-05
it# 279 , res = 2.0264310087555657e-05
it# 280 , res = 1.968063369145718e-05
it# 281 , res = 1.911376902662848e-05
it# 282 , res = 1.8563231861865738e-05
it# 283 , res = 1.802855191348807e-05
it# 284 , res = 1.7509272443225045e-05
it# 285 , res = 1.7004949868568575e-05
it# 286 , res = 1.6515153383502084e-05
it# 287 , res = 1.6039464590379097e-05
it# 288 , res = 1.5577477143179646e-05
it# 289 , res = 1.5128796399581022e-05
it# 290 , res = 1.4693039084396657e-05
it# 291 , res = 1.426983296184228e-05
it# 292 , res = 1.385881651752018e-05
it# 293 , res = 1.3459638650558788e-05
it# 294 , res = 1.3071958371886972e-05
it# 295 , res = 1.2695444514810404e-05
it# 296 , res = 1.232977545085969e-05
it# 297 , res = 1.1974638815751068e-05
it# 298 , res = 1.1629731241847108e-05
it# 299 , res = 1.1294758100014288e-05
it# 300 , res = 1.0969433246837716e-05
it# 301 , res = 1.0653478781279293e-05
it# 302 , res = 1.0346624806258583e-05
it# 303 , res = 1.0048609198970677e-05
it# 304 , res = 9.759177386197288e-06
it# 305 , res = 9.47808212752317e-06
it# 306 , res = 9.205083303808563e-06
it# 307 , res = 8.93994771193541e-06
it# 308 , res = 8.682448866006233e-06
it# 309 , res = 8.432366803544042e-06
it# 310 , res = 8.189487897434423e-06
it# 311 , res = 7.953604674120433e-06
it# 312 , res = 7.724515635677603e-06
it# 313 , res = 7.5020250880963685e-06
it# 314 , res = 7.28594297397562e-06
it# 315 , res = 7.0760847098981525e-06
it# 316 , res = 6.872271029318753e-06
it# 317 , res = 6.67432782927943e-06
it# 318 , res = 6.482086020977723e-06
it# 319 , res = 6.295381386642455e-06
it# 320 , res = 6.114054437644614e-06
it# 321 , res = 5.937950279925802e-06
it# 322 , res = 5.76691848053566e-06
it# 323 , res = 5.600812939242161e-06
it# 324 , res = 5.43949176438185e-06
it# 325 , res = 5.282817150910681e-06
it# 326 , res = 5.1306552632462264e-06
it# 327 , res = 4.982876120455496e-06
it# 328 , res = 4.839353485726994e-06
it# 329 , res = 4.699964757938274e-06
it# 330 , res = 4.564590867505201e-06
it# 331 , res = 4.433116174627793e-06
it# 332 , res = 4.305428369747959e-06
it# 333 , res = 4.1814183786555695e-06
it# 334 , res = 4.060980268562807e-06
it# 335 , res = 3.944011157935261e-06
it# 336 , res = 3.830411128625705e-06
it# 337 , res = 3.7200831405900385e-06
it# 338 , res = 3.6129329483024554e-06
it# 339 , res = 3.5088690214873126e-06
it# 340 , res = 3.407802465808122e-06
it# 341 , res = 3.30964694736853e-06
it# 342 , res = 3.2143186191879644e-06
it# 343 , res = 3.1217360491474662e-06
it# 344 , res = 3.0318201507096594e-06
it# 345 , res = 2.9444941153641406e-06
it# 346 , res = 2.859683346783644e-06
it# 347 , res = 2.7773153973727495e-06
it# 348 , res = 2.6973199060763634e-06
it# 349 , res = 2.6196285384985283e-06
it# 350 , res = 2.544174928831634e-06
it# 351 , res = 2.4708946223896937e-06
it# 352 , res = 2.399725021125423e-06
it# 353 , res = 2.330605330197754e-06
it# 354 , res = 2.2634765054031638e-06
it# 355 , res = 2.1982812037232047e-06
it# 356 , res = 2.1349637335304492e-06
it# 357 , res = 2.073470007303043e-06
it# 358 , res = 2.013747495334684e-06
it# 359 , res = 1.9557451811077795e-06
it# 360 , res = 1.8994135175396976e-06
it# 361 , res = 1.8447043844748387e-06
it# 362 , res = 1.791571047919463e-06
it# 363 , res = 1.7399681200718407e-06
it# 364 , res = 1.6898515201307419e-06
it# 365 , res = 1.6411784372724611e-06
it# 366 , res = 1.5939072934139005e-06
it# 367 , res = 1.5479977084625727e-06
it# 368 , res = 1.5034104652990243e-06
it# 369 , res = 1.460107475987965e-06
it# 370 , res = 1.4180517501337712e-06
it# 371 , res = 1.377207362576874e-06
it# 372 , res = 1.337539422962326e-06
it# 373 , res = 1.299014045705442e-06
it# 374 , res = 1.261598321594467e-06
it# 375 , res = 1.2252602888374641e-06
it# 376 , res = 1.1899689068495104e-06
it# 377 , res = 1.155694028462588e-06
it# 378 , res = 1.1224063752871304e-06
it# 379 , res = 1.0900775122129035e-06
it# 380 , res = 1.0586798227893458e-06
it# 381 , res = 1.0281864864196623e-06
it# 382 , res = 9.98571454936751e-07
it# 383 , res = 9.698094301057029e-07
it# 384 , res = 9.41875842865459e-07
it# 385 , res = 9.147468316638282e-07
it# 386 , res = 8.883992217879801e-07
it# 387 , res = 8.62810506821342e-07
it# 388 , res = 8.379588280324918e-07
it# 389 , res = 8.138229562025276e-07
it# 390 , res = 7.903822740934107e-07
it# 391 , res = 7.676167580597427e-07
it# 392 , res = 7.455069610100994e-07
it# 393 , res = 7.240339962682999e-07
it# 394 , res = 7.031795209779901e-07
it# 395 , res = 6.829257208412981e-07
it# 396 , res = 6.632552942156443e-07
it# 397 , res = 6.441514383404186e-07
it# 398 , res = 6.255978341156058e-07
it# 399 , res = 6.075786324080275e-07
it# 400 , res = 5.900784409938462e-07
it# 401 , res = 5.730823104426305e-07
it# 402 , res = 5.565757224933843e-07
it# 403 , res = 5.405445765776248e-07
it# 404 , res = 5.249751785342237e-07
it# 405 , res = 5.098542286040506e-07
it# 406 , res = 4.951688097197936e-07
it# 407 , res = 4.809063776942166e-07
it# 408 , res = 4.670547489353102e-07
it# 409 , res = 4.536020909158384e-07
it# 410 , res = 4.405369121519964e-07
it# 411 , res = 4.27848051932019e-07
it# 412 , res = 4.155246711722246e-07
it# 413 , res = 4.035562429372338e-07
it# 414 , res = 3.9193254331669837e-07
it# 415 , res = 3.806436430240624e-07
it# 416 , res = 3.6967989901723137e-07
it# 417 , res = 3.590319455903584e-07
it# 418 , res = 3.48690686791875e-07
it# 419 , res = 3.386472892449078e-07
it# 420 , res = 3.2889317342912394e-07
it# 421 , res = 3.194200068964767e-07
it# 422 , res = 3.1021969772915295e-07
it# 423 , res = 3.0128438662569666e-07
it# 424 , res = 2.9260644081936364e-07
it# 425 , res = 2.8417844734525875e-07
it# 426 , res = 2.759932068277077e-07
it# 427 , res = 2.6804372721340723e-07
it# 428 , res = 2.6032321786947923e-07
it# 429 , res = 2.528250835438541e-07
it# 430 , res = 2.455429192625459e-07
it# 431 , res = 2.384705045500921e-07
it# 432 , res = 2.3160179766467693e-07
it# 433 , res = 2.2493093150215811e-07
it# 434 , res = 2.1845220725729344e-07
it# 435 , res = 2.1216009095428013e-07
it# 436 , res = 2.06049207465847e-07
it# 437 , res = 2.001143368212774e-07
it# 438 , res = 1.9435040940013613e-07
it# 439 , res = 1.8875250121518675e-07
it# 440 , res = 1.8331583051036622e-07
it# 441 , res = 1.780357532459033e-07
it# 442 , res = 1.7290775888793225e-07
it# 443 , res = 1.6792746709421952e-07
it# 444 , res = 1.630906234740257e-07
it# 445 , res = 1.5839309627610792e-07
it# 446 , res = 1.5383087284536173e-07
it# 447 , res = 1.4940005581862628e-07
it# 448 , res = 1.450968604686929e-07
it# 449 , res = 1.4091761072212966e-07
it# 450 , res = 1.3685873658733965e-07
it# 451 , res = 1.3291677096719276e-07
it# 452 , res = 1.2908834643443196e-07
it# 453 , res = 1.2537019268676804e-07
it# 454 , res = 1.2175913347842722e-07
it# 455 , res = 1.182520843648559e-07
it# 456 , res = 1.1484604931170288e-07
it# 457 , res = 1.1153811884270967e-07
it# 458 , res = 1.0832546732694544e-07
it# 459 , res = 1.0520535028824282e-07
it# 460 , res = 1.0217510259304374e-07
it# 461 , res = 9.923213561313785e-08
it# 462 , res = 9.637393544003841e-08
it# 463 , res = 9.359806036314324e-08
it# 464 , res = 9.090213934831149e-08
it# 465 , res = 8.828386943033343e-08
it# 466 , res = 8.574101406605229e-08
it# 467 , res = 8.327140101228209e-08
it# 468 , res = 8.08729205629197e-08
it# 469 , res = 7.854352396630114e-08
it# 470 , res = 7.628122141877191e-08
it# 471 , res = 7.408408032738147e-08
it# 472 , res = 7.19502238635781e-08
it# 473 , res = 6.987782925465512e-08
it# 474 , res = 6.78651261140824e-08
it# 475 , res = 6.591039528299981e-08
it# 476 , res = 6.401196684682435e-08
it# 477 , res = 6.216821923553513e-08
it# 478 , res = 6.037757730425827e-08
it# 479 , res = 5.863851156660593e-08
it# 480 , res = 5.694953651983629e-08
it# 481 , res = 5.530920923224853e-08
it# 482 , res = 5.371612870335364e-08
it# 483 , res = 5.216893372971895e-08
it# 484 , res = 5.066630305022126e-08
it# 485 , res = 4.920695284076212e-08
it# 486 , res = 4.778963657567523e-08
it# 487 , res = 4.641314342335096e-08
it# 488 , res = 4.507629770521025e-08
it# 489 , res = 4.3777957383212005e-08
it# 490 , res = 4.251701339416773e-08
it# 491 , res = 4.1292388620330046e-08
it# 492 , res = 4.010303679407674e-08
it# 493 , res = 3.894794220887215e-08
it# 494 , res = 3.7826117927085566e-08
it# 495 , res = 3.673660572708499e-08
it# 496 , res = 3.5678475014837906e-08
it# 497 , res = 3.465082180803267e-08
it# 498 , res = 3.3652768215043675e-08
it# 499 , res = 3.2683461727862234e-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.12203689059859804
CG iteration 3, residual = 0.08630349626448036
CG iteration 4, residual = 0.05476506137338163
CG iteration 5, residual = 0.027039871611641366
CG iteration 6, residual = 0.014941028668260053
CG iteration 7, residual = 0.00669463091431055
CG iteration 8, residual = 0.002848398062743388
CG iteration 9, residual = 0.0009877015081514702
CG iteration 10, residual = 0.0003968124927947859
CG iteration 11, residual = 0.00012407544574560334
CG iteration 12, residual = 6.330209216165877e-05
CG iteration 13, residual = 2.9785257662869865e-05
CG iteration 14, residual = 1.3390306292216046e-05
CG iteration 15, residual = 5.363365741544917e-06
CG iteration 16, residual = 2.046266192471712e-06
CG iteration 17, residual = 8.231874195477879e-07
CG iteration 18, residual = 3.454456776090241e-07
CG iteration 19, residual = 1.2951230160193455e-07
CG iteration 20, residual = 5.3553595524320823e-08
CG iteration 21, residual = 1.761101316835354e-08
CG iteration 22, residual = 5.970141704993232e-09
CG iteration 23, residual = 2.4674252558032045e-09
CG iteration 24, residual = 7.865934453395222e-10
CG iteration 25, residual = 3.1722300643176456e-10
CG iteration 26, residual = 1.1006534206433609e-10
CG iteration 27, residual = 4.095450344654651e-11
CG iteration 28, residual = 1.5996345010343813e-11
CG iteration 29, residual = 5.326653784450928e-12
CG iteration 30, residual = 1.9887046841316592e-12
CG iteration 31, residual = 5.843515407393082e-13
CG iteration 32, residual = 2.2678569699607902e-13
CG iteration 33, residual = 6.065405903982064e-14

WebGLScene

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

20.3116176251635

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

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

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 0x7fbeb17a21b0>

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.