This page was generated from unit-2.1.3-bddc/bddc.ipynb.

2.1.3 Element-wise BDDC Preconditioner

The element-wise BDDC (Balancing Domain Decomposition preconditioner with Constraints) preconditioner in NGSolve is a good general purpose preconditioner that works well both in the shared memory parallel mode as well as in distributed memory mode. In this tutorial, we discuss this preconditioner, related built-in options, and customization from python.

Let us start with a simple description of the element-wise BDDC in the context of Lagrange finite element space \(V\). Here the BDDC preconditoner is constructed on an auxiliary space \(\widetilde{V}\) obtained by connecting only element vertices (leaving edge and face shape functions discontinuous). Although larger, the auxiliary space allows local elimination of edge and face variables. Hence an analogue of the original matrix \(A\) on this space, named \(\widetilde A\), is less expensive to invert. This inverse is used to construct a preconditioner for \(A\) as follows:

\[C_{BDDC}^{-1} = R {\,\widetilde{A}\,}^{-1}\, R^t\]

Here, \(R\) is the averaging operator for the discontinous edge and face variables.

To construct a general purpose BDDC preconditioner, NGSolve generalizes this idea to all its finite element spaces by a classification of degrees of freedom. Dofs are classified into (condensable) LOCAL_DOFs that we saw in 1.4 and a remainder that includes these types:

WIREBASKET_DOF
INTERFACE_DOF

The original finite element space \(V\) is obtained by requiring conformity of both the above types of dofs, while the auxiliary space \(\widetilde{V}\) is obtained by requiring conformity of WIREBASKET_DOFs only.

Here is a figure of a typical function in the default \(\widetilde{V}\) (and the code to generate this is at the end of this tutorial) when \(V\) is the Lagrange space:

title

[1]:

import netgen.gui
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner
from netgen.csg import unit_cube
from netgen.geom2d import unit_square
SetHeapSize(100*1000*1000)

[2]:

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.5))
# mesh = Mesh(unit_square.GenerateMesh(maxh=0.5))

Built-in options

Let us define a simple function to study how the spectrum of the preconditioned matrix changes with various options.

Effect of condensation

[3]:

def TestPreconditioner (p, condense=False, **args):
    fes = H1(mesh, order=p, **args)
    u,v = fes.TnT()
    a = BilinearForm(fes, eliminate_internal=condense)
    a += grad(u)*grad(v)*dx + u*v*dx
    c = Preconditioner(a, "bddc")
    a.Assemble()
    return EigenValues_Preconditioner(a.mat, c.mat)

[4]:

lams = TestPreconditioner(5)
print (lams[0:3], "...\n", lams[-3:])

 1.00408
 1.11168
 1.34769
 ...
  33.7436
 33.8948
 38.0012

Here is the effect of static condensation on the BDDC preconditioner.

[5]:

lams = TestPreconditioner(5, condense=True)
print (lams[0:3], "...\n", lams[-3:])

 1.00139
 1.04212
 1.12173
 ...
  6.95459
 7.03965
  7.3257

Tuning the auxiliary space

Next, let us study the effect of a few built-in flags for finite element spaces that are useful for tweaking the behavior of the BDDC preconditioner. The effect of these flags varies in two (2D) and three dimensions (3D), e.g.,

  • wb_fulledges=True: This option keeps all edge-dofs conforming (i.e. they are marked WIREBASKET_DOFs). This option is only meaningful in 3D. If used in 2D, the preconditioner becomes a direct solver.

  • wb_withedges=True: This option keeps only the first edge-dof conforming (i.e., the first edge-dof is marked WIREBASKET_DOF and the remaining edge-dofs are marked INTERFACE_DOFs).

The complete situation is a bit more complex due to the fact these options can take the three values True, False, or Undefined, the two options can be combined, and the space dimension can be 2 or 3. The default value of these flags in NGSolve is Undefined. Here is a table with the summary of the effect of these options:

wb_fulledges

wb_withedges

2D

3D

True

any value

all

all

False/Undefined

Undefined

none

first

False/Undefined

False

none

none

False/Undefined

True

first

first

An entry \(X \in\) {all, none, first} of the last two columns is to be read as follows: \(X\) of the edge-dofs is(are) WIREBASKET_DOF(s).

Here is an example of the effect of one of these flag values.

[6]:

lams = TestPreconditioner(5, condense=True,
                          wb_withedges=False)
print (lams[0:3], "...\n", lams[-3:])

 1.00468
 1.09359
 1.32949
 ...
  31.1877
 33.3524
 34.3632

Clearly, when moving from the default case (where the first of the edge dofs are wire basket dofs) to the case (where none of the edge dofs are wire basket dofs), the conditioning became less favorable.

Customize

From within python, we can change the types of degrees of freedom of finite element spaces, thus affecting the behavior of the BDDC preconditioner.

To depart from the default and commit the first two edge dofs to wire basket, we perform the next steps:

[7]:

fes = H1(mesh, order=10)
u,v = fes.TnT()

for ed in mesh.edges:
    dofs = fes.GetDofNrs(ed)
    for d in dofs:
        fes.SetCouplingType(d, COUPLING_TYPE.INTERFACE_DOF)

    # Set the first two edge dofs to be conforming
    fes.SetCouplingType(dofs[0], COUPLING_TYPE.WIREBASKET_DOF)
    fes.SetCouplingType(dofs[1], COUPLING_TYPE.WIREBASKET_DOF)

a = BilinearForm(fes, eliminate_internal=True)
a += grad(u)*grad(v)*dx + u*v*dx
c = Preconditioner(a, "bddc")
a.Assemble()

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

[7]:

18.310200831382026

This is a slight improvement from the default.

[8]:

lams = TestPreconditioner (10, condense=True)
max(lams)/min(lams)

[8]:

27.550548160831678

Combine BDDC with AMG for large problems

The coarse inverse \({\,\widetilde{A}\,}^{-1}\) of BDDC is expensive on fine meshes. Using the option coarsetype=h1amg flag, we can ask BDDC to replace \({\,\widetilde{A}\,}^{-1}\) by an Algebraic MultiGrid (AMG) preconditioner. Since NGSolve’s h1amg is well-suited
for the lowest order Lagrange space, we use the option wb_withedges=False to ensure that \(\widetilde{A}\) is made solely with vertex unknowns.
[9]:

p = 5
mesh = Mesh(unit_cube.GenerateMesh(maxh=0.05))
fes = H1(mesh, order=p, dirichlet="left|bottom|back",
         wb_withedges=False)
print("NDOF = ", fes.ndof)
u,v = fes.TnT()
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
f = LinearForm(fes)
f += v*dx

with TaskManager():
    pre = Preconditioner(a, "bddc", coarsetype="h1amg")
    a.Assemble()
    f.Assemble()

    gfu = GridFunction(fes)
    solvers.CG(mat=a.mat, rhs=f.vec, sol=gfu.vec,
               pre=pre, maxsteps=500)
Draw(gfu)

NDOF =  1165206

iteration 0 error = 0.7073079321639637
iteration 1 error = 0.29234899259626856
iteration 2 error = 0.2850423590749113
iteration 3 error = 0.2909985858220177
iteration 4 error = 0.3196033806324102
iteration 5 error = 0.2548899214031628
iteration 6 error = 0.19211663159074388
iteration 7 error = 0.15448092069205674
iteration 8 error = 0.1262939610830819
iteration 9 error = 0.10147830081475578
iteration 10 error = 0.07424583493548162
iteration 11 error = 0.056624349933461396
iteration 12 error = 0.04501661464420658
iteration 13 error = 0.03822244147240056
iteration 14 error = 0.03168341537618148
iteration 15 error = 0.02517295699223492
iteration 16 error = 0.019311927758283772
iteration 17 error = 0.01554005065209648
iteration 18 error = 0.01242698420061178
iteration 19 error = 0.009988699609349479
iteration 20 error = 0.0077720965317934895
iteration 21 error = 0.006098988979459672
iteration 22 error = 0.004845023443492067
iteration 23 error = 0.0038936689458989302
iteration 24 error = 0.0031211786508215595
iteration 25 error = 0.002514683062479327
iteration 26 error = 0.0020184587220798333
iteration 27 error = 0.001617829179391698
iteration 28 error = 0.001292397220162447
iteration 29 error = 0.0010424470017298738
iteration 30 error = 0.0008356482973016168
iteration 31 error = 0.0006714480065978774
iteration 32 error = 0.0005352623318995764
iteration 33 error = 0.0004215503628216195
iteration 34 error = 0.0003333305814681214
iteration 35 error = 0.0002673703978420636
iteration 36 error = 0.00021292667233358957
iteration 37 error = 0.00016781697413274704
iteration 38 error = 0.00013450312384210065
iteration 39 error = 0.00010875174584049092
iteration 40 error = 8.771907984040991e-05
iteration 41 error = 7.041354864255918e-05
iteration 42 error = 5.698975226821377e-05
iteration 43 error = 4.549968019506618e-05
iteration 44 error = 3.658807446970773e-05
iteration 45 error = 2.9152587032893793e-05
iteration 46 error = 2.3080793175006922e-05
iteration 47 error = 1.870094550457831e-05
iteration 48 error = 1.5030440049793567e-05
iteration 49 error = 1.2051856509027753e-05
iteration 50 error = 9.778231424356836e-06
iteration 51 error = 7.826482451346145e-06
iteration 52 error = 6.201639849639167e-06
iteration 53 error = 4.93963662333325e-06
iteration 54 error = 3.9148196875022585e-06
iteration 55 error = 3.1363011982838677e-06
iteration 56 error = 2.527140301843146e-06
iteration 57 error = 2.0410775677569746e-06
iteration 58 error = 1.6311701124297936e-06
iteration 59 error = 1.2845984617590747e-06
iteration 60 error = 1.0402375507476477e-06
iteration 61 error = 8.352505853166585e-07
iteration 62 error = 6.684146548764809e-07
iteration 63 error = 5.326854521179992e-07
iteration 64 error = 4.261471708895145e-07
iteration 65 error = 3.405062250303084e-07
iteration 66 error = 2.7074679827502577e-07
iteration 67 error = 2.1748059061703744e-07
iteration 68 error = 1.74868732837972e-07
iteration 69 error = 1.4059048256368237e-07
iteration 70 error = 1.1189128227516209e-07
iteration 71 error = 8.936904517096745e-08
iteration 72 error = 7.124604070123911e-08
iteration 73 error = 5.6448886926153085e-08
iteration 74 error = 4.559204075472101e-08
iteration 75 error = 3.6242751399217374e-08
iteration 76 error = 2.9297934878590947e-08
iteration 77 error = 2.3414802424909063e-08
iteration 78 error = 1.8471477201606494e-08
iteration 79 error = 1.4932632769238514e-08
iteration 80 error = 1.2092166506946097e-08
iteration 81 error = 9.765211792818168e-09
iteration 82 error = 8.010318939169762e-09
iteration 83 error = 6.529857520348664e-09
iteration 84 error = 5.2312634219840075e-09
iteration 85 error = 4.229971609082354e-09
iteration 86 error = 3.3963753192277273e-09
iteration 87 error = 2.755033461276114e-09
iteration 88 error = 2.174719486704714e-09
iteration 89 error = 1.7382346321335394e-09
iteration 90 error = 1.3897128952418485e-09
iteration 91 error = 1.1169138584896175e-09
iteration 92 error = 8.947100784381828e-10
iteration 93 error = 7.209509157574718e-10
iteration 94 error = 5.787644335896124e-10
iteration 95 error = 4.571069125557547e-10
iteration 96 error = 3.6641940320910413e-10
iteration 97 error = 2.945231786975309e-10
iteration 98 error = 2.34054976000148e-10
iteration 99 error = 1.8873870394297015e-10
iteration 100 error = 1.5096443829349653e-10
iteration 101 error = 1.2060518392866795e-10
iteration 102 error = 9.80834165396866e-11
iteration 103 error = 8.064042051881531e-11
iteration 104 error = 6.55627072133953e-11
iteration 105 error = 5.281077636414923e-11
iteration 106 error = 4.1934538685554554e-11
iteration 107 error = 3.3373902448792593e-11
iteration 108 error = 2.6432525787330204e-11
iteration 109 error = 2.131706035755686e-11
iteration 110 error = 1.7168188143119886e-11
iteration 111 error = 1.3753308386351179e-11
iteration 112 error = 1.1161014345903987e-11
iteration 113 error = 8.89884299169311e-12
iteration 114 error = 7.041278700003329e-12
iteration 115 error = 5.65236769157018e-12
iteration 116 error = 4.5006014246596205e-12
iteration 117 error = 3.5772226481261264e-12
iteration 118 error = 2.8852844436709226e-12
iteration 119 error = 2.303959741082989e-12
iteration 120 error = 1.8193225627423864e-12
iteration 121 error = 1.452934986533004e-12
iteration 122 error = 1.1756420347486093e-12
iteration 123 error = 9.443464678463255e-13
iteration 124 error = 7.582123510783575e-13
iteration 125 error = 6.097304021046931e-13

Postscript

By popular demand, here is the code to draw the figure found at the beginning of this tutorial:

[10]:

from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes_ho = Discontinuous(H1(mesh, order=10))
fes_lo = H1(mesh, order=1, dirichlet=".*")
fes_lam = Discontinuous(H1(mesh, order=1))
fes = FESpace([fes_ho, fes_lo, fes_lam])
uho, ulo, lam = fes.TrialFunction()

a = BilinearForm(fes)
a += Variation(0.5 * grad(uho)*grad(uho)*dx
               - 1*uho*dx
               + (uho-ulo)*lam*dx(element_vb=BBND))
gfu = GridFunction(fes)
solvers.Newton(a=a, u=gfu)
Draw(gfu.components[0])

Newton iteration  0
err =  0.3937351156371161
Newton iteration  1
err =  3.874027236775335e-15